matlab.txt

基于matlab 的 dft 的 fft 的 是吸纳

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;                                    %输出