几个月前写过一篇关于分析用傅里叶变换来计算频率与选取正弦波周期数的关系的博文，今天来讨论一下计算幅度的问题。

不过在此之前还是请让我重述一下下面这两个已经老生常谈的东西（写毕业论文不就是这样么，前面都是一大堆本科生都懂的常识性东西。。）：

一，对于离散序列的连续时间傅里叶变换DTFT，得到的结果是一个连续的，以2\(\pi\)为周期的函数。而FFT，是快速傅里叶变换的英文简称，它所做的，不是离散时间傅里叶变换DTFT，而是离散傅里叶变换DFT，也就是说，它把一个离散序列变换完，得到的是一个离散的序列，那么这个离散的序列和DTFT得到的连续周期函数有什么关系呢？我们可以认为它是后者的一个周期的采样的结果。

二，理论上来说正弦函数sin(ωt)的傅里叶变换结果应该是两个无穷细的冲激函数，但是FFT出来的结果应该是如下图所示，原因吉布斯(振铃)现象在之前那篇博文里面已经解释过了，这里就不累述了！然后，我们看看上图中下面那副频谱图，如果我们把这条曲线认为是DTFT得到的结果，也就是连续函数的那个，那么FFT的结果就是对它的采样！之前那篇博文中讨论的是从这个频谱中确定峰的频率的问题，今天我们讨论的是峰的幅度的问题！！

假设现在我们只对峰值的幅度感兴趣，但是，采样，就不可避免的会出现一个问题，FFT采样的结果不一定会采到峰值那个点，比如说下面这个图就知道了：如果得到的是上面那幅图的蓝线，那么我们频率幅度估计就错了，因为我们得到的是蓝线的最大点的值，顺带一提，这样连频率定位也会出错！！

我们理想的结果应该是下图这个样子的：其实我刚刚提到了，我们只对一个频率感兴趣，所以即便出现下面这个情况，我们也是无所谓的！！我们现在来分析一下产生偏差的来源，我们FFT得到的结果中，第一个点和最后一个点这两个点对应的频率分别是0和采样频率fs，或者也可以认为是-f_s/2到f_s/2，这个在matlab里面就是一个fftshift函数的差别而已，说白了就是：反正离散的结果都是2\(\pi\)为周期的，我们取0~2\(\pi\)还是-\(\pi\)到\(\pi\)都包含了全部信息，对吧！

然后我们假设我们做FFT用的是N点的变换，那么每两个点之间对应的频差就是f_s/N了(关于这个除以N还是除以N-1这个问题我们这里就简单处理吧~别在意了，这样理解起来简单一点)，那么显然的，FFT出来的结果第k个点对应的频率是\(f_s\dfrac{k-1}{N}\)；

那么现在我们可以理解了吧，如果不存在一个整数k，使得\(\dfrac{f_sk}{N}=f_c\)，其中f_c就是我们正弦波的频率的话，那么我们必然就得不到我们的正确的峰值的频率f_c了，那么也得不到f_c这个频率的幅度了，因为采样的时候根本没采到这个点！！

另外，如果你的信号是个宽带谱（说白了就是峰值那里很平坦的那种）的还好，错开一点点峰值的位置对幅度计算影响不大嘛，但是像一个正弦函数那样的信号，峰就是一个很陡的sinc函数，那么错开一点点对幅度影响是很大的！！

做个实验吧：有一个10k的正弦信号，然后做1024个点的FFT，然后采样频率是102.4k Hz，采样频率不仅满足奈奎斯特定理，而且还满足\(\dfrac{f_sk}{N}=f_c\)，这里k=100:

function FFTtest()
 test(102.4*1e3,1024)
function test(fs,N)
 t = 0:1/fs:(N-1)/fs;
 x = sin(2*pi*10000*t);
 y = fftshift(abs(fft(x,N)))/N*2;
 w = linspace(-fs/2,fs/2,length(y));
idx = find(w>=0);
 w = w(idx);
 y = y(idx);
plot(w,y,'b')
 string = ['fc = 10000 Hz    fs = ',num2str(fs),' Hz    N=',num2str(N)];
 title(string,'fontsize',24)
 xlim([0,fs/2])

幅度为1，结果完全没错。

但是我们换一个例子，变换点数N还是1024，但是把采样频率变成100k Hz，也就是不存在整数k满足：\(\dfrac{f_sk}{N}=f_c\)，这时候的结果就变成了：看到了吧，瞬间就变成了原来的80%以下了！！

就算我们狠狠地提高采样频率，现在换成2.3M Hz，我们发现峰值幅度还是误差很大。现在我们明白了，想要准确获得一个频率的幅度，不仅要考虑到这个频率大小f_c，还要以此设计适合的采样频率f_s和FFT变换点数N。绝对不要以为一味的提高采样频率就可以得到好的结果！！上图就是一个例子！！

