STFT (短时傅立叶变换)

短时傅立叶变换 (STFT) 详细介绍

1. 基本概念

傅立叶变换（Fourier Transform）用于将一个信号从时间域转换到频率域，然而它假设信号是平稳的，这意味着信号的频率成分在整个时间上是不变的。对于非平稳信号，傅立叶变换无法提供信号随时间变化的频率信息。

短时傅立叶变换通过将信号分割成多个短的、重叠的时间窗口，每个窗口内的信号近似为平稳的，然后对每个窗口分别进行傅立叶变换，从而获得信号在不同时刻的频率信息。

2. 数学表达

• STFT的基本思想是将信号分割成多个小的时间窗口，在每个窗口内假设信号是平稳的，然后对每个窗口内的信号进行傅立叶变换。其数学表达式为：

  $$X(t,f)={\int }_{-\infty }^{\infty }x(\tau )w(t-\tau )e^{-j2\pi f\tau }d\tau$$

  其中：

  • $x(\tau )$ 是原始信号。
  • $w(t-\tau )$ 是窗口函数。
  • $e^{-j2\pi f\tau }$ 是傅立叶变换的核函数。

STFT 结果

STFT 的计算结果通常是一个复数矩阵，这个矩阵描述了信号在不同时间和频率上的特性。具体来说，STFT 结果的每个元素都是一个复数，表示信号在某个时间段和频率上的振幅和相位。下面是对这些信息的详细解释：

1. 复数矩阵表示

STFT 结果是一个复数矩阵 (Z_{xx}(t, f))，其中 (t) 表示时间，(f) 表示频率。每个元素 (Z_{xx}(t, f)) 是一个复数，包含了信号在时间 (t) 和频率 (f) 上的振幅和相位信息。

2. 幅度 (Magnitude)

复数矩阵中的每个元素的模值表示信号在该时间点和频率上的振幅（或强度）。振幅的计算公式为：

$$|Z_{xx}(t,f)|=\sqrt{\text{Re}(Z_{xx}(t,f){)}^2+\text{Im}(Z_{xx}(t,f){)}^2}$$

其中，(\text{Re}(Z_{xx}(t, f))) 和 (\text{Im}(Z_{xx}(t, f))) 分别表示复数矩阵的实部和虚部。

3. 相位 (Phase)

复数矩阵中的每个元素的角度值表示信号在该时间点和频率上的相位。相位的计算公式为：

$$\angle Z_{xx}(t,f)=\text{atan2}(\text{Im}(Z_{xx}(t,f)),\text{Re}(Z_{xx}(t,f)))$$

1. 时频图 (Spectrogram)

时频图是通过短时傅立叶变换（STFT）得到的，显示信号在时间和频率上的变化情况。它的横轴表示时间，纵轴表示频率，颜色或亮度表示信号在不同时间和频率上的强度。

生成方法

时频图通过对信号进行STFT得到。STFT将信号分割成多个短时间片段，在每个片段上进行傅立叶变换，得到每个时间点上的频率成分。

示例代码

以下是使用Python生成时频图的示例代码：

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

# 生成一个示例信号：两个不同频率的正弦波
fs = 1000  # 采样频率
t = np.arange(0, 2, 1/fs)
x = np.sin(2 * np.pi * 50 * t) + np.sin(2 * np.pi * 120 * t)

# 计算STFT
f, t, Zxx = stft(x, fs, nperseg=256)

# 绘制时频图
plt.pcolormesh(t, f, np.abs(Zxx), shading='gouraud')
plt.title('STFT Magnitude')
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.colorbar(label='Magnitude')
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

2. 幅度谱 (Magnitude Spectrum)

幅度谱显示信号在每个时间点上的频率强度，即信号在频域的能量分布。幅度谱由STFT的结果的模值计算得到。

计算方法

幅度谱的计算公式为：

$$|X(t,f)|=\sqrt{\text{Re}(X(t,f){)}^2+\text{Im}(X(t,f){)}^2}$$

