多频窄带数字信号处理仿真系统设计--毕业论文

1、多频窄带数字信号处理仿真系统设计 多频窄带数字信号处理仿真系统设计摘 要本文主要是实现在TC2.0的环境下利用C语言实现多频窄带数字信号处理仿真系统。系统功能包括信号采集、信号时域及频域分析、FIR滤波器设计、数字信号滤波，在本文中主要设计四种滤波系统：低通滤波系统、高通滤波系统、带通滤波系统、带阻滤波系统。在FIR滤波器设计时使用的是窗函数设计的方法，涉及四种窗函数，分别是矩形窗、哈明窗、汉宁窗、布莱克曼窗。最后本文对设计的系统程序进一步改进，使之可以完成信号模拟到数字的转换功能、根据数字滤波器指标自动完成滤波器设计的功能、四种滤波功能以及个步骤的时域频域图形绘制。关键词：多频窄带数字信号处

2、理仿真系统、FIR滤波器、C语言、TC2.0前言: 本次数字信号处理小项目主要是利用数字信号的理论结合编程实现一个多频窄带数字信号处理仿真系统设计。MATLAB强大的数字图像处理能力为我们提供了很好的处理环境，但是由于其处理都是集成的，给我们的使用带来一定的不便性，对于我们理解信号处理过程和定制我们需要的信号处理软件没有参考价值。故我们选择利用tc2.0环境下的C语言进行系统设计与编写，预期实现信号的采集、信号时域及频域分析、FIR滤波器设计、数字信号滤波等功能的多频窄带数字信号处理软件仿真系统。正文1系统整体设计思想本报告主要是用C语言在tc2.0环境下实现信号采集、信号时域及频域分析、FI

3、R滤波器设计、数字信号滤波等功能的多频窄带数字信号处理软件仿真系统。由此可知，该报告设计需包括以下系统:A模拟信号采集模块：完成对模拟信号的采样、量化，使信号变成数字信号。B 显示模块：完成对信号的时域和频域显示，以便于对信号的频谱和时域分布进行分析。C 滤波器模块：根据模拟信号生成的数字信号的频域特点确定要完成的滤波指标，进而按照一定的原则选择合适的滤波器进行最后要完成的滤波过程。D 滤波模块：系统的最后一步，完成对要求滤除的信号频率的滤除和有用频率分量的保留。具体的系统整体设计流程见图1-1系统设计框图:图1-1系统设计框图2 各模块设计2.1 模拟信号采集模块本项目的研究的原始信号是模拟

4、信号，而计算机只能处理二进制编码形式的数字信号，故模拟信号与数字信号之间需有一个误差比较小或者无误差的转换过程，而这个过程就是通过我们的模拟信号采集模块实现的。由理论知识可知，这个转变过程至少需要满足两个条件，一是数字信号可以无失真的还原原模拟信号，二是数字信号的频域相对于模拟信号的频域发生的改变。无失真的还原原模拟信号就要满足采样频率fs大于等于模拟信号的最高频率以为例，采样后的数字信号为,。故编程实现即为：void Signal_low_high(int n)/*高通和低通信号滤波测试信号函数*/ int i;for ( i = 0; i n; i+)ARRAYIi = 0.0;ARRAY

