现代信号处理

1、 现代信号处理课程设计实验报告 实验课题：现代信号处理 专业班级： 学生姓名： 学生学号： 指导老师： 完成时间： 目录1 前言-22 课程设计内容要求及题目-33 设计思想和系统功能结构及功能说明-44 关键部分的详细描述和介绍，流程图描述关键模块和设计思想-75 问题分析及心得体会-206 参考文献-217 附录：程序源代码清单-21 一、前言 数字滤波在通信、图像编码、语音编码、雷达等许多领域中有着十分广泛的应用。目前，数字信号滤波器的设计在图像处理、数据压缩等方面的应用取得了令人瞩目的进展和成就。它是数字信号处理理论的一部分。数字信号处理主要是研究用数字或符号的序列来表示信号

2、波形，并用数字的方式去处理这些序列，以便估计信号的特征参量，或削弱信号中的多余分量和增强信号中的有用分量。具体来说，凡是用数字方式对信号进行滤波、变换、调制、解调、均衡、增强、压缩、固定、识别、产生等加工处理，都可纳入数字信号处理领域。 数字信号处理学科的一项重大进展是关于数字滤波器设计方法的研究。关于数字滤波器，早在上世纪40年代末期就有人讨论设计它的可能性问题，在50年代也有人讨论过数字滤波器，但直到60年代中期，才开始形成关于数字滤波器的一整套完整的正规理论。在这一时期，提出了各种各样的数字滤波器结构，有的以运算误差最小为特点，有的则以运算速度高见长，而有的则二者兼而有之。出现

3、了数字滤波器的各种实现方法，对递归和非递归两类滤波器作了全面的比较，统一了数字滤波器的基本概念和理论。  数字滤波器与模拟滤波器相比，具有精度高、稳定、体积小、重量轻、灵活、不要求阻抗匹配以及能实现模拟滤波器无法进行的特殊滤波等优点。上学期学习了数字信号处理这门课，这学期的课程设计使我更加形象具体的掌握这门课程，并且可以熟练的运用MATLAB进行编程，经过此次课程设计达到了其设计内容的基本要求，实现了基本功能，并对程序进行了优化，以达到最佳效果。在设计前制定一个设计方案，明确设计目标，方法和过程，参阅相关文献，咨询老师和同学，严格按照方案进行课程设计，最终解决问题，完成任务要求。 二

4、、课程设计内容要求及题目1、信号发生器用户根据测试需要，可任选以下两种方式之一生成测试信号：（1）直接输入（或从文件读取）测试序列；（2）输入由多个不同频率正弦信号叠加组合而成的模拟信号公式（如式所示）、采样频率（Hz）、采样点数，动态生成该信号的采样序列，作为测试信号。 px=100*sin(2*pi*15*t)+100*sin(2*pi*20*t)+100*sin(2*pi*30*t) 2、频谱分析 使用FFT对产生的测试信号进行频谱分析并展示其幅频特性与相频特性，指定需要滤除的频带，通过选择滤波器类型（IIR/FIR），确定对应的滤波器（低通、高通）技术指标。3、滤波器设计 根据以上技术

5、指标（通带截止频率、通带最大衰减、阻带截止频率、阻带最小衰减），设计数字滤波器，生成相应的滤波器系数，并画出对应的滤波器幅频特性与相频特性。（1）IIR DF设计：可选择滤波器基型（巴特沃斯或切比雪夫型）；（2）FIR DF设计：使用窗口法（可选择窗口类型，并比较分析基于不同窗口、不同阶数所设计数字滤波器的特点）。4、数字滤波 根据设计的滤波器系数，对测试信号进行数字滤波，展示滤波后信号的幅频特性与相频特性，分析是否满足滤波要求（对同一滤波要求，对比分析各类滤波器的差异）。（1）IIRDF：要求通过差分方程迭代实现滤波（未知初值置零处理）；（2）FIRDF：要求通过快速卷积实现滤波（对于长序列

