当前位置：首页 > ai >正文

Matplotlib 库来可视化频谱泄漏和加窗的效果

ai 2025/6/28 2:28:47

前言

很多朋友学习音频技术的时候，不理解这个频谱泄漏是什么，我们这次写个小代码直观地感受一下

代码演示：频谱泄漏与加窗

我们将生成一个简单的正弦波信号，然后分别用**不加窗（矩形窗）和加窗（汉明窗）**的方式对其进行傅里叶变换，并对比它们的频谱图。

你会清晰地看到加窗如何减少了频谱泄漏。

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
from scipy.signal import hamming# --- 1. 信号参数设置 ---
fs = 1000       # 采样频率 (Hz) - 每秒采集1000个点
T = 1           # 信号总时长 (秒)
N = int(fs * T) # 采样点数# 生成时间轴
t = np.linspace(0, T, N, endpoint=False)# 生成一个包含两个频率成分的信号
# 注意：这里我们故意让频率不是采样窗口长度的整数倍，以更容易观察泄漏
f1 = 50.5      # 第一个频率成分 (Hz)
f2 = 120.5     # 第二个频率成分 (Hz)
amplitude1 = 1.0
amplitude2 = 0.5# 原始信号 (纯粹的复合波)
signal_pure = amplitude1 * np.sin(2 * np.pi * f1 * t) + \amplitude2 * np.sin(2 * np.pi * f2 * t)# --- 2. 频谱分析 (不加窗 vs. 加汉明窗) ---# --- 不加窗 (等同于矩形窗) ---
# 直接对信号进行傅里叶变换
yf_no_window = fft(signal_pure)
# 计算对应的频率轴
xf_no_window = fftfreq(N, 1 / fs)# --- 加汉明窗 ---
# 生成汉明窗
window = hamming(N)
# 将信号与汉明窗相乘
signal_windowed = signal_pure * window# 对加窗后的信号进行傅里叶变换
yf_windowed = fft(signal_windowed)
# 计算对应的频率轴
xf_windowed = fftfreq(N, 1 / fs)# --- 3. 可视化 ---
# 设置中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号plt.figure(figsize=(15, 10))# 原始时域信号
plt.subplot(3, 1, 1)
plt.plot(t, signal_pure)
plt.title('原始时域信号 (Pure Signal)')
plt.xlabel('时间 (秒)')
plt.ylabel('幅值')
plt.grid(True)
plt.xlim(0, 0.2) # 只显示一部分，方便观察# 不加窗的频谱 (矩形窗)
plt.subplot(3, 1, 2)
# 我们只关注正频率部分，并取幅值的绝对值（能量）
plt.plot(xf_no_window[:N//2], 2.0/N * np.abs(yf_no_window[:N//2]))
plt.title('频谱 (不加窗/矩形窗) - 泄漏明显')
plt.xlabel('频率 (Hz)')
plt.ylabel('归一化幅值')
plt.grid(True)
plt.xlim(0, 150) # 限制显示范围，聚焦于有效频率
plt.axvline(f1, color='r', linestyle='--', label=f'实际频率 {f1} Hz')
plt.axvline(f2, color='g', linestyle='--', label=f'实际频率 {f2} Hz')
plt.legend()# 加汉明窗的频谱
plt.subplot(3, 1, 3)
plt.plot(xf_windowed[:N//2], 2.0/N * np.abs(yf_windowed[:N//2]))
plt.title('频谱 (加汉明窗) - 泄漏减轻')
plt.xlabel('频率 (Hz)')
plt.ylabel('归一化幅值')
plt.grid(True)
plt.xlim(0, 150) # 限制显示范围，聚焦于有效频率
plt.axvline(f1, color='r', linestyle='--', label=f'实际频率 {f1} Hz')
plt.axvline(f2, color='g', linestyle='--', label=f'实际频率 {f2} Hz')
plt.legend()plt.tight_layout() # 自动调整子图参数，使之填充整个图像区域
plt.show()# --- 4. 解释分析 ---
print(f"生成的信号频率为 {f1} Hz 和 {f2} Hz。")
print("观察上方频谱图：")
print("1. 不加窗（矩形窗）的频谱图：")
print("   - 在实际频率 {f1} Hz 和 {f2} Hz 附近，主峰虽然存在，但两侧有明显的拖尾（旁瓣）。")
print("   - 这些旁瓣使得能量扩散，可能模糊了相邻频率的区分，这就是频谱泄漏。")
print("2. 加汉明窗的频谱图：")
print("   - 相同频率 {f1} Hz 和 {f2} Hz 处的主峰变得更窄，能量更集中。")
print("   - 主峰两侧的旁瓣（拖尾）被大大抑制，几乎看不见了。")
print("   - 这表明加窗有效地减轻了频谱泄漏，使得频率成分的识别更加清晰和准确。")

代码解释：

信号生成：
- 我们创建了一个包含两个不同频率（f1 和 f2）的正弦波信号。
- 关键在于，f1 和 f2 都不是采样点数 N 的整数因子，这意味着在 T=1 秒的采样窗口内，它们不包含整数个周期。这是为了故意制造频谱泄漏，以便你观察。
傅里叶变换：
- scipy.fft.fft 用于执行快速傅里叶变换 (FFT)。
- scipy.fft.fftfreq 用于计算 FFT 结果对应的频率轴。
不加窗（矩形窗）：
- yf_no_window = fft(signal_pure)：直接对原始信号进行 FFT。这相当于给信号应用了一个矩形窗，因为所有采样点都被等权重地包含了，而边界是突然截断的。
加汉明窗：
- window = hamming(N)：生成一个长度为 N 的汉明窗。汉明窗的特点是两端平滑地衰减到零。
- signal_windowed = signal_pure * window：将原始信号与汉明窗相乘。这在时域上平滑了信号的起始和结束边界。
- yf_windowed = fft(signal_windowed)：对加窗后的信号进行 FFT。
可视化与对比：
- 代码绘制了三张图：原始时域信号、不加窗的频谱和加汉明窗的频谱。
- 请重点对比第二张图和第三张图。 你会发现在不加窗的频谱图中，每个频率峰值旁都有明显的“小山丘”或“拖尾”（旁瓣），这就是泄漏。而加汉明窗后，这些旁瓣被大大抑制了，频率峰值变得更尖锐，能量更集中。
  
  通过运行这段代码，你将能够直观地看到频谱泄漏是如何发生的，以及加窗在音频（或任何信号）分析中解决这个问题的实际效果。