5、Ri = (float)cos(2 * PI*i*f1 / fs) + (float)cos(2 * PI*i*f2 / fs); 程序中ARRAYR,ARRAYI数组表示采集数字信号的实部和虚部。2.2 显示模块为了实现信号的采样滤波过程的图形显示，本次项目选择用TC2.0编译环境实现对信号的时域和频域显示。考虑到信号包含实部和虚部，因此程序编写时为了方便，选择将一个信号用实部和虚部表示。因此显示的时候就出现了两种情况，时域显示和频域显示，信号不同域波形需要不同的求模函数实现对信号强度的表示，具体程序函数声明为：/*取模函数，主要用来对时域信号的取模*/void GetModu1(float

6、 a,float b,int ll)/*取模函数，主要用来对频域信号的取模*/void GetModu2(float a, float b,int NNN)在具体的DOS环境下的图形显示，考虑系统功能实现问题，设计两个绘图函数，GRAPH用于专门显示窗函数的衰减曲线，GRAPH1通用绘图函数为显示步骤中函数的时域和频域，具体见程序：/*画图程序，用以画窗函数的时域图*/void GRAPH(float xwMAX,int nn,int kk)/*xw表示窗函数，nn表示描点点数，kk表示滤波器类型*/*画图程序，用以画函数的时域和频域的图型*/void GRAPH1(float a,int n

7、,int ii,int high)/*n表示描点点数，ii表示步长，high表示所控制的画图函数的高度*/2.3 滤波器模块本项目选择使用FIR滤波器设计方法。FIR滤波器通常采用窗函数方法来设计。窗设计的基本思想是，首先选择一个适当的理想选频滤波器（它总是具有一个非因果，无限持续时间脉冲响应），然后截取（加窗）它的脉冲响应得到线性相位和因果FIR滤波器。因此这种方法的重点是选择一个合适的窗函数和理想滤波器。我们用表示理想的选频滤波器，它在通带上具有单位增益和线性相位，在阻带上具有零响应。一个带宽的低通滤波器由下式给定：为了从得到一个FIR滤波器，必须同时在两边截取。而要得到一个因果的线性相位

8、滤波器，它的长度为，必须有：这种操作叫做加窗，可以看作是与窗函数的乘积：其中 根据的不同定义，可以得到不同的窗结构。在频域中，因果FIR滤波器响应由和窗响应的周期卷积得到，即： 常用的窗函数有矩形窗、巴特利特（BARTLETT）窗、汉宁（HANNING）窗、海明（HAMMING）窗、布莱克曼（BLACKMAN）窗、凯泽（KAISER）窗等。窗函数设计的具体步骤见流程图2-1：图2-1 窗函数设计流程图图2-2 滤波器指标选择窗和理想响应流程本项目实现提供四个窗函数选择，分别是矩形窗、哈明窗、汉宁窗、布莱克曼窗，通过给定的数字滤波器指标可以求出wn(n)窗函数，wn(n)窗函数得到程序实现上见附

9、录程序中JXC(),HNC(),HMC(),BLKMC()四个窗函数子函数。具体选择的窗函数类型需根据as阻带最小衰减，ws阻带截止频率，wp通带截止频率,wsl、wsh、wpl、wph带通带阻滤波器的上下边带阻带通带截止频率确定，由指标选择窗函数的子程序见choose_window()子程序，具体见流程图2-2。理想单位脉冲函数Hd(n)获取上，因为我们设计的滤波器为低通、高通、带通、带阻四种，故需要四个能产生理想逼近单位响应函数的子函数。查阅资料可知，对于高通滤波器函数来说，其理想频率响应和理想单位脉冲响应函数分别为：对于带通滤波器函数来说，其理想频率响应和理想单位脉冲响应函数分别为：对于

10、带阻滤波器函数来说，其理想频率响应和理想单位脉冲响应函数分别为：所以，由公式可以编写具体的函数产生程序，具体函数的产生见附录子函数HDN(),HDN1(),HDN2(),HDN3()程序。滤波器的单位脉冲响应h(n)得到是窗函数设计滤波器的最后一步，又流程图3-1可知，h(n)=Wn(n)*Hd(n)，通过一个简单的for循环即可实现最后的相乘工作，即程序：for(i=0;in;i+)/*时域加窗处理*/ Hri=Hdi*wni; Hii=0.0; 2.4 滤波模块滤波模块是设计的多频窄带数字信号处理仿真系统的最后一个步骤。该步骤的设计流程图见图2-3。图2-3 滤波设计流程图流程图中的截断处

