【频域分析】频谱泄露、频率分辨率、栅栏效应

现实生活中的信号大部分是连续的，通过对连续的信号进行采样得到散时间信号，但是计算机所能处理的数据都是有限长的，因而我们可以对原始序列做加窗处理使其成为有限长序列。 以矩形窗为例，其时域表达式为： 式(1)中， N = M + 1 N=M+1 N=M+1，为矩形窗的长度。

对无限长序列进行加窗处理，就是对序列在时域上乘以一个窗函数。 由卷积定理可以得到，时域的相乘等于频域的卷积。

设仿真信号的时域表达式为： x ( t ) = A 0 ∗ c o s ( 2 π f 0 t ) + A 1 ∗ c o s ( 2 π f 1 t ) x(t)=A_{0}*cos(2πf_{0}t)+A_{1}*cos(2πf_{1}t) x(t)=A0​∗cos(2πf0​t)+A1​∗cos(2πf1​t) x ( t ) x(t) x(t)做傅里叶变换（FT）的频域表达式为： X ( j Ω ) = A 0 π δ ( Ω + Ω 0 ) + A 0 π δ ( Ω − Ω 0 ) + A 1 π δ ( Ω + Ω 0 ) + A 1 π δ ( Ω − Ω 0 ) X(jΩ)=A_{0}πδ(Ω+Ω_{0})+A_{0}πδ(Ω-Ω_{0})+A_{1}πδ(Ω+Ω_{0})+A_{1}πδ(Ω-Ω_{0}) X(jΩ)=A0​πδ(Ω+Ω0​)+A0​πδ(Ω−Ω0​)+A1​πδ(Ω+Ω0​)+A1​πδ(Ω−Ω0​) 连续信号 x ( t ) x(t) x(t)的波形及频谱如图1所示。

连续信号 x ( t ) x(t) x(t)经过采样后，得到的离散时间的表达式为： x [ n ] = x ( t ) ∣ t = n T s x[n]=x(t)|_{t=nT_{s}} x[n]=x(t)∣t=nTs​​ 离散序列 x [ n ] x[n] x[n]做离散时间傅里叶变换（DTFT）的频域表达式为： X ( e j w ) = 1 T s ∑ k = − ∞ ∞ X ( j w T s − j k 2 π T s ) X(e^{jw})=\frac{1}{T_{s}}\sum_{k=-∞}^{∞}X(j\frac{w}{T_{s}}-jk\frac{2π}{T_{s}}) X(ejw)=Ts​1​k=−∞∑∞​X(jTs​w​−jkTs​2π​) 离散序列 x [ n ] x[n] x[n]的波形及频谱如图2所示。 矩形窗函数 w [ n ] w[n] w[n]做离散时间傅里叶变换（DTFT）的频域表达式为： W ( e j w ) = e − j w ( N − 1 ) / 2 ∗ s i n ( w N / 2 ) s i n ( w / 2 ) W(e^{jw})=e^{-jw(N-1)/2} *\frac{sin(wN/2)}{sin(w/2)} W(ejw)=e−jw(N−1)/2∗sin(w/2)sin(wN/2)​

矩形窗函数 w [ n ] w[n] w[n]的波形及频谱如图3所示。 离散序列 x [ n ] x[n] x[n]与窗函数 w [ n ] w[n] w[n]的卷积为： V ( e j w ) = 1 2 π ∫ − π π X ( e j θ ) W ( e j ( w − θ ) ) d θ = A 0 2 W ( e j ( w + w 0 ) ) + A 0 2 W ( e j ( w − w 0 ) ) + A 1 2 W ( e j ( w + w 0 ) ) + A 1 2 W ( e j ( w − w 0 ) ) V(e^{jw})=\frac{1}{2π}\int_{-π}^{π}X(e^{jθ})W(e^{j(w-θ)})dθ=\frac{A_{0}}{2}W(e^{j(w+w_{0})})+\frac{A_{0}}{2}W(e^{j(w-w_{0})})+\frac{A_{1}}{2}W(e^{j(w+w_{0})})+\frac{A_{1}}{2}W(e^{j(w-w_{0})}) V(ejw)=2π1​∫−ππ​X(ejθ)W(ej(w−θ))dθ=2A0​​W(ej(w+w0​))+2A0​​W(ej(w−w0​))+2A1​​W(ej(w+w0​))+2A1​​W(ej(w−w0​)) 截断后的离散序列 v [ n ] v[n] v[n]的波形及频谱如图4所示。

信号的频率成分包括1MHz和1.05MHz，1MHz对应的幅值为1，但是1.05MHz的幅值减小了，且在其他频率点上都有不小的幅值。这就是出现了频谱泄露的现象。

由于计算机只能处理有限长的数据，所以需要对采集的信号进行截断，相当于对原始信号做了加窗处理。对信号加窗就是对信号在时域上乘以一个窗函数，时域的乘积对应频域的卷积，而窗函数的频域包括主瓣和旁瓣，旁瓣造成了信号频谱的泄漏。频域泄漏不可避免，只能减小。

可以取更长的数据点，与原始数据越接近越好，但缺点就是运算量加大； 可以选择窗谱的旁瓣能量较小的窗函数。 典型的窗函数中，矩形窗的频率分辨率最高，旁瓣泄露最大。

为了便于理解什么是频率分辨率，可以将频率分辨率划分为两种类型，一种是波形频率分辨率（Waveform Frequency Resolution，简记为波形分辨率），也称为视觉频率分辨率，另一种则为FFT频率分辨率（简记为FFT分辨率）。

