数字信号处理三级项目

1、燕山大学信息科学与工程学院燕 山 大 学多频窄带数字信号处理仿真系统设计 学生所在学院：信息学院 学生所在班级：电子信息工程2班 课程名称：数字信号处理 学生姓名： 姜升殿、李林、刘鹏、左园叶、陈敏、黄品程指导教师：陈书贞 日期：2014.12.202014年12月多频窄带数字信号处理仿真系统设计姜升殿、李林、黄品程、刘鹏、左园叶、陈敏摘 要本文主要是实现在TC2.0的环境下利用C语言实现多频窄带数字信号处理仿真系统。系统功能包括信号采集、信号时域及频域分析、FIR滤波器设计、数字信号滤波，在本文中主要设计四种滤波系统：低通滤波系统、高通滤波系统、带通滤波系统、带阻滤波系统。在FIR滤波器设计

2、时使用的是窗函数设计的方法，涉及四种窗函数，分别是矩形窗、哈明窗、汉宁窗、布莱克曼窗。最后本文对设计的系统程序进一步改进，使之可以完成信号模拟到数字的转换功能、根据数字滤波器指标自动完成滤波器设计的功能、四种滤波功能以及个步骤的时域频域图形绘制。关键词：多频窄带数字信号处理仿真系统、FIR滤波器、C语言、TC2.0前言: 本次数字信号三级项目是教改后用于帮助同学理解数字信号的理论和提高处理解决问题的能力而设立的，主要是利用数字信号的理论结合编程实现一个多频窄带数字信号处理仿真系统设计。目前的软件解决数字信号处理的有很多，比较好的是MATLAB，其强大的数字图像处理能力为我们提供了很好的处理环境

3、，但是由于其处理都是集成的，同时给我们的使用带来一定的不便性，对于我们理解信号处理过程和定制我们需要的信号处理软件没有参考价值。故我们选择利用tc2.0环境下的C语言进行系统设计与编写，预期实现信号的采集、信号时域及频域分析、FIR滤波器设计、数字信号滤波等功能的多频窄带数字信号处理软件仿真系统。项目组分工：组长姜升殿，负责程序的编写和调试以及报告和汇报ppt的部分制作；组员李林，负责报告主要编写；组员陈敏，负责汇报ppt的主要制作；组员黄品程，负责部分程序的编写；组员左园叶，负责报告的理论部分整理；组员刘鹏，负责报告的理论部分整理。研究报告正文1系统整体设计思想本报告主要是用C语言在tc2.

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

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

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