11、理是因为线性卷积过程中会在输出信号的两边形成过渡带，过渡带中包含的信号并不是我们想要的信号，而且过渡带中包含许多我们不想要的频率分量，因此为避免上述的过渡带带来的影响，我们对线性卷积的得到的结果进行截断处理，使得到的输出信号y(n)与原信号等长，这也满足我们对数字信号处理的基本要求。在程序的具体实现上，线性卷积程序为Cir_conv()子函数程序(具体见附录程序)，截断函数程序为for(I=0;INN;I+)/*对滤波信号进行截断处理*/ JieduanrI=DfCvARI+(N-1)/2; JieduaniI=DfCvAII+(N-1)/2; 在信号的分析过程中，对于时域频域之间的转换，我们

12、的设计思路是以DFT(程序见附录Dft子函数)变换为主，以FFT(程序见附录DFT_FFT子函数)变换为辅。因FFT程序运行速度快，但变换点数必须为2的幂次，因此主要用来完成窗函数的时域到频域的转换；其他时域频域之间的不定点转换则选择DFT实现。3 滤波实现3.1 低通滤波(程序见附录)测试信号函数为10Hz和45Hz的混合信号，采样频率为100Hz，采样时间为2.0s，测试信号DFT变换点数为200点。实现对45Hz高频信号的滤出，信号的时域波形和频域波形见图3-1-1与图3-1-2。图3-1-1 测试信号时域波形图3-1-2 测试信号频域波形运行结果图： 图3-1-3 滤波器指标输入图3-

13、1-4 窗函数时域波形图3-1-5 滤波未截断时域波形图3-1-6 滤波截断后时域波形图3-1-7 滤波截断后频域波形图3-1-8 窗函数衰减波形3.2 高通滤波(程序见附录)测试信号函数为10Hz和45Hz的混合信号，采样频率为100Hz，采样时间为2.0s，测试信号DFT变换点数为200点。实现对10Hz低频信号的滤出。信号的时域波形和频域波形见图3-1-1与图3-1-2。运行结果图为: 图3-2-1 滤波器指标输入图3-2-2滤波未截断时域波形图3-2-3滤波截断后时域波形图3-2-4滤波截断后频域波形图3-2-5 窗函数衰减波形3.3 带通滤波(程序见附录)测试信号函数为10Hz、30

14、Hz、45Hz的混合信号，采样频率为100Hz，采样时间为2.0s，测试信号DFT变换点数为200点。实现对10Hz和45Hz频率信号的滤出。信号的时域波形和频域波形见图3-3-1与图3-3-2。运行结果图为:图3-3-1 测试信号的时域波形图3-3-2 测试信号的频域波形 图3-3-3滤波器指标输入图3-3-4滤波未截断时域波形图3-3-5滤波截断后时域波形图3-3-6滤波截断后频域波形图3-3-7窗函数衰减波形3.4 带阻滤波(程序见附录)测试信号函数为10Hz、30Hz、45Hz的混合信号，采样频率为100Hz，采样时间为2.0s，测试信号DFT变换点数为200点。实现对30Hz频率信号

15、的滤出。信号的时域波形和频域波形见图3-3-1与图3-3-2。运行结果图为:图3-4-1滤波器指标输入图3-4-2滤波未截断时域波形图3-4-3滤波截断后时域波形图3-4-4滤波截断后频域波形图3-4-5窗函数衰减波形4 项目结论与感想 通过这次的项目研究，让我们加深了对C语言和数字信号处理理论的理解。同时，团队良好的配合是我们完成的这次项目完成的基础。我们不仅加深了基础知识的理解，更熟悉理解了一些我们在课堂上没有学习到的、没能理解的知识点。对于项目中遇到的问题，我们采取边查边改的方式，一边查询问题所在，一遍翻阅资料，查找改进方法，修改程序。附录：#include#include#includ

16、e#include#include#define PI 3.1415926#define MAX 600#define f1 10.0#define f2 45.0#define f3 30.0#define fs 100.0#define Tp 2int L,NN,k,N;float ARRAYRMAX, ARRAYIMAX,DfCvARMAX,DfCvAIMAX,MODULUSMAX,MODULUS1MAX,HrMAX,HiMAX,HdMAX,wnMAX;float wc,wl,wh;void Dft(float re, float im, int s,int N11);void GetM

