DFT的计算步骤是什么?

DFT的计算步骤如下：

离散傅里叶变换（Discrete Fourier Transform，DFT）傅里叶分析方法是信号分析的最基本方法，傅里叶变换是傅里叶分析的核心，通过它把信号从时间域变换到频率域，进而研究信号的频谱结构和变化规律。

基本性质

1.线性性质

如果X1（n）和X2(N)是两个有限长序列，长度分别为N1和N2，且Y(N)=AX1(N)+BX2(N)。

式中A,B为常数，取N=max[N1,N2]，则Y(N）的N点DFT为：

Y(K)=DFT[Y(N)]=AX1(K)+BX2(K), 0≤K≤N-1。

2.循环移位特性

设X(N)为有限长序列，长度为N，则X(N)地循环移位定义为：

Y(N)=X((N+M))下标nR（N）。

式中表明将X(N）以N为周期进行周期拓延得到新序列X'(N)=X((N))下标n，再将X'(N）左移M位，最后取主值序列得到循环移位序列Y(N)。

1. 用联立方程计算

已知时域中N个点的值，将正弦波加起来，可以列N个方程，求解N个未知数。

虽然实际中不这么做，但它说明了信号为什么能分解成N + 2个（其中两个幅值恒为0）个正弦波，为什么会刚好有解。而且这些基本函数必须线性不相关。

2. DFT标准算法——通过相关性计算

相关性算法成立的条件是基本函数正交。

3. FFT

#define PI 3.1415926

#define PI2 6.2831853

#include "stdio.h"

#include "math.h"

/* inv=1 forward transforminv=-1 inverse transform */

fft(float sr[],float sx[],int m0,int inv)

{

int i,j,lm,li,k,lmx,np,lix,mm2

float t1,t2,c,s,cv,st,ct

if(m0<0)

return

inv=-inv

lmx=1

for(i=1i<=m0++i)

lmx+=lmx //get 2**m0

cv=PI2/(double)lmx

ct=cos(cv)st=sin(cv)*inv

np=lmxmm2=m0-2

/* fft butterfly numeration */

for(i=1i<=mm2++i)

{

lix=lmxlmx/=2

c=cts=st

for(li=0li<npli+=lix)

{

j=lik=j+lmx

t1=sr[j]-sr[k]t2=sx[j]-sx[k]

sr[j]+=sr[k]sx[j]+=sx[k]sr[k]=t1sx[k]=t2

++j++k

t1=sr[j]-sr[k]t2=sx[j]-sx[k]

sr[j]+=sr[k]sx[j]+=sx[k]

sr[k]=c*t1-s*t2sx[k]=s*t1+c*t2

}

for(lm=2lm<lmx++lm)

{

cv=cc=ct*c-st*ss=st*cv+ct*s

for(li=0li<npli+=lix)

{

j=li+lmk=lmx+j

t1=sr[j]-sr[k]t2=sx[j]-sx[k]

sr[j]+=sr[k]sx[j]+=sx[k]

sr[k]=c*t1-s*t2sx[k]=s*t1+c*t2

}

}

cv=ctct=2.0*ct*ct-1.0st=2.0*st*cv

}

/* 4 points DFT */

if(m0>=2)

for(li=0li<npli+=4)

{

j=lik=j+2

t1=sr[j]-sr[k]t2=sx[j]-sx[k]

sr[j]+=sr[k]

sx[j]+=sx[k]

sr[k]=t1sx[k]=t2

++j++k

t1=sr[j]-sr[k]t2=sx[j]-sx[k]

sr[j]+=sr[k]sx[j]+=sx[k]

sr[k]=-inv*t2sx[k]=inv*t1

}

/* 2 points DFT */

for(li=0li<npli+=2)

{

j=lik=j+1

t1=sr[j]-sr[k]t2=sx[j]-sx[k]

sr[j]+=sr[k]sx[j]+=sx[k]

sr[k]=t1sx[k]=t2

}

/* sort according to bit reversal */

lmx=np/2j=0

for(i=1i<np-1++i)

{

k=lmx

while(k<=j)

{

j-=kk/=2

}

j+=k

if(i<j)

{

t1=sr[j]sr[j]=sr[i]sr[i]=t1

t1=sx[j]sx[j]=sx[i]sx[i]=t1

}

}

/* if Inverse FFT, multiply 1.0/np */

if(inv!=1.0)

return

t1=1.0/np

for(i=0i<np++i)

{

sr[i]*=t1sx[i]*=t1

}

}

main()

{

float xr[2048],xi[2048],xt[2048]

int i,np,nfft,k,nf

float t,dt,df,f,hf

float f0=10,f1=20,f2=30

FILE *fp1,*fp2

char fil1[60],fil2[60]

printf("ENTER TIME SIGNAL FILE\n")

scanf("%s",fil1)

printf("ENTER AMPLITUDE SPECTRUM FILE\n")

scanf("%s",fil2)

printf("SAMPLE POINTS?\n")

scanf("%d",&np)

printf("SAMPLE INTERVAL(S)?\n")

scanf("%f",&dt)

printf("HIGH CUTOFF FREQUENCY(Hz)?\n")

scanf("%f",&hf)

fp1=fopen(fil1,"w")

fp2=fopen(fil2,"w")

for(i=0i<npi++)

{ t=(i-np/2)*dt

if(t!=0.0)xt[i]=1.0/(PI*t)

else xt[i]=0.0

xr[i]=xt[i]

}

// calculate fft point

k=log(np*1.0)/log(2.0)

if(np>pow(2.0,k*1.0))k=k+1

nfft=pow(2.0,k*1.0)

df=1.0/(nfft*dt)

nf=hf/df+1

printf("nfft=%d k=%d\n",nfft,k)

printf("dt=%f df=%f nf=%d\n",dt,df,nf)

// fill zero

if(np<nfft)

{for(i=npi<nffti++)

xr[i]=0

}

for(i=0i<nffti++)

xi[i]=0.0

// calculate fft

fft(xr,xi,k,1)

// output amplitude and argument

for(i=0i<nfi++)

{ f=i*df

fprintf(fp2,"%d %8.2f %12.4f %12.4f\n",i,f,atan(xr[i]/xi[i]),sqrt(xr[i]*xr[i]+xi[i]*xi[i]))

}

fft(xr,xi,k,-1)

for(i=0i<npi++)

{ t=i*dt

fprintf(fp1,"%10d %10.4f %10.4f %10.4f\n",i+1,t,xt[i],xr[i])

}

fclose(fp1)

fclose(fp2)

}

---