用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);

grid on

m=nextpow2(x);N=2^m; % 求x的长度对应的2的最低幂次m

if length(x)

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芯片的存储量也将增大。