东莞建网站公司游戏代理是怎么赚钱的如何代理游戏

前言

最近一直在研究振动信号的最优解调带宽寻找方法，而快速谱峭度算法是最常见的一种，因此本篇文章准备对该算法做个深入解析，也相当于自己的学习记录了。

原理解析

谱峭度是一种信号处理中的时频分析工具，用于量化信号在不同频率成分上的非高斯性或冲击特性。它通过在频域中计算每个频率段上的峭度（统计学中衡量数据分布陡峭程度的指标，高峭度表示数据具有更尖锐的峰值和更重的尾部，如冲击信号），反映信号中瞬态成分或异常事件的分布情况。
快速谱峭度算法则是Antoni提出的一种高效计算谱峭度的方法，核心思想是通过多分辨率滤波器组快速分解信号频带，并优化峭度值的计算流程，从而显著提升效率。
1、信号分解
通过多级滤波器组将原始信号分解为多个子频带，目前较为常用的是1/3-二叉树的结合。

2、频带峭度计算
对每个子频带信号$x_i(t)$计算其经典峭度值$K_i$。
$$K_i=\frac{mean(|x_i(t){|}^4)}{mean(|x_i(t){|}^2{)}^2}$$
快速谱峭度算法里还提供了一种基于2阶统计量的鲁棒峭度值。
$$K_i=\frac{mean(|x_i(t){|}^2)}{mean(|x_i(t)|{)}^2}$$
3、生成峭度图
将各频带的峭度值映射到二维平面（频率 vs. 分解层数），标记峭度最大值对应的最优频带。
4、最优频带提取
根据峭度图的指示，选择峭度最大的频带进行带通滤波，提取冲击成分。

代码解析

快速谱峭度算法有一版Matlab程序，由Antoni老师编写的，而后笔者也找到一份Lecocq老师改写的Python版本，该代码同样通过树状滤波器来实现（也可以通过短时傅里叶变换、小波包分解等方式实现，核心就是将原信号分解到不同的频段，计算各频段的谱峭度指标，找到指标最大的频段）。
为了避免影响分析，在原有的函数基础上保留了最核心的逻辑，最完整的代码会在文末给出。

Fast_Kurtogram入口函数

