python 时频分析

### 使用Python实现时频分析的方法

对于不平稳信号而言，传统的傅里叶变换仅能提供全局的频率信息而忽略了局部的时间特性。为了克服这一局限性，多种时频分析技术被开发出来，其中包括但不限于短时傅里叶变换（STFT），小波变换（WT），以及经验模态分解加希尔伯特黄变换（EMD/HHT）。这些方法允许观察者同时获取时间和频率两个维度的信息。

#### 短时傅里叶变换 (STFT)

通过引入滑动窗口机制，STFT能够在保持一定时间分辨率的同时计算出不同时间段内的频谱特征[^2]。下面是一个简单的例子展示如何利用`matplotlib`库中的`spectrogram()`函数绘制语谱图：

import numpy as np
from matplotlib import pyplot as plt
from scipy.signal import get_window, stft

# 创建测试信号
fs = 10e3  # 采样率
t = np.linspace(0, .5, int(fs/2))
signal = np.sin(2*np.pi*10*t) + np.sin(2*np.pi*500*t)

frequencies, times, Sxx = stft(signal, fs=fs, window='hann', nperseg=1024)
plt.pcolormesh(times, frequencies, np.abs(Sxx), shading='gouraud')
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.show()

此段代码首先定义了一个复合正弦波作为输入数据集；接着应用汉宁窗(`'hann'`)进行了STFT运算，并最终以伪彩色图像的形式展示了所得的结果——即所谓的“语谱图”。

#### 小波变换 (Wavelet Transform)

相比于固定宽度的时间窗口，小波基的选择赋予了WT更灵活的时间-尺度表示能力。这里给出一段基于PyWavelets包的小波连续变换(CWT)实例:

import pywt
import matplotlib.pyplot as plt

waveletname = 'morl'
data = ... # 加载实际的数据序列
coeffs, freqs = pywt.cwt(data, scales=np.arange(1, 128), wavelet=waveletname)
plt.imshow(abs(coeffs), extent=[-1, 1, 1, 31], cmap='PRGn', aspect='auto',
           vmax=abs(coeffs).max(), vmin=-abs(coeffs).max())
plt.show()

上述脚本选择了Morlet小波作为基础函数执行CWT操作，并采用热力图的方式呈现了转换后的系数矩阵。

#### 经验模式分解与希尔伯特黄变换 (EMD & HHT)

当面对复杂的多分量非线性和非周期过程时，可以考虑使用EMD先将原始信号逐步分离成若干个内在模式函数(IMF)，再借助Hilbert变换求取瞬时相位角和幅值从而构建完整的时频分布图。以下是结合`pyhht`库完成整个流程的一个简单案例：

from pyhht.emd import EMD
from pyhht.visualization import plot_imfs

decomposer = EMD(s=data)
imfs = decomposer.decompose()

plot_imfs(imfs=imfs, time_samples=t, original_signal=s,
          title="Intrinsic Mode Functions",
          figsize=(12, 9))

for i in range(len(imfs)):
    analytic_signal = hilbert(imfs[i])
    instantaneous_phase = np.unwrap(np.angle(analytic_signal))
    inst_freq = (np.diff(instantaneous_phase)/(2.0*np.pi)*fs)
    
    fig, ax = plt.subplots(figsize=(12, 6))
    ax.plot(t[:-1], inst_freq)
    ax.set_xlabel("Time")
    ax.set_ylabel("Instantaneous Frequency")
    ax.grid(True)
    plt.title(f"Instantaneous frequency of IMF {i}")
    plt.show()

这段程序先是运用EMD算法提取出了所有的IMFs组件，之后针对每一个单独的IMF分别估计其对应的瞬时频率轨迹并绘制成图表形式输出。