大家好，欢迎来到IT知识分享网。

倒频谱分析也称为二次频谱分析，主要通过对信号频谱的对数进行傅里叶逆变换而实现，具有分析同族谐频、异族谐频和多成分边频的能力。工程上常用到的是功率倒频谱和幅值倒频谱，分别基于功率谱和幅值谱获得。倒频谱分析对信号频谱中较低幅值的分量给以较高的加权，有利于检测出弱周期信号，且能精确测量各分量的频率间隔，方便突出故障边频信息。比如在机械故障诊断中，信号频谱图上会出现难以识别的边频成分，采用倒频谱可以分解和识别这些边频成分包含的故障信息，并诊断产生故障的原因和位置。另外，倒频谱相比于频谱的主要优点在于对传感器安装位置、信号传输路径以及调制相位关系不敏感，保证了分析结果的稳定性。因此倒频谱分析在齿轮、轴承等故障诊断方面的应用较为广泛。但对于信号中故障特征较微弱时，倒频谱分析不能获得满意的效果，需要采用特定的手段对故障特征进行增强处理。

import numpy as np from scipy.fftpack import fft, rfft, irfft, fftfreq, rfftfreq import scipy.signal import matplotlib.pyplot as plt %matplotlib inline

# Create RC series dt = 0.001 # s, equiv. 1 ms T = 1.0 # s thick = 0.040 # s, equiv. 40 ms rc = 0.3 t = np.linspace(0, T, T/dt) rcs = np.zeros_like(t) rcs[1000*(T-thick)/2] = rc rcs[-1000*(T-thick)/2] = -rc ax1 = plt.subplot(111) ax1.plot(t, rcs) ax1.set_ylim(-0.4, 0.4) plt.show()

RCS = fft(rcs) freq = fftfreq(rcs.size, d=dt) keep = freq>=0 # only positive frequencies RCS = RCS[keep] freq = freq[keep]

import matplotlib.gridspec as gridspec gs = gridspec.GridSpec(1, 2,width_ratios=[1,2]) plt.figure(figsize=(15,5)) ax1 = plt.subplot(gs[0]) ax1.plot(t, rcs) ax1.set_ylim(-0.4, 0.4) ax2 = plt.subplot(gs[1]) ax2.plot(freq, np.abs(RCS)) ax2.set_xlim(0,200) plt.savefig("/home/matt/images/cepstrum.png") plt.show()

Basic signal and plot

def ricker(f, duration=0.512, dt=0.001): t = np.linspace(-duration/2, (duration-dt)/2, duration/dt) y = (1.0 - 2.0*(np.pi**2)*(f**2)*(t**2)) * np.exp(-(np.pi**2)*(f**2)*(t**2)) return t, y

dt = 0.001 # s, equiv. 1 ms T = 0.256 # s thick = 0.023 # s, equiv. 23 ms rc = 0.3 t = np.linspace(0, T, T/dt) rcs = np.zeros_like(t) rcs[1000*(T-thick)/2] = rc rcs[-1000*(T-thick)/2] = -rc tw, w = ricker(f=43., length=0.064, dt=dt) s = np.convolve(rcs, w, mode='same') ax1 = plt.subplot(111) ax1.plot(t, s) ax1.set_ylim(-0.4, 0.4) plt.show()

plt.plot(tw,w)

S = fft(s) freq = fftfreq(s.size, d=dt) keep = freq>=0 # only positive frequencies Sk = S[keep] freqk = freq[keep]

plt.figure(figsize=(15,5)) ax1 = plt.subplot(111) ax1.plot(freqk, np.abs(Sk)) ax1.set_xlim(0, 250) plt.show()

freqk = freqk[:] Sk = np.abs(Sk)

Q = fft(Sk) quef = fftfreq(Sk.size, d=freqk[1]-freqk[0]) keep = quef>=0 # only positive frequencies Qk = Q[keep] quefk = quef[keep]

quefk

