数字出版视野中的版权集体管理制度创新

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

2、数字滤波器设计（1）用脉冲响应不变法设计巴特沃斯数字滤波器。该实验所需M文件如下：、butterworth低通滤波器原型设计函数：function b,a=afd_butt(Wp,Ws,Rp,As)N=ceil(log10(10(Rp/10)-1)/(10(As/10)-1)/(2*log10(Wp/Ws)fprintf(n Butterworth Filter Order=%2.0fn,N)OmegaC=Wp/(10(Rp/10)-1)(1/(2*N) b,a=u_buttap(N,OmegaC)、非归一化Butterworth模拟低通滤波器设计函数：function b,a=u_butta

3、p(N,Omegac);z,p,k=buttap(N);p=p*Omegac; k=k*OmegacN;B=real(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(b); a=real(a);、频率响应函数freqz的修正：function db,mag,pha,w=freqz_m(b,a);H,w=freqz(b,a,1000,whole); H=(H(1:

4、501); w=(w(1:501); mag=abs(H); db=20*log10(mag+eps)/max(mag); pha=angle(H);本实验程序如下：wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;T=1;OmegaP=wp/T;OmegaS=ws/T;cs,ds=afd_butt(OmegaP,OmegaS,Rp,As);b,a=imp_invr(cs,ds,T)db,mag,pha,w=freqz_m(b,a);subplot(2,1,1);plot(w/pi,mag);title(digital filter Magnitude Response)axis(

5、0,1,0,1.1)subplot(2,1,2);plot(w/pi,db);title(digital filter Magnitude in DB)axis(0,1,-40,5);结果：N1 =5.8858N =6Butterworth Filter Order= 6OmegaC =0.7032b =0.0000 0.0006 0.0101 0.0161 0.0041 0.0001a =1.0000 -3.3635 5.0684 -4.2759 2.1066 -0.5706 0.0661本实验波形图如下：（2）用双线性变换法设计切比雪夫数字滤波器。本实验所需M文件如下：非归一化切比雪夫I型

6、模拟低通滤波器原型设计：function b,a=u_chb1ap(N,Rp,Omegac);z,p,k=cheb1ap(N,Rp); a=real(poly(p); aNn=a(N+1); p=p*Omegac;a=real(poly(p);aNu=a(N+1);k=k*aNu/aNn;B=real(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)-1);Ripple=sqrt(1/(1+ep*ep);

7、Attn=1/(10(As/20);A1=1/Attn;a1=sqrt(A1*A1-1)/ep;a2=OmegaS/OmegaP;N=ceil(logm(a1+sqrt(a1*a1-1)/logm(a2+sqrt(a2*a2-1);fprintf(n Chebyshev Filter Order=%2.0fn,N)cs,ds=u_chb1ap(N,Rp,OmegaP);b,a=bilinear(cs,ds,1/T)db,mag,pha,w=freqz_m(b,a);subplot(2,1,1);plot(w/pi,mag);title(digital filter Magnitude Resp

8、onse); 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);Chebyshev Filter Order= 4b = 0.0018 0.0073 0.0110 0.0073 0.0018a =1.0000 -3.0543 3.8290 -2.2925 0.5507本实验波形如下：（3）用双线性变换法设计巴特沃斯数字滤波器，并将直接型结构转换成级联型结构。将直接形式变为级联形式函数：function b0,B,A=dir2cas(b,a)b0=b

9、(1);b=b/b0;a0=a(1);a=a/a0;b0=b0/a0; M=length(b); N=length(a);if NM b=b zeros(1,N-M);elseif MN a=a zeros(1,M-N);else NM=0;endK=floor(N/2); B=zeros(K,3); A=zeros(K,3);if K*2=N b=b 0; a=a 0;end broots=cplxpair(roots(b); aroots=cplxpair(roots(a); for i=1:2:2*K Brow=broots(i:1:i+1,:); Brow=real(poly(Brow

10、); B(fix(i+1)/2,:)=Brow; Arow=aroots(i:1:i+1,:); Arow=real(poly(Arow); A(fix(i+1)/2,:)=Arow;end本实验程序如下：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);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);tit

11、le(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);b0,B,A=dir2cas(b,a)结果：N1 = 5.3044N = 6 Butterworth Filter Order= 6OmegaC = 0.7273b = 0.0006 0.0035 0.0087 0.0116 0.0087 0.0035 0.0006a =1.0000 -3.3143 4.9501 -4.1

12、433 2.0275 -0.5458 0.0628B = 1.0000 2.0335 1.0338 1.0000 1.9996 1.0000 1.0000 1.9669 0.9673A = 1.0000 -0.9459 0.2342 1.0000 -1.0541 0.37531.0000 -1.3143 0.7149本实验波形如下：2. FIR数字滤波器设计本实验所需M文件如下：function hd=ideal_lp(wc,M);alpha=(M-1)/2; n=0:(M-1); m=n-alpha+eps; hd=sin(wc*m)./(pi*m);本实验程序如下：f1=100;f2=20

13、0; fs=2000; m=(0.3*f1)/(fs/2); M=round(8/m); N=M-1; b=fir1(N,0.5*f2/(fs/2); figure(1)h,f=freqz(b,1,512); %H,W=freqz(B,A,N) plot(f*fs/(2*pi),20*log10(abs(h) xlabel(频率/赫兹);ylabel(增益/分贝);title(滤波器的增益响应);figure(2)subplot(211)t=0:1/fs:0.5;s=sin(2*pi*f1*t)+sin(2*pi*f2*t);plot(t,s);xlabel(时间/秒);ylabel(幅度);

14、title(信号滤波前时域图);subplot(212)Fs=fft(s,512);AFs=abs(Fs);f=(0:255)*fs/512;plot(f,AFs(1:256);xlabel(频率/赫兹);ylabel(幅度);title(信号滤波前频域图);figure(3)sf=filter(b,1,s); subplot(211)plot(t,sf) xlabel(时间/秒);ylabel(幅度);title(信号滤波后时域图);axis(0.2 0.5 -2 2); subplot(212)Fsf=fft(sf,512); AFsf=abs(Fsf);f=(0:255)*fs/512;

15、 plot(f,AFsf(1:256)xlabel(频率/赫兹);ylabel(幅度);title(信号滤波后频域图);本实验波形如图：3.利用数字滤波器实现DTMF信号的提取双音多频（Dual Tone Multi Frequency）信号是音频电话中的拨号信号，每一个数字（0-9）由两个不同频率单音组成（每个单音用正弦信号表示），所用频率分为高频带和低频带两组，每个数字由高、低频带中各一个频率组成，例如数字9使用852Hz和1477Hz两个频率。数字与符号对应频率关系见表3-1所示。表3-1 双频拨号的频率分配1209Hz1336Hz1477Hz1633Hz697Hz123A770Hz45

16、6B852Hz789C941Hz*0#D实验内容：电话中DTMF信号的产生于检测方法：在电话中，数字09的中每一个都用两个不同的单音频传输，所用的8个频率分成高频带和低频带两组，低频带有四个频率：679Hz,770Hz,852Hz和941Hz；高频带也有四个频率：1209Hz,1336Hz,1477Hz和1633Hz.。每一个数字均由高、低频带中各一个频率构成，例如1用697Hz和1209Hz两个频率，信号用表示，其中，。这样8个频率形成16种不同的双频信号。具体号码以及符号对应的频率如表3-2所示。表中最后一列在电话中暂时未用。列行1209Hz1336Hz1477Hz633Hz697Hz12

17、3A770Hz456B852Hz789C942Hz*0#D表3-2 双频拨号的频率分配DTMF信号在电话中有两种作用，一个是用拨号信号去控制交换机接通被叫的用户电话机，另一个作用是控制电话机的各种动作，如播放留言、语音信箱等。自己选择一个数字（0-9），通过数字方法产生该数字的双频信号；设采样频率为8000Hz。选择数字1，产生双频信号：N=205;K=18,20,22,24,31,34,38,42;f1=697,770,852,941;f2=1209,1336,1477,1633;n=0:1023; x = sin(2*pi*n*f1(1)/8000) + sin(2*pi*n*f2(1)/

18、8000);X=goertzel(x(1:N),K+1); val = abs(X);subplotstem(K,val,.);grid;xlabel(k);ylabel(|X(k)|); axis(10 50 0 120)本实验波形如图：自己设计实验：DTMF双频拨号信号的生成和检测程序tm=1,2,3,65;4,5,6,66;7,8,9,67;42,0,35,68; N=205;K=18,20,22,24,31,34,38,42;f1=697,770,852,941; f2=1209,1336,1477,1633; TN=input(键入6位电话号码= ); TNr=0; for m=1:6; d=fix(TN/10(6-m); TN=TN-d*10(6-m); for p=1:4; for q=1:4; if tm(p,q)=abs(d); break,end endif tm(p,q)=abs(d); break,end end n=0:1023; x = sin(2*pi*n*f1(p)/8000) + sin(2*pi*n*f2(q)/8000); sound(x,8000); pause(0.1