7、用来对时域信号的取模*/void GetModu1(float 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表示滤波器类型*/*画图程序，用以画函数的时域和频

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

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

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

11、波器函数来说，其理想频率响应和理想单位脉冲响应函数分别为：对于带阻滤波器函数来说，其理想频率响应和理想单位脉冲响应函数分别为：所以，由公式可以编写具体的函数产生程序，具体函数的产生见附录子函数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 滤波模块滤波模块是设计的多频窄带数字信号处理仿真系统的最后一个步骤。该步骤的

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

13、-1)/2; 在信号的分析过程中，对于时域频域之间的转换，我们的设计思路是以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 测

14、试信号频域波形运行结果图： 图3-1-3 滤波器指标输入图3-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 窗函数衰减

15、波形3.3 带通滤波(程序见附录)测试信号函数为10Hz、30Hz、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.0

16、s，测试信号DFT变换点数为200点。实现对30Hz频率信号的滤出。信号的时域波形和频域波形见图3-3-1与图3-3-2。运行结果图为:图3-4-1滤波器指标输入图3-4-2滤波未截断时域波形图3-4-3滤波截断后时域波形图3-4-4滤波截断后频域波形图3-4-5窗函数衰减波形4 项目结论与感想项目开展至结束，时间紧迫，但收获颇丰。 这次我们采用分组合作来开展我们的课题，所谓“男女搭配干活不累”，在这次分组开展的研究中体现的淋漓尽致。男生热爱钻研，冲劲十足，是小组攻坚的信心来源和决断力来源；女生心思细腻，做事谨慎，在发现问题、解决问题以及后期材料整理方面贡献了重要力量。 通过这次的项目研究，让

17、我们加深了对C语言和数字信号处理理论的理解，使得之后的运用更加熟练。除此之外，讨论的过程让我们学到了别人知道但自己还不是很熟悉的知识，还纠正了一些理论知识学习中自己没有发现的错误。同时感谢老师的帮助，才能使我们在问题不能解决的情况下走出困境。参考文献1 郑君里.信号与系统M.北京：高等教育出版社，2000.2 高西全数字信号处理第三版 M西安：西安电子科技大学出版社，20013 林建秋 韩静萍. C语言程序人设计M . 北京：机械工业出版社，2004.1附录：/*本程序可以实现信号的低通、高通、带通、带阻的滤波*/#include#include#include#include#include

18、#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 2 /*信号截取时长*/int L,NN,k,N;/*信号频域采样点L，原信号DFT变换的点数为NN,k表示选择窗函数的类型，N表示窗函数的长度*/float ARRAYRMAX, ARRAYIMAX,DfCvARMAX,DfCvAIMAX, MODULUSMAX,MODULUS1MAX,HrMAX,HiMAX,HdMAX,wnMAX;/*ARRAYR,ARRAYI数

19、组表示原信号的实部和虚部，DfCvAR,DfCvAI表示DFT变换后的实部和虚部；MODULUS，MODULUS1表示时域的模值和频域的模值*/*Hr，Hi数组表示原信号的实部和虚部，Hd表示理想逼近的单位脉冲响应；wn表示窗函数数组*/float wc,wl,wh;/*过渡带中心频点*/void Dft(float re, float im, int s,int N11);/*DFT函数和IDFT函数,N11表示DFT和IDFT的点数*/void GetModu1(float a, float b,int ll);/*求原信号的幅值*/void GetModu2(float a, float

20、 b,int NNN);/*求DFT变换后的信号的幅值*/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;/*as表示阻带最小衰减,lll表示循环卷积的点数*/ void HDN(),JXC(),HNC(),HMC(),BLKMC(),HDN1(),HDN2(),HDN3();/*HDN,HDN1,

21、HDN2,HDN3分别表示低通,高通,带通,带阻理想逼近函数*/ float JieduanrMAX,JieduaniMAX,HkMAX,HkxwMAX,arrayrMAX,arrayiMAX;/*Hk表示设计的函数的频域模值，Hkxw表示转换为db形式的频域模值*/ 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 c

22、hoose what style of filter you wantn);/*选择设计滤波器类型*/ printf(1-low-pass filter 2- high-pass filter 3-band pass filter,4-band stop filter*/;n);/*1-低通；2-高通；3-带通；4-带阻*/ scanf(%d,&fff);if (fff=2)/*选择测试滤波器函数*/ Signal_low_high(L);/*低通高通信号函数*/else Signal_ee(L);/*带通带阻信号函数*/printf(your signal styles are %dn,L)

23、;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_window(fff,n); for(i=0;in;i+) Hri=Hdi*wni; Hii=0.0; /*时域加窗处理*/ GRAPH1(Hr,n,6,300);/

24、*画加窗后的时域图*/ 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,120);/*画滤波

25、后时域的波形*/ 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,in

26、t nn,int kk)/*xw表示窗函数，nn表示描点点数，kk表示滤波器类型*/ int gdriver,mode,i; gdriver=DETECT,mode; initgraph(&gdriver,&mode,C:TC20BGI); 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)

27、line(30,30+i,25,30+i); for(i=30+i;i0)/*选择窗函数的长度的估计指标*/ wp=wph;ws=wsh;elsewp=wpl;ws=wsl; else /*设计低通和高通滤波器的指标*/printf(please 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:);/*衰减的DB绝对值*/scanf(%d,&as);if(as=21)/*根据阻带最小衰减选

28、择窗函数类型*/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(CHOSE 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

29、)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) case 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(e

30、rrorn); 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)/*n表示描点点数，ii表示步长，high表示所控制的画图函数的高度*/ int gdriver,mode,x,y,m; gdriver=DETECT,mode; initgraph(&gdrive

31、r,&mode,C:TC20BGI); 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();/*低通理想逼近函数，n表示fft变换的点数，N表示窗函数的长度*/void HDN(int n,int N) int i; float a; a=(N-1)/2; for(i=0;iN;i+) if(i=a) Hdi=

32、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;void HDN2(int n,int N)/*

33、带通理想逼近函数*/ 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); else Hd

34、i=sin(PI*wl)*(i-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;/*n表示fft变换的点数，N表示窗函数的长度*/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/(

35、N-1); else wni=0.0; void HMC(int n,int N) int i; for(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;fo

36、r (i = 0; i ll; i+)if (bi = 0)MODULUS1i = ai;elseMODULUS1i = (float)sqrt(pow(ai, 2) + 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