6、，可以选择使用重叠相加或重叠保留法进行卷积运算）。5、选做内容将一段语音作为测试信号，通过频谱展示和语音播放，对比分析滤波前后语音信号的变化，进一步加深对数字信号处理的理解。三、设计思想和系统功能结构及功能说明1、设计思想 首先设定由多个不同频率正弦信号叠加组合而成的模拟信号函数、决定采样频率（Hz）和采样点数，然后使用FFT对产生的测试信号进行频谱分析并展示其幅频特性与相频特性，指定需要滤除的频带，接下来使用汉宁窗低通滤波器与巴特沃斯高通滤波器进行滤波，分析滤波后的信号的幅频特性以及相频特性，最后设计语音滤波，检查是否满足滤波要求。2、系统结构功能说明（1）信号发生器 定义一个由4个不同频率

7、叠加的正弦函数，输出该函数，观察其相频以及幅频特性。（2）频谱分析 使用FFT对产生的测试信号进行频谱分析并展示其幅频特性与相频特性，指定需要滤除的频带，通过选择滤波器类型汉宁窗和巴特沃斯滤波器，确定对应的滤波器（低通、高通）技术指标Fp、Fc、Rp、Rs。然后对其滤波后信号进行频谱分析。（3）数字滤波器设计1、FIR 数字滤波器设计原理 FIR 数字滤波器可以用下面的差分方程来描述其中x（n）是输入序列，y（n）是输出序列，h（n）是系统的单位脉冲响应。由于其采用卷积的方法得到输出，故可以采用快速算法FFT进行卷积计算，提高信号处理速度。FIR 数字滤波器满足线性相位的条件是其乘法器的系数关

8、于中点对称，又因为其结构无反馈的回路，所以是无条件的稳定系统。2、 窗函数法设计及MATLAB 实现FIR 滤波器目前常用的设计方法有窗函数法和频率采样法，窗函数法是从时域进行设计，而频率采样法是从频域进行设计。窗函数法由于简单、物理意义清晰，因而得到了较为广泛的应用。窗函数法设计的基本思想是：首先根据技术指标要求，选取合适的阶数N和窗函数的类型，使其幅频特性逼近理想滤波器幅频特性。其次，因为理想滤波器的是无限长的，所以需要对进行截断，数学上称这种方法为窗函数法1。其中是希望逼近的滤波器的单位脉冲响应，是希望逼近的滤波器的幅频特性，是所设计的滤波器的单位脉冲响应，是窗函数。在截断的过程中要保证

9、FIR滤波器的稳定性与线性相位的特性。由于截断效应使所设计滤波器的幅频特性存在误差，称之为吉布斯效应。最后，验证所设计的滤波器是否满足要求，若不满足，重新设计。常用的窗函数有矩形窗、三角窗、汉宁窗和凯塞窗等。窗函数的选择一般要满足以下三个条件：（1）具有较低的旁瓣幅度，尤其是第一旁瓣幅度。（2）旁瓣幅度下降速度要快，以利于增加阻带衰减。（3）主瓣宽度要窄，以获得较陡的过渡带。FIR 滤波器的设计任务是选择有限长度的 h(n), 使函数H(ej)满足技术要求。我选择的是汉宁窗低通滤波。3、 IIR数字滤波器的设计IIR数字滤波器冲激响应无限宽，与模拟滤波器相匹配。因此可采取在模拟滤波器设计的基础

10、上进一步变换的方法设计IIR滤波器，即通过已知的模拟滤波器系统的系统函数H(s)来设计数字滤波器的系统函数H(z)，主要是通过脉冲响应不变法或双线性变换法完成s平面到Z平面的变换。通过典型的模拟滤波器（如巴特沃斯滤波器、切比雪夫滤波器等）可实现一定参数要求的数字滤波器。IIR数字滤波器的设计具体步骤：<1> 按照既定规则将数字滤波器技术指标转为模拟低通滤波器技术指标。 <2> 以转换后的技术指标为依据设计模拟低通滤波器H(s)（传递函数) 。<3> 再按照一定规则将H(s)转换成H(z)（数字滤波器的传递函数）。若设计的是高通、带通或带阻滤波器，还需进行以下

