matlab 基2时间抽取（DIT）快速傅里叶变换

clear;%信号Fs=1000;%?采样频率?T=1/Fs;%? 采样时间间隔??N=1024;%? 信号长度??t=(0:N-1)*T;%?时间点向量x=sin(2*pi*50*t)+0.7*sin(2*pi*200*t);figure(1)plot(t,x)%系数矩阵M=log2(N);xishu=zeros(N,M+1);%计算旋转因子k=0:1:N-1;W=exp(-1i*2*pi/N);Wk=W.^k;%先码位倒置得到序号xuhao=zeros(N,2);xuhao(:,1)=0:1:N-1;for i=1:1:N    xuhao(i,2)=bin2dec(fliplr(dec2bin(xuhao(i,1),M)));end%得到序号后变换位置x_=zeros(1,N);for i=1:1:N    x_(i)=x(xuhao(i,2)+1);endxishu(:,1)=x_;%循环M次，计算得到结果for i=0:1:M-1%i为级数     step=2^i;    %计算每级组数    zushu=N/2^(i+1);        %每次计算之前对系数进行处理，也是分组，每组的后半组*Wk    for j=1:1:zushu%一组一组的来        for k=(N/zushu/2+1):1:N/zushu            xishu((j-1)*(N/zushu)+k,i+1)=xishu((j-1)*(N/zushu)+k,i+1)*Wk(1+2^(M-i-1)*(k-((N/zushu/2+1))));        end    end        for j=1:1:zushu%一组一组的来        for k=1:1:N/zushu/2            %计算向下的            xishu((j-1)*(N/zushu)+k,i+2)=xishu((j-1)*(N/zushu)+k,i+1)+xishu((j-1)*(N/zushu)+k+step,i+1);        end        for k=(N/zushu/2+1):1:N/zushu            %计算向上的            xishu((j-1)*(N/zushu)+k,i+2)=xishu((j-1)*(N/zushu)+k-step,i+1)-xishu((j-1)*(N/zushu)+k,i+1);        end    endendfigure(2)f=1/(N*T)*(0:N/2);plot(f,2*abs(xishu(1:N/2+1,M+1))/N)