ditfft.m

来自「fft自编源码」· M 代码 · 共 41 行

%该程序实现快速傅立叶变换
%输入参数可以是行或列向量
%函数返回值为x的傅立叶变换值及其做变换所用的点数N
function [y1,N]=ditfft(x)

for p=1:inf  %找x的长度对应的2的最低幂次
    if 2^p>=length(x);break;end
end 

N=2^p;

if size(x,1)>1, x=x';end %若为列向量则转换为行向量

if length(x<N)
  x1=[x,zeros(1,N-length(x))];%补零
end

handle1=@flip1;%调用flip1倒序函数
y1=feval(handle1,x1,N);

for p1=1:p     %将DFT作P次分解，依次计算
 
 nm=2^p1;u1=1;
 wn=exp(-i*2*pi/nm); %本次分解的基本因子wn
  
 for j=1:nm/2       %本次跨越间隔内的各碟形运算
       
     for k=j:nm:N  %本次碟形运算的跨越间隔为nm=2^p1
           kp=k+nm/2; %k,kp为碟形运算的对应单元下标
           t=y1(kp)*u1;
           y1(kp)=y1(k)-t;
           y1(k)=y1(k)+t;
       end
    
       u1=u1*wn;
       
   end
   
end