实验五数字信号处理综合设计一、实验目的1．掌握在Windows环境下语音信号采集的方法；2．掌握MATLAB设计FIR和IIR数字滤波器的方法；3．学会用MATLAB 对信号进行分析和处理。

二、实验内容1．语音信号的采集要求利用windows下的录音机或其他软件，录制一段自己的话音，时间定为10秒。

然后在MATLAB软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。

通过wavread函数的使用，要求理解采样频率、采样位数等概念。

wavread函数调用格式：y=wavread(file)，读取file所规定的wav文件，返回采样值放在向量y中。

[y,fs,nbits]=wavread(file)，采样值放在向量y 中，fs 表示采样频率（Hz），nbits 表示采样位数。

y=wavread(file,N)，读取前N 点的采样值放在向量y中。

y=wavread(file,[N1,N2])，读取从N1点到N2点的采样值放在向量y中。

例如：x1=wavread(h:\课程设计2\shuzi.wav); %读取语音信号的数据，赋给变量x12．语音信号的频谱分析要求首先画出语音信号的时域波形,然后对语音信号进行频谱分析，在MATLAB中，可以利用函数fft对信号进行快速付立叶变换，得到信号的频谱特性；从而加深对频谱特性的理解。

解析：clear;clc;clf;%语音信号的频谱分析y=wavread('2.wav');[y,Fs,nbits]=wavread('2.wav');N=2048;Y=fft(y,N);Y1=fftshift(Y);plot(abs(Y));title('语音信号的幅度谱');f=0:1/Fs;(size(y)-1)/Fs;%将所加噪声信号的点数调整到与原信号相同050100150200250300350400450500246805010015020025030035040045050002468103．设计数字滤波器和画出频率响应根据语音信号的特点给出有关滤波器的性能指标：1）低通滤波器性能指标，fp=1000Hz ，fc=1200 Hz ，As=50dB ，Ap=1dB ； 2）高通滤波器性能指标，fc=4800 Hz ，fp=5000 Hz As=50dB ，Ap=1dB ； 3）带通滤波器性能指标，fp1=1200 Hz ，fp2=3000 Hz ，fc1=1000 Hz ，fc2=3200 Hz ， As=50dB ，Ap=1dB 。

首先用窗函数法设计上面要求的三种滤波器（FIR 滤波器）；然后在用双线性变换法设计上面要求的三种滤波器，可以设计巴特沃斯、切比雪夫Ⅰ、切比雪夫Ⅱ、椭圆型IIR 滤波器；最后，利用 MATLAB 中的函数 freqz 画出各滤波器的频率响应。

