时频分析法——短时傅里叶变换(STFT)

时频分析法是对时序信号的一种时间-频率联合分析手段，也是将一维时间序列信号转化为二维图像的一种方法。其中，短时傅里叶变换、连续小波变换、希尔伯特-黄变换是时频分析法中最常用的三种方法。本文将以心电信号(ECGs)为例子来着重介绍STFT，相关的代码和数据也将在下面提供。

所使用的数据来自于公开的心电信号数据库，获取网址如下：PTB Diagnostic ECG Database v1.0.0 (physionet.org)

本文所有插图均由其中同一份经过预处理的数据得到，有想复现的小伙伴可以下载：文章顶部的"processed_ecg segment_csdn"

STFT的原理解释：

为了更直观的了解STFT，我们先画出时序心电信号图，以及它的由傅里叶变换(FFT)得到的频谱图。

时序信号图及绘制代码：

import matplotlib.pyplot as plt
import json
import matplotlib

# 全局字体和大小
matplotlib.rcParams['font.family'] = 'Times New Roman'
matplotlib.rcParams['font.size'] = 26

#加载ECG数据
with open(r"F:\ecg segment_csdn.json", "r") as file:   #将地址"F:\ecg segment_csdn.json"换成你的本地储存地址
    ecg_data = json.load(file)
signal_data = ecg_data['data']
fs = ecg_data['fs']  #加载采样频率 hz=1000

# 获取v1导联的数据(有12个导联数据，其中一个叫V1)
signal = signal_data['V1']

# 创建图表并绘制V1导联
plt.figure(figsize=(15, 5))
plt.plot(signal)
plt.title('V1 Lead ECG Signal')
plt.xlabel('Time (samples)')
plt.ylabel('Amplitude')
plt.show()

频谱图及绘制代码： 

import numpy as np
import matplotlib.pyplot as plt
import json
import matplotlib

matplotlib.rcParams['font.family'] = 'Times New Roman'
matplotlib.rcParams['font.size'] = 26

#加载ECG数据
with open(r"F:\ecg segment_csdn.json", "r") as file:
    ecg_data = json.load(file)
signal_data = ecg_data['data']
fs = ecg_data['fs']  # 加载采样频率

#获取V1导联的数据
signal = signal_data['V1']

#傅里叶变换（FFT）
fft_result = np.fft.fft(signal)
frequencies = np.fft.fftfreq(len(signal), 1/fs)

#获取频率
half_n = len(fft_result) // 2
fft_magnitude = np.abs(fft_result[:half_n]) * 2 / len(fft_result)
frequencies = frequencies[:half_n]

# 创建图表并绘制频谱图
plt.figure(figsize=(15, 5))
plt.plot(frequencies, fft_magnitude)
plt.title('Frequency Spectrum of V1 Lead ECG Signal')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.xlim(0, fs/2)
plt.show()

 通过上面的操作我们对整个信号进行FFT得到它的频谱图，假设信号为x(t)，则FFT的计算公式：

既然FFT能够得到信号的频率，那么我们使用类似于积分的办法对每一小段信号进行FFT是不是就可以视为得到该信号每一个时刻的频率分布呢？没错，你已经get到STFT的核心思想了。PS：重要的事再说一遍“STFT就是对一小段一小段的信号进行FFT，然后再全部连起来”。                        

于是为了将整条信号分成多个小段，引入了一个窗函数，使得我们可以在对窗口内的信号段进行FFT，并且通过移动实现分窗，如下图所示：                                                                                                  

 在这里我们使用流行的汉明窗口，M是窗口长度，它的公式如下：  

假设窗口中心在，那么将要在该窗口下进行FFT的信号则可表示为y(t)； 

最终我们对窗口信号做FFT得到时频图S(t, f)：

 

详细的实现STFT的代码如下，代码解释可以看注释 ：                                                                                                          

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import spectrogram
import json
import matplotlib

matplotlib.rcParams['font.family'] = 'Times New Roman'
matplotlib.rcParams['font.size'] = 26

#加载ECG数据
with open(r"F:\ecg segment_csdn.json", "r") as file:
    ecg_data = json.load(file)
fs = ecg_data['fs']  #加载采样频率

#短时傅里叶函数实现     其中，neperseg是窗口长度，noverlap是窗口移动长度，这里设置了一半重叠
def compute_stft(signal, fs):
    signal = np.array(signal)  
    frequencies, times, Sxx = spectrogram(signal, fs=fs, window='hamming', nperseg=500, noverlap=250)  
    return frequencies, Sxx 

lead = 'V1'
signal = np.array(ecg_data['data'][lead])
frequencies, Sxx = compute_stft(signal, fs)

#限制频率范围至0到80Hz进行绘图
idx_max = np.max(np.where(frequencies <= 80)[0])
frequencies_clipped = frequencies[:idx_max+1]
Sxx_clipped = Sxx[:idx_max+1, :]

#绘制STFT频谱图
plt.figure(figsize=(15, 5))
plt.imshow(Sxx_clipped, aspect='auto', cmap='jet', origin='lower', extent=[0, 4, 0, 80])
plt.title(f'Spectrogram of Lead {lead}')
plt.xlabel('Time (s)')
plt.ylabel('Frequency (Hz)')
plt.colorbar(label='Magnitude (dB)')
plt.show()

最终运行得到的时频图如下：