c=rand(1,512);% 定制512點抽樣%
%%%%求反序，實現反序輸入，正序輸出%%%%
c(1)=0;c(2)=256;
for i=1:8   % log(512)-1=8%
     c((2^i+1):2^(i+1))=c(1:2^i)+2^(8-i);
end
d=rand(1,512);
t=linspace(1,100,512)/100;b=[t;sin(2*pi*t*5)];  %輸入5HZ正弦波，
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 取均分512點的正弦值作為fft的初始值
for i=1:512
     d(i)=b(2,(c(i)+1));
end%完成初始值的反序%
%%%%求出所有的蝶算指數因子%%%%
x=rand(512,10);a=rand(1,256);q=1;w=rand(1,256);
for o=1:256
     w(o)=exp(-j*2*pi*(o-1)/512);
end  
%%%%fft變換%%%%
x(:,1)=d';
for k=2:10
     i=1;p=1;
     for m=1:2^(k-2)%求出每級所需的蝶算指數因子%
         q=(m-1)*512/(2^(k-1));
         a(p)=w(q+1);p=p+1;
     end
     for n=1:(512/(2^(k-1)))%對偶節點個數的確定及變換計算%
         p=1;
         for h=i:(i+(2^(k-2)-1))
             x(h,k)=x(h,k-1)+a(p)*x(h+2^(k-2),k-1);
             x(h+2^(k-2),k)=x(h,k-1)-a(p)*x(h+2^(k-2),k-1);
             p=p+1;
         end
         i=(2^(k-1))*n+1;
     end
end
%%%%取出fft變換后的點數值，取1到50點的值進行作圖%%%%
S=abs(x(:,10));
plot((0:49)/50*50,abs(S(1:50)),'linewidth',2)























我是根據一個DITFFT（時域抽選）程序進行修改的DIFFFT（頻域抽選），但是畫圖的時候發現信號和原先他的程序差異很大，我懷疑是我的程序編寫的有問題，請大家幫忙改改！
% DIF-FFT

N=256;
x1=[1 1 1 1];
x=[x1 zeros(1,N-length(x1))];


L=1;
M=log2(N);
while L<=M
LE=2^L;
LE1=LE/2;
U=1;
W=exp(-j*2*pi/LE);
J=0;
while J<=LE1-1;
I=J;
while I<=N-1
IP=I+LE1;
T=x(I+1)+x(IP+1);
x(IP+1)=(x(I+1)-x(IP+1))*U;
x(I+1)=T;
I=I+LE;
end
U=U*W;
J=J+1;
end
L=L+1;
end
x;

NV2=N/2;
NM1=N-1;
I=0;
J=0;
while I<NM1
if I<J
T=x(J+1);
x(J+1)=x(I+1);
x(I+1)=T;
end
K=NV2;

while K<=J
J=J-K;
K=K/2;
end
J=J+K;
I=I+1;
end
x;

n=0:N-1;
q=n*2*pi/N;
stem(q,angle(x),'.')

---------------------------------------------
這是DIT-FFT的程序，我就是根據這個來修改的，測試通過，而DIF-FFT則波形不對。
N=256;
x1=[1 1 1 1];
x=[x1 zeros(1,N-length(x1))];

NV2=N/2;
NM1=N-1;
I=0;
J=0;
while I<NM1
if I<J
T=x(J+1);
x(J+1)=x(I+1);
x(I+1)=T;
end
K=NV2;

while K<=J
J=J-K;
K=K/2;
end
J=J+K;
I=I+1;
end
x;

L=1;
M=log2(N);
while L<=M
LE=2^L;
LE1=LE/2;
U=1;
W=exp(-j*2*pi/LE);
J=0;
while J<=LE1-1;
I=J;
while I<=N-1
IP=I+LE1;
T=x(IP+1)*U;
x(IP+1)=x(I+1)-T;
x(I+1)=x(I+1)+T;
I=I+LE;
end
U=U*W;
J=J+1;
end
L=L+1;
end
x;
n=0:N-1;
q=n*2*pi/N;
stem(q,abs(x),'.') 

 
作者：fang5566      發表時間：2007-4-29 15:09:00

　第1樓  
哪位兄弟幫我改改啊，急需啊！！！！！！！！

 
作者：zhoujb      發表時間：2007-5-2 13:45:00

　第2樓  
你好，不知道你還在不在？

 
作者：fang5566      發表時間：2007-5-4 10:53:00

　第3樓  
多謝樓上關心 我已經自行解決了

就是LE=2^L;改成LE=2^(M-L+1) 















基于Matlab的時間抽取基2FFT算法
function y=myditfft(x)
%本程序對輸入序列實現DIT-FFT基2算法，點數取大于等于長度的2的冪次
%------------------------------------
%            Leo's fft program(改編網上的一個程序)
%------------------------------------
m=log2(2^nextpow2(length(x)));                %求的x長度對應的2的最低冪次m
N=2^m;
if length(x)<N
    x=[x,zeros(1,N-length(x))];               %若長度不是2的冪，補0到2的整數冪
end
x;
%--------------------------------------------------------------------------
%對輸入序列進行倒序
%如果輸入序列的自然順序號I用二進制數（例如n2n1n0）表示
%則其倒位序J對應的二進制數就是（n0n1n2），這樣，在原來自然順序時應該放x(I)的
%單元，現在倒位序后應放x(J)。
%--------------------------------------------------------------------------
%以下程序相當于以下程序:
%nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1;    %求1：2^m數列的倒序
%y=x(nxd);                                     %將倒序排列作為初始值
%--------------------------------------------------------------------------
NV2=N/2;
NM1=N-1;
I=0;
J=0;
while I<NM1
    if I<J
        T=x(J+1);
        x(J+1)=x(I+1);
        x(I+1)=T;
    end
    K=NV2;
   
    while K<=J
        J=J-K;
        K=K/2;
    end
    J=J+K;
    I=I+1;
end              
x;
%--------------------------------------------------------------------------
%以下程序解釋:
%第一級從x(0)開始,跨接一階蝶形,再取每條對稱
%第二級從x(0)開始,跨接兩階蝶形,再取每條對稱
%第m級從x(0)開始,跨接2^(m-1)階蝶形,再取每條對稱....
%--------------------------------------------------------------------------
for mm=1:m                                    %將DFT做m次基2分解，從左到右，對每次分解作DFT運算
    Nmr=2^mm;
    u=1;                                      %旋轉因子u初始化
    WN=exp(-j*2*pi/Nmr);                      %本次分解的基本DFT因子WN＝exp(-i*2*pi/Nmr)
      for n=1:Nmr/2                           %本次跨越間隔內的各次碟形運算
          for k=n:Nmr:N                       %本次碟形運算的跨越間隔為Nmr=2^mm
              kp=k+Nmr/2;                     %確定碟形運算的對應單元下標(對稱性)
              t=x(kp)*u;                      %碟形運算的乘積項
              x(kp)=x(k)-t;                   %碟形運算的加法項
              x(k)=x(k)+t;
          end
          u=u*WN;                             %修改旋轉因子，多乘一個基本DFT因子WN
      end
end
y=x;                                    %輸出