数字信号处理课程设计

目录前言第一章离散傅里叶变换DFT (3)1.1 DFT定义 (3)1.2 DFT的快速算法——快速傅里叶变换(FFT) (3)第二章基2 DIT-FFT算法 (4)2.1 按时域抽取的基2 DIT-FFT算法 (4)第三章基于C语言设计16点基2DIT-FFT程序及运行结果 (6)3.1 按时间抽取的基2DIT-FFT程序 (6)3.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 算法21j 0()DFT[()]()e , 0, 1, , 1N k n N N n X k x n x n k N π--====-∑ 2j e N N W π-=10()DFT[()](), 0, 1, , 1N k n N N n X k x n x n W k N -====-∑ m N W (2)()0,1,1(21)2x r N x n r x r ⎧→=-⎨+⎩ (2)()0,1,1(21)2X m N X k m X m ⎧→=-⎨+⎩第二章基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 ，即10()DET[()]() 0,1,,1N k n N N n X k x n x n W k N -====-∑/21/212(21)00()()()(2)(21)k n k n N N n n N N k l k l N N l l X k x n W x n W x l W x l W --+===+=++∑∑∑∑ 2/2k l k l N N W W = 21/211/22/200()()() 0,1,,1M N k l k k l N N N l l X k x l W W x l W k N --===+=-∑∑ /2111/21/20()DFT[()]() 0,1,,12N k l N N l N X k x l x l W k -====-∑ /2122/22/2()DFT[()]() 0,1,,12N k l N N N X k x l x l W k -===-∑有上述公式，及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 程序及运行结果1212()()() 0,1,,12()()2k N k N X k X k W X k N k N X k X k W X k ⎫=+⎪=-⎬⎛⎫+=- ⎪⎪⎝⎭⎭设计参数及要求：本次试验要求利用C语言设计16点基2DIT-FFT程序，并对如下程序进行分析，给出输出结果并画出信号的幅值图和相位谱。

3.1 按时间抽取的基2DFT-FFT程序：/*时间抽选基2 FFT算法C语言实现*/#include <stdio.h>#include <math.h>#include <stdlib.h>#include <conio.h>#define N 1000#include<graphics.h>/*定义复数类型*/typedef struct{double real;double img;}complex;complex x[N], *W; /*输入序列,变换核*/int size_x=0; /*输入序列的大小，在本程序中仅限2的次幂*/double PI; /*圆周率*/int z;float w[N],p[N];char str[N];/*初始化变换核*/void initW(){int i;W=(complex *)malloc(sizeof(complex) * size_x);for(i=0;i<size_x;i++){W[i].real=cos(2*PI/size_x*i);W[i].img=-1*sin(2*PI/size_x*i);}}/*加法*/void add(complex a,complex b,complex *c){c->real=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;i<size_x;i++){k=i;j=0;t=(log(size_x)/log(2));while( (t--)>0 ){j=j<<1;j|=(k & 1);k=k>>1;}if(j>i){temp=x[i];x[i]=x[j];x[j]=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=1<<i;for(j=0;j<size_x;j+= 2*l ){ /*一组蝶形运算*/for(k=0;k<l;k++){ /*一个蝶形运算*/mul(x[j+k+l],W[size_x*k/2/l],&product); add(x[j+k],product,&up);sub(x[j+k],product,&down);x[j+k]=up;x[j+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;j<size_x;j+= 2*l ){ /*一组蝶形运算*/for(k=0;k<l;k++){ /*一个蝶形运算*/ add(x[j+k],x[j+k+l],&up);up.real/=2;up.img/=2;sub(x[j+k],x[j+k+l],&down);down.real/=2;down.img/=2;divi(down,W[size_x*k/2/l],&down);x[j+k]=up;x[j+k+l]=down;}}}change();}/*输出傅里叶变换的结果和幅值结果*/void output(){int i;clrscr() ;printf("The result are as follows\n\n");for(i=0;i<size_x;i++){printf("%.4f",x[i].real);if(x[i].img>=0.0001)printf("+%.4fj\n",x[i].img);else if(fabs(x[i].img)<0.0001)printf("\n");else printf("%.4fj\n",x[i].img);for( z=0;z<16;z++){w[z]=sqrt((x[z].real)*(x[z].real)+(x[z].img)*(x[z].img)); }}printf("\n\nThe Amplitudes are as follows\n\n");for(z=0;z<size_x;z++){printf("%.4f\n",w[z]);}}/*输出傅里叶变换相位的结果*/void output2(){int i;setbkcolor(0);for(i=0;i<size_x;i++){for( z=0;z<16;z++){p[z]=atan((x[z].img)/(x[z].real));}}printf("\n\n\nThe Phase are as follows\n\n");for(z=0;z<size_x;z++){setcolor(1);printf("%.4f\n",p[z]);}}void draw(int n) /* 作幅值图 */{int k,driver,mode,step,x,y;int i,j,xstep=512/n,ystep=55;char buffer[10],buffera[N];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,"x[n]");/*横坐标标注*//* 绘制频谱图 */for(k=0;k<n;k++){setcolor(4);setlinestyle(0,0,3);line(x,y,x,(y-(w[k])*11));x+=step;sprintf(buffera,"%.4f",w[k]);if(w[k]>=0)outtextxy(x-step+2,(y-(w[k])*12)-10,buffera); elseouttextxy(x-step+2,(y-(w[k])*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 buffer[10],bufferp[N],bufferkd[N];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,"x[n]");/*横坐标标注*//* 绘制相位图 */for(k=0;k<n;k++){setcolor(4);setlinestyle(0,0,3);line(x,y,x,(y-(p[k])*60));x+=step;sprintf(bufferp,"%.4f",p[k]);if(p[k]>=0)outtextxy(x-step+2,(y-(p[k])*60)-10,bufferp);elseouttextxy(x-step+2,(y-(p[k])*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 x[N]:\n");for(i=0;i<size_x;i++)scanf("%lf%lf",&x[i].real,&x[i].img);initW();printf("Use FFT(0) or IFFT(1)?\n");scanf("%d",&method);if(method==0)fft();elseifft();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 X[N]：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 04.49118667428684 - 22.5787201341516i-0.863137084989848 + 2.08379725676967i -0.708894177324789 + 1.06093511090429i -0.660000000000000 + 0.660000000000000i -0.667400913065437 + 0.445943033038161i -0.636862915010152 + 0.263797256769669i -0.634891583896613 + 0.126287787982296i -0.640000000000000-0.634891583896613 - 0.126287787982296i -0.636862915010152 - 0.263797256769669i -0.667400913065437 - 0.445943033038161i -0.660000000000000 - 0.660000000000000i -0.708894177324789 - 1.06093511090429i -0.863137084989848 - 2.08379725676967i 4.49118667428684 + 22.5787201341516i序列X(n):傅里叶变换后：第四章课程设计的总结通过课程设计，熟悉了DFT的性质。