前言
最近一直在研究振动信号的最优解调带宽寻找方法，而快速谱峭度算法是最常见的一种，因此本篇文章准备对该算法做个深入解析，也相当于自己的学习记录了。
原理解析
谱峭度是一种信号处理中的时频分析工具，用于量化信号在不同频率成分上的非高斯性或冲击特性。它通过在频域中计算每个频率段上的峭度（统计学中衡量数据分布陡峭程度的指标，高峭度表示数据具有更尖锐的峰值和更重的尾部，如冲击信号），反映信号中瞬态成分或异常事件的分布情况。
快速谱峭度算法则是Antoni提出的一种高效计算谱峭度的方法，核心思想是通过多分辨率滤波器组快速分解信号频带，并优化峭度值的计算流程，从而显著提升效率。
1、信号分解
通过多级滤波器组将原始信号分解为多个子频带，目前较为常用的是1/3-二叉树的结合。
2、频带峭度计算
对每个子频带信号 x i ( t ) x_i(t) xi(t)计算其经典峭度值 K i K_i Ki。
K i = m e a n ( ∣ x i ( t ) ∣ 4 ) m e a n ( ∣ x i ( t ) ∣ 2 ) 2 K_i=\frac{mean(|x_i(t)|^4)}{mean(|x_i(t)|^2)^2} Ki=mean(∣xi(t)∣2)2mean(∣xi(t)∣4)
快速谱峭度算法里还提供了一种基于2阶统计量的鲁棒峭度值。
K i = m e a n ( ∣ x i ( t ) ∣ 2 ) m e a n ( ∣ x i ( t ) ∣ ) 2 K_i=\frac{mean(|x_i(t)|^2)}{mean(|x_i(t)|)^2} Ki=mean(∣xi(t)∣)2mean(∣xi(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 vector
freq_w, the complex envelope of the signal c and the extreme
frequencies 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) - 7
if 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 envelope
Kwav = K_wpQ(x, h, g, h1, h2, h3, nlevel, verbose, 'kurt2')
else:
# variance of the envelope magnitude
Kwav = K_wpQ(x, h, g, h1, h2, h3, nlevel, verbose, 'kurt1')
# keep positive values only!
# 计算最后一层各频段的频率，最大频率是采样率的一半
Kwav[Kwav <= 0] = 0
Level_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 = 1
lev = l1
while 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 bandpass
f_lower = Fs*np.round((fc-Bw/2.)*10**3)/10**3
f_upper = Fs*np.round((fc+Bw/2.)*10**3)/10**3
return 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 = nlevel
x = 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是三分段的峭度值，组合到一起，形成最终的K
K[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，以及高频系数d
a, 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 = 0
Ka2 = 0
Ka3 = 0
Kd1 = 0
Kd2 = 0
Kd3 = 0
if 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的长度乘2
K1 = 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 filter
a = si.lfilter(h, 1, x)
a = a[1::2]
a = a.ravel()
# highpass filter
d = 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 filter
a1 = si.lfilter(h1, 1, x)
a1 = a1[2::3]
a1 = a1.ravel()
# passband filter
a2 = si.lfilter(h2, 1, x)
a2 = a2[2::3]
a2 = a2.ravel()
# highpass filter
a3 = 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
总结
本篇文章对快速谱峭度算法做了一个较为详细的介绍，实际上就是通过多组滤波器将信号分解到不同的频段，再计算每个频段的峭度，从而找到值最大的频段，进一步地我们也可以选取其他指标作为频段的选择，以便更好地找到振动信号中的故障信息。