其中，$\text{Re}(X(t,f))$ 和 $\text{Im}(X(t,f))$ 分别表示STFT结果的实部和虚部。

示例代码

下面是从STFT结果中提取并绘制幅度谱的示例代码：

# 计算幅度谱
magnitude_spectrum = np.abs(Zxx)

# 绘制幅度谱
plt.pcolormesh(t, f, magnitude_spectrum, shading='gouraud')
plt.title('Magnitude Spectrum')
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.colorbar(label='Magnitude')
plt.show()
1
2
3
4
5
6
7
8
9
10

3. 相位谱 (Phase Spectrum)

相位谱显示信号在每个时间点上的相位信息，相位谱表示信号在不同频率上的相位变化。相位谱由STFT的结果的角度值计算得到。

计算方法

相位谱的计算公式为：

$$\angle X(t,f)=\text{atan2}(\text{Im}(X(t,f)),\text{Re}(X(t,f)))$$

示例代码

以下是从STFT结果中提取并绘制相位谱的示例代码：

# 计算相位谱
phase_spectrum = np.angle(Zxx)

# 绘制相位谱
plt.pcolormesh(t, f, phase_spectrum, shading='gouraud')
plt.title('Phase Spectrum')
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.colorbar(label='Phase (radians)')
plt.show()
1
2
3
4
5
6
7
8
9
10

4. 窗口函数

常用的窗口函数包括：

• 矩形窗：简单但频谱泄漏严重。
• 汉宁窗：减少频谱泄漏，适用于大多数应用。
• 海明窗：与汉宁窗相似，但具有稍微不同的频谱特性。
• 高斯窗：具有良好的时频局部化特性。

窗口长度的选择通常是一个折中：窗口过长会导致频率分辨率高但时间分辨率低；窗口过短则时间分辨率高但频率分辨率低。

4.1 矩形窗 (Rectangular Window)

矩形窗是最简单的窗口函数，定义为在窗口长度内取1，其余为0。其数学表达式为：