def Fast_Kurtogram(x, nlevel, verbose=False, Fs=1, NFIR=16, fcut=0.4,opt1=1, opt2=1):"""计算信号x分解到第nlevel层的快速谱峭度，最大的分解层数是log2(length(x))，但是推荐低于该最大层数的1/8倍，返回Also returns the vector of k-levels Level_w, the frequency vectorfreq_w, the complex envelope of the signal c and the extremefrequencies of the "best" bandpass f_lower and f_upper.:param x: 待分析的信号，numpy数组:param nlevel: 分解的层数，整数:param verbose: 是否输出调试信息:param Fs: 待分析信号的采样率:param NFIR: FIR滤波器的长度，整数:param fcut: 滤波器的截止频率，浮点数:param opt1: [1 | 2]:* opt1 = 1: 基于4阶统计量的经典峭度* opt1 = 2: 基于2阶统计量的鲁棒峭度（如果这两个峭度有差异，则是因为冲击噪声的存在）:param opt2: [1 | 2]:* opt2=1: 基于树状滤波器实现* opt2=2: 基于短时傅里叶算法实现，暂未实现（方法1更快，并且在滤波器的设计上更灵活：一个短的滤波器和方法2可以实现近乎一样的效果）:returns: Kwav, Level_w, freq_w, c, f_lower, f_upper* Kwav: kurtogram* Level_w: vector of levels* freq_w: frequency vector* c: complex envelope of the signal filtered in the frequency band that maximizes the kurtogram* f_lower: lower frequency of the band pass* f_upper: upper frequency of the band pass"""N = len(x)N2 = np.log2(N) - 7if nlevel > N2:logging.error('Please enter a smaller number of decomposition levels')# Fast computation of the kurtogram####################################if opt1 == 1:# 1) Filterbank-based kurtogram############################# Analytic generating filters# h是二分段低通滤波器，g是二分段高通滤波器，h1是三分段第一段滤波器，h2是三分段第二段滤波器，h3是三分段第三段滤波器h, g, h1, h2, h3 = get_h_parameters(NFIR, fcut)# 获取所有分解频段的谱峭度if opt2 == 1:# kurtosis of the complex envelopeKwav = K_wpQ(x, h, g, h1, h2, h3, nlevel, verbose, 'kurt2')else:# variance of the envelope magnitudeKwav = K_wpQ(x, h, g, h1, h2, h3, nlevel, verbose, 'kurt1')# keep positive values only!# 计算最后一层各频段的频率，最大频率是采样率的一半Kwav[Kwav <= 0] = 0Level_w = np.arange(1, nlevel+1)Level_w = np.array([Level_w, Level_w + np.log2(3.)-1])Level_w = sorted(Level_w.ravel())Level_w = np.append(0, Level_w[0:2*nlevel-1])freq_w = Fs*(np.arange(0, 3*2.0**nlevel)/(3*2**(nlevel+1)) +1.0/(3.*2.**(2+nlevel)))M, index = get_GridMax(Kwav)level_index = index[0]freq_index = index[1]bw, fc, fi, l1 = getBandwidthAndFrequency(nlevel, Fs, Level_w, freq_w,level_index, freq_index)# 是否绘制快速谱峭度图if verbose:logging.info("max kur:{}".format(M))plotKurtogram(Kwav, freq_w, nlevel, Level_w, Fs, fi, level_index)else:logging.error('stft-based is not implemented')# Signal filtering !c = []test = 1lev = l1while test == 1:test = 0# 找出峭度最大的频段c, s, threshold, Bw, fc = Find_wav_kurt(x, h, g, h1, h2, h3,nlevel, lev, fi, Fs=Fs,verbose=verbose)# Determine the lowest and the uppest frequencies of the bandpassf_lower = Fs*np.round((fc-Bw/2.)*10**3)/10**3f_upper = Fs*np.round((fc+Bw/2.)*10**3)/10**3return Kwav, Level_w, freq_w, c, f_lower, f_upper

整个函数

get_h_parameters构造滤波器

该函数的作用就是构造二分段滤波器和三分段滤波器。

def get_h_parameters(NFIR, fcut):"""构造FIR滤波器:param NFIR: FIR滤波器长度，整数:param fcut: 滤波器截止频率，浮点数:returns: h-parameters: h（二分段低通滤波器）, g（二分段高通滤波器）, h1（三分段第一段滤波器）, h2（三分段第二段滤波器）, h3（三分段第三段滤波器）"""h = si.firwin(NFIR+1, fcut) * np.exp(2*1j*np.pi*np.arange(NFIR+1) * 0.125)n = np.arange(2, NFIR+2)g = h[(1-n) % NFIR]*(-1.)**(1-n)NFIR = int(np.fix((3./2.*NFIR)))h1 = si.firwin(NFIR+1, 2./3*fcut)*np.exp(2j*np.pi*np.arange(NFIR+1) *0.25/3.)h2 = h1*np.exp(2j*np.pi*np.arange(NFIR+1)/6.)h3 = h1*np.exp(2j*np.pi*np.arange(NFIR+1)/3.)return (h, g, h1, h2, h3)

由于笔者对FIR滤波器的设计并不熟悉，此处并未深入研究过其形式为何如此，感兴趣的朋友们可以自行搜索一下。

K_wpQ谱峭度计算函数

在获取到FIR滤波器参数后，下面就可以进行谱峭度的计算了，主要的函数就是K_wpQ。

