信号STFT处理

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)