python 之窗函数
python 之窗函数
在数字信号处理中，由于只能使用有限个时域数据进行傅立叶变换，因此需要对时域数据进行截断处理。
信号截断分为周期截断和非周期截断。
- 周期截断是指截断后的信号为周期信号，这种截断不会导致频谱泄漏。
- 非周期截断是指截断后的信号不再是周期信号，哪怕原始信号本身是周期信号，这中截断会导致频谱泄漏。
为了减少截断导致的频谱泄漏，我们需要使用窗函数对截断的信号进行处理。
窗函数相关总结：
- 主瓣越窄，频率分辨越准确
- 旁瓣越高，说明能量泄漏越严重
- 旁瓣衰减越慢，频谱拖尾越严重
周期截断
当对周期时域数据进行截断时，只要保证傅立叶变换的 N 个时域点包含完整的周期即可保证频谱不泄漏。 实验方法如下：
- 使用采样频率为 8000 Hz, 使用 N = 512 个采样点进行傅立叶变换（采样分辨率为 15.625 Hz)
- 令原始信号包含 156.25 和 234.375 两个正弦信号
- 不加窗（默认为矩形窗）对原始进行傅立叶变换，查看 FFT 变换是否有频谱泄漏
%matplotlib notebook
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
sample_rate = 8000
fft_size = 512
t = np.arange(0, 1.1, 1.0/sample_rate)
## 两个正弦信号，频率分辨为 156.25 和 234.375
x = np.sin(2*np.pi*156.25*t) + 2*np.sin(2*np.pi*234.375*t)
## 取 N 个采样点进行 FFT 变换
# rfft 函数的返回 N/2+1 个复数，分别表示从 0(Hz) ~ fs/2 频率分量
xs = x[:fft_size]
xf = np.fft.rfft(xs) / fft_size
## 生成频率和振幅(dB)
freqs = np.linspace(0, sample_rate/2, int(fft_size/2 + 1))
xfp = 20 * np.log10(np.clip(np.abs(xf), 1e-20, 1e1000))
## 画出波形及其频谱图
plt.figure(figsize=(8,4))
plt.subplot(211)
plt.plot(t[:fft_size], xs)
plt.xlabel("时间(秒)")
plt.ylabel("振幅")
plt.title("振幅和频率（156.25Hz 和 234.375Hz）")
plt.subplot(212)
plt.plot(freqs, xfp)
plt.xlabel("频率(Hz)")
plt.ylabel("振幅(dB)")
plt.subplots_adjust(hspace=0.4)
plt.show()
从上面的图形可以看出，FFT 变换可以准确的计算出 156.25 和 234.375 两个频率点，未发生频谱泄漏。也就是说只有当输入数据序列中的分析频率是频率分辨率（$f_x / N$) 整数倍上包含能量时，FFT 才能得出正确的结果。
非周期截断
当输入信号包含其他非分辨率整数倍信号时，则使用 FFT 计算会出现频谱泄漏。
下面我们用同样的方法试一下 100Hz （非 15.625Hz 倍数）时，FFT 变化出来的频谱图。
## 100Hz 正弦信号，不是 15.625 倍数时，出现频谱泄漏
x = np.sin(2*np.pi*100*t)
## 取 N 个采样点进行 FFT 变换
# rfft 函数的返回 N/2+1 个复数，分别表示从 0(Hz) ~ fs/2 频率分量
xs = x[:fft_size]
xf = np.fft.rfft(xs) / fft_size
## 生成频率和振幅(dB)
freqs = np.linspace(0, sample_rate/2, int(fft_size/2+1))
xfp = 20 * np.log10(np.clip(np.abs(xf), 1e-20, 1e1000))
## 画出波形及其频谱图
plt.figure(figsize=(8,4))
plt.subplot(211)
plt.plot(t[:fft_size], xs)
plt.xlabel("时间(秒)")
plt.ylabel("振幅")
plt.title("振幅和频率（100Hz）")
plt.subplot(212)
plt.plot(freqs, xfp)
plt.xlabel("频率(Hz)")
plt.ylabel("振幅(dB)")
plt.subplots_adjust(hspace=0.4)
plt.show()
上面频谱图看出，100Hz 信号使用 512 点 FFT 变换后频谱泄漏严重。
100Hz 峰值处的幅值减小，泄漏到整个频带的其他谱线上。频谱图已经不能准确反映原始信号的频率了。
泄漏的主要原因时信号截断导致，截断信号等价与给信号加了一个矩形窗，其幅频响应为 sinc 函数。
上图 FFT 计算出来的频谱图之所以不时关于 100Hz 为中心的 sinc 对称函数，是因为频谱周期性缠绕和 N/2 对称缠绕导致的。
频谱泄漏时无法完全避免的，但是可以通过对输入信号加窗来减小它的不良影响。
窗函数
窗函数通过最小化 sinc 的旁瓣来降低频谱泄漏。这项工作通过同时将输入信号在采样间隔的起始和末尾点的振幅强制平滑为一个单一的共同值来完成（使其周期化，会降低高频分量的振幅）。下面用 python 代码画出矩形窗、三角窗和汉宁窗的振幅及其幅频响应。
%matplotlib notebook
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
from scipy.fftpack import fft, fftshift
# 主瓣和旁瓣个数
win_size = 11
fft_size = 2048
# 生成矩形窗、三角窗、含宁窗数据
rect_win = np.ones(win_size)
triang_win = signal.triang(win_size)
hanning_win = signal.hanning(win_size)
# 计算三种窗的频率响应
rect_h = fft(rect_win, fft_size)
triang_h = fft(triang_win, fft_size)
hanning_h = fft(hanning_win, fft_size)
freqs = np.linspace(-0.5, 0.5, fft_size)
def a2db(amp):
return 20 * np.log10(np.clip(np.abs(fftshift(amp / abs(amp).max())), 1e-5, 1e100))
plt.figure(figsize=(10,4))
plt.subplot(121)
plt.xlabel('采样点数')
plt.ylabel('振幅')
plt.title('窗函数')
plt.plot(rect_win, color='b', label='矩形窗')
plt.plot(triang_win, color='g', label='三角窗')
plt.plot(hanning_win, color='r', label='汉宁窗')
plt.legend()
plt.subplot(122)
plt.title('幅频响应')
plt.xlabel('归一化频率')
plt.ylabel('归一化增益(dB)')
plt.axis([-0.5, 0.5, -120, 0])
plt.plot(freqs, a2db(rect_h), color='b', label='矩形窗')
plt.plot(freqs, a2db(triang_h), color='g', label='三角窗')
plt.plot(freqs, a2db(hanning_h), color='r', label='汉宁窗')
plt.legend()
plt.show()
从上图可以看出，矩形窗（蓝色）的主瓣比较窄更有利于识别出指定频率，但是旁瓣增益较高，泄漏严重导致振幅信息不准。三角窗的主瓣比矩形窗宽一些，旁瓣泄漏也相对小一些。汉宁窗主瓣更宽，旁瓣衰减更快。
一般来说，我们希望主瓣越窄越好，旁瓣衰减越快越好，但是这两点不可同时获得，需要根据实际情况，折中选择。
下面我们再看一下汉宁窗的幅频及其相频响应，并从图中看出汉宁窗是线性相位的。
%matplotlib notebook
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
win_size = 11
## 矩形窗
hann_window = signal.hanning(win_size)
hann_w, hann_h = signal.freqz(rect_window)
db = 20 * np.log10(np.clip(np.abs(hann_h/abs(hann_h).max()), 1e-5, 1e100))
plt.figure()
plt.title('汉宁窗幅频、相频响应')
plt.plot(hann_w, db, color='b', label='振幅')
plt.ylabel('振幅 [dB]', color='b')
plt.xlabel('频率 [rad/sample]')
plt.legend()
plt.twinx()
rect_angles = np.unwrap(np.angle(rect_h))
plt.plot(rect_w, rect_angles, color='g', label='相位')
plt.ylabel('相位 (弧度)', color='g')
plt.grid()
plt.axis('tight')
plt.legend()
plt.show()