array([ 49.0 +0.j , 38.-26.j, 14.-39.j , -7.-34.j, -19.-18.j, -18. -3.j, -9. +4.j, -1. +4.j, 3. -1.j, 1. -7.j, -5.-10.j, -12. -7.j, -16. +1.j, -12. +9.j , -4.+14.0j])

plt.figure(figsize=(15,5)) ax1 = plt.subplot(111) ax1.plot(quefk, np.abs(Qk)) plt.show()

import agilegeo.wavelet as ag

dt = 0.001 # s, equiv. 1 ms T = 0.256 # s thick = 0.021 # s, equiv. 21 ms rc = 0.3 # RC SERIES t = np.linspace(0, T, T/dt) rcs = np.zeros_like(t) rcs[1000*(T-3*thick)/2] = rc rcs[-1000*(T-3*thick)/2] = -rc rcs[1000*(T-2*thick)/2] = -rc rcs[-1000*(T-2*thick)/2] = rc rcs[1000*(T-thick)/2] = rc rcs[-1000*(T-thick)/2] = -rc # WAVELET duration = 0.256 f = f=[10,16,64,80] #tw, w = ricker(f=32., duration=0.128, dt=dt) w = ag.ormsby(f=f, duration=duration, dt=dt) tw = np.linspace(-duration/2, (duration-dt)/2, duration/dt) # CONVOLVE s = np.convolve(rcs, w, mode='same') # Plot it plt.figure(figsize=(20,3)) ax1 = plt.subplot(131) ax1.plot(tw, w) ax1.set_ylim(-0.5, 1.1) ax1 = plt.subplot(132) ax1.plot(t, rcs) ax1.set_xlim(0, 0.256) ax1.set_ylim(-0.4, 0.4) ax1 = plt.subplot(133) ax1.plot(t, s) ax1.set_xlim(0, 0.256) ax1.set_ylim(-0.4, 0.4) plt.show()

W = fft(w) freqw = fftfreq(w.size, d=dt) RCS = fft(rcs) S = fft(s) freq = fftfreq(s.size, d=dt) keep = freq>=0 keepw = freqw>=0 Wk = W[keepw] freqwk = freqw[keepw] RCSk = RCS[keep] Sk = S[keep] freqk = freq[keep] plt.figure(figsize=(20,3)) ax1 = plt.subplot(131) ax1.plot(freqwk, np.abs(Wk)) ax1.set_xlim(0,200) ax2 = plt.subplot(132) ax2.plot(freqk, np.abs(RCSk)) ax2.set_xlim(0,200) ax3 = plt.subplot(133) ax3.plot(freqk, np.abs(Sk)) ax3.set_xlim(0, 200) plt.show()

Wk = np.abs(Wk) RCSk = np.abs(RCSk) Sk = np.abs(Sk) Qw = fft(Wk) Qr = fft(RCSk) Qs = fft(Sk) quefw = fftfreq(Wk.size, d=freqwk[1]-freqwk[0]) quefr = fftfreq(RCSk.size, d=freqk[1]-freqk[0]) quefs = fftfreq(Sk.size, d=freqk[1]-freqk[0]) keepw = quefw>=0 # only positive frequencies Qwk = Qw[keepw] quefwk = quefw[keepw] keep = quefr>=0 # only positive frequencies Qrk = Qr[keep] Qsk = Qs[keep] quefk = quef[keep] plt.figure(figsize=(15,3)) ax1 = plt.subplot(131) ax1.plot(quefwk, np.abs(Qwk)) ax2 = plt.subplot(132) ax2.plot(quefk, np.abs(Qrk)) ax3 = plt.subplot(133) ax3.plot(quefk, np.abs(Qsk)) plt.show()

知乎学术咨询： https://www.zhihu.com/consult/people/?isMe=1 工学博士，担任《Mechanical System and Signal Processing》《中国电机工程学报》《控制与决策》等期刊审稿专家，擅长领域：现代信号处理，机器学习，深度学习，数字孪生，时间序列分析，设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。