用matlab实现fft算法

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);gridonm=nextpow2(x);N=2^m;%求x的长度对应的2的最低幂次miflength(x)Nx=[x,zeros(1,N-length(x))];%若x的长度不是2的幂，补零到2的整数幂endnxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1;%求1:2^m数列序号的倒序y=x(nxd);%将x倒序排列作为y的初始值forL=1:m;%将DFT作m次基2分解,从左到右，对每次分解作DFT运算,共做m级蝶形运算，每一级都有2^(L-1)个蝶形结B=2^(L-1);%两个输入数据相距2^(L-1)forj=0:B-1;%J代表了不同的旋转因子p=2^(m-L)*j;fork=(j+1):2^L:N;%本次蝶形运算的跨越间隔为2^LWN=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地址，注意必须先进行复数加法运算endendendaxes(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芯片的存储量也将增大。