我这里画了一个对于f_c=10k Hz的信号在不同采样频率下做1024点FFT的结果的幅度变化曲线，可以发现只有一些比较特定的值才可以让计算出来的幅度很接近于1，而取得不好的话，得到的幅值可能只有理论的65%！！我们前面说了，要存在整数k使得下式成立的采样频率f_s和FFT点数N可以得到较好的幅度计算：

\(k=\dfrac{f_c\times N}{f_s}\)

为了佐证我上面的说法，看下图，蓝线和上图意义相同，红线画的是上式k的小数部分减去0.5的绝对值，也就是说这个红线越小，说明这个k值小数部分越接近0.5，而越大，就越接近0或者1，也就是说k离整数越近，就是我们想要的理想的结果；可以看出，确实当k的小数部分见0.5越大的时候，得到的幅度更加接近于1，也就是更加的精确！！但是也不完全正确，因为我们只说了k为整数的时候效果很好，但是没说k小数部分是0.1的时候比0.2的时候好，但是基本可以这么认为了！！
其实上面讨论的问题，本质上还是一个频谱泄露问题，只是我换了一种解释方法而已。

不信的话，可以简易证明如下：

FFT能准确计算频率的假设是：将截断信号无限延拓之后可以得到原来的信号。在这个例子中，原来的信号指的就是无限长的一个正弦函数，而截断信号就是我们编程中一直用的那个部分正弦信号，也就是用矩形窗截断无限长的正弦波得到的东西。

那么什么情况下上面那个无限延拓的假设可以成立呢？我们截断的信号t是从0开始，到\(t=\dfrac{N}{f_s}\)为止的【其实是\(\dfrac{N}{f_s}\)，但是还是别在意,为了方便说明而已】，起始时刻t=0的时候相位是0，那么终止时刻\(t=\dfrac{N}{f_s}\)的时候相位也要是0，那么无限延拓的话才不会出现相位突变，这样就可以完美的重现正弦波了！！

\(sin(2\pi f_c t)\)在终止时刻\(t=\dfrac{N}{f_s}\)的相位是：\(2\pi f_c \dfrac{N}{f_s}\)，这个相位应该等于2k\(\pi\)，化简一下就会得到之前的式子:

\(k=\dfrac{f_c\times N}{f_s}\)

Q.E.D

数字信号处理里面我们学过，加窗，可以弱化频谱泄露的影响，所以我们可以预知到，加窗，也可以增加频谱幅度的计算准确度！！

作为实验，我们设置这个窗函数：

\(1-1.93cos(\dfrac{2\pi n}{N})+1.29cos(\dfrac{4\pi n}{N})-0.388cos(\dfrac{6\pi n}{N})+0.0322cos(\dfrac{8\pi n}{N})\)

这个窗函数的形状是这个样子的，江湖人称——平坦窗！
我们用这个窗来实践一下上面那个100k Hz采样的情况，和之前的区别就在于加了个窗，代码是：

function FFTtest()
test(2.3*1e6,1024)

function test2(fs,N)
t = 0:1/fs:(N-1)/fs;
x = sin(2*pi*10000*t);
y = fftshift(abs(fft(x,N)))/N*2;
w = linspace(-fs/2,fs/2,length(y));

%窗函数
n = 0:N-1;
windows = 1-1.93*cos(2*pi*n/N)+1.29*cos(4*pi*n/N)-0.388*cos(6*pi*n/N)+0.0322*cos(8*pi*n/N);
%加窗后的频谱
y_window = fftshift(abs(fft(x.*windows,N)))/N*2;

%取正频率部分
idx = find(w>=0);
w = w(idx);
y = y(idx);
y_window = y_window(idx);

hold on;
plot(w,y_window,'r','linewid',2)
plot(w,y,'b','linewid',2)
hold off

加窗效果前后对比如下，可以看到频谱幅度的估计马上就被矫正了！！
我们还可以实践一下上面对102.4k Hz采样的结果，发现同样保持住了水准！！再对上面那个2.3M Hz的实验进行加窗，效果也是很显著的！！
我们可以按照上面的方法，再画一次加窗后，不同采样频率情况下对幅度计算的影响，可以发现，呵呵，不管采样频率怎么变(当然，至少奈奎斯特条件还是要满足的)，幅度结果基本完全不变了！！
所以说，通过加窗可以很好地稳定住因为采样频率或者说频谱泄露引起的误差。

但是这里说一句，虽然有点废话，加窗不是随便加什么窗都行的，怎么选择加窗，窗的参数怎么取，这个暂时超出了本文的研究范畴了，有空我再写一篇讲这个的吧。。。

好吧，今天暂时讲到这。。

【完】

本文内容遵从CC版权协议，转载请注明出自http://www.kylen314.com