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后的頻譜圖