MATLAB中fft函数用法、性质、特性、缺陷全面深入解析(含程序)

MATLAB中fft函数用法、性质、特性、缺陷全面深入解析（含程序）前言目前网络上在信号处理方面对fft函数使用泛滥，但大多数人根本就不了解它的性能，使用方法也是千奇百怪，各种程序代码鱼龙混杂，乱七八糟，误导了很多人，所以我在这里进行一下全面的演示说明，希望后来人能够少走弯路。fft信号处理中获得时域信号的频域的方法叫做傅里叶变换，但教材上的变换都是理论推导的理想、连续、无限情况，在实际的matlab计算中，必须是实际、离散、有限的计算，所以实现起来比较麻烦，需要各种转换。fft函数就是为了实现快速傅里叶变换而设计的计算机算法，与传统傅里叶变换的步骤不完全一样，但通过适当的变换也能够获得时域信号的频域图形。fft只是数值计算，不是理论推导，所以频域图形肯定有误差，但我们能够通过适当调整变换参数来提高频域图形的精确性。建议大家能够理论推导的傅里叶变换就直接在纸上推导好了再用matlab直接画出来就是了。所以不是因为fft精确，而是因为fft能够对大数据复杂信号变换计算快，自由度可控制性高，所以使用的人才多。fft(x)图示fft的一些定性问题解析时域信号x有多少个点，fft(x)就有多少个点，这时的fft(x)之所以是对称的，是因为fft的方法是进行周期计算从，所以后面的一半要翻转到负半轴，是负频率。所以有用的是fft(x)的前一半的点。我是经过各项测试实验，才完全明白fft的性能。利用程序解释如下：%这里的fft变换有补零clearall;Fs=1e3;%这是原始时域信号x的采样频率multiple=0.6;%x进行fft前的补零倍数，可以增加频率分辨率time=2;t=0:1/Fs:time-1/Fs;%这个会影响x的采样点数，当然点数越多，频域越精确越平滑%x=1.5*cos(2*pi*100*t)+3*sin(2*pi*202.5*t)+2*sin(2*pi*1250*t);x=1.5*cos(2*pi*100*t)+3*sin(2*pi*201.25*t);%频率最小分辨率————————————————-x=[x,zeros(1,length(x)*multiple)];%{注意可以通过对x进行补零来增加频率分辨率，只是因为fft的周期性计算方法，需要在后面对fft(x)后的数据进行幅度值重新计算时，乘以多的零点的倍数，见后面的幅度值计算。%}%freqres=Fs/length(x);freqres=1/time;%{x的频率必须是freqres的整数倍时，x的频率才能被正确分辨，也就是相应频域值才能够正常，比如freqres=1,在频率202.5处的幅度就不正常，当freqres=2时，在频率1001处的幅度就不正常。%}%——————————————————%{x有多少个点，fft(x)就有多少个点，这时的fft(x)之所以是对称的，-是因为fft的方法是进行周期计算从[-2π,2π]，所以后面的一半要翻转到负半轴，是负频率。所以有用的是fft(x)的前一半的点再加上0点。%}xdft=fft(x);xdft=xdft(1:length(x)/2+1);%————这与fft方法的周期计算方法相关，也就是幅度值要重新计算一下——————xdft=1/length(x).*xdft;xdft(2:end-1)=(1+multiple)*2*xdft(2:end-1);%{幅度值首先由于周期对称性，需要乘以2，然后由于补零的倍数是multiple，所以是(1+multiple)*2%}%重新计算频率横坐标—————————————————————————————freq=0:Fs/length(x):Fs/2;%{Fs/2就是最高能够显示出来的频率，但要注意实际上，fft(x)实际是计算出了任何频率，但是在[Fs/2,Fs]这个区间的频率，是关于Fs/2这里对称的，也就是正负的频率开始混合或者混叠在一起。而超出了Fs的频率则会由于周期性，它的值先减去Fs后，重新从0点开始。所以这个程序并不是说画出了在[0,Fs/2]的一定正确的频率，而是可能会出现频率的混合和混叠，一定要先预判一下频率的区间范围后，再使用fft这个函数。%}%——————————————————————————plot(freq,abs(xdft));%画频域幅度谱title('频域图');010020030040050000.511.522.533.5频域图fft的一些定量详细解析与测试利用程序解释如下：%为了更好的展示，这里信号没有补零clearall;clc;Fs=1e3;%这是原始时域信号x的采样频率multiple=0;%x进行fft前的补零倍数，可以增加频率分辨率，也会影响信号频域值time=10;t=0:1/Fs:time-1/Fs;%这个会影响x的采样点数，当然点数越多，频域越精确越平滑%x=1.5*cos(2*pi*100*t)+3*sin(2*pi*202.5*t)+2*sin(2*pi*1450*t);%x=1.5*cos(2*pi*100*t)+3*sin(2*pi*200*t)+2*cos(2*pi*100000.1*t);x=1.5*cos(2*pi*100*t)+3*sin(2*pi*200*t)+2*cos(2*pi*0.1*t);%频率最小分辨率————————————————-x=[x,zeros(1,length(x)*multiple)];%{注意可以通过对x进行补零来增加频率分辨率，只是因为fft的周期性计算方法，需要在后面对fft(x)后的数据进行幅度值重新计算时，乘以多的零点的倍数，见后面的幅度值计算。%}freqres=Fs/length(x);%{x的频率必须是freqres的整数倍时，x的频率才能被正确分辨，也就是幅度值才能够正常，比如freqres=1,在频率202.5处的幅度就不正常，当freqres=2时，在频率1001处的幅度就不正常。%}%——————————————————%{有多少个点，fft(x)就有多少个点，这时的fft(x)之所以是对称的，-是因为fft的方法是进行周期计算从[-2π,2π]，所以后面的一半要翻转到负半轴，是负频率。所以有用的是fft(x)的前一半的点再加上0点。%}xdft=fft(x);subplot(3,1,1);plot(0:length(xdft)-1,abs(xdft));title('abs');title('频域图');subplot(3,1,2);plot(0:length(xdft)-1,real(xdft));title('real');subplot(3,1,3);plot(0:length(xdft)-1,imag(xdft));title('imag');back=ifft(xdft);error=x-back;fori=2:length(xdft)/2store(i)=abs(xdft(i))-abs(xdft(length(xdft)-i+2));%这时的误差为0,%store(i)=real(xdft(i+1))-real(xdft(length(xdft)-i+1));%这时的误差为0,%store(i)=imag(xdft(i+1))+imag(xdft(length(xdft)-i+1));%这时的误差为0,注意这时是加号endfigure;plot(store);title('正负频率误差');0200040006000800010000012x104频域图0200040006000800010000-202x104real0200040006000800010000-202x104imag010002000300040005000-1-0.8-0.6-0.4-0.200.20.40.60.81正负频率误差首先如上面程序所示，xdft=fft(x)，频率分辨率是freqres，设L=length(xdft)=length(x)，显示的频率范围是，点数范围是。xdft的幅度谱在频率范围，也就是点数范围内，是与右边点数内的频率的实部是完全对称的，虚部是正负对称的，所以幅度谱是完全对称的。xdft的幅度谱在频率范围，也就是点数范围，是与左边边点数内的频率的实部是完全对称的，虚部是正负对称的，所以幅度谱是完全对称的，但是这时候左右的正负频率是混叠的，这就是奈奎斯特准则的限制，采样率要大于信号的最高频率才不会出现混叠，因此只有的频率是准确可用的。然后是缺的几个特征点，在频率为0，也就是第1个点的地方，频率虚部为0，实部为实际振幅值的两倍；在频率为，也就是点数为的地方，频率虚部为0，实部为实际振幅值的两倍；在频率为的地方，也就是超出了频率范围，那么这个点的频率值变为，与频率为0时的情况一致；同理，在频率为的地方，那么这个点的频率值变为，与前面的情况一致，又是一个循环。横轴必须是以Hz单位进行转换。0246810x104012x105abs0246810x104-202x105real0246810x104-202x105imagfft变换存在的缺陷与问题有人问，为什么我画的频域图与书上的不一样啊，这就是你在对参数进行设置时，没有考虑到它们对于fft变换后的结果的影响，或者根本就没设置正确。不要直接就进行fft变换，变换前必须先估计信号的频率范围，和你想要获得的图形精确程度，然后适当改变参数，确定不会发生下面这些问题。1.混叠效应如果x(t)的频谱是带限的，X(f)=0,|f|fm，则由抽样定理，抽样间隔最大满足Ts=1/2fm如果f(t)的频谱不是带限的，则抽样后频谱总要发生混叠，减小抽样间隔Ts，fs增大，可减小混叠，但工作量增加。解决办法：预滤波，再抽样，或者一般选择Ts1/(3~5)fm。2.泄漏(leakage)若X(f)为有限带宽频谱，则x(t)为时间无限的。为利用FFT分析x(t)的频谱，必须截取x(t)有限范围，即加窗。频域卷积后，使原频带受限的频谱扩展开来(有限带宽拖了尾巴)，这种现象称为泄漏。解决方法：改善窗的形状或长度。3.栅栏效应（PicketFenceEffect）N点序列的DFT只能在有限的N个频点上观察频谱，这相当于从栅栏的缝隙中观察景色，对于了解信号在整个频域上的特性是不够的。为了观察到其他频率的信息，需要对原信号x[n]做一些处理，以便在不同的频点上采样。解决方法：增加采样点数，信号补零。4.补零与fft分辨率通过补零观察到更多的频点，但是这并不意味着补零能够提高真正的频谱分辨率。这是因为x[n]实际上是x(t)采样的主值序列，而将x[n]补零得到的x'[n]周期延拓之后与原来的序列并不相同，也不是x(t)的采样。因此与是不同离散信号的频谱。对于补零至M点的x'的DFT，只能说它的分辨率2π/M仅具有计算上的意义，并不是真正的、物理意义上的频谱。物理频谱分辨率的提高只能在满足采样定理的条件下增加时域采样长度来实现。5.物理分辨率与计算分辨率频率分辨率是指所用的算法能将信号中两个考得很近的谱峰保持分开的能力。若信号x(t)的长度为T秒，通过傅里叶变换后得到X（jw)，那么X（jw)的频率分辨率是：Δf=1/T(Hz)而计算分辨率是Δf=fs/M(Hz)当M长度不够时，可以通过对x(t)后面补零，使补零后的序列达到所需要的长度N，但这并没有增加频谱的分辨率，因为有效的数据长度M并没有增加。计算分辨率是靠计算得出的，并不能反映真实的频率分辨能力，与“物理分辨率”是不相同的。2015年5月22日BYXUHAO百度账号：xuhaobestone2