用matlab实现fft算法.docVIP

A1=str2double(get(handles.edit8,String)); A2=str2double(get(handles.edit9,String)); F1=str2double(get(handles.edit10,String)); F2=str2double(get(handles.edit11,String)); Fs=str2double(get(handles.edit12,String)); N=str2double(get(handles.edit13,String)); t=[0:1/Fs:(N-1)/Fs]; x=A1*sin(2*pi*F1*t)+A2*sin(2*pi*F2*t); %信号x的离散值 axes(handles.axes1) %在axes1中作原始信号图 plot(x); grid on m=nextpow2(x);N=2^m; % 求x的长度对应的2的最低幂次m if length(x)N x=[x,zeros(1,N-length(x))];% 若x的长度不是2的幂，补零到2的整数幂 end nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1;% 求1:2^m数列序号的倒序 y=x(nxd); %将x倒序排列作为y的初始值 for L=1:m; %将DFT作m次基2分解,从左到右，对每次分解作DFT运算,共做m级蝶形运算，每一级都有2^(L-1)个蝶形结 B=2^(L-1);%两个输入数据相距2^(L-1) for j=0:B-1 ;%J代表了不同的旋转因子 p=2^(m-L)*j; for k=(j+1):2^L:N ;%本次蝶形运算的跨越间隔为2^L WN=exp(-i*2*pi*p/N);%计算本次运算的旋转因子 T=y(k)+y(k+B)*WN ;%计算K地址上蝶形项 y(k+B)=y(k)-y(k+B)*WN ;%计算（K+B）地址上的蝶形项并仍然放回（K+B） y(k)=T ;%将原来计算的K地址蝶形项放回K地址，注意必须先进行复数加法运算 end end end axes(handles.axes2) %在axes2中作自编fft运算后的幅频特性图 Ayy = (abs(y)); %取模 Ayy=Ayy/(N/2); %换算成实际的幅度 F=([1:N]-1)*Fs/N; %换算成实际的频率值 stem(F(1:N/2),Ayy(1:N/2)); %显示换算后的FFT模值结果 axes(handles.axes3) %在axes3中作系统fft运算后的幅频特性图 G=fft(x,N); %利用系统作N点FFT变换，作对比 Ayy1= (abs(G)); %取模 Ayy1=Ayy1/(N/2); %换算成实际的幅度 F=([1:N]-1)*Fs/N; %换算成实际的频率值 stem(F(1:N/2),Ayy1(1:N/2)); %显示换算后的FFT模值结果 GUI界面 输入A1=4，A2=3，f1=50Hz，f2=80Hz，采样频率Fs=256Hz，采样点数N=512 所得频谱图： 图1 由图1可知，我所编写的fft算法与系统自带的fft算法结果一致，是正确的 倘若我未知信号频率f1，f2，仍然设f1=50Hz，f2=80Hz，在保证Fs大于两倍信号频率的前提下（Fs仍然为256Hz） 改变采样点数，取N=64 图2 由图2可知，50Hz频率点处产生误差 取N=40 图3 由图3可知，50Hz频率点和80Hz频率点处均产生误差 再取N=8 图4 误差更大了！ 取N=1000 图5 由图五可知，当N=1000时，所得图形接近理想值 再取N=1024 图6 与理想值符合，但与图1比，分辨率更高 结论：当我们已知一个周期信号的周期时，我们对其进行频谱分析的时候，（由图1和图6可知）对其进行fft运算所取的点数N必须满足分辨率（Fs/N）能够被信号频率整除才能保证没有误差，但是在实际的情况中，我们往往不知道信号的周期是多少，所以我们在对这类未知信号进行频谱分析的时候，在保证时域采样频率满足乃奎斯特条件的情况下，对其作FFT运算的时候采样点数N尽量地取大一点，这样才能保证减小误差，由图3与图5对比可知，当N的值较大时，所得到的频谱图误差较小，反之误差很大。但是，N的值不能无限制地取大，因为当N增大时，意味着运算量的增大，也意味着所占DSP芯片的存储量也将增大。