def K_wpQ(x, h, g, h1, h2, h3, nlevel, verbose, opt, level=0):"""计算2维矩阵形式的谱峭度:param x: 待分析的信号:param h: 二分段低通滤波器:param g: 二分段高通滤波器:param h1: 三分段第一段滤波器:param h2: 三分段第二段滤波器:param h3: 三分段第三段滤波器:param nlevel: 分解层数:param verbose: 是否输出调试信息:param opt: ['kurt1' | 'kurt2']* 'kurt1' = 基于4阶统计量的经典峭度* 'kurt2' = 基于2阶统计量的鲁棒峭度:param level: 当前分解层数:returns: 二维矩阵的谱峭度"""L = np.floor(np.log2(len(x)))if level == 0:if nlevel >= L:logging.error('nlevel must be smaller')level = nlevelx = x.ravel()# 计算谱峭度KD, KQ = K_wpQ_local(x, h, g, h1, h2, h3, nlevel, verbose, opt, level)K = np.zeros((2*nlevel, 3*2**nlevel))# KD是二分段的峭度值，KQ是三分段的峭度值，组合到一起，形成最终的KK[0, :] = KD[0, :]for i in range(1, nlevel):K[2*i-1, :] = KD[i, :]K[2*i, :] = KQ[i-1, :]K[2*nlevel-1, :] = KD[nlevel, :]return K

可以看到该函数的核心功能就是K_wpQ_local函数实现的。

def K_wpQ_local(x, h, g, h1, h2, h3, nlevel, verbose, opt, level):"""递归计算二维矩阵的谱峭度:param x: 待分析的信号:param h: 二分段低通滤波器:param g: 二分段高通滤波器:param h1: 三分段第一段滤波器:param h2: 三分段第二段滤波器:param h3: 三分段第三段滤波器:param nlevel: 分解层数:param verbose: 是否输出调试信息:param opt: ['kurt1' | 'kurt2']* 'kurt1' = 基于4阶统计量的经典峭度* 'kurt2' = 基于2阶统计量的鲁棒峭度:param level: 当前分解层数:returns: K: 二分段的谱峭度值KQ: 三分段的谱峭度值"""# 计算低频系数a，以及高频系数da, d = DBFB(x, h, g)N = len(a)# 保持频带的频率顺序，使高通序列转换为低通序列d = d*np.power(-1., np.arange(1, N+1)) # indices pairs multipliés par -1# 计算低频部分峭度K1 = kurt(a[len(h)-1:], opt)# 计算高频部分峭度K2 = kurt(d[len(g)-1:], opt)if level > 2:# 计算低频段三分后的系数a1, a2, a3 = TBFB(a, h1, h2, h3)# 计算高频段三分后的系数d1, d2, d3 = TBFB(d, h1, h2, h3)# 计算各频段的峭度Ka1 = kurt(a1[len(h)-1:], opt)Ka2 = kurt(a2[len(h)-1:], opt)Ka3 = kurt(a3[len(h)-1:], opt)Kd1 = kurt(d1[len(h)-1:], opt)Kd2 = kurt(d2[len(h)-1:], opt)Kd3 = kurt(d3[len(h)-1:], opt)else:Ka1 = 0Ka2 = 0Ka3 = 0Kd1 = 0Kd2 = 0Kd3 = 0if level == 1:# 分解到最小单元时，K保持和KQ一样的长度，直接返回这两个值# K = [K1, K1, K1, K2, K2, K2]K = np.array([K1*np.ones(3), K2*np.ones(3)]).flatten()KQ = np.array([Ka1, Ka2, Ka3, Kd1, Kd2, Kd3])if level > 1:# 递归分解# level为2的时候，Ka和KaQ的长度均为6，Ka是两个频段的峭度，KaQ是将Ka的两个频段分解到六个频段(三分段滤波器)的峭度Ka, KaQ = K_wpQ_local(a, h, g, h1, h2, h3, nlevel, verbose, opt,level-1)Kd, KdQ = K_wpQ_local(d, h, g, h1, h2, h3, nlevel, verbose, opt,level-1)# K1表示Ka整个频带的峭度，K2表示Kd整个频带的峭度，每向上递归一层，K的长度乘2K1 = K1*np.ones(np.max(Ka.shape))K2 = K2*np.ones(np.max(Kd.shape))K12 = np.append(K1, K2)Kad = np.hstack((Ka, Kd))K = np.vstack((K12, Kad))# KQ的长度和K保持一致Long = int(2./6*np.max(KaQ.shape))Ka1 = Ka1*np.ones(Long)Ka2 = Ka2*np.ones(Long)Ka3 = Ka3*np.ones(Long)Kd1 = Kd1*np.ones(Long)Kd2 = Kd2*np.ones(Long)Kd3 = Kd3*np.ones(Long)# tmp是将低频部分a通过二分段滤波器、三分段滤波器后的6个频段值与高频部分同样得到的6个频段峭度值组合到一起tmp = np.hstack((KaQ, KdQ))# KQ是当前信号分解到低频部分a通过三分段滤波器后的3个频段峭度值和高频部分同样得到的3个频段峭度值组合到一起KQ = np.concatenate((Ka1, Ka2, Ka3, Kd1, Kd2, Kd3))KQ = np.vstack((KQ, tmp))if level == nlevel:K1 = kurt(x, opt)K = np.vstack((K1*np.ones(np.max(K.shape)), K))# 只计算了低频部分的峭度a1, a2, a3 = TBFB(x, h1, h2, h3)Ka1 = kurt(a1[len(h)-1:], opt)Ka2 = kurt(a2[len(h)-1:], opt)Ka3 = kurt(a3[len(h)-1:], opt)Long = int(1./3*np.max(KQ.shape))Ka1 = Ka1*np.ones(Long)Ka2 = Ka2*np.ones(Long)Ka3 = Ka3*np.ones(Long)tmp = np.array(KQ[0:-2])KQ = np.concatenate((Ka1, Ka2, Ka3))KQ = np.vstack((KQ, tmp))return K, KQ

