dit.m

按时间抽选（DIT）的基-2 FFT算法

clear,format compact
[y,Fs,bits]=wavread('joy.wav');%读出信号，采样率和采样位
y=y(:,1);%取单声道
y=y';
m=nextpow2(y);N=2^m; % 求y的长度对应的2的最低幂次m
if size(y) 
    y=[y,zeros(1,N-length(y))]; % 若y的长度不是2的幂，补零到2的整数幂
end
tt=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1; % 求1:2^m数列的倒序
t=y(tt); % 将y倒序排列作为t的初始值
for mm=1:m % 将DFT作m次基2分解,从左到右，对每次分解作DFT运算,mm为层数
u=1; % 旋转因子u初始化为WN^0=1

for j=1:2^mm:N-(2^mm-1) % 本次跨越间隔内的各次蝶形运算
for k=j:1:j+2^(mm-1)-1 % 本次蝶形运算的跨越间隔为1
kp=k+2^(mm-1); % 确定蝶形运算的对应单元下标
p=t(kp)*u; % 蝶形运算的乘积项
t(kp)=t(k)-p; % 蝶形运算
t(k)=t(k)+p; % 蝶形运算
u=u*exp(-i*2*2^(m-mm)*pi/N);
end
u=1; % 旋转因子恢复初值
end
end
subplot(2,1,1),stem(y,'.'),title('原波形图');%画出原波形图
subplot(2,1,2),stem(abs(t),'.'),title('DIT后的频谱图');%画出DIT后的频谱图