数字信号处理数字滤波器论文.docx

1、数字信号处理数字滤波器论文1 绪 论1.1 课程设计目的1. 在学习了数字信号处理课程之后能把所学的理论知识与实际的工程应用联系起来。2. 在熟练掌握数字信号处理的基本原理基础上，能结合工程实际学习更多的数字信号处理实现技术及其应用技能。3. 通过课程设计树立正确的设计思想，提高学生分析问题、解决问题的能力。1.2 课程设计要求1. 掌握数字信号处理的基本概念，基本理论和基本方法。2. 熟悉离散信号和系统的时域特性。3. 掌握序列快速傅里叶变换方法。4. 学会MATLAB的使用，掌握MATLAB的程序设计方法。5. 掌握利用MATLAB对语音信号进行频谱分析。6. 掌握滤波器的网络结构7. 掌

2、握MATLAB设计IIR数字滤波器的方法和对信号进行滤波的方法。1.3 课程设计内容1. 设计一卷积演示程序： 可输入任意2待卷积序列x1(n)、x2(n)，长度不做限定。例如x1(n)=1,1,1,1,0,0,1,1,1,1,0,0，x2(n)=0,1,2,1,0,0,0,1,2,1,0,0；请采用自己设计的数据。 分别动态演示2序列进行线性卷积x1(n)x2(n)和圆周卷积x1(n)x2 (n)的过程；要求分别动态演示翻转、移位、乘积、求和的过程； 圆周卷积默认使用2序列中的最大长度，但卷积前可以指定卷积长度N用以进行混叠分析； 根据实验结果分析2类卷积的关系。 产生时域序列x1(n)、x

3、2(n)，长度不小于10000，序列内容自定义。利用FFT实现快速卷积，验证时域卷积定理，并与直接卷积进行效率对比（验证时采用matlab子函2. 编写程序演示采样定理（时域采样、频谱周期延拓），同时演示采样频率小于2fc时，产生的混叠效应：1 对下面连续信号采样， A为幅度因子，a为衰减因子， 为模拟角频率。要求输入采样频率fs（根据程序处理需要指定范围）后，在时域演示信号波形、采样脉冲及采样后信号；在频域演示对应的信号频谱、采样脉冲及频域周期拓延。3. 滤波设计及其应用 利用Windows下的录音机或其他软件，进行语音信号的采集（*.wav）； 语音信号的频谱分析，画出采样后语音信号的时域

4、波形和频谱图； 产生噪声信号并加到语音信号中，得到被污染的语音信号，并回放语音信号； 污染信号的频谱分析，画出被污染的语音信号时域波形和频谱； 根据有关的频谱特性，采用直接法设计IIR模拟、数字滤波器，并画出相应滤波器的幅频、相频图（设计6个IIR滤波器）；a. 滤波器类型：椭圆滤波器（低通、高通、带通）b 总体要求：Matlab原程序+仿真波形+技术指标 根据有关的频谱特性，采用直接法设计FIR数字滤波器，并画出相应滤波器的幅频；（设计3个FIR滤波器）a.滤波器类型：汉明窗（低通、高通、带阻）b总体要求：Matlab原程序+仿真波形+技术指标+窗函数 用自己设计的3种滤波器分别对被不同噪声

5、污染的信号进行滤波； 分析得到信号的频谱，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化； 回放语音信号。2 设计的详细原理2.1 卷积线性卷积定理的公式是： 式中的符号“*”代表卷积运算，此式表示线性时不变系统的输出等于输入序列和该系统的单位脉冲响应的卷积。观察此式，计算卷积的基本运算是翻转、移位、相乘和相加，这类卷积称为序列的线性卷积。如果两个序列的长度分别为N和M，那么卷积结果的长度为NM1。循环卷积定理：有限长序列x1(n)和x2(n)， 长度分别为N1和N2， N=max N1, N2 。 x1(n)和x2(n)的N点DFT分别为： X1(k)=DFTx1

