数字信号处理课程设计-DFT的快速算法——快速傅里叶变换(FFT).doc

燕 山 大 学 课 程 设 计 说 明 书目录前言第一章 离散傅里叶变换dft31.1 dft定义31.2 dft的快速算法快速傅里叶变换(fft) 3第二章 基2 dit-fft算法42.1 按时域抽取的基2 dit-fft算法4第三章 基于c语言设计16点基2dit-fft程序及运行结果 63.1 按时间抽取的基2dit-fft程序 63.2 程序运行结果 15第四章 课程设计的总结 17参考文献 17 前 言信号(signal)是一种物理体现，或是传递信息的函数。而信息是信号的具体内容。模拟信号(analog signal)：指时间连续、幅度连续的信号。数字信号(digital signal)：时间和幅度上都是离散(量化)的信号。 数字信号可用一序列的数表示，而每个数又可表示为二制码的形式，适合计算机处理。数字信号处理是将信号以数字方式表示并处理的理论和技术。数字信号处理与模拟信号处理是信号处理的子集。 其目的是对真实世界的连续模拟信号进行测量或滤波。因此在进行数字信号处理之前需要将信号从模拟域转换到数字域，这通常通过模数转换器实现。而数字信号处理的输出经常也要变换到模拟域，这是通过数模转换器实现的。数字信号处理的核心算法是离散傅立叶变换(dft)，是dft使信号在数字域和频域都实现了离散化，从而可以用通用计算机处理离散信号。而使数字信号处理从理论走向实用的是快速傅立叶变换(fft)，fft的出现大大减少了dft的运算量，使实时的数字信号处理成为可能、极大促进了该学科的发展。随着大规模集成电路以及数字计算机的飞速发展，加之从60年代末以来数字信号处理理论和技术的成熟和完善，用数字方法来处理信号，即数字信号处理，已逐渐取代模拟信号处理。 随着信息时代、数字世界的到来，数字信号处理已成为一门极其重要的学科和技术领域。第一章 离散傅里叶变换dft离散傅立叶变换（dft）实现了信号首次在频域表示的离散化，使得频域也能够用计算机进行处理。并且这种dft变换可以有多种实用的快速算法。使信号处理在时、频域的处理和转换均可离散化和快速化，因而具有重要的理论意义和应用价值。1.1 dft定义设序列x(n)长度为m，定义x(n)的n点dft为式中，n称为离散傅里叶变换区间长度，要求n m。为书写简单，令 ，因此通常将n点dft表示为1.2 dft的快速算法快速傅里叶变换(fft)dft使计算机在频域处理信号成为可能，但是当n很大时，直接计算n点dft的计算量非常大。快速傅里叶变换（fft，fast fourier transform）可使实现dft的运算量下降几个数量级，从而使数字信号处理的速度大大提高，工程应用成为可能。人们已经研究出多种fft算法，它们的复杂度和运算效率各不相同。快速傅里叶变换就是不断地将长序列的dft分解为短序列的dft，并利用 的周期性和对称性及其一些特殊值来减少dft运算量的快速算法。fft算法分类：1.时间抽选法 dit：decimation-in-time2.频率抽选法 dif：decimation-in-frequency时间域抽取：基2时间抽取(decimation in time)dit-fft算法频率域抽取：基2频率抽取(decimation in frequency)dif-fft算法第二章 基2 dit-fft算法2.1 按时间抽取的基2 dit-fft算法：1、按时间抽取的基2 dit-fft算法原理先设序列点数为n=2m，m为整数。如果不满足这个条件，可以人为地加上若干零值点，使之达到这一要求。这种n为2的整数幂的fft称基-2 fft。设输入序列长度为n=2m (m为正整数) ，将该序列按时间顺序的奇偶分解为越来越短的子序列，称为按时间抽取(dit )的fft算法。 序列x(n)的n点dft为 将上面的和式按n的奇偶性分解为 令x1(l)=x(2l),x2(l)=x(2l+1)，因为 ，所以上式可写成上式说明，按n的奇偶性将x(n)分解为两个n/2长的序列x1(l)和x2(l)，则n点dft可分解为两个n/2点dft来计算。用x1(k)和x2(k)分别表示x1(l)和x2(l)的n/2点dft，即有上述公式，及x1（k）、x2(k)的隐含周期性得到：这样，就将n点dft的计算分解为计算两个n/2点离散傅里叶变换x1(k)和x2(k)，再计算上式。蝶形图：2、 按时间抽选的fft算法的特点：（1）原位运算由图4.2.4可以看出，dit-fft的运算过程很有规律。n=2m点的fft共进行m级运算，每级由n/2个蝶形运算组成。同一级中，每个蝶形的两个输入数据只对计算本蝶形有用，而且每个蝶形的输入、输出数据结点又同在一条水平线上，这就意味着计算完一个蝶形后，所得输出数据可立即存入原输入数据所占用的存储单元。这样，经过m级运算后，原来存放输入序列数据的n个存储单元中便依次存放 的n个值。这种利用同一存储单元存储蝶形计算输入、输出数据的方法称为原位计算。原位计算可节省大量内存，从而使设备成本降低。（2）倒序规律按原位计算时，fft的输出 是按正常顺序排列在存储单元中，但是这时输入x(n)却不是、按自然顺序存储的，而是按x(0)，x(4)，x(7)的顺序存入存储单元，看起来好象是“混乱无序”的，实际上是有规律的，我们称之为倒序。第三章 基于c语言设计16点基2dft-fft程序及运行结果设计参数及要求：本次试验要求利用c语言设计16点基2dit-fft程序，并对如下程序进行分析，给出输出结果并画出信号的幅值图和相位谱。n012345678910x(n)01.222.232.85 2.982.601.760.62-0.62-1.76-2.601112131415-2.98-2.85-2.23-1.2203.1 按时间抽取的基2dft-fft程序：/*时间抽选基2 fft算法c语言实现*/#include #include #include #include #define n 1000#include/*定义复数类型*/typedef struct double real; double img;complex;complex xn, *w; /*输入序列,变换核*/int size_x=0; /*输入序列的大小，在本程序中仅限2的次幂*/double pi; /*圆周率*/int z;float wn,pn;char strn;/*初始化变换核*/void initw() int i; w=(complex *)malloc(sizeof(complex) * size_x); for(i=0;ireal=a.real+b.real; c-img=a.img+b.img; /*乘法*/void mul(complex a,complex b,complex *c) c-real=a.real*b.real - a.img*b.img; c-img=a.real*b.img + a.img*b.real; /*减法*/void sub(complex a,complex b,complex *c) c-real=a.real-b.real; c-img=a.img-b.img; /*除法*/void divi(complex a,complex b,complex *c) c-real=( a.real*b.real+a.img*b.img )/( b.real*b.real+b.img*b.img); c-img=( a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img); /*变址计算，将x(n)码位倒置*/void change() complex temp; unsigned short i=0,j=0,k=0; double t; for(i=0;i0 ) j=j1; if(ji) temp=xi; xi=xj; xj=temp; /*快速傅里叶变换*/void fft() int i=0,j=0,k=0,l=0; complex up,down,product; change(); for(i=0;i log(size_x)/log(2) ;i+) /*一级蝶形运算*/ l=1i; for(j=0;jsize_x;j+= 2*l ) /*一组蝶形运算*/ for(k=0;kl;k+) /*一个蝶形运算*/ mul(xj+k+l,wsize_x*k/2/l,&product); add(xj+k,product,&up); sub(xj+k,product,&down); xj+k=up; xj+k+l=down; /*快速傅里叶逆变换*/void ifft() int i=0,j=0,k=0,l=size_x; complex up,down; for(i=0;i (int)( log(size_x)/log(2) );i+) /*一级蝶形运算*/ l/=2; for(j=0;jsize_x;j+= 2*l ) /*一组蝶形运算*/ for(k=0;kl;k+) /*一个蝶形运算*/ add(xj+k,xj+k+l,&up); up.real/=2;up.img/=2; sub(xj+k,xj+k+l,&down); down.real/=2;down.img/=2; divi(down,wsize_x*k/2/l,&down); xj+k=up; xj+k+l=down; change();/*输出傅里叶变换的结果和幅值结果*/void output() int i; clrscr() ; printf(the result are as followsnn); for(i=0;i=0.0001)printf(+%.4fjn,xi.img); else if(fabs(xi.img)0.0001)printf(n); else printf(%.4fjn,xi.img); for( z=0;z16;z+) wz=sqrt(xz.real)*(xz.real)+(xz.img)*(xz.img); printf(nnthe amplitudes are as followsnn); for(z=0;zsize_x;z+) printf(%.4fn,wz); /*输出傅里叶变换相位的结果*/void output2() int i; setbkcolor(0); for(i=0;isize_x;i+) for( z=0;z16;z+) pz=atan(xz.img)/(xz.real); printf(nnnthe phase are as followsnn); for(z=0;zsize_x;z+) setcolor(1); printf(%.4fn,pz); void draw(int n) /* 作幅值图 */ int k,driver,mode,step,x,y; int i,j,xstep=512/n,ystep=55; char buffer10,bufferan; step=512/n; driver=vga; mode=vgahi; x=64; y=400; initgraph(&driver,&mode,); setbkcolor(15); /* 绘制坐标及刻度 */ setcolor(1); moveto(200,30); outtext(the figure of amplitude is as follow!); line(x,80,x,400); /* 画纵坐标*/ line(x-4,100,x,80); line(x,80,x+4,100); line(40,y,600,y); /*画横坐标*/ line(600,y,580,y+4); line(600,y,580,y-4); j=400; settextjustify(center_text,center_text); /*纵坐标刻度*/ for(i=0;i=25;i+=5) line(64,j,54,j); itoa(i,buffer,10); outtextxy(30,j,buffer); j-=ystep; setcolor(1); outtextxy(x,70,amplitude);/*纵坐标标注*/ j=64; settextjustify(center_text,top_text); /*画横坐标刻度*/ for(i=0;i=15;i+) line(j,400,j,410); itoa(i,buffer,10); outtextxy(j,415,buffer); j+=xstep; outtextxy(600,y+10,xn);/*横坐标标注*/ /* 绘制频谱图 */ for(k=0;k=0) outtextxy(x-step+2,(y-(wk)*12)-10,buffera); else outtextxy(x-step+2,(y-(wk)*12)+10,buffera); getch(); cleardevice(); void draw2(int n) /* 作相位图 */ int k,driver,mode,step,x,y; float kd; int i,j,xstep=512/n,ystep=30; char buffer10,bufferpn,bufferkdn; step=512/n; driver=vga; mode=vgahi; x=64; y=250; initgraph(&driver,&mode,); setbkcolor(15); /* 绘制坐标及刻度 */ setcolor(1); moveto(200,30); outtext(the figure of phase is as follow!); line(x,50,x,450); /* 画纵坐标*/ line(x-4,70,x,50); line(x,50,x+4,70); line(40,y,600,y); /*画横坐标*/ line(600,y,580,y+4); line(600,y,580,y-4); j=250; settextjustify(center_text,center_text); /*纵坐标刻度*/ outtextxy(15,j,0); for(kd=0.5;kd=3.0;kd+=0.5) line(64,j-ystep,54,j-ystep); sprintf(bufferkd,%.4f,kd); outtextxy(30,j-ystep,bufferkd); j-=ystep; j=250; for(kd=0.5;kd=3.0;kd+=0.5) line(64,j+ystep,54,j+ystep); sprintf(bufferkd,%.4f,-kd); outtextxy(30,j+ystep,bufferkd); j+=ystep; setcolor(1); outtextxy(x,40,phase);/*纵坐标标注*/ j=64; settextjustify(center_text,top_text); /*画横坐标刻度*/ for(i=0;i=15;i+) line(j,240,j,250); itoa(i,buffer,10); outtextxy(j+10,265,buffer); j+=xstep; outtextxy(600,y+10,xn);/*横坐标标注*/ /* 绘制相位图 */ for(k=0;k=0) outtextxy(x-step+2,(y-(pk)*60)-10,bufferp); else outtextxy(x-step+2,(y-(pk)*60)+5,bufferp); int main() int i,method; void fft(); /*快速傅里叶变换*/ void ifft(); void initw(); /*初始化变换核*/ void change(); /*变址*/ void add(complex ,complex ,complex *); /*复数加法*/ void mul(complex ,complex ,complex *); /*复数乘法*/ void sub(complex ,complex ,complex *); /*复数减法*/ void divi(complex ,complex ,complex *);/*复数除法*/ void output(); void output2(); /*输出结果*/ system(cls); pi=atan(1)*4; printf(please input the size of x:n); scanf(%d,&size_x); printf(please input the data in xn:n); for(i=0;isize_x;i+) scanf(%lf%lf,&xi.real,&xi.img); initw(); printf(use fft(0) or ifft(1)?n); scanf(%d,&method); if(method=0) fft(); else ifft(); output(); getch();/* 画频谱图 */ draw(size_x); output2(); getch(); draw2(size_x); getch(); closegraph();3.2 程序运行结果：please input the size of x:16please input the data in xn：xn = 0 1.2200 2.2300 2.8500 2.9800 2.6000 1.7600 0.6200 -0.6200 -1.7600 -2.6000 -2.9800 -2.8500 -2.2300 -1.2200 004.49118667428684 - 22.5787201341516i-0.863137084989848 +