其实这是比较久之前就抱有一点疑问的问题了，除去项目各种相关知识后的问题起源描述是这样子的，有一个正弦波，它只有2~5个周期左右，然后请确定出它的频率来。

不过有两点需要说明一下：

第一点，在实际项目中采回来的数据，不可能是完美的正弦信号，有可能有相位噪声或者幅度噪声什么的，因为我们微波光子学研究的频率是GHz~THz这种数量级的，所以畸变是必然的，但是！！！！我们这里不涉及这方面，讨论的就是完美的正弦信号。

第二点，实际上当一个信号只有2~5个周期这种时候，是不应该用傅里叶变换的，因为傅里叶对瞬时信号或者短信号相对之下是比较”无能”的，但是我们还是不管，我们这里分析的就是福利叶变换对完美的短周期数正弦信号的频率定位问题。

首先，我们知道一个正弦信号的傅里叶变换应该如下图所示的，下图是一个频率为100Hz的正弦波在1kHz的采样频率下得到的数字信号及其傅里叶变换。
必然如预期所想，在100Hz出有个峰，但是同时还有很多小的旁瓣，这个原因很简单，我们知道正弦信号的傅里叶变换是个冲激函数，只在100Hz出有值，但是这是有个前提的，就是这个正弦信号是无限长，我们现在这个只有10个周期的正弦信号可以看成一个无限长的正弦信号在时域上乘上了一个窗，然后就相当于只在10个周期有值，其余的值都是0，时域上的相乘作用到频域上是卷积，也就是说十个周期的正弦波的傅里叶变换（频谱）相当于无限长的正弦波的傅里叶变换(冲激函数)卷积上窗函数的福利叶变换(sinc函数)，而任何函数卷积冲激函数相当于把那个函数频率为0的位置搬移到冲积函数的位置上，所以我们看到的10个周期的正弦波的傅里叶变换才会像上图所说，冲激的函数在100Hz处变成了sinc函数。【这也就是吉布斯效应相关的现象，或者叫振铃效应(ringing effect)】。

不过扯远了，我们现在假设有10个周期的正弦波，然后要怎么得到它的频率呢？在已知道采样频率的情况下，就可以通过福利叶变换得到频谱，然后找峰值就可以知道频率了。

但是，问题在于，虽然数字信号的频谱在理论上是连续的，以2\(\pi\)为周期的，但是我们在计算机中使用的是DFT，或者FFT，得到的频谱是数字的，说白了，就是对那个连续的X(jω)的采样，这样我们通过找采样后的数字信号的峰值就不一定可以准确的找到原本连续的DTFT的峰值，甚至可以说是不可能的，所以这样就产生了频率计算的误差了。扯了这么多，我这里就是分析一下这个误差相关的乱七八糟的东西的。。。。【哎哟，总算入正题了。。】

之前之所以会想到这个问题就在于，我们仪器采到的信号大概是3.8个周期左右，老板跟我说，让我做傅里叶变换的时候先要找出整数个周期，比如3个，然后再对那一段做傅里叶变换，然后当时我的直觉是，怎么会？难道3个会比3.8个要好？老板凭他多年的经验告诉我，是的！！但是没解释。

我的想法很简单，既然我们的信号可以看成无限长的正弦波乘窗，那么取得周期越多，就代表窗宽越大，时域越宽的信号，频域必然越窄，所以那个sinc函数就会越窄，那么自然就越接近冲激函数啊，那么也就是越精确，所以怎么回事整数周期要好呢？我做了一下简单的实验，如下图所示，显然和理论分析一样，3.8个周期的峰要比3个的窄。
但是，这里要说的是，我完全忽略了一件事情，我们所感兴趣的并不是波形是否接近于完美的冲激函数的形状，而是，频率定位准不准确，也就是找到峰值对应的频率是不是等于理想的波形的频率，这个就不一定和形状相关了。

