频谱泄露与加窗

wave1,wave2实际在时域中是一样的正弦波，表达式y=sin(2πft)，其中频率f=20Hz,但是经过快速傅里叶FFT变换之后，wave2的频域图中除了有冲激响应的峰之外，还出现了额外的频率元素，这个就是频谱泄露(Spectrum Leakage)。

很显然，频谱泄露和信号泄露无关，从根本上说，它就是算法引起的。那同样的快速傅里叶算法,怎么会有这个差异？

简单说来，造成这个问题的原因是：采样信号的相位在始端和终端不连续。

从上面图中，可以发现，wave1的时域图形是完整的周期，而wave2不是。当我们知道频率的时候，当然可以很容易算出周期，但是我们实际应用中，FFT的目的就是为了得到复杂信号的各个频率分量，所以周期也是未知的，也无法同时满足所有信号频率。

这个时候，我们需要采用加窗的方式来减轻频谱泄露。

加窗，就是通过算法让采样信号首尾不连续的相位趋向连续。这种窗函数有很多，不同的窗函数对应不同的应用场景会有比较好的效果，其中汉明(Hanmming)窗是一个比较好的窗函数。

我们可以调用了Python中的numpy的hamming窗函数来实现。把窗函数和原信号函数相乘，就得到了处理过的时域采样信号，然后再进行快速傅里叶转换。

 

import numpy as np
import matplotlib.pyplot as plt

# ======================
# 参数设置
# ======================
sample_rate = 1000       # 采样率 (Hz)
signal_freq = 20         # 信号频率 (Hz)
duration_perfect = 1.0   # 完整周期时长 (整数倍周期)
duration_leak = 1.02     # 非完整周期时长 (制造频谱泄漏)

# ======================
# 信号生成函数
# ======================
def generate_signal(duration):
    """生成正弦波信号"""
    n_samples = int(sample_rate * duration)  # 关键修正：确保样本数正确
    t = np.linspace(0, duration, n_samples)
    return t, np.sin(2 * np.pi * signal_freq * t)

# ======================
# 计算FFT函数
# ======================
def compute_fft(signal, sample_rate):
    """计算信号的FFT"""
    n = len(signal)
    freq = np.fft.fftfreq(n, 1/sample_rate)
    fft_values = np.fft.fft(signal)
    magnitude = np.abs(fft_values) / n * 2   # 幅度归一化
    return freq[:n//2], magnitude[:n//2]     # 仅取正频率部分

# ======================
# 生成原始信号
# ======================
# 完整周期的信号 (无泄漏)
t1, y1 = generate_signal(duration_perfect)
freq1, mag1 = compute_fft(y1, sample_rate)

# 非完整周期的信号 (有泄漏)
t2, y2 = generate_signal(duration_leak)
freq2, mag2 = compute_fft(y2, sample_rate)

# ======================
# 加窗处理
# ======================
hamming_window = np.hamming(len(y2))          # 使用正确样本数生成汉明窗
y3 = y2 * hamming_window
freq3, mag3 = compute_fft(y3, sample_rate)

# ======================
# 可视化
# ======================
plt.figure(figsize=(12, 10))

# ---- 原始信号对比 ----
# 时域波形
plt.subplot(3, 2, 1)
plt.plot(t1, y1)
plt.title('Perfect Signal (Time Domain)\n{}周期'.format(int(signal_freq*duration_perfect)))
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.grid(True)
plt.ylim(-1.5, 1.5)

plt.subplot(3, 2, 3)
plt.plot(t2, y2, color='orange')
plt.title('Leaky Signal (Time Domain)\n{:.2f}周期'.format(signal_freq*duration_leak))
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.grid(True)
plt.ylim(-1.5, 1.5)

# 频域分析
plt.subplot(3, 2, 2)
plt.plot(freq1, mag1)
plt.title('Frequency Spectrum (No Leakage)')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.xlim(0, 50)  # 聚焦关注频段
plt.grid(True)

plt.subplot(3, 2, 4)
plt.plot(freq2, mag2, color='orange')
plt.title('Frequency Spectrum (With Leakage)')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.xlim(0, 50)
plt.grid(True)

# ---- 加窗处理效果 ----
# 时域窗函数
plt.subplot(3, 2, 5)
plt.plot(t2, hamming_window, color='green')
plt.title('Hamming Window (Time Domain)')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.grid(True)

# 加窗后频域
plt.subplot(3, 2, 6)
plt.plot(freq3, mag3, color='red')
plt.title('Windowed Signal Spectrum')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.xlim(0, 50)
plt.grid(True)

plt.tight_layout()
plt.show()

进一步优化建议：

添加交互

from ipywidgets import interact

@interact(freq=(10, 50, 1), duration=(0.5, 2.0, 0.1))
def interactive_plot(freq=20, duration=1.0):
    # 更新信号生成和绘图逻辑
    ...

添加信噪比指标

def calculate_snr(magnitude, target_freq, bw=5):
    main_lobe = np.where((freq > target_freq-bw) & (freq < target_freq+bw))[0]
    noise_floor = np.delete(magnitude, main_lobe)
    return 10*np.log10(np.max(magnitude)**2 / np.mean(noise_floor**2))

支持多窗

windows = {
    'Hamming': np.hamming,
    'Hann': np.hanning,
    'Blackman': np.blackman,
    'Flattop': lambda n: np.flatTop(n)
}