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

什么是频谱泄露？

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

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

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

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

$$x_r(t)=Acos({\omega }_{0}n)w(n)$$
那么根据频域卷积定理我们得到：

如果该窗是一个矩形窗的话,即
$$w(n)=\{\begin{array}{rl}& =1,-L/2\le n\le L/2-1\\ & =0,others\end{array}$$
它的离散时间傅里叶变换(DTFT)为
$$W(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}$$
下面以 N = 16为例给出其矩形窗函数幅度谱:

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

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

包含2个完整的周期。由对 $cos(10\times 2\pi t){|}_{t=\frac{n}{f_s}}$ 采样得到 $x[n]=cos(\frac{20\pi }{100}n)$ ,那么我们得到:
$$X_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}$$
接下来就是对 $X_r(e^{j\omega })$ 进行N=64点采样得到DFT了，但是在这之前思考一个问题，我们在采样点$\omega =\frac{2\pi k}{N}$ 处得到的频率点对应的模拟角频率是多少？
这个问题可以根据模拟角频率和数字角频率对应的关系得到:$\Omega =\frac{\omega }{T_s}=\omega f_s$ ,那么模拟角频率对应DFT采样点的频率为： $\Omega =\frac{2\pi k}{N}fs,-\frac{N}{2}\le k\le \frac{N}{2}-1$ ，至此得到模拟域对应点处的频率。
继续进行我们的DFT采样，我们得到：
$$X_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}$$
采样图如下：

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

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

可以消除频谱泄露吗？

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

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

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

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

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

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

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

X = fft(x)/SampleNum;

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