最近看论文，看到窗函数的内容，之前也了解过，但是感觉不具体，所以这里再整理一下，方便以后查阅。
“泄漏”
在做信号处理时，经常涉及到“泄漏”。那泄漏是什么，是什么原因造成了泄漏呢？
信号截断
一次FFT分析截取1帧长度的时域信号，这1帧的长度总是有限的，因为FFT分析一次只能分析有限长度的时域信号。
而实际采集的时域信号总时间很长，因此，需要将采样时间很长的时域信号截断成一帧一帧长度的数据块。
这个截取过程叫做信号截断。
信号截断分为周期截断和非周期截断。周期截断是指截断后的信号为周期信号，而非周期截断是指截断后的信号不再是周期信号，哪怕原始信号本身是周期信号。
周期截断
周期信号最明显的特征是信号的起始和结束时刻的幅值相等，哪怕是一个周期。
在这假设采样时间很长的信号为单频正弦波（周期信号），若1帧的时间长度等于这个正弦波周期的整数倍，那么，截断后的信号仍为周期信号。
将这个截断后的信号再重构，可以得到原始的正弦波，如下图所示。
对截断的这一帧信号做FFT分析，得到它的频谱如下图所示。从图中可以看出，得到的频率成分为原始信号的真实频率，并且幅值与原始信号的幅值相等（100%幅值）。
非周期截断
倘若信号截断的长度不为原始正弦信号周期的整数倍，那么，截断后的信号则不为周期信号，哪怕原始信号是周期信号。并且现实世界中，我们进行FFT分析时，绝大多数情况都是非周期截断。
对之前的正弦信号进行非周期截断，如下图所示。截断后的信号起始时刻和结束时刻的幅值明显不等，将这个信号再进行重构，在连接处信号的幅值不连续，出现跳跃，如图中黑色圆圈区域所示。
对截断后的信号做FFT分析，得到的频谱如下图所示。
对比周期截断的频谱，可以看出，此时频谱在整个频带上发生“拖尾”现象。
峰值处的频率与原始信号的频率相近，但并不相等。
另一方面，峰值处的幅值已不再等于原始信号的幅值，为原始信号幅值的64%（矩形窗的影响）。而幅值的其他部分（36%幅值）则分布在整个频带的其他谱线上。
泄漏
由于信号的非周期截断，导致频谱在整个频带内发生了拖尾现象。这是非常严重的误差，称为泄漏，是数字信号处理所遭遇的最严重误差。
但是为什么会出现这种误差呢？原始实际信号为一条单频正弦波，它的频谱怎么会变得如此失真？这个问题很容易解释。这是因为截断后的信号不再是周期信号。
对比一下正确的频谱与发生泄漏的频谱，如下图所示，可以看出，泄漏后的频谱的幅值更小，频谱拖尾更严重。当截断后的信号不为周期信号时，就会发生泄漏。而现实世界中，在做FFT分析时，很难保证截断的信号为周期信号，因此，泄漏不可避免。
现在返回到非周期截断的正弦波，可以看出在截断时间长度内没有捕捉到整数倍个周期正弦波，导致波形发生了失真，似乎在信号周期的末端波形出现了不连续。这就解释了为什么FFT会在整个频带上发生拖尾现象了。
本质上，这需要多个傅立叶展开项（多条谱线）去近似这个明显不连续的信号，因此，频谱出现了拖尾现象。
窗函数
为了将上述·泄漏·误差减少到最小程度（注意是减少，而不是消除），我们需要使用·加权函数·，也叫·窗函数·。加窗主要是为了使时域信号似乎更好地满足FFT处理的周期性要求，减少泄漏。
为什么加窗函数
如下图所示，若周期截断，则FFT频谱为单一谱线。若为非周期截断，则频谱出现拖尾，如图中部所示，可以看出泄漏很严重。
为了减少泄漏，给信号施加一个窗函数（如图中上部红色曲线所示），原始截断后的信号与这个窗函数相乘之后得到的信号为上面右侧的信号。可以看出，此时，信号的起始时刻和结束时刻幅值都为0，也就是说在这个时间长度内，信号为周期信号，但是只有一个周期。
对这个信号做FFT分析，得到的频谱如下部右侧所示。相比较之前未加窗的频谱，可以看出，泄漏已明显改善，但并没有完全消除。因此，窗函数只能减少泄漏，不能消除泄漏。
窗函数的定义
信号截断时，只能截取一定长度，哪怕原始信号是无限长的，因此，好像是用一个“窗”（确切地说更像个“框”）去作这样的截取了。
如下图所示，原始信号是周期信号，时间很长，截取时用红色的“窗”去截取这个周期信号，截取得到的信号如图中下部所示。
上图中用于截取信号的时域截取函数（就是上图中红色的那个“窗”）就称为·窗函数·，它是一种·加权函数·，不同的窗函数加权是不一样的。
也就是说，可以用不同的截取函数（窗函数）来做信号截取。到底用何种窗函数基于信号类型和分析目的。常用的窗函数有矩形窗、汉宁窗、平顶窗、指数窗等。
窗函数的典型频谱特征
为了减少泄漏，可采用不同的窗函数来进行信号截取，因而，泄漏与窗函数的频谱特征相关的。窗函数的典型频谱特征如下图所示：
各种窗函数频谱特征的主要差别在于：主瓣宽度（也称为有效噪声带宽，ENBW）、幅值失真度、最高旁瓣高度和旁瓣衰减速率等参数。
加窗的主要想法是用比较光滑的窗函数代替截取信号样本的矩形窗函数，也就是对截断后的时域信号进行特定的不等加权，使被截断后的时域波形两端突变变得平滑些，以此压低谱窗的旁瓣。
因为旁瓣泄露量最大，旁瓣小了泄露也相应减少了。不同的窗函数具有不同的频谱特征，下表列出了一些常用窗函数的特征。
主瓣宽度主要影响信号能量分布和频率分辨能力。
频率的实际分辨能力为有效噪声带宽乘以频率分辨率，因此，主瓣越宽，有效噪声带宽越宽，在频率分辨率相同的情况下，频率的分辨能力越差。
如下图所示，红色为平顶窗（3.77∆f），黑色为汉宁窗（1.5∆f），蓝色为信号频率，可以明显地看出，主瓣越窄，频率分辨越准确。
对于窗函数宽的主瓣而言，如果有邻近的小峰值频率，则越难辨别出来。
旁瓣高低及其衰减率影响能量泄漏程度（频谱拖尾效应）。旁瓣越高，说明能量泄漏越严重，衰减越慢，频谱拖尾越严重。
对50.5Hz（频率分辨率为1Hz）的信号分别施加矩形窗（红色）、汉宁窗（绿色）和平顶窗（蓝色），用对数显示幅值，加窗后的结果如下图所示。从图中可以看出，矩形窗的频谱拖尾更严重。
相对而言，如果旁瓣能量较小，高度趋于零，使得信号能量相对集中于主瓣，则较为接近真实的频谱。
不同的窗函数对信号频谱的影响是不一样的，这主要是因为不同的窗函数，产生泄漏的大小不一样，频率分辨能力也不一样。
加窗函数的原则
加窗函数时，应使窗函数频谱的主瓣宽度应尽量窄，以获得高的频率分辨能力；旁瓣衰减应尽量大，以减少频谱拖尾，但通常都不能同时满足这两个要求。各种窗的差别主要在于集中于主瓣的能量和分散在所有旁瓣的能量之比。
窗的选择取决于分析的目标和被分析信号的类型。一般说，有效噪声频带越宽，频率分辨能力越差，越难于分清有相同幅值的邻近频率。选择性（即分辨出强分量频率邻近的弱分量的能力）的提高与旁瓣的衰减率有关。通常，有效噪声带宽窄的窗，其旁瓣的衰减率较低，因此窗的选择是在二者中取折衷。
因而，窗函数的选择一般原则如下：
-
如果截断的信号仍为周期信号，则不存在泄漏，无须加窗，相当于加矩形窗。
-
如果信号是随机信号或者未知信号，或者有多个频率分量，测试关注的是频率点而非能量大小，建议选择汉宁窗，像LMS Test.Lab中默认加的就是汉宁窗。
-
对于校准目的，则要求幅值精确，平顶窗是个不错的选择。
-
如果同时要求幅值精度和频率精度，可选择凯塞窗。
-
如果检测两个频率相近、幅值不同的信号，建议用布莱克曼窗。
窗函数带来的影响
窗函数会使信号幅值失真，那么窗函数对计算RMS值是否有影响呢？由于加窗使得频率峰值失真，因此，如果计算峰值处的RMS值，必然也是有影响的。如下图所示，由于峰值高低不一样，则对应的RMS也不一样。但如果计算窄带RMS或整个频带的总RMS值呢？
从上图可以看出，不同的窗函数下，计算19-87Hz内的总有效值都为0.71，因此，对于不同的窗函数下，计算总有效值是没有影响的。因为能量虽然泄漏到旁瓣上，但总的能量是不变的。
一些常用窗函数
窗函数有很多很多种，Scipy中的get_window提供了多种窗函数的实现，具体请参看scipy.signal.get_window
rectangular window
API reference
Examples
Plot the window and its frequency response:
import numpy as np
from scipy import signal
from scipy.fft import fft, fftshift
import matplotlib.pyplot as plt
plt.figure()
plt.subplot(211)
window = signal.windows.boxcar(51)
plt.plot(window)
plt.title("Boxcar window")
plt.ylabel("Amplitude")
plt.xlabel("Sample")
plt.subplot(212)
A = fft(window, 2048) / (len(window)/2.0)
freq = np.linspace(-0.5, 0.5, len(A))
response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
plt.plot(freq, response)
plt.axis([-0.5, 0.5, -120, 0])
plt.title("Frequency response of the boxcar window")
plt.ylabel("Normalized magnitude [dB]")
plt.xlabel("Normalized frequency [cycles per sample]")
plt.show()
hanning window
import numpy as np
from scipy import signal
from scipy.fft import fft, fftshift
import matplotlib.pyplot as plt
plt.figure()
plt.subplot(211)
window = signal.windows.hanning(51)
plt.plot(window)
plt.title("hanning window")
plt.ylabel("Amplitude")
plt.xlabel("Sample")
plt.subplot(212)
A = fft(window, 2048) / (len(window)/2.0)
freq = np.linspace(-0.5, 0.5, len(A))
response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
plt.plot(freq, response)
plt.axis([-0.5, 0.5, -120, 0])
plt.title("Frequency response of the hanning window")
plt.ylabel("Normalized magnitude [dB]")
plt.xlabel("Normalized frequency [cycles per sample]")
plt.show()
hamming window
API reference
The Hamming window is defined as
Examples
Plot the window and its frequency response:
import numpy as np
from scipy import signal
from scipy.fft import fft, fftshift
import matplotlib.pyplot as plt
plt.figure()
plt.subplot(211)
window = signal.windows.hamming(51)
plt.plot(window)
plt.title("hamming window")
plt.ylabel("Amplitude")
plt.xlabel("Sample")
plt.subplot(212)
A = fft(window, 2048) / (len(window)/2.0)
freq = np.linspace(-0.5, 0.5, len(A))
response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
plt.plot(freq, response)
plt.axis([-0.5, 0.5, -120, 0])
plt.title("Frequency response of the hamming window")
plt.ylabel("Normalized magnitude [dB]")
plt.xlabel("Normalized frequency [cycles per sample]")
plt.show()
Hann window
API reference
The Hann window is defined as
Examples
Plot the window and its frequency response:
import numpy as np
from scipy import signal
from scipy.fft import fft, fftshift
import matplotlib.pyplot as plt
plt.figure()
plt.subplot(211)
window = signal.windows.hann(51)
plt.plot(window)
plt.title("hann window")
plt.ylabel("Amplitude")
plt.xlabel("Sample")
plt.subplot(212)
A = fft(window, 2048) / (len(window)/2.0)
freq = np.linspace(-0.5, 0.5, len(A))
response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
plt.plot(freq, response)
plt.axis([-0.5, 0.5, -120, 0])
plt.title("Frequency response of the hann window")
plt.ylabel("Normalized magnitude [dB]")
plt.xlabel("Normalized frequency [cycles per sample]")
plt.show()
flat top window
API reference
Examples
Plot the window and its frequency response:
import numpy as np
from scipy import signal
from scipy.fft import fft, fftshift
import matplotlib.pyplot as plt
plt.figure()
plt.subplot(211)
window = signal.windows.flattop(51)
plt.plot(window)
plt.title("flattop window")
plt.ylabel("Amplitude")
plt.xlabel("Sample")
plt.subplot(212)
A = fft(window, 2048) / (len(window)/2.0)
freq = np.linspace(-0.5, 0.5, len(A))
response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
plt.plot(freq, response)
plt.axis([-0.5, 0.5, -120, 0])
plt.title("Frequency response of the flattop window")
plt.ylabel("Normalized magnitude [dB]")
plt.xlabel("Normalized frequency [cycles per sample]")
plt.show()
Blackman window
API reference
The Blackman window is defined as
Examples
Plot the window and its frequency response:
import numpy as np
from scipy import signal
from scipy.fft import fft, fftshift
import matplotlib.pyplot as plt
plt.figure()
plt.subplot(211)
window = signal.windows.blackman(51)
plt.plot(window)
plt.title("blackman window")
plt.ylabel("Amplitude")
plt.xlabel("Sample")
plt.subplot(212)
A = fft(window, 2048) / (len(window)/2.0)
freq = np.linspace(-0.5, 0.5, len(A))
response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
plt.plot(freq, response)
plt.axis([-0.5, 0.5, -120, 0])
plt.title("Frequency response of the blackman window")
plt.ylabel("Normalized magnitude [dB]")
plt.xlabel("Normalized frequency [cycles per sample]")
plt.show()
Blackman-Harris window
API reference
Examples
Plot the window and its frequency response:
import numpy as np
from scipy import signal
from scipy.fft import fft, fftshift
import matplotlib.pyplot as plt
plt.figure()
plt.subplot(211)
window = signal.windows.blackmanharris(51)
plt.plot(window)
plt.title("blackmanharris window")
plt.ylabel("Amplitude")
plt.xlabel("Sample")
plt.subplot(212)
A = fft(window, 2048) / (len(window)/2.0)
freq = np.linspace(-0.5, 0.5, len(A))
response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
plt.plot(freq, response)
plt.axis([-0.5, 0.5, -120, 0])
plt.title("Frequency response of the blackmanharris window")
plt.ylabel("Normalized magnitude [dB]")
plt.xlabel("Normalized frequency [cycles per sample]")
plt.show()
2144
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
