快速傅里叶算法实现

function xk=fft_my(xn,N)% *********************Declaration***************************% File name:        fft_my                                                             % Author:           @Harry                                              % Date:             2016/8/25 13:42:14                                  % Version Number:   1.0                                                                                                                                                                                                              % Modification history:[including time, version, author and abstract]   %     2016/8/25       version 1.0     xxx                                  % Abstract: %     用时间抽取的基2FFT-DFT进行快速傅里叶变换                                                                                                                          % *********************************end******************************* %计算FFT的级数M M=log2(N); %对输入序列进行码位倒序 xn=bit_reverse(xn); tempdata=xn; for m=1:M   wcount=2^(m-1); %计算此级旋转因子的个数wcount      %计算这一级中wcount个旋转因子的值wvalue   for w=1:wcount      wvalue(w)=exp(-1i*pi*(w-1)/wcount);   end     %按照此级中旋转因子的个数分组   for num=1:wcount     group=(N/2)/(2^(m-1)); %同一旋转因子可分为group个组分别计算     for g=1:group       p=num+(g-1)*(2^m);       q=p+2^(m-1);       xk(p)=tempdata(p)+wvalue(num)*tempdata(q);       xk(q)=tempdata(p)-wvalue(num)*tempdata(q);     end   end   tempdata=xk; endfunction xk=bit_reverse(xn)% *********************Declaration***************************% File name:        bit_reverse                                                            % Author:           @Harry                                              % Date:             2016/8/25 8:47:16                                   % Version Number:   1.0                                                                                                                                                                                                              % Modification history:[including time, version, author and abstract]   %     2016/8/25       version 1.0     xxx                                  % Abstract: %     计算一个序列的码位倒序(序列长度为2^n)。把自然顺序的十进制转换成二进制数，% 然后将这些二进制的首末尾倒序再重新转换成十进制，那么，这时的十进制% 的排列就是码位倒序排列。(自然顺序必须从0开始)                                                                                                                          % *********************************end***********************%检查输入是否符合要求[row,column]=size(xn);if row~=1   error('输入必须是一维序列');endlength=column; %序列长度xk=zeros(1,length);for i=0:length-1 %自然顺序必须从0开始    tempbinary=dec2bin(i,log2(length));    m=log2(length); %码位宽度    %开始对二进制的码位倒序    for j=1:m        tempdata(j)=tempbinary(m-j+1);    end    n=bin2dec(tempdata);    xk(i+1)=xn(n+1); %数组下标必须从1开始end