用matlab编写实现fft的程序。

function y=myditfft(x)

%本程序对输入序列实现DIT-FFT基2算法，点数取大于等于长度的2的幂次

%------------------------------------

%

myditfft.c

%------------------------------------

m=nextpow2(x)

%求的x长度对应的2的最低幂次m

N=2^m

if length(x)<N

x=[x,zeros(1,N-length(x))]

%若的长度不是2的幂，补0到2的整数幂

end

nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1

%求1：2^m数列的倒序

y=x(nxd)

%将倒序排列作为的初始值

for mm=1:m

%将DFT做m次基2分解，从左到右，对每次分解作DFT运算

Nmr=2^mm

u=1

%旋转因子u初始化

WN=exp(-i*2*pi/Nmr)

%本次分解的基本DFT因子WN＝exp(-i*2*pi/Nmr)

for j=1:Nmr/2

%本次跨越间隔内的各次碟形运算

for k=j:Nmr:N

%本次碟形运算的跨越间隔为Nmr=2^mm

kp=k+Nmr/2

%确定碟形运算的对应单元下标

t=y(kp)*u

%碟形运算的乘积项

y(kp)=y(k)-t

%碟形运算的加法项

y(k)=y(k)+t

end

u=u*WN

%修改旋转因子，多乘一个基本DFT因子WN

end

end

% 基于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 %输出

---