python实现时间序列信号的频谱、倒频谱以及功率谱

以振动信号为例

---

认识傅里叶变换

这里我就不多说了 百度谷歌一大堆说的明明白白。

一、频谱

频谱，在复杂的振动中又称振动谱 。反映振动现象最基本的物理量就是频率，例如简单周期振动只有一个频率。
将信号通过快速傅里叶变换会得到幅频谱和相频谱两种图，不过常用的的是幅频谱。若对加速度进行变换，得到的振动幅频谱中，横坐标表示分振动的圆频率***f（Hz）***，纵坐标则表示分振动的加速度振幅***A(m/s^2）***

1.引入库

两种方式：用numpy库 或 scipy库（示例）：

#基于numpy的方法
import numpy.fft as fft
#基于scipy的方法
from scipy.fftpack import fft

2.频谱函数封装

def get_fft_values(y_values, N, f_s):
    f_values = np.linspace(0.0, f_s/2.0, N//2)
    fft_values_ = fft.fft(y_values)
    fft_values = 2.0/N * np.abs(fft_values_[0:N//2])
    return f_values, fft_values

其中y_values是时间序列下的振动信号值，N为采样点数，f_s是采样频率，返回f_values希望的频率区间, fft_values真实幅值，代码运行结果。

二、功率谱

一开始我常常把功率谱和能量谱搞混。
能量谱是原信号傅立叶变换的平方。
功率谱是原信号傅立叶变换的平方并除以采样点数N，称功率谱密度函数，它定义为单位频带内的信号功率。它表示了信号功率随着频率的变化情况，即信号功率在频域的分布状况。
此外维纳-辛钦定理指出：一个信号的功率谱密度就是该信号自相关函数的傅里叶变换。

功率谱谱函数封装

代码如下：

def get_fft_power_spectrum(y_values, N, f_s, fs_n):
    f_values = np.linspace(0.0, f_s/fs_n, N//fs_n)
    fft_values_ = np.abs(fft(y_values))
    fft_values = 2.0/N * (fft_values_[0:N//fs_n])    #频率真实幅值分布
    # power spectrum 直接周期法
    ps_values = fft_values**2 / N
    # power spectrum using correlate
    cor_x = np.correlate(y_values, y_values, 'same')    #自相关
    cor_X = fft(cor_x, Sampling_points)
    ps_cor = np.abs(cor_X)
    ps_cor_values = 10*np.log10(ps_cor[0:N//fs_n] / np.max(ps_cor))
    
    return f_values, fft_values, ps_values, ps_cor_values

其中返回ps_values是直接周期法功率, ps_cor_values是自相关下的对数功率。

从图中可以看出直接周期法功率的方差大，频率分辨率差，会将一些低幅值的频率给淹没，可以采用对数的方式解决，做对数处理的目的是使低振幅成分得以拉高，便于观察噪声中的周期信号。如自相关下的对数功率谱，虽然有明显的变化，但波形较为杂乱。
直接周期法和自相关法都是经典的功率谱估计方法，为了克服以上的缺点出现了现代参数法功率谱估计，比较有效且实用的是AR模型法，Burg谱估计法

三、倒频谱

倒频谱是指信号的对数功率谱的逆它具有时间因次,为了与通常频率的频谱相区别,有时称它为时谱。

倒频谱与对数功率谱是一对傅里叶变换，正如上述自相关函数与功率谱是一对傅里叶变换。其差异在于倒频谱是功率谱在对数坐标下的傅里叶逆变换,而自相关函数是功率谱在线性坐标下的傅里叶逆变换

该方法方便提取、分析原频谱图上肉眼难以识别的周期性信号，能将原来频谱图上成族的边频带谱线简化为单根谱线，受传感器的测点位置及传输途径的影响小

倒频谱谱函数封装

def get_fft_cepstrum(y_values, N, f_s, fs_n):
    f_values = np.linspace(0.0, f_s/fs_n, N//fs_n)
    fft_values_ = np.abs(fft(y_values))
    fft_values = 2.0/N * np.abs(fft_values_[0:N//fs_n])    #频率真实幅值分布
    
    spectrum = fft(y, N)
    ceps_values = np.fft.ifft(np.log(np.abs(spectrum))).real #信号→傅里叶→对数→傅里叶逆变换
    
    return f_values, fft_values, ceps_values