DBFB函数实际上就是信号通过低通、高通滤波器后降采样，得到低频系数和高频系数。

def DBFB(x, h, g):"""双频带滤波，计算低频近似系数a贺高频细节系数d:param x: 待分析的信号:param h: 二分段低通滤波器:param g: 二分段高通滤波器:returns: a, d"""# lowpass filtera = si.lfilter(h, 1, x)a = a[1::2]a = a.ravel()# highpass filterd = si.lfilter(g, 1, x)d = d[1::2]d = d.ravel()return (a, d)

TBFB函数实际上是计算三段频段的系数，同样是滤波后降采样，只不过这里是3倍降采样。

def TBFB(x, h1, h2, h3):"""三频带滤波器:param x: 待分析的信号:param h1: 三分段第一段滤波器:param h2: 三分段第二段滤波器:param h3: 三分段第三段滤波器:returns: a1, a2, a3"""# lowpass filtera1 = si.lfilter(h1, 1, x)a1 = a1[2::3]a1 = a1.ravel()# passband filtera2 = si.lfilter(h2, 1, x)a2 = a2[2::3]a2 = a2.ravel()# highpass filtera3 = si.lfilter(h3, 1, x)a3 = a3[2::3]a3 = a3.ravel()return (a1, a2, a3)

下面我们再通过图示的方式，更为清晰地展示K_wpQ函数运算的逻辑。
如果设定level为7，那么当递归到level为1的时候，我们可以得到如下图所示的K和KQ。

当level为2的时候，可以看出K实际上就是本次待分解信号的低频峭度K1和高频峭度K2，在长度扩充后，再与低频部分a分解得到的Ka、高频部分d分解得到的Kd组合而得。此时level2层的K长度为12，即level1层K长度的一倍。

而KQ则是本次待分解信号的Ka1、Ka2、Ka3、Kd1、Kd2、Kd3，在长度扩充后，再与低频部分a分解得到的KaQ、高频部分d分解得到的KdQ组合而得。此时level2层的KQ长度为12，与K一致。

实际上后续的level都是按如此规律生成：K是本level层信号的K1、K2与level-1层的Ka、Kd组合而成，KQ是本level层信号的Ka1、Ka2、Ka3、Kd1、Kd2、Kd3、低频系数a在level-1层的KaQ、高频系数d在level-1层的KdQ组合而成。
计算得到K和KQ后，我们就可以按照自上而下一层用K系数填充，一层用KQ系数填充的方式得到最终的快速谱峭度图：

完整的代码可查看：https://download.csdn.net/download/qq1234534215/90471995

总结

本篇文章对快速谱峭度算法做了一个较为详细的介绍，实际上就是通过多组滤波器将信号分解到不同的频段，再计算每个频段的峭度，从而找到值最大的频段，进一步地我们也可以选取其他指标作为频段的选择，以便更好地找到振动信号中的故障信息。