6、(n) X2(k)=DFTx2(b) 如果 X(k)=X1(k)X2(k) 则2.2 采样定理对连续信号进行等间隔采样形成采样信号，采样信号的频谱是原连续信号的频谱以采样频率为周期进行周期性的延拓形成的，用公式(3-1)表示。 （3-1）设连续信号xa(t)属带限信号，最高截止频率为c，如果采样角频率s2c，那么让采样信号xa(t)通过一个增益为T，截止频率为s/2的理想低通滤波器，可以唯一地恢复出原连续信号xa(t)。否则s2c会造成采样信号中的频谱混叠现象，不可能无失真地恢复原连续信号。 2.3 椭圆滤波器的设计 椭圆（Elliptic）滤波器在通带和阻带内都具有等波纹幅频响应特性。由于其

7、极点位置与经典场论中的椭圆函数有关，所以由此取名为椭圆滤波器。又因为在1931年考尔（Cauer）首先对这种滤波器进行了理论证明，所以其另一个通用名字为考尔（Cauer）滤波器。椭圆滤波器的典型幅频响应特性曲线如图2.1所示。由图2.1(a)可见，椭圆滤波器通带和阻带波纹幅度固定时，阶数越高, 过渡带越窄； 由图2.1(b)可见，当椭圆滤波器阶数固定时，通带和阻带波纹幅度越小, 过渡带就越宽。所以椭圆滤波器的阶数N由通带边界频率p、阻带边界频率s、通带最大衰减 p和阻带最小衰减 s共同决定。后面对五种滤波器的比较将证实，椭圆滤波器可以获得对理想滤波器幅频响应的最好逼近，是一种性能价格比最高的滤

8、波器，所以应用非常广泛。 MATLAB信号处理工具箱提供椭圆滤波器设计函数ellipap、 ellipord和ellip。其调用格式如下: (1)z, p, k= ellipap(N, Rp, As) (2) N, wpo= ellipord(wp, ws, Rp, As) (3) N, wpo= ellipord(wp, ws, Rp, As, s) (4) B, A= ellip(N, Rp, wpo, ftype) (5) B, A= ellip(N, Rp, wpo, ftype, s) 2.4用汉明窗设计数字滤波器汉明窗能量更加集中在主瓣中，主瓣的能量约占99.96%，瓣峰值幅度为4

9、0 dB，但其主瓣宽度为8/N。此次设计利用MATLAB软件完成。fir1是用窗函数法设计线性相位FIR数字滤波器的工具箱函数，以实现线性相位FIR数字滤波器的标准窗函数法设计。这里的所谓“标准”，是指在设计低通、高通、带通和带阻FIR滤波器时， Hd(ej)分别表示相应的线性相位理想低通、高通、带通和带阻滤波器的频率响应函数。因而将所设计的滤波器的频率响应称为标准频率响应。Fir1的调用格式及功能如下： hn=fir1(M, wc)，返回6 dB截止频率为wc的M阶（单位脉冲响应h(n)长度N=M+1）FIR低通（wc为标量）滤波器系数向量hn，默认选用哈明窗。滤波器单位脉冲响应h(n)与向

10、量hn的关系为h(n)=hn(n+1) n=0, 1, 2, , M而且满足线性相位条件: h(n)=h(N1n)。其中wc为对归一化的数字频率，0wc1。当wc=wcl, wcu时，得到的是带通滤波器，其6 dB通带为wclwcu。hn=fir1(M, wc, ftype)，可设计高通和带阻FIR滤波器。当ftype=high时，设计高通FIR滤波器；当ftype=stop，且wc=wcl, wcu时，设计带阻FIR滤波器。3设计的过程和程序运行结果3.1 卷积3.1.1 线性卷积首先定义两个序列：x1 =2 0 0 7 8 4 2 5 0 1 0 1，x2=1 9 8 9 0 5 2 7，