11、步骤：将高通、带通或带阻数字滤波器的技术指标转换为低通模拟滤波器的技术指标, 然后设计出低通H(s), 再将H(s)转换为H(z)。我选择的是巴特沃斯高通滤波器。（4）语音滤波1、语音信号的采集 MATLAB软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。通过wavread函数的使用，我们很快理解了采样频率、采样位数等概念。 这里我直接采用了一段现成的wav格式的语音信号。2、语音信号的频谱分析 首先画出语音信号的时域波形；然后对语音号进行快速傅里叶变换，得到信号的频谱特性。四、关键部分的详细描述和介绍，流程图描述关键模块和设计思想1、信号发生器（1）输入信号程序

12、如下：N=512;T=0.01;Fs=1/T;t=0:pi/256:2*pi;f1=7;f2=37;f3=39;f4=43;x=100*sin(2*pi*f1*t)+100*sin(2*pi*f2*t)+100*sin(2*pi*f3*t)+100*sin(2*pi*f4*t);plot(t,x,'g');grid;title('原始输入信号');xlabel('时间/t');ylabel('振幅');X=fft(x,N);plot(X);grid;title('输入信号频谱');设计输入信号并对其进行傅里叶变换。

13、输入信号如下图（2）输入信号频谱分析幅频特性Fs=100;len=512;f=Fs*(0:len/2-1)/len;plot(f,abs(X(1:len/2);gridtitle('原始信号幅频特性曲线');相频特性Fs=100;len=512;f=Fs*(0:len/2-1)/len;plot(f,angle(X(1:len/2);grid;title('相频特性曲线');幅频特性相频特性2、巴特沃斯高通滤波器 高通滤波器采用Butterworth滤波器, 因为它具有较平坦的幅频特性, 而且其元件值比较合乎实际、对精度要求稍宽。用N,Wn=buttord(Wp

14、,Ws,Rp,Rs,'s')语句实现Butterworth滤波器阶数和边缘频率的确定；b,a=butter(N,Wn,'s')语句来实现Butterworth滤波器的设计。其中Wp为通带截止频率；Rp为通带截止频率处的衰减；Ws为阻带截止频率；Rs为阻带截止频率处的衰减； b,a是滤波器传递函数多项式的系数；n为滤波器的阶数；Wn为截止频率。该通带截止频率为14HZ。（1）程序如下：T=0.01;Fs=1/T;wp=0.4*pi;ws=0.3*pi;Rp=1;Rs=15;Nn=512;%采样间隔Wp=wp/T;Ws=ws/T;N,Wn=buttord(Wp,Ws

15、,Rp,Rs,'s');%模拟滤波器阶数z,p,k=buttap(N);%数字滤波器设计Bap,Aap=zp2tf(z,p,k);b,a=lp2hp(Bap,Aap,Wn);bz,az=bilinear(b,a,1/T);%脉冲响应不变法figurefreqz(bz,az,Nn);H,F=freqz(bz,az,Nn,1/T);plot(F,20*log10(abs(H),'r');grid;title('巴特沃斯滤波器幅频特性');xlabel('频率/HZ');ylabel('振幅/dB');巴特沃斯幅频特性（

16、2） 滤波后信号巴特沃斯高通滤波器滤除了7HZ的信号，保留了高于14HZ的信号。滤波后信号幅频特性相频特性3、汉宁窗低通滤波器通带截止频率为16HZ，滤除16HZ以上的信号，只保留7HZ的信号。程序如下：wp=0.3*pi;ws=0.45*pi;wdelta=ws-wp;N=ceil(8*pi/wdelta);if rem(N,2)=0 N=N+1;end;Nw=N;wc=(wp+ws)/2;n=0:N-1;alpha=(N-1)/2;%相位延迟m=n-alpha+0.00001;%0.00001为MATALB系统的精度hd=sin(wc*m)./(pi*m);%理性滤波器脉冲响应win=ha

17、nn(Nw);h=hd.*win'%时域乘积对应频域卷积b=h;H,f=freqz(b,1,512,100);%100HZ的采样频率plot(f,20*log10(abs(H),'r');grid;xlabel('频率/HZ');ylabel('振幅/dB');title('汉宁窗幅频特性曲线');汉宁窗幅频特性滤波后信号幅频特性相频特性4、 语音滤波器利用函数wavread对语音信号进行采样，利用巴特沃斯滤波器进行滤波。程序如下：fs=32768;x1=wavread('C:UsersloryDesktopwo.

18、wav');t=0:1/32768:(size(x1)-1)/32768;Au=0.03;d=Au*cos(2*pi*10000*t)'x2=x1+d;wp=0.25*pi;ws=0.3*pi;Rp=1;Rs=15;Fs=32768;Ts=1/Fs;wp1=2/Ts*tan(wp/2); %将模拟指标转换成数字指标ws1=2/Ts*tan(ws/2); N,Wn=buttord(wp1,ws1,Rp,Rs,'s');%选择滤波器的最小阶数Z,P,K=buttap(N); %创建butterworth模拟滤波器Bap,Aap=zp2tf(Z,P,K);b,a=lp

19、2lp(Bap,Aap,Wn); bz,az=bilinear(b,a,Fs); %用双线性变换法实现模拟滤波器到数字滤波器的转换H,W=freqz(bz,az); %绘制频率响应曲线figure(1)plot(W*Fs/(2*pi),abs(H)gridxlabel('频率Hz')ylabel('频率响应幅度')title('Butterworth')f1=filter(bz,az,x2);figure(2)subplot(2,1,1)plot(t,x2) %画出滤波前的时域图title('滤波前的时域波形');subplot(

20、2,1,2)plot(t,f1); %画出滤波后的时域图title('滤波后的时域波形');sound(f1,32768); %播放滤波后的信号F0=fft(f1,1024);f=fs*(0:511)/1024;figure(3)y2=fft(x2,1024);subplot(2,1,1);plot(f,abs(y2(1:512); %画出滤波前的频谱图title('滤波前的频谱')xlabel('Hz');ylabel('fuzhi');subplot(2,1,2)F1=plot(f,abs(F0(1:512); %画出滤波后的

21、频谱图title('滤波后的频谱')xlabel('Hz');ylabel('fuzhi');巴特沃斯幅频特性 五、问题分析及心得体会 经过这次课程设计，我的收获良多，也遇到了很多问题。在这之前我们并没有接触过MATLAB，所以在课程设计时查阅了大量有关MATLAB的书籍以及与课程相关的设计案例。但是编程过程中总是出现一些错误，变量的定义，参数的传递，函数的调用等各方面大大小小的问题层出不穷，为此，只有仔细认真、不急不躁才能解决这些小问题，这次课程设计我分成了四个模块进行编程，信号发生器，巴特沃斯高通滤波器，汉宁窗低通滤波器，还有语音滤波部分。在

22、巴特沃斯沃斯滤波时由于参数取的不好，有滤波残留，改一下参数，滤波效果很好。开始时幅频特性与相频特性是对称分布的，通过加入几行代码，让它只显示一半。做语音滤波时，读不到我的语音文件，后来通过音频转换器把它的格式转换成了wav格式的文件后，就可以正常运行了。最后，在做GUI界面时，因为所需要出的图太多，有很多个按钮，后来我把各个模块分成了四部分，再建立一个总的调用函数来调出各个模块。形成了主界面与子界面的模式。画面看着直观，清晰。这次的课程设计，使我对MATLAB软件运用能力的有所提高，并且将所学的数字信号处理知识应用在课程设计之中，提升了对数字信号处理技术的掌握。参数的使用，函数的调用，信号的产

23、生，频谱分析，滤波器参数的设定及不同滤波器的设计都要注意，避免细节错误。经过这次课程设计，提高了我对问题的处理能力，有了问题，在经过不断的努力和老师与同学的帮助下成功的解决了。我不仅在数字信号处理知识的理解、设计思想的构建有了进一步的提升，也对MATLAB数学处理软件有了一定的认识和掌握，并认识到编程习惯的重要性。锻炼了我的动手能力，丰富了我的思维模式。我们要尽最大努力完成任务。 六、参考文献1 S.K.Mitra.Digital Signal Processing: A Computer Based Approach, 3rdEditionM,New York, USA:McGraw-Hil

24、l,20002 R.G.Lyons.Understanding Digital Signal Processing,2nd EditionM.New Jersey, USA:Prentice Hall,20053 程佩青.数字信号处理教程,第二版M.北京:清华大学出版社,20014 赵树杰等.数字信号处理M.西安:西安电子科技大学出版社,19975 丁玉美等.数字信号处理时域离散随机信号处理M.西安:西安电子科技大学出版社,20026 陈怀琛等.MATLAB及在电子信息课程中的应用M,北京:电子工业出版社出版,2002 七、程序源代码清单N=512;T=0.01;Fs=1/T;t=0:pi/2

25、56:2*pi;f1=7;f2=37;f3=39;f4=43;x=100*sin(2*pi*f1*t)+100*sin(2*pi*f2*t)+100*sin(2*pi*f3*t)+100*sin(2*pi*f4*t);figure;subplot(411)plot(t,x,'g');grid;title('原始输入信号');xlabel('时间/t');ylabel('振幅');X=fft(x,N);subplot(412)plot(X);grid;title('输入信号频谱');Fs=100;len=512;f=

26、Fs*(0:len/2-1)/len;subplot(413)plot(f,angle(X(1:len/2);grid;title('相频特性曲线');Fs=100;len=512;f=Fs*(0:len/2-1)/len;subplot(414)plot(f,abs(X(1:len/2);gridtitle('原始信号幅频特性曲线');T=0.01;Fs=1/T;wp=0.4*pi;ws=0.3*pi;Rp=1;Rs=15;Nn=512;%采样间隔Wp=wp/T;Ws=ws/T;N,Wn=buttord(Wp,Ws,Rp,Rs,'s');%模拟

27、滤波器阶数z,p,k=buttap(N);%数字滤波器设计Bap,Aap=zp2tf(z,p,k);b,a=lp2hp(Bap,Aap,Wn);bz,az=bilinear(b,a,1/T);%脉冲响应不变法figurefreqz(bz,az,Nn);H,F=freqz(bz,az,Nn,1/T);figure;subplot(411)plot(F,20*log10(abs(H),'r');grid;title('巴特沃斯滤波器幅频特性');xlabel('频率/HZ');ylabel('振幅/dB');N=512;T=0.01;

28、t=0:pi/256:2*pi;f1=7;f2=37;f3=39;f4=43;x=100*sin(2*pi*f1*t)+100*sin(2*pi*f2*t)+100*sin(2*pi*f3*t)+100*sin(2*pi*f4*t);y2=filtfilt(bz,az,x);subplot(412)plot(t,y2,'r');title('滤波后信号');xlabel('时间/s');grid;legend('filtfilt');Fs=100;len=200;Y1=fft(y2,len);f=Fs*(0:len/2-1)/le

29、n;subplot(413)plot(f,abs(Y1(1:len/2);grid;title('滤波后信号幅频特性曲线');Fs=100;len=512;f=Fs*(0:len/2-1)/len;subplot(414)plot(f,abs(X(1:len/2);gridtitle('原始信号幅频特性曲线');wp=0.3*pi;ws=0.45*pi;wdelta=ws-wp;N=ceil(8*pi/wdelta);if rem(N,2)=0 N=N+1;end;Nw=N;wc=(wp+ws)/2;n=0:N-1;alpha=(N-1)/2;%相位延迟m=n-

30、alpha+0.00001;%0.00001为MATALB系统的精度hd=sin(wc*m)./(pi*m);%理性滤波器脉冲响应win=hann(Nw);h=hd.*win'%时域乘积对应频域卷积b=h;H,f=freqz(b,1,512,100);%100HZ的采样频率figure;subplot(411)plot(f,20*log10(abs(H),'r');grid;xlabel('频率/HZ');ylabel('振幅/dB');title('汉宁窗幅频特性曲线');N=512;T=0.01;t=0:pi/256:

31、2*pi;f1=7;f2=37;f3=39;f4=43;x=100*sin(2*pi*f1*t)+100*sin(2*pi*f2*t)+100*sin(2*pi*f3*t)+100*sin(2*pi*f4*t);y4=filtfilt(b,1,x);subplot(412)plot(t,y4,'r');grid;title('FIR滤波后输出信号')legend('filtfilt');xlabel('时间/t');ylabel('振幅');Fs=100;len=200;Y2=fft(y4,len);f=Fs*(0:len/2-1)/len;subplot(413)plot(f,abs(Y2(1:len/2);grid;title('滤波后信号幅频特性曲线');global Y2subplot(414)plot(angle(Y2);grid