用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

fft是快速傅立叶变换,可直接调用,如fft(A)

离散傅立叶变换是dft

这是我自己做的dft

function

X=dft(x)

N=length(x)

W=exp(-2i*pi/N)

X=zeros(1,N)

for

k=1:N

X(k)=sum(x.*W.^((0:N-1)*(k-1)))

end

% 下面的程序里Pn 存的就是基波相位 如果求的是谐波相位，稍微修改即可

x = load('data.dat')%load 数据

fs=10000% 采样频率，自己根据实际情况设置

N=length(x)% x 是待分析的数据

n=1:N

%1-FFT

X=fft(x)% FFT

X=X(1:N/2)

Xabs=abs(X)

Xabs(1) = 0%直流分量置0

[Amax,index]=max(Xabs)

if(Xabs(index-1) >Xabs(index+1))

a1 = Xabs(index-1) / Xabs(index)

r1 = 1/(1+a1)

k01 = index -1

else

a1 = Xabs(index) / Xabs(index+1)

r1 = 1/(1+a1)

k01 = index

end

Fn = (k01+r1-1)*fs/N%基波频率

An = 2*pi*r1*Xabs(k01)/(N*sin(r1*pi))%基波幅值

Pn = phase(X(k01))-pi*r1%基波相位 单位弧度

Pn = mod(Pn(1),pi)

---