信号处理 | 短时傅里叶变换实战

短时傅里叶变换(STFT)原理

短时傅里叶变换（Short-Time Fourier Transform, STFT）是一种分析时变信号频率特性的方法。它通过将长时间的信号分割成较短的时间片段，然后对每个时间片段进行傅里叶变换，从而克服了传统傅里叶变换无法同时提供时间和频率信息的限制。
原理

1. 分割信号：STFT首先将连续的信号分割成较短的时间片段。这通常通过乘以一个滑动窗口函数来实现，窗口函数在特定的时间区间内非零，并随着时间滑动。

2. 窗口函数：窗口函数的选择对STFT的结果有重要影响。常用的窗口函数包括矩形窗、汉宁窗、汉明窗等。窗口函数的宽度（或称为窗口长度）决定了时间分辨率和频率分辨率的平衡：窗口越宽，频率分辨率越高，时间分辨率越低；窗口越窄，时间分辨率越高，频率分辨率越低。

3. 傅里叶变换：对每个时间片段应用傅里叶变换，计算该时间片段内信号的频率成分。这样，每个时间片段都对应一个频谱。

4. 时间-频率表示：将所有时间片段的傅里叶变换结果组合起来，就可以得到信号的时间-频率表示，即STFT的结果。这个结果通常表示为一个二维数组，其中一个维度表示时间，另一个维度表示频率。
   数学表达式

STFT的数学表达式为：
S T F T { x ( t ) } ( τ , ω ) = ∫ − ∞ + ∞ x ( t ) ⋅ w ( t − τ ) ⋅ e − j ω t d t STFT\{x(t)\}(τ, ω) = \int_{-\infty}^{+\infty} x(t) \cdot w(t-τ) \cdot e^{-jωt} dt STFT{ x(t)}(τ,ω)=∫−∞+∞​x(t)⋅w(t−τ)⋅e−jωtdt
其中，$x(t)$是原始信号，$w(t-\tau )$是窗口函数，$\tau$是时间变量，表示当前窗口的中心位置，$\omega$是频率变量。
应用

STFT广泛应用于信号处理领域，如语音分析、音乐处理、地震数据分析等，它能够提供信号随时间变化的频率信息，对于非平稳信号分析尤为重要。

示例代码

生成模拟信号

import numpy as np
import scipy.io as scio

from matplotlib import pyplot as plt
from matplotlib import rcParams

绘制原始时域信号

def plt_time_domain(arr, fs=1600, ylabel='Amp(mg)', title='原始数据时域图', img_save_path=None, x_vline=None, y_hline=None):
    """
    :fun: 绘制时域图模板
    :param arr: 输入一维数组数据
    :param fs: 采样频率
    :param ylabel: y轴标签
    :param title: 图标题
    :return: None
    """
    import matplotlib.pyplot as plt
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 显示中文
    plt.rcParams['axes.unicode_minus'] = False  # 显示负号
    font = {
   'family': 'Times New Roman', 'size': '20', 'color': '0.5', 'weight': 'bold'}
    
    plt.figure(figsize=(12,4))
    length = len(arr)
    t = np.linspace(0, length/fs, length)
    plt.plot(t, arr, c='g')
    plt.xlabel('t(s)')
    plt.ylabel(ylabel)
    plt.title(title)
    if x_vline:
        plt.vlines(x=x_vline, ymin=np.min(arr), ymax=np.max(arr), linestyle='--', colors='r')
    if y_hline:
        plt.hlines(y=0.2, xmin=np.min(t), xmax=np.max(t), linestyle=':', colors='y')
    #===保存图片====#
    if img_save_path:
        plt.savefig(img_save_path, dpi=500, bbox_inches = 'tight')
    plt.show()

fs = 100  # 采样频率
f = 200    # 模拟正弦信号频率
time = 5  # 采样时长
t = np.linspace(0, time, time*fs)
data = 1*np.sin(2*np.pi*f*t) + np.random.normal(0, 0.1, time*fs)
plt_time_domain(data, fs=fs)

绘制STFT图

import scipy.signal as signal
import numpy as np
import matplotlib.pyplot as plt
 
f, t, nd = signal.stft(data, fs=fs, window='hann', nperseg=128, noverlap=64,nfft=None,
                       detrend=False, return_onesided=True, boundary='odd', padded=False, axis=-1)
#  fs:时间序列的采样频率,  nperseg:每个段的长度，默认为256(2^n)   noverlap:段之间重叠的点数。如果没有则noverlap=nperseg/2
 
#window ： 字符串或元组或数组，可选需要使用的窗。
# #如果window是一个字符串或元组，则传递给它window是数组类型，直接以其为窗，其长度必须是nperseg。
# 常用的窗函数有boxcar，triang，hamming， hann等，默认为Hann窗。
 
#nfft ： int，可选。如果需要零填充FFT，则为使用FFT的长度。如果为 None，则FFT长度为nperseg。默认为无
 
# detrend ： str或function或False，可选
# 指定如何去除每个段的趋势。如果类型参数传递给False，则不进行去除趋势。默认为False。
 
# return_onesided ： bool，可选
# 如果为True，则返回实际数据的单侧频谱。如果 False返回双侧频谱。默认为 True。请注意，对于复杂数据，始终返回双侧频谱。
 
# boundary ： str或None，可选
# 指定输入信号是否在两端扩展，以及如何生成新值，以使第一个窗口段在第一个输入点上居中。
# 这具有当所采用的窗函数从零开始时能够重建第一输入点的益处。
# 有效选项是['even', 'odd', 'constant', 'zeros', None].
# 默认为‘zeros’,对于补零操作[1, 2, 3, 4]变成[0, 1, 2, 3, 4, 0] 当nperseg=3.
 
# padded： bool，可选
# 指定输入信号在末尾是否填充零以使信号精确地拟合为整数个窗口段，以便所有信号都包含在输出中。默认为True。
# 填充发生在边界扩展之后，如果边界不是None，则填充为True，默认情况下也是如此。
 
# axis ： int，可选
# 计算STFT的轴; 默认值超过最后一个轴(即axis=-1)。

plt.figure(figsize=(12,4))
plt.pcolormesh(t, f, np.abs(nd), vmin = np.min(np.abs(nd)), vmax = np.max(np.abs(nd)))
plt.title('STFT')
plt.ylabel('frequency')
plt.xlabel('time')
plt.show()

可见在整个时间上存在10Hz的信号。

保存图片

plt.save_fig(file_path, bbox_inches='tight')