目錄
本教程為腦機學習者Rose原創(轉載請聯系作者授權)發表於公眾號：腦機接口社區(微信號：Brain_Computer).QQ交流群：903290195
倒頻譜定義
倒頻譜可以分析復雜頻譜圖上的周期結構，分離和提取在密集調頻信號中的周期成分，對於具有同族諧頻、異族諧頻和多成分邊頻等復雜信號的分析非常有效。倒頻譜變換是頻域信號的傅立葉積分變換的再變換。時域信號經過傅立葉積分變換可轉換為頻率函數或功率譜密度函數，如果頻譜圖上呈現出復雜的周期結構而難以分辨時，對功率譜密度取對數再進行一次傅立葉積分變換，可以使周期結構呈便於識別的譜線形式。第二次傅立葉變換的平方就是倒功率譜，即“對數功率譜的功率譜”。倒功率譜的開方即稱幅值倒頻譜，簡稱倒頻譜。
簡言之，倒頻譜分析技術是將時域振動信號的功率譜對數化，然后進行逆傅里葉變化后得到的。倒頻譜的水平軸為“倒頻率”的偽時間，垂直軸為對應倒頻率的幅值，其計算公式為：
倒頻譜python案例
實現如下：
from scipy.fftpack import fft, fftshift, ifft
from scipy.fftpack import fftfreq
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")
fs = 1000
#采樣點數
num_fft = 1024
"""
生成原始信號序列
在原始信號中加上噪聲
np.random.randn(t.size)
其中y1是主頻為5/10/20Hz的低頻信號+噪聲信號;
y2是主頻為50、100/200Hz的高頻信號+噪聲信號;
y是y1和y2的調制結果
"""
t = np.arange(0, 5, 1/fs)
y1 = 10*np.cos(2*np.pi*5*t) + 7*np.cos(2*np.pi*10*t) + 5*np.cos(2*np.pi*20*t) + np.random.randn(t.size)
y2 = 20*np.cos(2*np.pi*50*t) + 15*np.cos(2*np.pi*100*t) + 25*np.cos(2*np.pi*200*t) + np.random.randn(t.size)
y = y1*y2
plt.figure(figsize=(20, 12))
ax=plt.subplot(331)
ax.set_title('y1')
plt.plot(y1)
ax=plt.subplot(332)
ax.set_title('y2')
plt.plot(y2)
ax=plt.subplot(333)
ax.set_title('y=y1*y2')
plt.plot(y)
"""
對低頻信號y1進行 FFT(Fast Fourier Transformation)快速傅里葉變換
"""
Y1 = fft(y1, num_fft)
Y1 = np.abs(Y1)
ax=plt.subplot(334)
ax.set_title('y1 fft')
plt.plot(Y1[:num_fft//2])
"""
對高頻信號y2進行 FFT
"""
Y2 = fft(y2, num_fft)
Y2 = np.abs(Y2)
ax=plt.subplot(335)
ax.set_title('y2 fft')
plt.plot(Y2[:num_fft//2])
"""
對信號y進行 FFT
"""
Y = fft(y, num_fft)
Y = np.abs(Y)
ax=plt.subplot(336)
ax.set_title('y fft')
plt.plot(Y[:num_fft//2])
plt.tight_layout()
plt.show()
"""
倒頻譜的定義表述為：信號→功率譜→對數→傅里葉逆變換
"""
spectrum = np.fft.fft(y, n=num_fft)
ceps = np.fft.ifft(np.log(np.abs(spectrum))).real
plt.figure(figsize=(10, 5))
plt.plot(np.abs(ceps)[:num_fft//2])
plt.title('y->spectrum->log->ifft')
plt.ylim([0, 0.2])
plt.show()
代碼來源於網絡，本文對代碼進行注釋並整理
腦機學習者Rose筆記分享,QQ交流群：903290195
更多分享，請關注公眾號