《用matlab实现fft算法.pdf》由会员分享，可在线阅读，更多相关《用matlab实现fft算法.pdf（5页珍藏版）》请在三一文库上搜索。
1、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*
2、F2*t);%信号 x 的离散值 axes(handles.axes1)%在 axes1中作原始信号图 plot(x); grid on m=nextpow2(x);N=2m;% 求 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:2m 数列序号的倒序 y=x(nxd); %将 x 倒序排列作为y 的初始值 for L=1:m; %将 DFT 作 m 次基 2 分解 ,从左到右，对每次分解作
3、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):2L:N;%本次蝶形运算的跨越间隔为2L 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
4、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;%换算成
5、实际的频率值 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 频率
6、点和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芯片的存储量也将增大。