11、然后进行x2序列的翻转、平移。最终得到x1、x2两序列线性卷积的结果。图3.1.1 （a）线性卷积图3.1.1 （b）线性卷积图3.1(a)、(b)表示了序列x1、x2以及序列x2的翻转、移位过程，从图3.1（b）可看出线性卷积的最后结果。3.1.2 循环卷积首先定义两个序列：序列x=2 0 0 7 8 4 2 5，y=7 5 5 6, 对y序列进行翻转，然后向右平移，再把两序列进行循环卷积，如下图3.2所示。图3.1.2 循环卷积3.2 采样对下面连续信号采样： ，A为幅度因子，a为衰减因子，为模拟角频率。输入采样频率fs（根据程序处理需要指定范围）后，在时域演示信号波形、采样脉冲及采样后信

12、号；在频域演示对应的信号频谱、采样脉冲及频域周期拓延。 分别对1500Hz,300Hz,200Hz的采样信号进行采样，得到不同采样频率下的时域图和频域图，如下图3.2所示。图3.2 不同频率下的信号采样3.3 椭圆模拟滤波器的设计 将采集到的语音信号进行频谱分析，画出采样后语音信号的时域波形和频谱图，然后产生噪声信号并加到语音信号中，得到被污染的语音信号，并回放语音信号。之后进行污染信号的频谱分析，画出被污染的语音信号时域波形和频谱，再利用自己设计滤波器分别对被不同噪声污染的信号进行滤波。图3.3.1 原始信号的波形及频谱图3.3.2 椭圆模拟带通滤波器滤波图3.3.3 椭圆模拟高通滤波器滤波

13、图3.3.1为采集的原始信号的波形及频谱图，加入噪声后，利用椭圆滤波器对其进行滤波，以椭圆模拟高通、带通滤波器为例，分别见图3.3.2、3.3.3。对于椭圆模拟滤波器，这里就不再一一介绍了。3.4 椭圆数字滤波器的设计3.4.1 椭圆数字低通滤波器图3.4.1 椭圆数字低通滤波器滤波3.4.2 椭圆数字高通滤波器图3.4.2 椭圆数字高通滤波器滤波3.4.3 椭圆数字带通滤波器图3.4.3 椭圆数字带通滤波器滤波3.5 利用汉明窗设计FIR滤波器 （3-2）其频谱函数WHm(ej)为 （3-3）其幅度函数WHmg()为 （3-4）3.5.1 汉明窗数字带通滤波器图3.5.1 汉明窗数字带通滤波

14、器滤波3.5.2 汉明窗数字带阻滤波器 图3.5.2 汉明窗数字带阻滤波器滤波3.5.3 汉明窗数字高通滤波器图3.5.3 汉明窗数字高通滤波器滤波4 心得体会在本次课程设计过程中，我发现自己的基础知识和应用能力都非常不够。刚接过任务书，就发现之前学的matlab知识都忘得差不多了，以至于这次在设计的额时候花费了较多时间来学习如何熟练操作matlab。随着课程设计的进行，更加发现自己在数字信号处理课程学的远远不够，不能很好地根据仿真出来的信号频谱对滤波器的参数进行调整以达到最佳的滤波效果。然而，经过参考数字信号处理及MATLAB程序设计与应用书籍，再通过利用参考文献与网络，在老师的指导下，终于

15、完成了本次课程设计。这次课程设计的任务较多，但是通过这次学习，我们树立了正确的设计思想，锻炼了自己理解问题、分析问题、解决问题的能力，加深了对课堂抽象概念的理解，同时巩固了课堂上所学的理论知识。能够更好地理解与掌握了数字信号处理的基本概念、基本原理、基本分析方法。熟悉了matlab的运行环境，初步掌握了matlab语言在数字信号处理中一些基本函数的调用和编写基本程序等应用；熟悉了滤波器设计的一般原理，对滤波器有了更加深刻的认识和理解；学会了滤波器设计的一般步骤；加深了对滤波器设计产生误差的原因的理解和认识；也在学习了数字信号处理课程之后能把所学的理论知识与实际的工程应用联系起来。参考文献1 数