17、odu1(float a, float b,int ll);void GetModu2(float a, float b,int NNN);void Signal_low_high(int n);void Signal_ee(int n);void Cir_conv(float h,float x,int h_size,int x_size);main()int i,m,n,fff,nn,kk,I,LL,lll;void HDN(),JXC(),HNC(),HMC(),BLKMC(),HDN1(),HDN2(),HDN3();float JieduanrMAX,JieduaniMAX,HkMA

18、X,HkxwMAX,arrayrMAX,arrayiMAX;void DFT_FFT();void Convolution_LLL();void GRAPH();void GRAPH1();void choose_window();m=7;n=(int)pow(2,m);/*DFT_FFT LENGTH*/L=Tp*fs;printf(n);printf(please choose what style of filter you wantn);printf(1-low-pass filter 2- high-pass filter 3-band pass filter,4-band stop

19、 filter*/;n);scanf(%d,&fff);if (fff=2) Signal_low_high(L);else Signal_ee(L);printf(your signal styles are %dn,L);printf(Please Input x(n)-DFT NN:n);scanf(%d, &NN);GetModu1(ARRAYR,ARRAYI,NN);GRAPH1(MODULUS1,NN,3,20);Dft(ARRAYR, ARRAYI, 1,NN);GetModu2(DfCvAR, DfCvAI,NN);GRAPH1(MODULUS,NN,3,5);choose_w