w ( t ) = { 1 if  ∣ t ∣ ≤ T 2 0 otherwise w(t) =

$$\begin{cases} 1 & \text{if } |t| \leq \frac{T}{2} \\ 0 & \text{otherwise} \end{cases}$$ w(t)={10​if ∣t∣≤2T​otherwise​

其中 $T$ 是窗口长度。

优点：

• 计算简单。

缺点：

• 频谱泄漏严重，频率分辨率低。

4.2 汉宁窗 (Hanning Window)

汉宁窗是常用的余弦窗之一，具有良好的频谱特性。其数学表达式为：

$$w(t)=0.5\left(1-\cos \left(\frac{2\pi t}{T}\right)\right),t\in [0,T]$$

优点：

• 频谱泄漏较少，适用于大多数应用。

缺点：

• 有一定的主瓣宽度，影响频率分辨率。

4.3 海明窗 (Hamming Window)

海明窗与汉宁窗类似，但其加权系数略有不同。其数学表达式为：

$$w(t)=0.54-0.46\cos \left(\frac{2\pi t}{T}\right),t\in [0,T]$$

优点：

• 减少了频谱泄漏，性能与汉宁窗相近。

缺点：

• 主瓣宽度稍大，频率分辨率略低。

4.4 高斯窗 (Gaussian Window)

高斯窗是一种具有良好时频局部化特性的窗函数。其数学表达式为：

$$w(t)=e^{-\frac{1}{2}{\left(\frac{t}{\sigma }\right)}^2},t\in [-T/2,T/2]$$

其中 $\sigma$ 是标准差，决定了窗口的宽度。

优点：

• 优异的时频局部化特性，适用于需要高时频分辨率的应用。

缺点：

• 计算复杂度较高。

4.5 布莱克曼窗 (Blackman Window)

布莱克曼窗是另一种常用的余弦窗，具有更低的频谱泄漏。其数学表达式为：

$$w(t)=0.42-0.5\cos \left(\frac{2\pi t}{T}\right)+0.08\cos \left(\frac{4\pi t}{T}\right),t\in [0,T]$$

优点：

• 较低的频谱泄漏。

缺点：

• 主瓣宽度较大，频率分辨率较低。

5. 窗口长度选择

窗口长度的选择通常是一个折中问题：窗口越长，频率分辨率越高，但时间分辨率越低；反之，窗口越短，时间分辨率高但频率分辨率低。因此，在实际应用中，需根据信号的特性和具体分析需求选择合适的窗口长度。

例如，在语音信号处理中，通常选择20-40ms的窗口长度，因为语音信号在这个范围内可以认为是平稳的。而在快速变化的生物医学信号处理中，可能需要更短的窗口以捕捉快速变化的特征。

频谱泄漏、主瓣与频率分辨率

1. 频谱泄漏 (Spectral Leakage)

频谱泄漏是指信号的频率成分在进行傅立叶变换后，其能量扩散到邻近频率的现象。频谱泄漏的产生是因为在有限时间窗口内截取信号，导致截取的信号边界不连续，从而在频率域上产生不期望的能量分布。

原因

• 窗口函数的选择：不同的窗口函数会对频谱泄漏有不同程度的影响。
• 窗口长度：较短的窗口可能导致更严重的频谱泄漏。

影响

• 频率分辨率降低，频率成分难以区分。
• 频谱图中出现伪频率成分，增加了解读的难度。

解决方法

• 选择适当的窗口函数（如汉宁窗、海明窗）以减少频谱泄漏。
• 适当增加窗口长度。

2. 主瓣 (Main Lobe) 和 旁瓣 (Side Lobe)

在使用窗口函数对信号进行傅立叶变换时，频谱的幅度响应通常包含一个主瓣和多个旁瓣。

主瓣

主瓣是窗口函数频谱响应的中心部分，包含信号主要的频率成分。主瓣宽度越窄，频率分辨率越高。

旁瓣

旁瓣是主瓣两侧较小的频率响应部分，旁瓣越高，频谱泄漏越严重。

关系

• 主瓣宽度：决定了信号的频率分辨率。
• 旁瓣高度：影响频谱泄漏的程度，旁瓣越低，频谱泄漏越少。

3. 频率分辨率 (Frequency Resolution)

频率分辨率是指在频谱分析中能够区分出两个相邻频率成分的能力。频率分辨率由窗口函数的主瓣宽度决定，主瓣越窄，频率分辨率越高。

计算

频率分辨率通常表示为频谱中的频率间隔，可以通过以下公式计算：

$$\Delta f=\frac{1}{T}$$

其中，$T$ 是时间窗口的长度。

4. 举例说明

考虑一个时间长度为 $T$ 的信号 $x(t)$，进行傅立叶变换时使用不同的窗口函数：

• 矩形窗：主瓣较窄，但旁瓣较高，频谱泄漏严重。
• 汉宁窗：主瓣适中，旁瓣较低，频谱泄漏相对较少。
• 高斯窗：主瓣较窄且旁瓣较低，频谱泄漏最少，但计算复杂度高。

5. 图示说明

下面是使用Python和Matplotlib绘制的示例图，用于展示不同窗口函数的频谱响应：

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

# 定义窗口函数
windows = ['rectangular', 'hann', 'hamming', 'gaussian']
window_lengths = 256
fs = 1000  # 采样频率

# 绘制不同窗口函数的频谱响应
plt.figure(figsize=(10, 6))
for window in windows:
    if window == 'gaussian':
        w = get_window(('gaussian', 7), window_lengths)
    else:
        w = get_window(window, window_lengths)
    W = np.fft.fft(w, 1024)
    W = np.fft.fftshift(W)
    freq = np.linspace(-fs/2, fs/2, len(W))
    plt.plot(freq, 20 * np.log10(np.abs(W) / np.max(np.abs(W))), label=window)

plt.title('Window Functions Frequency Response')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude (dB)')
plt.legend()
plt.grid()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27