我们知道Fourier分析是信号处理里很重要的技术,matlab提供了强大的信号处理能力,但是有一些细节部分需要我们注意。
记信号f(t)的起始时间为t_start, 终止时间为t_end, 采样周期为t_s, 可以计算信号的持续时间Duration为 t_end – t_start, 信号离散化造成的采样点数 N = Duration/t_s + 1;
根据Fourier分析的相关结论,我们知道时域的采样将会造成频域的周期化,该周期为采样频率f_s(著名的香农采样定理基于此).
于是, 经过matlab的fft函数处理后,得到数据的横坐标为0:f_s/(N-1):f_s。相关代码如下所示:
%matlab fft 测试代码
t_s = 0.01;
t_start = 0.5; t_end = 5;
t = t_start:t_s:t_end;
y = 0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);
y_f = fft(y);subplot(3,1,1);
plot(t,y); title('original signal');Duration = t_end - t_start;
Sampling_points = Duration/t_s + 1;f_s = 1/t_s;
f_x = 0:f_s/(Sampling_points-1):f_s;subplot(3,1,2);
plot(f_x,abs(y_f)); title('fft transform');subplot(3,1,3);
plot(f_x-f_s/2,abs(fftshift(y_f))); title('shift fft transform');也就是说,如果我们不使用fftshift,其变换后的横坐标为0:f_s/(N-1):f_s, 如果使用fftshift命令,0频率分量将会移到坐标中心,这也正是matlab中帮助中心给出的意思:对fft的坐标进行了处理。实际上由于频谱的周期性,我们这样做是合理的,可以接受的。
请读者特别要注意横坐标的差别。另外,根据函数的特性,频谱应当只有在15Hz,40Hz出现峰值,但是fft变换后在60Hz,及85Hz处同样出现了峰值,应当可以从fft的计算过程中得到相应的解释。
事实上,如果我们用15Hz,60Hz来测试fft变换,也即是 y = 0.5*sin(2*pi*15*t)+2*sin(2*pi*60*t);图像如下所示,没有任何变化。
这种现象提醒我们,频率在f_s以内,即 0<f<f_s,f 以及 f_s – f 都有可能是测试信号的频率谱,这就给我们带来了歧义。并且从第三个子图也可以看出,这时候的fftshift会给我们带来错误的引导,也就是说,如果我们试图采样fft或者fftshift来分析信号的频率谱显得不那么靠谱了,matlab的fft谱线与信号的实际频率并不是一一对应的映射关系。这当然不是我们所期望看到的结果,所以实际分析信号时,有关这个问题需要额外的注意。
实际上,这也就间接地证明了Nyquist采样定理的合理性:采样频率要高于截止频率的两倍,上面的处理中我们所使用的采样频率为100Hz,于是当截止频率超过50Hz时,就会出现混叠效应,特殊情况就如上图所示:完全一样。于是,这也就告诉我们若要正确的显示频谱,需要仔细地考量采样频率与截止频率的关系,若太小,则有可能出现混叠,若太大,则计算代价过高。