数字滤波器设计及应用综合实验

1、一、实验目的数字信号处理实验实验四、数字滤波器设计及应用综合实验淫院：信息工程序院班级：电了 101班姓名：学号：1. 熟悉hr数字滤波器的设计原理及方法。2. 熟悉ftr数字滤波器的设计原理及方法。3. 掌握利用matlab实现数字滤波器的方法4. 掌握利用数字滤波器进行信号处理的方法。5. 了解基于simulink的动态仿真实现信号滤波的基本方法。二、实验内容及要求实验内容：综合运川数字滤波器设计的相关知识，根据给定设计方法要求，川脉冲响应不变法和双 线性变换法设计iir数字滤波器；利用窗函数设计法设计fir数字滤波器。根据实际信号的频 谱特性，分析、确定滤波器设计技术指标，实现对信号的滤

2、波。1. iir数字滤波器设计(1) 用脉冲响应不变法设计巴特沃斯数字滤波器。该实验所需m文件如下： 、butterworth低通滤波器原型设计函数： function b, a=afd_butt (wp, ws, rp, as)n 二 ce 订(logl0(107rp/10)-l)/(107as/10)-l)/(2*logl0(wp/ws)fpri ntf n butterworth filter order二2. of n', n)0megac=wp/(10" (rp/10)t)" (1/(2*n)b, a=u_buttap (n, omegac) 、非归一化b

3、utterworth模拟低通滤波器设计函数：function b, a=u_buttap(n, omegac);z, p, k=buttap (n);p二p*0megac;k二k*0megac"n;b=rcal(poly(z);b二k*b;a=real(poly(p); 、利用脉冲响应不变法从模拟到数字滤波器变换函数：function b, a=imp_invr (c, d, t)r, p, k=residue(c, d);p二exp(p*t);b, a=residuez (r, p, k);b=real (bj);a=real (a，); 、频率响应函数freqz的修正：funct

4、ion db, mag, pha, w二freqz m(b, a);h, w二freqz (b, a, 1000, ' whole');i匸(1:501)'w=(w(l:501)'mag=abs(h);db二20*log10(mag+cps)/max(mag);pha二angle(h);本实验程序如下：wp=0. 2*pi;ws二0. 3*pi;rp=l;as=15;t=l;omegap=wp/t;omegas=ws/t;cs, ds=afd butt (omegap, omegas, rp, as);b, a=imp_invr (cs, ds, t)db, m

5、ag, pha, w=freqz_m(b, a);subplot (2, 1, 1) ; plot (w/pi, mag);title(? digital filter magnitude response，)axis(0, 1, 0, 1. 1)subplot (2, 1, 2) ;plot (w/pi, db);title(j digital filter magnitude in db')axis(0, 1, -40, 5);结果：n1 =5. 8858n =6butterworth filter 0rder= 6omegac =0. 7032b -0. 0000 0. 0006

6、 0.0101 0.0161 0. 0041 0. 0001a =1. 0000 -3. 3635 5. 0684 -4. 2759 2. 1066-0. 5706 0. 0661本实验波形图如下：（2）用双线性变换法设计切比雪夫数字滤波器。本实验所需m文件如下：非归一化切比雪夫t型模拟低通滤波器原型设计: function b,a=u_chblap(n, rp, omegac);z, p, k=cheblap (n, rp); areal(poly(p);ann=a(n+l);p=p*0megac; a二real(poly(p); anu二a(n+l);k二k*dnu/ann;boreal(

7、poly (z); b二k*b;本实验程序如下：wp二0. 2*pi;ws=0. 3*pi;rp二1;as二15;t=1;omegap二(2/t)*tan(wp/2);omegas=(2/t)*tan(ws/2);ep=sqrt (10' (rp/10)t);ripple二sqrt(l/(l+ep*ep);attn=l/(1(t (as/20);al = l/attn;al=sqrt (al*alt)/op;a2=0megas/0megap;n=ceil (logm(al+sqrt (al*al-l) )/logm(a2+sqrt (a2*a2t); fprintf c n cheby

8、shev filter 0rder=%2. ofn，, n)in db');0. 00730. 0018-2. 29250. 5507response') ; axis(0, 1, 0, 1. 1)cs, ds=u_chblap (n, rp, omegap); b, a=bilinear (cs, ds, 1/t) db, mag, pha, w=frcqz_m(b, a); subplot (2, 1, 1) ;plot (w/pi, mag); title (j digital f i 1 ter magnitude subplot (2, 1, 2); plot (w/

9、pi, db); titlec1 digital filter magnitude axis(0, 1,-40, 5);chebyshev filter 0rder= 4b =0.00180.00730.0110a =1. 0000-3. 05433. 8290本实验波形如下：digital filter magnitude response(3) 用双线性变换法设计巴特沃斯数字滤波器，并将直接型结构转换成级联型结构。 将直接形式变为级联形式函数：function bo, b, a=dir2cas (b, a)bo=b(1);b=b/bo;ao=a(1);a=a/ao;bo=bo/ao;m=l

10、ength(b); n=length(a);if n>mb=b zeros (1, nm);elseif m>na= a zeros (1, m-n);elsenm 二 0;endk=floor (n/2); b二zeros (k, 3) ; a=zeros (k, 3);if k*2=nb二b 0; a=a 0;endbroots=cplxpair(roots (b);aroots=cplxpair(roots (a);for i=l:2:2*kbrow=broots(i:1:i+1,:);brow-real(poly(brow);b(fix(i+l)/2, :) =brow;a

11、row=aroots(i:1:i+1,:);arow=real(poly(arow);a(fix(i+l)/2, :)=arow;end本实验程序如下：wp=o. 2*pi;ws=o. 3*pi;rp=l;as=15;t=1;omegap二(2/t)*tan(wp/2);omegas=(2/t)*tan(ws/2);cs, ds=afd butt (omegap, omegas, rp, as);b, a=bilinear (cs, ds, t)db, mag, pha, w=freqz_m(b, a);subplot (2, 1, 1) ; plot (w/pi, mag);title(?

12、digital filter magnitude response，); axis(0, 1, 0, 1. 1)subplot (2, 1, 2) ;plot (w/pi, db);title(，digital filter magnitude in db'); axis(0, 1,-40, 5);bo, b, a =dir2cas (b, a)结果：n1 =5. 3044n =6butterworth filter 0rder= 6omegac =0. 7273b 二0. 00060. 00350. 00870.01160. 00870. 0035d 一1.0000-3.31434.

13、 9501-4. 14332. 0275-0. 5458b =1.00002. 03351.03381.00001.99961. 00001.00001.96690. 9673a =1.0000-0. 94590. 23421.0000-1.05410. 37531.0000-1. 31430. 7149本实验波形如下：0. 00060. 0628digital filter magnitude response2. fir数字滤波器设计 本实验所需m文件如下： function hd=ideal_lp(wc, m); alpha=(m-l)/2;n=0:(m-l);m二n-alpha+eps

14、; hd=sin(wc*m). /(pi*m); 本实验程序如下：f1=100;f2=200;fs=2000;m=(0 3*fl)/(fs/2);m二round(8/m);n=m-1;b=firl(n, 0. 5*f2/(fs/2);figure(l)h, f=freqz (b, 1, 512);%h, w=freqz(b, a, n)plot(f*fs/(2*pi), 20*1oglo(abs(h) xlabel ('频率/赫兹');ylabelc增益/分贝')； titlec滤波器的增益响应')； figure(2) subplot (211) t=0:l/

15、fs:0.5;s=sin(2*pi*fl*t)+sin(2*pi *f2*t);plot (t, s);xlabel ('时间/秒');ylabel ('幅度);titlec信号滤波前时域图')；subplot(212)fs=fft(s, 512);afs=abs(fs);f=(0:255)*fs/512;plot(f, afs (1:256);xlabel ('频率/赫兹');ylabel ('幅度')； titlec信号滤波前频域图')；figure(3)sf=filter (b, 1, s) ; subplot (2

16、11)plot (t, sf) xlabel c 时间/秒);ylabel c 幅度);titlec信号滤波后时域图')； axis(0. 2 0. 5 -2 2); subplot(212) fsf=fft(sf, 512); afsf=abs(fsf);f= (0:255) * fs/512; plot(f, afsfd： 256) xlabel c频率/赫兹')；ylabel c 幅度')；titlec信号滤波后频域图')；本实验波形如图：滤波器的增益响应200-20406081:-100-120010020030040050060070080090010

17、00频率脇兹信号滤波前时域图信号滤波前频域图频率脇兹信号滤波后时域图时间砂 信号滤波后频域图频率侨赤兹3. 利用数字滤波器实现dtmf信号的提取双音多频（dual tone multi frequency）信号是音频电话中的拨号信号，每一个数字（0-9） iii两个不同频率单音组成（毎个单音用正弦信号表示），所用频率分为高频带和低频 带两组,每个数字由高、低频带中各一个频率组成,例如数字9使用852i1z和1477hz两个频率。数字与符号对应频率关系见表3-1所示。表3-1双频拨号的频率分配1209hz1336hz1477hz1633hz697hz123a770hz456b852hz789c9

18、41hz*0d实验内容：电话中dtmf信号的产生于检测方法：在电话中，数字09的屮每一个都用两个不同的单音频传输，所用的8个频率分成高频 带和低频带两组，低频带有-四个频率：679hz, 770hz, 852hz和941h/高频带也冇阻个频率: 1209hz, 1336hz, 1477hz和1633hz.。每一个数字均由高、低频带中各一个频率构成,例如1 用697hz和1209hz两个频率,信号用sin(2对')+ sin(2对/)表示,其屮.齐=679hz， f2 = 1209hz。这样8个频率形成16种不同的双频信号。具体号码以及符号对应的频率如 表3-2所示。表中最后一列在电话屮

19、暂时未用。x1209hz1336hz1477hz633hz697hz123a770hz456b852hz789c942hz*0#d表3-2双频拨号的频率分配dtmf信号在电话屮有两种作用，一个是用拨号信号去控制交换机接通被叫的用户电话 机，另一个作用是控制电话机的各种动作，如播放留言、语音信箱等。自己选择一个数字(0-9),通过数字方法产生该数字的双频信号；设采样频率为8000hzo 选择数字1,产生双频信号：n=205;k=l8, 20, 22, 24, 31, 34, 38, 42;fl=697, 770, 852,941;f2=1209, 1336, 1477, 1633;n=0:102

20、3;x = sin (2*pi*n*f1(1)/8000) + sin(2*pi*n*f2(1)/8000);x=goertzel(x(1:n),k+1);val = abs(x);subplotstem(k, val,'');grid;xlabel (' k') ;y label ('| x(k) |');axis(10 50 0 120)本实验波形如图：10080f 6040i t t.200101520253035404550kh己设计实验：dtmf双频拨号信号的生成和检测程序tm=l, 2, 3, 65;4, 5, 6, 66;7, 8, 9, 67;42, 0, 35, 68;n=205;k=l8, 20, 22, 24, 31, 34, 38, 42;fl=697, 770, 852,941;f2 二1209, 1336, 1477, 1633;t