方案一：首先用窗函数法设计上面要求的三种滤波器（FIR 滤波器） 由于，所以选择哈明窗设计三种滤波器 A 、FIR 低通滤波器 代码：%FIR 低通滤波器 clear;clc;clf;fp=1000;fs=1200;Wp = 2*pi*fp; Ws = 2*pi*fs;[y,Fs,nbits]=wavread('2.wav'); wp=Wp/Fs;ws=Ws/Fs; wc = (ws+wp)/2; Rp = 1;As = 50; tr_width=ws-wp;M= ceil(6.6*pi/tr_width) + 1; n=[0:1:M-1];hd=ideal_lp(wc,M);w_ham= (hamming(M))'; h = hd .* w_ham ;[db,mag,pha,grd,w]= reqz_m(h,[1]); figure(1);subplot(311); plot(w/pi,mag); title('低通加哈明窗幅度谱'); grid on ;subplot(312); plot(w/pi,pha);title('低通加哈明窗窗相位谱');grid on ; subplot(313); plot(w/pi,db);title('低通加哈明窗对数幅度响应');grid on ;0.10.20.30.40.50.60.70.80.9100.511.5低通加哈明窗幅度谱00.10.20.30.40.50.60.70.80.91-505低通加哈明窗窗相位谱00.10.20.30.40.50.60.70.80.91-2000200低通加哈明窗对数幅度响应B 、FIR 高通滤波器 %FIR 高通滤波器 clear;clc;clf;fp=1000;fs=1200;Wp = 2*pi*fp; Ws = 2*pi*fs; [y,Fs,nbits]= wavread('2.wav'); sound(y, Fs);wp=Wp/Fs;ws=Ws/Fs; wc = (ws+wp)/2; Rp = 1;As = 50; tr_width=ws - wp;M = ceil(6.6*pi/tr_width) + 1; n=[0:1:M-1];hd = ideal_lp(wc,M);w_ham = (hamming(M))'; h = hd .* w_ham ; fph=5000;fsh=4800;Wph = 2*pi*fph; Wsh = 2*pi*fsh; wph=Wph/Fs;wsh=Wsh/Fs;alpha = -(cos((wp+wph)/2))/(cos((wp-wph)/2)) Nz = -[alpha,1]; Dz = [1,alpha]; [bhp,ahp] = zmapping(h,1,Nz,Dz);[db,mag,pha,grd,w]= freqz_m(bhp,ahp); figure(1)subplot(311); plot(w/pi,mag);title('高通加哈明窗幅度谱');grid on ; subplot(312); plot(w/pi,pha);title('高通加哈明窗窗相位谱');grid on ; subplot(313); plot(w/pi,db);title('高通对数幅度响应');grid on ;00.10.20.30.40.50.60.70.80.905x 10210高通加哈明窗幅度谱00.10.20.30.40.50.60.70.80.9-505高通加哈明窗窗相位谱00.10.20.30.40.50.60.70.80.9-400-2000高通对数幅度响应C 、FIR 带通滤波器%FIR 带通滤波器clear;clc;clf;fp=1000;fs=1200; Wp = 2*pi*fp; Ws = 2*pi*fs; [y,Fs,nbits]=wavread('2.wav'); sound(y, Fs);wpl=Wp/Fs;wsl=Ws/Fs;wc = (wsl+wpl)/2; Rp = 1;As = 50; tr_width = wsl - wpl; M = ceil(6.6*pi/tr_width) + 1;n=[0:1:M-1];hd = ideal_lp(wc,M); w_ham = (hamming(M))'; h = hd .* w_ham ; fp1=1200;fp2=3000; fc1=1000;fc2=3200; Wp1 = 2*pi*fp1;Wp2 =2*pi*fp2;Ws1 = 2*pi*fc1;Ws2 = 2*pi*fc2;wp1=Wp1/Fs;ws1=Ws1/Fs;wp2=Wp2/F s;ws2=Ws2/Fs; K=cot(wp2-wp1)*ta n(wpl/2);beta=(cos((wpl+wp2)/2))/(cos((wp2-wp1)/2)); alpha1 = -2*beta*K/(K+1); alpha2=(K-1)/(K+1); Nz = -[alpha2,-alpha1,1]; Dz = [1,-alpha1,alpha2]; [bdp,adp] = zmapping(h,1,Nz,Dz );[db,mag,pha,grd,w]= freqz_m(bdp,adp); figure(1) subplot(311); plot(w/pi,mag);title('带通滤波器加哈明窗幅度谱');grid on ;subplot(312);plot(w/pi,pha);title('带通滤波器加哈明窗窗相位谱');grid on ; subplot(313);plot(w/pi,db);title('带通滤波器加哈明窗对数幅度响应');grid on ;0.10.20.30.40.50.60.70.80.9105x 10150带通滤波器加哈明窗幅度谱0.10.20.30.40.50.60.70.80.91-505带通滤波器加哈明窗窗相位谱00.10.20.30.40.50.60.70.80.91-400-2000带通滤波器加哈明窗对数幅度响应方案二：用双线性变换法设计三种滤波器（IIR 滤波器），在MATLAB 中，采用巴特沃斯、切比雪夫Ⅰ、切比雪夫Ⅱ、椭圆型滤波器设计IIR 滤波器A 、IIR 低通滤波器：%IIR 低通滤波器clear;clc;close; fp=1000;fc=1200;W p = 2*pi*fp; Ws = 2*pi*fc; [y,Fs,nbits]=wavread('2.wav'); sound(y, Fs);T=1/Fs; wpl=Wp/Fs;wsl=Ws/Fs;Rp = 1;As = 50; Omega = (2/T)*tan(wpl/2);Om egaS = (2/T)*tan(wsl/2); [cs1,ds1]=afd_butt(OmegaP,O megaS,Rp,As); [cs2,ds2]=afd_chb1(OmegaP,O megaS,Rp,As); [cs3,ds3]=afd_chb1(OmegaP,O megaS,Rp,As); [cs4,ds4]=afd_elip(OmegaP,O megaS,Rp,As); [blp1,alp1] = bilinear(cs1,ds1,Fs);[blp2,alp2] = bilinear(cs2,ds2,Fs); [blp3,alp3] = bilinear(cs3,ds3,Fs); [blp4,alp4] = bilinear(cs4,ds4,Fs); [db1,mag1,pha1,grd1,w1]=freqz_m(blp1,alp1); [db2,mag2,pha2,grd2,w2]=freqz_m(blp2,alp2); [db3,mag3,pha3,grd3,w3]=freqz_m(blp3,alp3); [db4,mag4,pha4,grd4,w4]=freqz_m(blp4,alp4); figure(1);subplot(211);plot(w1/pi,mag1);title('低通巴特沃斯幅度响应');grid;subplot(212);plot(w1/pi,db1);title('低通巴特沃斯对数幅度响应');grid; figure(2)subplot(211);plot(w2/pi,mag2);title('低通切比雪夫Ⅰ型幅度响应');grid;subplot(212);plot(w2/pi,db2);title('低通切比雪夫Ⅰ型对数幅度响应');grid; figure(3)subplot(211);plot(w3/pi,mag3);title('低通切比雪夫Ⅱ型幅度响应');grid;subplot(212);plot(w3/pi,db3);title('低通切比雪夫Ⅱ型对数幅度响应');grid; figure(4)subplot(211);plot(w4/pi,mag4);title('低通椭圆型幅度响应');grid;subplot(212);plot(w4/pi,db4);title('低通椭圆型对数幅度响应');grid;0.5102468低通巴特沃斯幅度响应00.51510低通切比雪夫Ⅰ型幅度响应00.510.511.5低通切比雪夫Ⅱ型幅度响应00.510.51低通椭圆型幅度响应B 、IIR 高通滤波器： %IIR 高通滤波器 clear;clc;clf;fp1=1000;fc1=1200; Wp1 = 2*pi*fp1; Ws1 = 2*pi*fc1;Rp = 1;As = 50; [y,Fs,nbits]=wavread('2.wav'); sound(y, Fs);T=1/Fs; wpl=Wp1/Fs;wsl=W s1/Fs;fp2=5000;fc2=4800; Wp2= 2*pi*fp2; Ws2 = 2*pi*fc2; wph=Wp2/Fs;wsh=Ws2/Fs; OmegaP = (2/T)*tan(wpl/2); OmegaS = (2/T)*tan(wsl/2); [cs1,ds1]=afd_butt(OmegaP,O megaS,Rp,As) [cs2,ds2]=afd_chb1(OmegaP,O megaS,Rp,As); [cs3,ds3]=afd_chb2(OmegaP,O megaS,Rp,As); [cs4,ds4]=afd_elip(OmegaP,OmegaS,Rp,As); [blp1,alp1] = bilinear(cs1,ds1,Fs); [blp2,alp2] = bilinear(cs2,ds2,Fs); [blp3,alp3] = bilinear(cs3,ds3,Fs); [blp4,alp4] = bilinear(cs4,ds4,Fs); alpha=-(cos((wpl+wph)/2))/(cos((wpl-wph)/2)); Nz = -[alpha,1]; Dz = [1,alpha]; [bhp1,ahp1]=zmapping(blp1,alp1,Nz,Dz);[bhp2,ahp2]=zmapping(blp2,alp2,Nz,Dz);[bhp3,ahp3]=zmapping(blp3,alp3,Nz,Dz);[bhp4,ahp4]=zmapping(blp4,alp4,Nz,Dz); [cs11,ds11]=afd_butt(OmegaP,O megaS,Rp,As)[db1,mag1,pha1,grd1,w1]=freqz_m(bhp1,ahp1); [db2,mag2,pha2,grd2,w2]=freqz_m(bhp2,ahp2); [db3,mag3,pha3,grd3,w3]=freqz_m(bhp3,ahp3); [db4,mag4,pha4,grd4,w4] = freqz_m(bhp4,ahp4); figure(1)subplot(211);plot(w1/pi,mag1);title('高通巴特沃斯幅度响应');grid;subplot(212);plot(w1/pi,db1);title('高通巴特沃斯对数幅度响应');grid; figure(2)subplot(211);plot(w2/pi,mag2);title('高通切比雪夫Ⅰ型幅度响应');grid;subplot(212);plot(w2/pi,db2);title('高通切比雪夫Ⅰ型对数幅度响应');grid; figure(3)subplot(211);plot(w3/pi,mag3);title('高通切比雪夫Ⅱ型幅度响应');grid;subplot(212);plot(w3/pi,db3);title('高通切比雪夫Ⅱ型对数幅度响应');grid; figure(4)subplot(211);plot(w4/pi,mag4);title('高通椭圆型幅度响应');grid;subplot(212);plot(w4/pi,db4);title('高通椭圆型对数幅度响应');grid;0.5100.511.52高通巴特沃斯幅度响应00.511234高通切比雪夫Ⅰ型幅度响应0.5100.511.52高通切比雪夫Ⅱ型幅度响应00.5100.51高通椭圆型幅度响应B 、IIR 带通滤波器： %IIR 带通滤波器clear;clc;clf;fp=1000;fc =1200;Wp = 2*pi*fp; Ws = 2*pi*fc;[y,Fs,nbits]=wavread('2.wav ’');sound(y, Fs);T=1/Fs; wpl=Wp/Fs;wsl=Ws/Fs; fp1=1200;fp2=3000;fc1=1000;fc2=3200;Wp1 = 2*pi*fp1;Wp2 = 2*pi*fp2;Ws1 = 2*pi*fc1;Ws2= 2*pi*fc2;wp1=Wp1/Fs;ws1=Ws1/Fs;wp2=Wp2/Fs;ws2=Ws2/Fs;Rp = 1;As = 50; OmegaP = (2/T)*tan(wpl/2);OmegaS = (2/T)*tan(wsl/2); [cs1,ds1]=afd_butt(OmegaP,Omeg aS,Rp,As) [cs2,ds2]=afd_chb1(OmegaP,Ome gaS,Rp,As); [cs3,ds3]=afd_chb1(OmegaP,Ome gaS,Rp,As); [cs4,ds4]=afd_elip(OmegaP,Omeg aS,Rp,As); [blp1,alp1] = bilinear(cs1,ds1,Fs); [blp2,alp2] = bilinear(cs2,ds2,Fs); [blp3,alp3] = bilinear(cs3,ds3,Fs); [blp4,alp4] = bilinear(cs4,ds4,Fs); K=cot(wp2-wp1)*tan(wpl/2);beta=(cos((wpl+w p2)/2))/(cos((wp2-wp1)/2));alpha1=-2*beta*K/(K+1);alpha2 =(K-1)/(K+1); Nz=-[alpha2,-alpha1,1]; Dz = [1,-alpha1,alpha2]; [bhp1,ahp1]=zmapping(blp1,alp1,Nz,Dz);[bhp2,ahp2]=zmapping(blp2,alp2,Nz,Dz);[bhp3,ahp3]=zmapping(blp3,alp3,Nz,Dz);[bhp4,ahp4]=zmapping(blp4,alp4,Nz,Dz);[cs11,ds11]=afd_butt(OmegaP,OmegaS,Rp,As)[db1,mag1,pha1,grd1,w 1]=freqz_m(bhp1,ahp1); [db2,mag2,pha2,grd2,w 2]=freqz_m(bhp2,ahp2); [db3,mag3,pha3,grd3,w 3]=freqz_m(bhp3,ahp3); [db4,mag4,pha4,grd4,w 4]=freqz_m(bhp4,ahp4); figure(1)subplot(221);plot(w1/pi,mag1);title('带通巴特沃斯型幅度不同响应');grid;subplot(222);plot(w2/pi,mag2);title('带通切比雪夫Ⅰ型型幅度响应');grid;subplot(223);plot(w3/pi,mag3);title('带通切比雪夫Ⅱ型型幅度响应');grid;subplot(224);plot(w4/pi,mag4);title('带通椭圆型幅度响应');grid;0.5101234带通巴特沃斯幅度响应00.510.0050.010.015带通切比雪夫Ⅰ型幅度响应00.510.0050.010.015带通切比雪夫Ⅱ型幅度响应00.510.51带通椭圆型幅度响应4.用滤波器对信号进行滤波 （1）、FIR 低通滤波器滤波-4-224x 10401234x 10-3滤波前的频谱-4-224x 10401234x 10-3滤波后的频谱00.51 1.52x 106-1-0.500.51滤波前的时域波形00.51 1.52x 106-2-1012滤波后的时域波形（2）、FIR 高通滤波器滤波-4-2024x 10400.0050.010.0150.02滤波前的频谱-4-224x 104051015x 10225滤波后的频谱00.51 1.52x 106-1-0.500.51滤波前的时域波形00.51 1.52x 106-4-224x 10228滤波后的时域波形（3）、FIR 带通滤波器滤波-4-224x 10401234x 10-3滤波前的频谱-4-224x 1040246x 10158滤波后的频谱00.51 1.52x 106-1-0.500.51滤波前的时域波形00.51 1.52x 106-4-224x 10162滤波后的时域波形（4）、IIR 低通滤波器滤波-4-224x 10401234x 10-3滤波前的频谱-4-2024x 10401234x 10-3椭圆低通滤波后的频谱0.51 1.52x 106-1-0.500.51滤波前的时域波形00.51 1.52x 106-2-1012滤波前的频谱（5）、IIR 高通滤波器滤波-4-224x 10401234x 10-3滤波前的频谱-4-2024x 10402468x 10-5椭圆高通滤波后的频谱00.51 1.52x 106-1-0.500.51滤波前的时域波形00.51 1.52x 106-0.4-0.200.20.4滤波前的频谱（6）、IIR 带通滤波器滤波-4-224x 10401234x 10-3滤波前的频谱-4-2024x 10400.511.5x 10-5椭圆带通滤波后的频谱0.51 1.52x 106-1-0.500.51滤波前的时域波形00.51 1.52x 106-0.0200.020.04滤波前的频谱分析：由于巴特沃斯型、切比雪夫型滤波器的阶数太高，所以使用椭圆型IIR 滤波器进行滤波。