import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import stft
#STFT 短时傅里叶变换
fs = 1000
N = 1000
k = np.arange(2000)
frq = k*fs/N
frq1 = 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))/N
data_f1 = data_f[range(int(N / 2))]
plt.plot(frq1, data_f1)
data_ff = abs(np.fft.fft(y2)) / N
data_f2 = data_ff[range(int(N / 2))]
plt.plot(frq1, data_f2, 'red')
data_fff = abs(np.fft.fft(y3)) / N
data_f3 = data_fff[range(int(N / 2))]
plt.plot(frq1, data_f3, 'k')
data_ffff = abs(np.fft.fft(y4)) / N
data_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
# STFT 整个是与过程分解成无数个等长的小过程,每个小过程近似平稳,再进行短时傅里叶变换
# scipy.signal.stft(x, fs, window, nperseg, noverlap, nfft, detrend, return_oneside, boundary, padded, axis)
'''
常用参数:
x : 传入STFT变换的时域信号
fs : 时域信号的采样频率
window: 时域信号分割需要的窗函数
noverlap:窗函数的重叠数,默认为窗长的一半
nfft: FFT的长度
return_oneside : True返回复数实部,None返回复数,默认为False
注意:窗长(即为nperseg)设置的太短,窗内截取的信号太短,会导致频率分析不够精确,频率分辨率差;窗长设置太长,时域不够精确,时间分辨率低。
'''
f, 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.colorbar()
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.tight_layout()
plt.show()
def sigh_stft(x):
f,t,z = stft(x,fs=20,window='hann',nperseg=128)
plt.pcolormesh(t, f, np.abs(z))
plt.colorbar()
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.tight_layout()
plt.show()
if __name__ == '__main__':
path = ""
contents = np.load(path, allow_pickle=True)
content = contents.tolist()
sigh_stft(heart_data)
信号STFT处理
于 2022-08-05 16:33:43 首次发布