什么是频谱泄露？

一些人可能会说，频谱泄露就是频谱中除了本来该有的主瓣之外，还会出现本不该有的旁瓣；一下子主瓣、旁瓣这些概念都来了，这些一些很抽象的概念，搞得我们很难去理解。其实频谱泄露就一句话：原来信号的频谱中出现了本不该有的频率分量！

那频谱泄露的原因是什么？

原来的信号的序列我们认为是无限长的，但我们要分析某个信号的频谱的时候只能对它有限长度的序列进行分析，这么一来，就相当于时域上对它进行了加窗，但这并不能代表我们的信号啊，于是就导致了频谱泄露。
归根结底，频谱泄露的原因就是加窗导致的序列长度有限。
下来我们看加窗是如何导致信号的频谱发生泄露的？

时域加窗对正弦信号的影响

xr(t)=Acos(ω0n)w(n)x_r(t)=Acos(\omega _0n) w(n) xr(t)=Acos(ω0n)w(n)
那么根据频域卷积定理我们得到：

如果该窗是一个矩形窗的话,即
w(n)={=1,−L/2≤n≤L/2−1=0,othersw(n)=\left\{ \begin{aligned} & = 1 ,-L/2\leq n\leq L/2-1\\ & = 0,others \\ \end{aligned} \right. w(n)={=1,−L/2≤n≤L/2−1=0,others
它的离散时间傅里叶变换(DTFT)为
W(jΩ)=sin(ΩT/2)Ω/2e−jΩT/2W(ejω)=ejω(1/2)sinωL2sinω/2W(j\Omega) = \frac{sin(\Omega T/2)}{\Omega /2}e^{-j\Omega T/2} W(e^{j\omega }) = e^{j\omega(1/2) } \frac{sin\frac{ \omega L}{2}}{sin\omega/2} W(jΩ)=Ω/2sin(ΩT/2)e−jΩT/2W(ejω)=ejω(1/2)sinω/2sin2ωL
下面以 N = 16为例给出其矩形窗函数幅度谱:

如果直接调用工具包，会给出以分贝(dB)表示的频域特性：

接下来我们在做的是对Xr(ejω)X_r(e^{j\omega})Xr(ejω)进行的的频域抽样(即离散傅里叶变换-DFT)，而不是 X(ejω)X(e^{j\omega})X(ejω)的频域抽样。对Xr(ejω)X_r(e^{j\omega})Xr(ejω)的N点频域抽样(这里的N和时域加窗的L没有必然的联系，事实上只要N>=L,不发生时域混叠，都是可以的）为：
Xr(k)=A2W(e(j2πkN−ω0))+A2W(e(j2πkN+ω0))X_r(k)=\frac{A}{2}W(e^{(j\frac{2\pi k}{N}-\omega_0)})+\frac{A}{2}W(e^{(j\frac{2\pi k}{N}+\omega_0)}) Xr(k)=2AW(e(jN2πk−ω0))+2AW(e(jN2πk+ω0))
可能你会觉得到了这里不是很理解我上面所说的内容，没关系，下来我们给出相应的一个例子：
假设矩形窗的长度为L=20，余弦信号频率为10Hz，采样频率设置为fs = 100Hz，采样点数N=64。即对 cos(10×2πt)∣t=nTs=cos(10×2πt)∣t=nfscos(10\times 2\pi t)|_{t = nT_s}=cos(10\times 2\pi t)|_{t = \frac{n}{f_s}}cos(10×2πt)∣t=nTs=cos(10×2πt)∣t=fsn 进行采样，且采样点总长度L = 20。时域采样图如下：

包含2个完整的周期。由对 cos(10×2πt)∣t=nfscos(10\times 2\pi t)|_{t = \frac{n}{f_s}}cos(10×2πt)∣t=fsn 采样得到 x[n]=cos(20π100n)x[n] = cos(\frac{20 \pi}{100}n)x[n]=cos(10020πn) ,那么我们得到:
Xr(ejω)∣ω=2πk/N=12W(ej(ω−20π100))+12W(ej(ω+20π100))∣ω=2πk/NX_r(e^{j\omega})|_{\omega = 2\pi k/N} =\frac{1}{2}W(e^{j(\omega-\frac{20 \pi}{100})})+\frac{1}{2}W(e^{j(\omega+\frac{20 \pi}{100})}) |_{\omega = 2\pi k/N} Xr(ejω)∣ω=2πk/N=21W(ej(ω−10020π))+21W(ej(ω+10020π))∣ω=2πk/N
接下来就是对 Xr(ejω)X_r(e^{j\omega})Xr(ejω) 进行N=64点采样得到DFT了，但是在这之前思考一个问题，我们在采样点ω=2πkN\omega = \frac{ 2 \pi k}{N}ω=N2πk 处得到的频率点对应的模拟角频率是多少？
这个问题可以根据模拟角频率和数字角频率对应的关系得到:Ω=ωTs=ωfs\Omega = \frac{\omega}{T_s}=\omega f_sΩ=Tsω=ωfs ,那么模拟角频率对应DFT采样点的频率为： Ω=2πkNfs,−N2≤k≤N2−1\Omega =\frac{ 2 \pi k}{N}fs,-\frac{N}{2}\leq k\leq\frac{N}{2}-1Ω=N2πkfs,−2N≤k≤2N−1 ，至此得到模拟域对应点处的频率。
继续进行我们的DFT采样，我们得到：
Xr(ejω)∣ω=2πk/N=12W(ej(ω−20π100))+12W(ej(ω+20π100))∣ω=2πk/NX_r(e^{j\omega})|_{\omega = 2\pi k/N} =\frac{1}{2}W(e^{j(\omega-\frac{20 \pi}{100})})+\frac{1}{2}W(e^{j(\omega+\frac{20 \pi}{100})}) |_{\omega = 2\pi k/N} Xr(ejω)∣ω=2πk/N=21W(ej(ω−10020π))+21W(ej(ω+10020π))∣ω=2πk/N
采样图如下：

我们从中可以看到，DFT采样点出现了很多的频率分量，而不只是出现在了单根谱线10Hz所对应的数字角频率 ω=20π/100=π/5\omega =20\pi/100= \pi/5ω=20π/100=π/5 处，这就是频谱泄露：在除去±π/5\pm \pi/5±π/5 的频率分量处都出现了不该有的频率分量。为了更好理解，我们还原到模拟角频率看一下它对应的频谱图：

我们可以看到，它的频谱图像不光出现在了10Hz处，其他频率分量处对应幅值不是0，这就是频谱泄露。而且我们可以看到，它在10Hz的频率分量处不是最大值，最大值出现在了它附近，这是由于刚才DFT对它进行采样没有采到 π/5\pi/5π/5 处，而是采到了它的附近。
现在是不是感觉清楚了一点，我们再来回看频谱泄露的原因：因为时域的有限长度(加窗)导致了我们是在对窗函数平移后的频谱进行的采样，可能会采样到“主瓣”“旁瓣”等等，取决于你的DFT采样点数N。

可以消除频谱泄露吗？

既然频谱泄露是不好的，那我们就想要消除它。但问题是频谱泄露能够消除吗？
答：不能。因为时域样点的长度L(即窗长度)始终是有限的，如果我们进行DFT采样，那么始终就会采到“主瓣”“旁瓣”，所以频谱泄露是不能消除的。但有一种情况例外：
还是上面那个例子，其余参数保持不变，但我们对 Xr(ejω)X_r(e^{j\omega})Xr(ejω) 进行N=20点采样，我们能得到什么呢？

采样点只有在主瓣峰值 π/5\pi/5π/5 处不为0，其余采样点都是0，这时候恰好采样点把0点采样到了，就相当于把其余导致频谱泄露的点漏掉了，所以这时候可以认为我们得到的频谱就真真切切是 cos(10×2πt)cos(10\times 2\pi t)cos(10×2πt) 的频谱，对应的模拟域频谱图如下：

但这种情况概率极小，我们一般在做fft的时候是不会专门计算这些东西的，所以一般情况下频谱泄露都会存在。但你可能会问，为什么我平时做频谱，看不到频谱泄露现象呢？
其实是有的。我们平时时域的采样点取多少，L = 30000或者更高，也就意味着我们时域的窗的长度是L =30000,而且我们在使用fft的函数的时候只传入了x[n],并没有传入第二个参数N，即DFT的采样点数，而我们说N和L没有必然的联系，只要求N>L,matlab默认给我们选定了。而30000的窗函数的离散傅里叶变换(DTFT)是一个什么概念？
这是L=20的矩形窗的 W(ejω)W(e^{j\omega})W(ejω)

L=100的矩形窗的 W(ejω)W(e^{j\omega})W(ejω)

L=1000的矩形窗的 W(ejω)W(e^{j\omega})W(ejω)

接着往后，如果矩形窗的长度到了L =30000，那么它的主瓣峰值是30000，旁瓣值要更小，所以我们采样的时候一般采取的采样点很大，让它的频率间距很小，看起来像是我们的频谱。但实际上旁瓣仍然存在，这一点可以从L=1000长的频谱局部中看出：

旁瓣在w=0附近还是非常大的，只是因为它的范围太小导致映射回模拟域的时候都在中心频率附近导致我们认为它是一根谱线罢了！我们做频谱的时候有时候会很疑惑，明明余弦函数的频谱是一个冲激，无穷大的，为什么fft是一个有限值呢？其实就是我们矩形窗长度的限制：当矩形窗的长度L越大，时域的采样点就越多，那么它的DTFT也就越大，那采样得到的DFT也就越大，对应的幅值越也就大。如果时域采样点无穷多，自然而然就得到了一个冲激了；我们也经常能看到有些程序中有下面的形式：

X = fft(x)/SampleNum;

这里除以了时域的采样点数，就是因为它的峰值是在和我们的时域的采样点数是有关系的，所以在这里要进行归一化。
通过以上讨论，我们发现“一般情况下”频谱泄露都是会存在的。
但是怎么样解决呢？换窗，改变它的主瓣旁瓣的相对大小。
换窗的本质是一种“加权”，
未完待续…

MATLAB中的FFT函数以及频谱泄露相关推荐

MATLAB中的常用函数小结
1. MATLAB中的常用函数小结文章目录 1. MATLAB中的常用函数小结 1. MATLAB图像处理工具箱 1.1 图像显示 1.2 图像文件输入/输出 1.3. 图像像素值及其统计 1.4 ...
利用卷积定理进行信道估计 + 深入探究 DFT 与 matlab 中的 fft 用法
文章目录 1 背景 2 DFT & FFT 公式级别解析 2.1 DFT / IDFT 2.2 FFT 3 利用卷积定理进行信道估计 1 背景最近尝试用时域卷积定理来进行信道估计假设 TX ...
matlab ifft频率分辨率,[FFT] matlab中关于FFT的使用(理解频率分辨率、补零问题)
[FFT] matlab中关于FFT的使用(理解频率分辨率.补零问题).txt我这人从不记仇,一般有仇当场我就报了.没什么事不要找我,有事更不用找我!就算是believe中间也藏了一个lie! 我那 ...
fftw3/gsl/kissfft/OouraFFT库中傅里叶变换/反傅里叶变换函数和Matlab中的fft/ifft的对应关系
先分析一维度的一.fftw_plan_dft_1d 正变换: fftw_complex *in = fftw_malloc ( sizeof ( fftw_complex ) * n ); fftw ...
Matlab中常见实用函数（敲代码碰到的）
目录 1.norm函数 2.varargin函数(varargout) 3.nargout函数(nargin) 4.ndgrid函数 5.ndims函数 6.surface函数 7.gcbf函数 8. ...
Matlab中的lsqcurvefit函数的使用
Matlab中的lsqcurvefit函数的使用 lsqcurvefit函数调用示例 lsqcurvefit函数非线性曲线拟合是已知输入向量xdata和输出向量ydata,并且知道输入与输出的函数 ...
Matlab：Matlab中常用的函数、案例详细攻略
Matlab:Matlab中常用的函数.案例详细攻略目录常用函数 1.与文件相关 2.MATLAB GUI不同控件函数间变量传递方法常用函数 Matlab中的bwmorph函数解释 bwmorp ...
matlab作动态函数曲线图,[转载]Matlab中使用Plot函数动态画图方法总结
本帖最后由 sonictl 于 2012-12-31 12:18 编辑请删除我清楚超靠靠靠没办法,一会儿限制这不能发表,那不能发表的.... [转载]Matlab中使用Plot函数动态画图方法总 ...
python实现Matlab中的circshift函数
circshift是Matlab中矩阵循环移位函数,具体使用参照该链接. 但是python中并没有封装好的该函数,因此需要自己实现. 思路:将矩阵分为两部分,然后按照自己的需要堆叠在一起就可以了. n ...
matlab的数学函数,matlab中常见数学函数的使用
matlab中常见数学函数的使用 MATLAB 基本知识 Matlab 的内部常数 pi 圆周率 exp(1) 自然对数的底数 e i 或 j 虚数单位 Inf 或 inf 无穷大 Matlab 的常 ...

MATLAB中的FFT函数以及频谱泄露

什么是频谱泄露？

那频谱泄露的原因是什么？

时域加窗对正弦信号的影响

可以消除频谱泄露吗？

MATLAB中的FFT函数以及频谱泄露相关推荐

最新文章

热门文章