20、indow(fff,n);for(i=0;in;i+)Hri=Hdi*wni;Hii=0.0;GRAPH1(Hr,n,6,300);for(I = 0; I n; I+)arrayrI = HrI;arrayiI = HiI;DFT_FFT(Hr,Hi,m,n);for(i=0;in;i+)Hki=sqrt(pow(Hri,2)+pow(Hii,2);for(i=0;i=%d):n,(N+NN-1);scanf(%d, &lll);Cir_conv(arrayr,ARRAYR,N,NN);GetModu1(DfCvAR, DfCvAI,lll);GRAPH1(MODULUS1,lll,2,12

21、0);for(I=0;INN;I+)JieduanrI=DfCvARI+(N-1)/2;JieduaniI=DfCvAII+(N-1)/2;GRAPH1(Jieduanr,NN,2,80);Dft(Jieduanr,Jieduani,1,NN);GetModu2(DfCvAR, DfCvAI,NN);GRAPH1(MODULUS,NN,2,2);getch();void GRAPH(float xwMAX,int nn,int kk) int gdriver,mode,i;gdriver=DETECT;mode = 0;registerbgidriver(EGAVGA_driver);init

22、graph(&gdriver,&mode,C:Progra1WINYESTC20HBGI);setbkcolor(42);setcolor(30);moveto(30,30);line(30,30,619,30);line(30,30,30,430);for(i=0;inn;i+)xwi=abs(xwi);for(i=0;inn/2;i+)lineto(30+i*3,30+xwi*3);for(i=30;i240;i+=30)line(30,30+i,25,30+i);for(i=30+i;i0) wp=wph;ws=wsh;elsewp=wpl;ws=wsl;else printf(plea

23、se input ws:PI*);scanf(%f,&ws);printf(please input wp:PI*);scanf(%f,&wp);wc=(wp+ws)/2;printf(wc=%f*PIn,wc);printf(please input as:);scanf(%d,&as);if(as=21)k=1;printf(CHOSE Rectangle WINDOWn);N=1.8/(fabs(ws-wp);if(N%2=0)N=N+1;else N=N+2;printf(window size N=%dn,N);getch();else if(as=44)k=2;printf(CHO

24、SE hanning WINDOWn);N=6.2/(fabs(ws-wp);if(N%2=0)N=N+1;else N=N+2;printf(window size N=%dn,N);getch();else if(as=53)k=3;printf(CHOSE hanming WINDOWn);N=6.6/(fabs(ws-wp);if(N%2=0)N=N+1;else N=N+2;printf(window size N=%dn,N);getch();else if(as74)printf(can not realize,please quit.n);getch();switch(k)ca

25、se 1:JXC(n,N);break;case 2:HNC(n,N);break;case 3:HMC(n,N);break;case 4:BLKMC(n,N);break;default:printf(errorn);switch(fff)case 1:HDN(n,N);break;case 2:HDN1(n,N);break;case 3:HDN2(n,N);break;case 4:HDN3(n,N);break;default:printf(errorn);void GRAPH1(float a,int n,int ii,int high)int gdriver,mode,x,y,m

26、;gdriver=DETECT;mode = 0;registerbgidriver(EGAVGA_driver);initgraph(&gdriver,&mode,C:Progra1WINYESTC20HBGI);setbkcolor(42);setcolor(30);cleardevice();x=20;y=200;line(30,200,609,200);for(m=0;mn;m+)line(x+m*ii,y,x+m*ii,y-am*high);getch();setbkcolor(0);setcolor(15);closegraph();void HDN(int n,int N)int

27、 i;float a;a=(N-1)/2;for(i=0;iN;i+)if(i=a) Hdi=wc;else Hdi=sin(PI*wc)*(i-a)/(float)(PI*(i-a);for(i=N;in;i+)Hdi=0.0;void HDN1(int n,int N)int i;float a;a=(N-1)/2;for(i=0;iN;i+)if(i=a) Hdi=(float)(1.000-wc);else Hdi=sin(PI*(i-a)/(float)(PI*(i-a)-sin(PI*wc)*(i-a)/(float)(PI*(i-a);for(i=N;in;i+)Hdi=0.0;

28、void HDN2(int n,int N)int i;float a;a=(N-1)/2;for(i=0;iN;i+)if(i=a) Hdi=(float)(wh-wl);else Hdi=sin(PI*wh)*(i-a)/(float)(PI*(i-a)-sin(PI*wl)*(i-a)/(float)(PI*(i-a);for(i=N;in;i+)Hdi=0.0;void HDN3(int n,int N)int i;float a;a=(N-1)/2;for(i=0;iN;i+)if(i=a) Hdi=(float)(1.000+wl-wh);elseHdi=sin(PI*wl)*(i

29、-a)/(float)(PI*(i-a)+sin(PI*(i-a)/(float)(PI*(i-a)-sin(PI*wh)*(i-a)/(float)(PI*(i-a); for(i=N;in;i+)Hdi=0.0;void JXC(int n,int N)int i;for(i=0;in;i+)if(iN) wni=1.0;else wni=0.0;void HNC(int n,int N)int i;for(i=0;in;i+)if(iN) wni=0.54-0.46*cos(2*PI*i/(N-1);else wni=0.0;void HMC(int n,int N)int i;for(

30、i=0;in;i+)if(iN) wni=0.5*(1-cos(2*PI*i/(N-1);else wni=0.0;void BLKMC(int n,int N)int i;for(i=0;in;i+)if(iN) wni=0.42-0.5*cos(2*PI*i/(N-1)+0.08*cos(4*PI*i/(N-1);else wni=0.0;void GetModu1(float a,float b,int ll)int i;for (i = 0; i ll; i+)if (bi = 0)MODULUS1i = ai;elseMODULUS1i = (float)sqrt(pow(ai, 2

31、) + pow(bi, 2);void GetModu2(float a, float b,int NNN)int i;for (i = 0; i NNN; i+)MODULUSi = (float)sqrt(pow(ai, 2) + pow(bi, 2);void Signal_low_high(int n) int i;for ( i = 0; i n; i+)ARRAYIi = 0.0;ARRAYRi = (float)cos(2 * PI*i*f1 / fs) + (float)cos(2 * PI*i*f2 / fs);void Signal_ee(int n)int i;for ( i = 0; i n; i+)ARRAYIi = 0.0;ARRAYRi = (float)cos(2 * PI*i*f1 / fs) + (float)cos(2 * PI*i*f2 / f