时间序列信号处理（四）——傅里叶变换和短时傅里叶变换python实现

一、傅里叶变换（FFT）

一维时间序列信号为典型的离散信号，需要使用离散傅里叶变换，FFT变换就是将时域信号转到频域阶段分析，分析更深入，了解时域信号的频域特性。注意傅里叶变换只适用于平稳信号。

F(u)为离散傅里叶变换值，f(n)为逆变换值。

python实现案例：

import matplotlib.pyplot as pltimport numpy as npfrom scipy.fftpack import fft, ifftFs = 2000Ts = 1.0/FsN = 2000t = np.linspace(0, N*2000, N)k = np.arange(2000)frq = k*Fs/Nfrq1 = frq[range(int(N/2))]data = 2*np.sin(4*np.pi*50*t) + 4*np.sin(4*np.pi*120*t)plt.plot(data, 'grey')plt.xlabel('time')plt.ylabel('amplitude')plt.title("data")plt.show()# 原始图像单边谱data_f = abs(np.fft.fft(data))/Ndata_f1 = data_f[range(int(N/2))]plt.plot(frq1, data_f1, 'red')plt.xlabel('pinlv(hz)')plt.ylabel('amplitude')plt.title("pinputu")plt.show()

可以明显看出来这个信号的主要频率为100和250hz的样子，对上述时域信号又有了不同了解了。

去噪示例

# 去噪data1 = 3*np.cos(3*np.pi*2*t)+4*np.sin(3*np.pi*4*t)data1 = data1 + 0.3 * np.random.normal(size=data1.size)plt.plot(data1, 'grey')plt.xlabel('time ')plt.ylabel('amplitude')plt.title("data")plt.show()y = fft(data1)threadhold = 50y[threadhold:(N-threadhold)] = 0 # 滤波器data_c = ifft(y)plt.plot(data_c)plt.title('signal after filtering')plt.xlabel('time')plt.ylabel('amplitude')plt.show()

可以看出去噪效果还是可以，这样就有利于信号的重构，利用FFT和Ifft变换实现去噪。

二、短时傅里叶变换STFT

短时傅里叶变换就是在傅里叶变换的基础上进行了加窗处理，这也就提升了对非平稳信号的处理能力，强化了特征提取能力。

式中， x(m)是输入信号；w(m)是窗函数；它在时间上翻转并且有 n个样本的偏移量；X(n,w)是定义在样本（时间）和频率上的二维函数。

注意：FFT也不能辨别信号不同。

案例：

import matplotlib.pyplot as pltimport numpy as npfrom scipy.fftpack import fft, ifftfrom scipy.signal import stftfs = 1000N = 1000k = np.arange(2000)frq = k*fs/Nfrq1 = frq[range(int(N/2))]t = np.linspace(0, 1-1/fs, fs)x = np.linspace(0, 4000, 4000)y1 = 100 * np.cos(2 * np.pi * 100 * t)y2 = 200 * np.cos(2 * np.pi * 200 * t)y3 = 300 * np.cos(2 * np.pi * 300 * t)y4 = 400 * np.cos(2 * np.pi * 400 * t)y = np.append(y1, y2)yy = np.append(y3, y4)yyy = np.append(y, yy)# 画原始图plt.plot(x, yyy)plt.xlabel('time')plt.ylabel('amplitude')plt.title("data")plt.show()data_f = abs(np.fft.fft(y1))/Ndata_f1 = data_f[range(int(N/2))]plt.plot(frq1, data_f1)data_ff = abs(np.fft.fft(y2))/Ndata_f2 = data_ff[range(int(N/2))]plt.plot(frq1, data_f2, 'red')data_fff = abs(np.fft.fft(y3))/Ndata_f3 = data_fff[range(int(N/2))]plt.plot(frq1, data_f3, 'k')data_ffff = abs(np.fft.fft(y4))/Ndata_f4 = data_ffff[range(int(N/2))]plt.plot(frq1, data_f4, 'b')plt.xlabel('pinlv(hz)')plt.ylabel('amplitude')plt.title("pinputu")plt.show()window = 'hann'# frame长度n = 256# STFTf, t, Z = stft(yyy, fs=fs, window=window, nperseg=n)# 求幅值Z = np.abs(Z)# 如下图所示plt.pcolormesh(t, f, Z, vmin=0, vmax=Z.mean()*10)plt.show()# 反着画k1 = np.append(y4, y3)k2 = np.append(y2, y1)k = np.append(k1, k2)plt.plot(x, k)plt.xlabel('time')plt.ylabel('amplitude')plt.title("data")plt.show()data_f = abs(np.fft.fft(y1))/Ndata_f1 = data_f[range(int(N/2))]plt.plot(frq1, data_f1)data_ff = abs(np.fft.fft(y2))/Ndata_f2 = data_ff[range(int(N/2))]plt.plot(frq1, data_f2, 'red')data_fff = abs(np.fft.fft(y3))/Ndata_f3 = data_fff[range(int(N/2))]plt.plot(frq1, data_f3, 'k')data_ffff = abs(np.fft.fft(y4))/Ndata_f4 = data_ffff[range(int(N/2))]plt.plot(frq1, data_f4, 'b')plt.xlabel('pinlv(hz)')plt.ylabel('amplitude')plt.title("pinputu")plt.show()window = 'hann'# frame长度n = 256# STFTf, t, Z = stft(k, fs=fs, window=window, nperseg=n)# 求幅值Z = np.abs(Z)# 如下图所示plt.pcolormesh(t, f, Z, vmin=0, vmax=Z.mean()*10)plt.show()

正着画图，其原始图、频谱图及时频图如下：

反着画图，其原始图、频谱图及时频图如下：

由上述对比可以知道，FFT只能看出两个信号数据的频率分布，并不能分别两个信号，而时频图可以，说明stft对信号又有更深的理解。