16、字信号处理（第二版），丁玉美等，西安电子科技大学出版社；2 数字信号处理试验指导书王创新、文卉编 长沙理工大学印刷（内部使用）3 数字信号处理及其MATLAB实现，陈怀琛等译，电子工业出版社；4 MATLAB及在电子信息课程中的应用，陈怀琛等，电子工业出版社。5 数字信号处理A.V.奥本海姆，R.W.谢弗著，北京：科学出版社6 数字信号处理理论、算法与实现（第二版）胡广书编著，北京：电子工业出版社附录 MATLAB程序代码%线性卷积x1=0,0,2,0,0,7,8,4,2,5,0,1,0,7;subplot(5,1,1);stem(x1);axis(0 16 0 9);x2=1,9,8,9,0

17、,5,2,7;x2=fliplr(x2);subplot(5,1,2);stem(x2);axis(0 16 0 9);x2(1,16)=0;XL=0,0,0,0,0,0,0,0,0,0,0,0,0,0;sum=0,0,0,0,0,0,0,0,0,0,0,0,0,0;for i=0:13 x3=circshift(x2,-1,i); disp(x3); subplot(5,1,3); stem(x3); axis(0 16 0 16); for j=1:12 XL(1,j+2)=x1(1,j+2)*x3(1,j+2); sum(1,i+1)=sum(1,i+1)+XL(1,j); end di

18、sp(sum); disp(XL); subplot(5,1,4); stem(XL); axis(0 16 0 100); pause(2);end subplot(5,1,5) stem(sum); axis(0 16 0 200)；循环卷积:x=2 0 0 7 8 4 2 5 0 1 0 7; %定义两个序列，其中序列x为自己的学号y=7 5 5 6 ; N=12;length1 = length(x); %序列x的长度length2 = length(y); %序列y的长度if (length1=N&length2=N)x=x zeros(1,N-length1); %长度小于12的序

19、列其后补零y=y zeros(1,N-length2);endv=zeros(N,N); % zeros(m,n)表示生成m行n列的0矩阵 u=zeros(N,1);for k = N:-1:1 A(N-k+1)=y(k); %序列y翻转后用序列A来表示endfor i=1:N R=A(N); for k=N-1:-1:1 A(k+1)=A(k); %翻转后的序列A向右平移 end A(1)=R; for j=1:N v(i,j)=A(j); end v pause(1) figure; subplot(2,1,1) stem(v(i,:);title(序列y的平移动态分解图-刘偲);for

20、i=1:N u(i,1)=x(i); end u subplot(2,1,2) stem(u);title(序列x：-刘偲); w=v*u;end w pause(1) figure; stem(w);title(循环卷积-刘偲);时域采样:Tp=128/1500;%产生M长采样序列x(n)% Fs=1500;T=1/Fs;Fs=1500;T=1/Fs;M=Tp*Fs;n=0:M-1;A=444.128;alph=pi*50*sqrt(2);omega=pi*50*sqrt(2);xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M); %M点

21、FFTxnt)yn=xa(nT);X=fft(xnt);k=0:M-1;fk=k/Tp;figure(1);subplot(2,1,1);stem(xnt);title(波形-刘偲);subplot(2,1,2);plot(fk,abs(X);title(幅频-刘偲);figure(2);subplot(3,2,1);tstem(xnt,yn);box on;title(a) s=1500Hz);k=0:M-1;fk=k/Tp;subplot(3,2,2);plot(fk,abs(Xk);title(a) T*FFTxa(nT),Fs=1500Hz);xlabel(f(Hz);ylabel(幅

22、频);axis(0,Fs,0,1.2*max(abs(Xk)% Fs=300;T=1/Fs;Fs=300;T=1/Fs;M=ceil(Tp*Fs);n=0:M-1;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M); %M点FFTxnt)yn=xa(nT);subplot(3,2,3);tstem(xnt,yn);box on;title(b) Fs=300Hz);k=0:M-1;fk=k/Tp;subplot(3,2,4);plot(fk,abs(Xk);title(b) T*FFTxa(nT),Fs=300Hz);xlabel(f(Hz

23、);ylabel(幅频);axis(0,Fs,0,1.2*max(abs(Xk)% Fs=200;T=1/Fs;Fs=200;T=1/Fs;M=ceil(Tp*Fs);n=0:M-1;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M); %M点FFTxnt)yn=xa(nT);subplot(3,2,5);tstem(xnt,yn);box on;title(c) Fs=200Hz);k=0:M-1;fk=k/Tp;subplot(3,2,6);plot(fk,abs(Xk);title(c) T*FFTxa(nT),Fs=200Hz);x