波形分辨率：在频谱图中，两个频率可以被分辨率的最小间隔，与原始信号的时间长度有关。 △ R w = 1 T △R_{w}=\frac{1}{T} △Rw​=T1​ 其中， T T T为原始数据的时长。

FFT分辨率：在频谱图中的数据点数，跟信号做FFT计算时的点数有关。 △ R f f t = F s N f f t △R_{fft}=\frac{F_{s}}{N_{fft}} △Rfft​=Nfft​Fs​​ 其中， F s F_{s} Fs​为采样频率， N f f t N_{fft} Nfft​为信号做FFT计算时的点数。

例如，有一个复合信号的时域表达式为 x ( t ) = c o s ( 2 π f 1 t ) + c o s ( 2 π f 2 t ) x(t)=cos(2πf_{1}t)+cos(2πf_{2}t) x(t)=cos(2πf1​t)+cos(2πf2​t) 其中 f 1 = 1 M H z f_{1}=1MHz f1​=1MHz， f 2 = 1.05 M H z f_{2}=1.05MHz f2​=1.05MHz。 x ( t ) x(t) x(t)的波形及频谱如图6所示。 由上图中的频谱可以发现，在1MHz附近两个频率出现混叠，无法有效区分1MHz和1.05MHz，说明频率分辨率不够。

对时域数据进行补零，能否改变频率分辨率呢？例如在原始数据点后面再补充6000个数值为0的点，对信号本身数据没有影响，只是增加了参与FFT计算的数据点数，得到信号的波形及频谱如图7所示。 通过对原始数据进行补零操作，在频谱图中的数据点更密集了。但是仍然无法区分1MHz和1.05MHz，由此证明信号的波形频率分辨率与参与FFT计算的数据点数 N f f t N_{fft} Nfft​无关，只与原始数据的时长 T T T有关。在时域上补零等价于频域上进行插值，由此增加了频率的点数，使得频谱曲线变得更光滑，即增加了FFT频率分辨率。

为了有效区分1MHz和1.05MHz，必须延长原始数据的时长以提高波形分辨率。以相同的采样频率对原始信号进行采样，采集7000个数据点。得到信号的波形及频谱如图8所示。 有图8可见，1MHz和1.05MHz可以被区分开，但是也出现了频谱泄露现象。此时信号的波形分辨率为： △ R w = 1 70 u s ≈ 14 K H z △R_{w}=\frac{1}{70us}≈14KHz △Rw​=70us1​≈14KHz，小于1MHz和1.05MHz之间的距离 50 K H z 50KHz 50KHz。

为了减小频谱泄露，对原始信号取更长的数据点，采集8000个数据点。得到信号的波形及频谱如图9所示。 由图9可知，1MHz和1.05MHz对应的幅值均为1，且可以有效的区分开，此时信号的FFT分辨率为： △ R f f t = F s 8000 = 12.5 K H z △R_{fft}=\frac{F_{s}}{8000}=12.5KHz △Rfft​=8000Fs​​=12.5KHz，刚好是1MHz和1.05MHz的公约数，即 1 M H z = 12.5 K H z × 80 1MHz=12.5KHz×80 1MHz=12.5KHz×80， 1.05 M H z = 12.5 K H z × 84 1.05MHz=12.5KHz×84 1.05MHz=12.5KHz×84。 所以取合适长度的数据点可以减轻频谱泄露。

周期序列的DFS的系数 X ( k ) X(k) X(k)与 x ( n ) x(n) x(n)的一个周期的 Z Z Z变换在单位圆的 N N N个均匀点上的抽样值相等，这就是频域采样。

频域采样定理： 当 N ≥ L N≥L N≥L，即DFT计算的频域采样点数大于等于信号的长度时，频域采样不会造成时域混叠。 有限长序列 x ( n ) x(n) x(n)的 Z Z Z变换为： X ( Z ) = ∑ n = 0 N − 1 x ( n ) Z − n = ∑ n = 0 N − 1 [ 1 N ∑ k = 0 N − 1 X ( k ) W N − k n ] Z − n = 1 − Z − N N ∑ k = 0 N − 1 X ( k ) 1 − W N − k Z − 1 X(Z)=\sum_{n=0}^{N-1}x(n)Z^{-n}=\sum_{n=0}^{N-1}[\frac{1}{N}\sum_{k=0}^{N-1}X(k)W_{N}^{-kn}]Z^{-n}=\frac{1-Z^{-N}}{N}\sum_{k=0}^{N-1}\frac{X(k)}{1-W_{N}^{-k}Z^{-1}} X(Z)=n=0∑N−1​x(n)Z−n=n=0∑N−1​[N1​k=0∑N−1​X(k)WN−kn​]Z−n=N1−Z−N​k=0∑N−1​1−WN−k​Z−1X(k)​ 其中， W N − k n = e j 2 π k n N W_{N}^{-kn}=e^{j\frac{2πkn}{N}} WN−kn​=ejN2πkn​， W N − k = e j 2 π k N W_{N}^{-k}=e^{j\frac{2πk}{N}} WN−k​=ejN2πk​。

在进行DFT计算时需要对信号的频域进行采样，由于采样间隔为 △ w = 2 π N △w=\frac{2π}{N} △w=N2π​，得到的频谱图都是由一根根离散的谱线组成，就像透过栅栏观看外景。

增加频域采样点数N（不改变时域数据的情况下，在时域数据末端添加一些零值点，使得谱线更密），可缩小谱线间距，减轻栅栏效应。

[1] 数字信号处理 [2] 傅里叶变换的波形分辨率与频率分辨率 [3] 补零、频谱泄露、栅栏效应的关系？