如何用matlab定义DFT

1、首先双击matlab软件图标，打开matlab软件，可以看到matlab软件的界面。

2、使用syms命令，创建四个符号变量a、b、c、x、t。

3、使用符号变量a，创建代数式A，其中A=7*sin(a)。

4、使用函数fourier(A,a,t)，对代数式A进行DFT变换。得到的结果中diract(t-1)是狄拉克函数。

5、使用符号变量c，创建代数式B，其中A=3*c^2。

6、使用函数fourier(B,c,t)，对代数式B进行DFT变换。得到的结果中dirac(2,t)是对狄拉克函数的二阶导数。

#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)

}

求解过程如下:

clc,clear,close all

N=128

n = 0:1:N

m = 0:1:N

x=cos(2.*pi.*m.*n./N)

figure,plot(n,x)

MATLAB（矩阵实验室）是MATrix LABoratory的缩写，是一款由美国The MathWorks公司出品的商业数学软件。MATLAB是一种用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。除了矩阵运算、绘制函数/数据图像等常用功能外，MATLAB还可以用来创建用户界面及与调用其它语言（包括C，C++和FORTRAN）编写的程序。

---