24、label(f(Hz);ylabel(幅频);axis(0,Fs,0,1.2*max(abs(Xk);原始信号设计：clf;%原始音频抽样Fs=10000;bits=16;T=1;ts=1/Fs;N=38;x,Fs,bits=wavread(LS.wav); %利用Windows下的录音机，录制一段自己的语音信号，时间为1s x=x(:,1); m=length(x);sound(x,Fs,bits);y=fft(x,m); %对录制的信号进行傅立叶变换f=(Fs/m)*1:m;t=1:m/Fs;subplot(211);plot(t,x);title(原始信号波形);xlabel(time(

25、s);subplot(212);plot(f,abs(y);title(原始信号频谱);xlabel(frequency(hz);椭圆模拟带通滤波器：clear;z1,fs,bits=wavread(E:数字信号处理椭圆LS.wav)y1=z1(1:8192,1);Y1=fft(y1);fp1=1200 ;fp2=3000; fc1=1000 ;fc2=3200 ;As=30;Ap=1; Fs=22050;wc=2*fc1/Fs,2* fc2/Fs;wb=2*fp1/Fs,2*fp2/Fs; n,wpo=ellipord(wc,wb,Ap,As,s);b,a=ellip(n,Ap,As,wpo

26、,s);figure(1);freqs(b,a);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1);axis(0,4000,0,100);title(滤波前信号频谱-刘偲);subplot(2,2,2);plot(abs(X);axis(0,4000,0,2.5);title(滤波后信号频谱-刘偲);subplot(2,2,3);plot(z1);title(滤波前信号波形-刘偲);subplot(2,2,4);plot(x);title(滤波后信号波形-刘偲);sound(x,fs,bits);椭圆模拟高

27、通滤波器：clear;close allz1,fs,bits=wavread(LS.wav)y1=z1(1:8192);Y1=fft(y1);fc=2800 ;fp=3000 ;As=100;Ap=1; Fs=22050;wc=2*fc/Fs;wb=2*fp/Fs;n,wpo=ellipord(wc,wb,Ap,As);b,a=ellip(n,Ap,As,wpo,high);figure(1);freqz(b,a);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1);axis(0,4000,0,100);ti

28、tle(滤波前信号频谱-刘偲);subplot(2,2,2);plot(abs(X);axis(0,4000,0,30);title(滤波后信号频谱-刘偲);subplot(2,2,3);plot(z1);title(滤波前信号波形-刘偲);subplot(2,2,4);plot(x);title(滤波后信号波形-刘偲);sound(x,fs,bits);椭圆数字低通滤波器：clf;close all;Fs=8000;bits=8;T=1;ts=1/Fs;N=38;x,Fs,bits=wavread(LS.wav); %利用Windows下的录音机，录制一段自己的语音信号，时间为1s x=x(

29、:,1); m=length(x);sound(x,Fs,bits);% 语音信号的频谱分析，画出采样后语音信号的时域波形和频谱图；y=fft(x,m); %对录制的信号进行傅立叶变换f=(Fs/m)*1:m;t=1:m/Fs;figure(1);subplot(221);plot(t,x);title(a)原始信号波形图-刘偲);xlabel(time(s);axis(0,1.5,-0.5,0.5);subplot(222);plot(f,abs(y);title(b)原始信号频谱图-刘偲);xlabel(frequency(hz);axis(0,4000,0,300);% 产生噪声信号并加到语音信号中，得到被污染的语音信号，并回放语音信号；% 噪声信号noise1=0.01*sin(2*pi*3000*t); %高频噪声noise2=0.01*sin(2*pi*200*t); %低频噪声noise3=0.01*sin(2*pi*4000*t)+0.01*sin(2*pi*5000*t); %带阻噪声noise4=0.01*sin(2*pi*5000*t)+0.01*sin(2*pi*200*t); %带通噪声% 椭圆数字低通滤波 fp