欢迎访问个人博客：blog.spursgo.com

matlab中fft函数易忽略的地方

数字信号处理中快速傅里叶算法可以说是非常伟大的一个算法，当然也就免不了被matlab青睐了。

数字信号处理实验课中多次使用fft这个函数，之前只是简单的调用它，直到这次做数字信号处理大作业的时候，再次涉及到fft函数，于是就去看了下matlab的官方文档，却发现了一个值得探究的问题。

下面是matlab官方文档中给出的一个例子：

使用傅立叶变换来查找被噪声掩埋的信号的频率分量。
指定采样频率为1赫的信号的参数, 信号持续时间为1.5 秒。

Fs = 1000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector

形成一个信号, 其中包含振幅为0.7 的50赫兹正弦波和振幅1的120赫兹正弦波。

S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

用零平均白噪声破坏信号, 方差为4。

X = S + 2*randn(size(t));

在时间域中绘制噪声信号。通过查看信号X(t) , 很难识别频率分量。

plot(1000*t(1:50),X(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('t (milliseconds)')
ylabel('X(t)')

时域信号

计算信号的傅里叶变换。

Y = fft(X);

计算双侧频谱P2。然后计算基于P2和偶数值信号长度L的单面频谱P1.

P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

定义频率域f并绘制单面振幅频谱P1。因为增加的噪声, 振幅不是确切地在0.7 和 1, 如期望。平均来说, 较长的信号产生更好的频率近似。

f = Fs*(0:(L/2))/L;
plot(f,P1) 
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')

频域信号
注：L是信号的长度，官方文档没有给出。

问题来了

实验课上，我每次都是使用P2 = abs(Y);这种方式得到信号的幅度，可是官方文档却给出了这样的表述P2 = abs(Y/L);。

当时非常疑惑，想了很久没有想明白，在网上没有得到理想的答案。于是，和同学学长学姐讨论了这个问题。最后得到的解释还是不能让我满意。

这个问题并没有被搁置。结合和同学学长学姐讨论的，然后结合课本，找到了答案。
注：fft是DFT的快速算法，我们这里通过探究DFT来解释上述现象。

书本上是这样定义DFT的：

教材定义

听课的时候，没有去想为什么逆变换的时候一定会除以N，而正变换里却没有。而官方文档中的L和教材里扽N正是同一个东西。
我们将教材中的正变换记为(1)式，逆变换记为(2)式。

观察逆变换：变换的对象是X(k)，它和正变换的形式非常相似，最大的区别是逆变换还除以了N，也就是说二者的形式是对称的，这里我们做一个假设，如果我们去掉逆变换中的除以N，将此时的逆变换记为(3)式，从形式上将二者变成了真正意义上的对称。但是，此时的逆变换的结果却不是我们想要的了，x(n)已经变成了另外的序列。

接下来，我们把除以N放到(1)式的正变换中，将此时的正变换记为(4)式，再来看一下会有什么效果：X(k)已经变化，我们将此时的X(k)带到(3)式的逆变换中，对比(2)式，二者形式相同。解释一下前面的做法：我们将除以N拿到正变换中进行，仍然能够得到正确的结果，说明一个什么问题呢？

从x(n)到X(k)，定义式的这种处理并没有得到真正的频谱值，相比于真实值，它被扩大了N倍，而逆变换中的除以N正是将之前放大的效果抵消，这样才能得到原序列。

那么就可以解释为什么我们平时我们做实验的时候用的是P2 = abs(Y)，而官方文档却用P2 = abs(Y/L)。

因为我们平时并没有在意频谱图像中纵轴的值具体是多少，我们的关注点在它们的相对大小和变化趋势，因此N倍的处理并不重要。但是matlab官方文档力求严谨，将真实的频谱展现给我们。

虽然这个N到底除不除并不重要，但是这个探究过程对我们理解DFT还是很有帮助的。