于是我做了另一个实验，就是先构造一个频率为 f1 的数字正弦信号，然后截取不同的周期数，1到10个周期，步进为0.005个周期，然后计算它的傅里叶变换，找到峰值对应的频率f2，这个就是我们算出来的频率，然后计算算出来的频率的相对误差|(f1-f2)/f1|，画出相对误差关于周期数的图，结果如下图所示：

可以看出来，额。。我老板说的没错，确实在整数点的时候取得误差的极小值，虽然在2.5 , 3.5 , 4.5这些也可以取得极小值，但是这个也完全可以理解，也就是周期数取N\(\pi\)的时候，可以在频率分析的结果中取得较小误差。

而且和一开始的分析也相符，也就是周期数越多，误差确实会逐渐减小的，只是在周期数变化不大的情况下(小于1)，取0.5k个周期算出来的频点是比较准的。虽然其他的算出来相对误差也差不多1‰这个数量级，不是很大，但是对于高精度的分析这也是不允许的！！

虽然一时理解上没那么直观，但是我突然想到，我刚刚做的傅里叶变换的程序用的是fft(y)这种形式，matlab中这么写就是说傅里叶变换后的点数和你时序序列的点数相同，这样就出现了一个问题，我周期数取得越多，那么点数就越多啊，这样傅里叶变换后的点数就越多，傅里叶变换的点数就是对DTFT频域[-\(\pi\)，\(\pi\))这段范围的采样数，采样数越多自然采样频率越大，这样不同周期数FFT出来的结果实际上不是在同一水平下衡量的，于是，我算了一下统一傅里叶变换长度的实验，结果如下：

单看趋势的话，还真是比较符合预期所想，就是有损震荡的那种形状，但是仔细看看，我们这次就发现了，在整数周期的时候，误差取得是局部最大值，这和不固定傅里叶长度时的情况完全相反！！现象先放在这里，我写完这篇再做做推导，虽然推出来也不一定贴上来，因为敲公式麻烦，就算贴上来估计也是贴图片。

下图是放在一起的两幅图，可以看得出来相对误差相差还是有点大的，当然也是因为傅里叶变换点数取得有点多的缘故。

之后就继续做了一下别的研究；

一个是固定傅里叶变换后，点数选取造成的差异，对于大小我倒是完全不感兴趣，但是对于奇偶我倒想验证一个想法，因为我们知道傅里叶变换完-\(\pi\)和\(\pi\)分别代表的是-f_s/2和f_s/2，但是当我们取奇数个点，那么就可以取出后一半(N+1)/2个点来表征0到f_s/2，但是如果取得是偶数个点，这就办不到了，因为这种时候没有点表征的是0这个频率，然后我试了一下计算相对误差：

一定程度上如预计的一样，也就是说底限，奇数个点的情况会出现更加接近0误差的情况，而偶数的最小值一直在奇数情况以上一些。(如果图看不清，可以点击看大图。)

然后发现，为什么总是取不到误差为0的那种情况呢？然后我调大了原本正弦信号的频率：

继续，还有一些参数可以调，比如。。采样频率。。

Code:

f = 100;
fs = 1000;
N_range = 1: 0.005 : 10;
error_m(length(N_range)) = 0;
for i = 1:length(N_range)
 [x y] = GetCos(f*2*pi,N_range(i),fs);
 w = GetFFT_w(y,fs,1024*16+1);
 error_m(i) = abs(w-f)/f;
end
plot(N_range,error_m,‘b’,‘linewidth’,2);
xlim([N_range(1) N_range(end)]);
set(gcf,‘color’,‘w’)
hold on;

function [x y] = GetCos(w,PeriodNum,fs)
x = 0 : 1/fs : 2*pi/w*PeriodNum;
y = cos(w*x);

function w = GetFFT_w(y,fs,Point)
if nargin == 2
 Y = abs(fftshift(fft(y)));
else
 Y = abs(fftshift(fft(y,Point)));
end
w = linspace(-fs/2,fs/2,length(Y));
plot(w,Y)
[m idx] = max(Y);
w = abs(w(idx(1)));

【完】

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