语音信号处理与FIR IIR滤波器

1、数字信号处理课程设计10电子信息科学与技术 赵祎 10380051一、实验目的1学会MATLAB的使用，掌握MATLAB的程序设计方法；2掌握在Windows环境下语音信号采集的方法；3掌握数字信号处理的基本概念、基本理论和基本方法；4掌握MATLAB设计FIR和IIR数字滤波器的方法；5学会用MATLAB对信号进行分析和处理。二、实验原理参考数字信号处理教材，MATLAB的signal processing toolbox.三、主要实验仪器及材料微型计算机、Matlab软件。四、实验内容（一）语音信号的采集及频谱分析使用windows的录音功能或其他软件录制一段话音，时间控制在5秒左右。然后

2、在MATLAB软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。通过wavread函数的使用，要求理解采样频率、采样位数等概念。在录制的语音信号中加入高斯白噪声。对于无高斯白噪声影响和有高斯白噪声影响的语音信号：画出语音信号的时域波形；然后对语音信号进行频谱分析，在MATLAB中，可以利用函数FFT对信号进行快速傅里叶变换，得到信号的频谱特性。（二）设计FIR和IIR数字滤波器并滤波基于FIR（滤波器窗函数法，包括Hamming窗、Kaiser窗）和IIR（双线性变换法设计巴特沃兹滤波器、切比雪夫滤波器I型滤波器）的设计方法，分别设计以下三种数字滤波器，根据语音信号

3、的特点给出有关滤波器的性能指标：1）低通滤波器性能指标，fp=1000Hz，fc=1200 Hz，As=100dB，Ap=1dB；2）高通滤波器性能指标，fc=2800 Hz，fp=3000 Hz As=100dB，Ap=1dB；3）带通滤波器性能指标，fp1=1200 Hz，fp2=3000 Hz，fc1=1000 Hz，fc2=3200 Hz，As=100dB，Ap=1dB。利用MATLAB软件工具，设计上面要求的三种滤波器：第一步、设计出各数字滤波器并画出滤波器的频率特性曲线；第二步、比较FIR和IIR两种滤波器的性能，在FIR滤波器和IIR滤波器中分别选择性能好的一项，然后用各滤波器分

4、别对采集的语音信号进行滤波。第三步、比较滤波前后语音信号的波形及频谱，并回放语音信号。y, fs, bits = wavread( 'E:ringsmiao.wav' );五、实验代码1. 语音信号的采集与噪声的加入代码如下：%打开语音信号并绘制波形和频谱（见右图）y, fs, bits = wavread( 'E:ringsmiao.wav' );sound(y, fs, bits); n = length (y) ; Y = fft(y, n); subplot(2, 1, 1);plot(y);title( '语音信号波形' );subpl

5、ot(2, 1, 2);plot( abs(Y) );title( '语音信号频谱' );%在原语音信号中加入加性噪声并画出含噪信号波形和频谱（见右图）z = awgn(y, 10);Y2 = fft(z, n); Y3 = fftshift( Y2 ); f = 0 : fs / n : fs * (n - 1) / n;subplot(2, 1, 1);plot(z);title( '含噪信号波形' );subplot(2, 1, 2);plot( f, abs(Y3) ); title( '含噪信号频谱' );sound(z, fs, bi

6、ts);wavwrite(z,'E:miao2.wav');2. 滤波器的设计(1)低通滤波器 FIR低通（凯撒窗）% FIR凯撒窗口低通滤波器%参数Ft = 8000; Fp = 1000; Fc = 1200;wp = 2 * Fp / Ft; wc = 2 * Fc / Ft; p = 1-10.( -1 / 20 );s = 10.( -100 / 20 );%凯撒窗滤波fpts = wp wc;mag = 1 0;dev = p s;n21, wc21, beta, ftype = kaiserord(fpts, mag, dev);b21 = fir1( n21,

7、wc21, kaiser (n21+1, beta) );h, w = freqz(b21, 1);%画图plot(w * 8000 * 0.5 / pi, abs(h);title('FIR低通滤波器（凯撒窗口）');xlabel( 'Hz' ); ylabel( '幅度' );grid; FIR低通（hamming窗）%FIR hamming窗口低通滤波器%参数Ft = 8000; Fp = 1000; Fs = 1200;wp = 2 * pi * Fp / Ft;ws = 2 * pi * Fs / Ft;wc = (wp + ws) /

8、 2;B = ws - wp;%hamming窗滤波N = ceil(6.6 * pi / B) + 1;n = 0 : 1 : N-1;hd = ideal_lp(wc, N);w_ham = ( hamming(N) )' hn = hd.* w_ham;hw, w = freqz(hn, 1);%画图plot(w * 8000 * 0.5 / pi, abs(hw);title('FIR低通滤波器（hamming）');xlabel( 'Hz' ); ylabel( '幅度' );grid;两个窗函数低通滤波器的频率特性图如下：由图

9、中可以看出，凯撒窗低通滤波器的过渡带更窄，且阻带波纹更平缓，故选用FIR凯撒窗低通滤波器对信号进行滤波。 IIR低通（巴特沃斯）%IIR巴特沃斯低通滤波器%参数Ft = 8000; Fp = 1000; Fc = 1200;wp = 2 * pi * Fp / Ft;wc = 2 * pi * Fc / Ft;fp = tan(wp / 2);fc = tan(wc / 2);%巴特沃斯滤波n11, wn11 = buttord(fp, fc, 1, 50, 's');b11, a11 = butter(n11, wn11, 's'); num11, den11

10、 = bilinear(b11, a11, 0.5); h, w = freqz(num11, den11); %画图plot(w * 8000 * 0.5 / pi, abs(h);title('IIR低通滤波器（巴特沃斯）');xlabel( 'Hz' ); ylabel( '幅度' );grid; IIR低通（契比雪夫）%IIR契比雪夫低通滤波器%参数Ft = 8000; Fp = 1000; Fs = 1200;wp1 = tan(pi * Fp / Ft);ws1 = tan(pi * Fs / Ft);%契比雪夫滤波n12, wn12

11、 = cheb1ord(wp1, ws1, 1, 100, 's');b12, a12 = cheby1(n12, 1, wn12, 's');num12, den12 = bilinear(b12, a12, 0.5); h, w = freqz(num12, den12);%画图plot(w * 8000 * 0.5 / pi, abs(h);title('IIR高通滤波器（契比雪夫）');xlabel('Hz'); ylabel('幅度');grid;两个双线性变换低通滤波器的频率特性图如下：由图中可以看出，

12、虽然切比雪夫滤波器的过渡带更窄，但是阻带波纹的波形波动过大，故选用IIR巴特沃夫低通滤波器对信号进行滤波。(2) 高通滤波器 FIR高通（凯撒窗）% FIR凯撒窗口高通滤波器%参数Ft = 8000; Fp = 3000; Fc = 2800; wp = 2 * Fp / Ft; wc = 2 * Fc / Ft; p = 1 - 10.( -1 / 20); s = 10.( -100 / 20);%凯撒窗滤波fpts = wc wp;mag = 0 1;dev = p s;n22, wc22, beta2, ftype2 = kaiserord(fpts, mag, dev);b22 =

13、fir1(n22, wc22, ftype2, kaiser(n22+1, beta2); h, w = freqz(b22, 1); %画图 plot(w * 8000 * 0.5 / pi, abs(h);title('FIR高通滤波器（凯撒窗口）');xlabel( 'Hz' ); ylabel( '幅度' );grid; FIR高通（hamming）%FIR hamming窗口高通滤波器%参数Ft = 8000; Fp = 3000; Fs = 2800;wp = 2 * pi * Fp / Ft;ws = 2 * pi * Fs / F

14、t;wc = (ws + B / 2) / pi;B = wp - ws;%hamming窗滤波N = ceil(8 * pi / B);window = hamming(N+1);hn = fir1(N, wc,'high', window);n = 0 : N-1;hw, w = freqz(hn, 1);%画图plot(w * 8000 * 0.5 / pi, abs(hw);title('FIR高通滤波器（hamming）');xlabel( 'Hz' ); ylabel( '幅度' );grid;两个窗函数高通滤波器的频

15、率特性图如下：由图中可以看出，凯撒窗高通滤波器的过渡带更窄，且阻带波纹更平缓，故选用FIR凯撒窗低通滤波器对信号进行滤波。 IIR高通（巴特沃斯）%IIR巴特沃斯高通滤波器%参数Ft = 8000; Fp = 3000; Fs = 2800;wp = 2 * pi * Fp / Ft;ws = 2 * pi * Fs / Ft;fp = tan(wp / 2);fs = tan(ws / 2);%巴特沃斯滤波n22, wn22 = buttord(fp, fs, 1, 50, 's');b22, a22 = butter(n22, wn22, 's'); %低通

16、原型BT, AT = lp2hp(b22, a22, 2 * ws);num, den = bilinear(BT, AT, 0.5); %双线性变换H,W = freqz(num, den);%画图plot(W * 8000 * 0.5 / pi, abs(H);title('IIR高通滤波器（巴特沃斯）');xlabel( 'Hz' ); ylabel( '幅度' );grid; IIR高通（契比雪夫）%IIR契比雪夫高通滤波器%参数Ft = 8000; Fp = 3000; Fs = 2800;wp1 = tan(pi * Fp / Ft)

17、; ws1 = tan(pi * Fs / Ft);%契比雪夫滤波n12, wn12 = cheb1ord(wp1, ws1, 1, 100, 's'); b12, a12 = cheby1(n12, 1, wn12, 's'); num, den = lp2hp(b12, a12, 2 * ws1);num12, den12 = bilinear(num, den, 0.5); h, w = freqz(num12, den12);%画图plot(w * 8000 * 0.5 / pi, abs(h);title( 'IIR高通滤波器（契比雪夫）

18、9; );xlabel('Hz'); ylabel('幅度');grid;两个双线性高通滤波器的频率特性图如下：由图中可以看出，虽然切比雪夫滤波器的过渡带更窄，但是阻带波纹的波形波动过大，故选用IIR巴特沃夫高通滤波器对信号进行滤波。(3) 带通滤波器 FIR带通（凯撒窗）% FIR凯撒窗口带通滤波器%参数Fp1 = 1200; Fc1 = 1000;Fp2 = 3000; Fc2 = 3200;wp1 = 2 * Fp1 / Ft;wc1 = 2 * Fc1 / Ft;wp2 = 2 * Fp2 / Ft;wc2 = 2 * Fc2 / Ft;p = 1 -

19、10.( -1 / 20); s = 10.( -100 / 20);%凯撒窗滤波fpts2 = wc1 wp1 wp2 wc2;mag1 = 0 1 0;dev1 = s p s;n23, wc23, beta3, ftype3 = kaiserord(fpts2, mag1, dev1);b23 = fir1(n23, wc23, ftype3, kaiser(n23+1, beta3);h, w = freqz(b23, 1);%画图plot(w * 8000 * 0.5 / pi, abs(h);title('FIR带通滤波器（凯撒窗口）');xlabel( '

20、Hz' ); ylabel( '幅度' );grid; FIR带通（hamming）%FIR hamming窗口带通滤波器%参数Ft = 8000;Fp1 = 1200; Fs1 = 1000;Fp2 = 3000; Fs2 = 3200;wp1 = 2 * pi * Fp1 / Ft;ws1 = 2 * pi * Fs1 / Ft;wp2 = 2 * pi * Fp2 / Ft;ws2 = 2 * pi * Fs2 / Ft;wc1 = (wp1 + ws1) / 2;wc2 = (wp2 + ws2) / 2;As = 100;B = wp1 - ws1;%hamm

21、ing窗滤波N = ceil(6.6 * pi / B) + 1;n = 0 : 1 : N-1;hd = ideal_lp(wc2, N) - ideal_lp(wc1, N);w_ham = ( hamming(N) )'hn = hd.* w_ham; hw, w = freqz(hn, 1);%画图plot(w * 8000 * 0.5 / pi, abs(hw);title('FIR带通滤波器（hamming）');xlabel( 'Hz' ); ylabel( '幅度' ); grid;两个窗函数带通滤波器的频率特性图如下：由

22、图中可以看出，凯撒窗高通滤波器的过渡带更窄，且阻带波纹更平缓，故选用FIR凯撒窗带通滤波器对信号进行滤波。 IIR带通（巴特沃斯）%IIR巴特沃斯带通滤波器%参数Ft = 8000;Fp1 = 1200; Fs1 = 1000;Fp2 = 3000; Fs2 = 3200;wp1 = tan(pi * Fp1 / Ft); wp2 = tan(pi * Fp2 / Ft);ws1 = tan(pi * Fs1 / Ft);ws2 = tan(pi * Fs2 / Ft);w = wp1 * wp2 / ws2;bw = wp2 - wp1;wp = 1;ws = (wp1 * wp2 - w.

23、2)/(bw * w);%巴特沃斯滤波n12, wn12 = buttord(wp, ws, 1, 50, 's');b12, a12 = butter(n12, wn12, 's');num2, den2 = lp2bp(b12, a12, sqrt(wp1 * wp2), bw);num12, den12 = bilinear(num2, den2, 0.5);h, w = freqz(num12, den12);%画图plot(w * 8000 * 0.5 / pi, abs(h);title('IIR带通滤波器（巴特沃斯）');xlabe

24、l( 'Hz' ); ylabel( '幅度' );grid; IIR带通（契比雪夫）%IIR契比雪夫带通滤波器%参数Ft = 8000;Fp1 = 1200; Fs1 = 1000;Fp2 = 3000; Fs2 = 3200;wp1 = tan(pi * Fp1 / Ft);wp2 = tan(pi * Fp2 / Ft);ws1 = tan(pi * Fs1 / Ft);ws2 = tan(pi * Fs2 / Ft);w = wp1 * wp2 / ws2;bw = wp2 - wp1;ws = (wp1 * wp2 - w.2) / (bw * w);

25、%契比雪夫滤波n12, wn12 = cheb1ord(1, ws, 1, 100, 's');b12, a12 = cheby1(n12, 1, wn12, 's');num2, den2 = lp2bp(b12, a12, sqrt(wp1 * wp2), bw);num12, den12 = bilinear(num2, den2, 0.5);h, w = freqz(num12, den12);%画图plot(w * 8000 * 0.5 / pi, abs(h);axis(0 4000 0 1.5);title('IIR带通滤波器（契比雪夫）&

26、#39;);xlabel('Hz'); ylabel('幅度');grid;两个双线性带通滤波器的频率特性图如下：由图中可以看出，虽然切比雪夫滤波器的过渡带更窄，但是阻带波纹的波形波动过大，故选用IIR巴特沃夫带通滤波器对信号进行滤波。3.使用选择的滤波器进行滤波（1）FIR低通滤波%打开含噪语音信号y, fs, bits = wavread( 'E:miao2.wav' );sound(y, fs, bits); n = length (y) ; Y = fft(y, n);%凯撒窗滤波器Ft = 8000; Fp = 1000; Fc = 1

27、200;wp = 2 * Fp / Ft; wc = 2 * Fc / Ft; p = 1-10.( -1 / 20 );s = 10.( -100 / 20 );fpts = wp wc;mag = 1 0;dev = p s;n21, wc21, beta, ftype = kaiserord(fpts, mag, dev);b21 = fir1( n21, wc21, kaiser (n21+1, beta) );h, w = freqz(b21, 1);%滤波前后的图像比较z11 = fftfilt(b21, y);sound( z11 );m11 = fft( z11 ); %求滤波

28、后的信号subplot(2, 2, 1);plot( abs(Y);title( '滤波前信号频谱' );grid;subplot(2, 2, 2);plot( abs(m11) );title( '滤波后信号频谱' );grid;subplot(2, 2, 3);plot(y);title( '滤波前信号波形' );grid;subplot(2, 2, 4);plot( z11 );title( '滤波后的信号波形' );grid;wavwrite(y, 'E:FIR低通.wav');（2）IIR低通滤波%打开含

29、噪语音信号y, fs, bits = wavread( 'E:miao2.wav' );sound(y, fs, bits); n = length (y) ; Y = fft(y, n);%巴特沃斯滤波Ft = 8000; Fp = 1000; Fc = 1200;wp = 2 * pi * Fp / Ft;wc = 2 * pi * Fc / Ft;fp = tan(wp / 2);fc = tan(wc / 2);n11, wn11 = buttord(fp, fc, 1, 50, 's');b11, a11 = butter(n11, wn11, 

30、9;s'); num11, den11 = bilinear(b11, a11, 0.5); h, w = freqz(num11, den11); %滤波前后的图像比较z11= filter(num11, den11, y);sound( z11 );m11 = fft( z11 ); %求滤波后的信号subplot(2, 2, 1);plot( abs(Y);title( '滤波前信号频谱' );grid;subplot(2, 2, 2);plot( abs(m11) );title( '滤波后信号频谱' );grid;subplot(2, 2, 3

31、);plot(y);title( '滤波前信号波形' );grid;subplot(2, 2, 4);plot( z11 );title( '滤波后的信号波形' );grid;wavwrite(y, 'E:IIR低通.wav');（3）FIR高通滤波%打开含噪语音信号y, fs, bits = wavread( 'E:miao2.wav' );sound(y, fs, bits); n = length (y) ; Y = fft(y, n);% FIR凯撒窗口高通滤波器Ft = 8000; Fp = 3000; Fc = 280

32、0; wp = 2 * Fp / Ft; wc = 2 * Fc / Ft; p = 1 - 10.( -1 / 20); s = 10.( -100 / 20);fpts = wc wp;mag = 0 1;dev = p s;n22, wc22, beta2, ftype2 = kaiserord(fpts, mag, dev);b22 = fir1(n22, wc22, ftype2, kaiser(n22+1, beta2); h, w = freqz(b22, 1); 滤波前后的图像比较z22 = fftfilt(b22, y);%sound( z22 );m22 = fft( z2

33、2 ); %求滤波后的信号subplot(2, 2, 1);plot( abs(Y);title( '滤波前信号频谱' );grid;subplot(2, 2, 2);plot( abs(m22) );title( '滤波后信号频谱' );grid;subplot(2, 2, 3);plot(y);title( '滤波前信号波形' );grid;subplot(2, 2, 4);plot( z22 );title( '滤波后的信号波形' );grid;wavwrite(y, 'E:FIR高通.wav');（4）II

34、R高通滤波%打开含噪语音信号y, fs, bits = wavread( 'E:miao2.wav' );sound(y, fs, bits); n = length (y) ; Y = fft(y, n);%IIR巴特沃斯高通滤波器Ft = 8000; Fp = 3000; Fs = 2800;wp = 2 * pi * Fp / Ft;ws = 2 * pi * Fs / Ft;fp = tan(wp / 2);fs = tan(ws / 2);n22, wn22 = buttord(fp, fs, 1, 50, 's');b22, a22 = butter

35、(n22, wn22, 's'); %低通原型BT, AT = lp2hp(b22, a22, 2 * ws);num, den = bilinear(BT, AT, 0.5); %双线性变换H,W = freqz(num, den);%滤波前后的图像比较z22 =filter(num, den, y);sound( z22 );m22 = fft( z22 ); %求滤波后的信号subplot(2, 2, 1);plot( abs(Y);title( '滤波前信号频谱' );grid;subplot(2, 2, 2);plot( abs(m22) );titl

36、e( '滤波后信号频谱' );grid;subplot(2, 2, 3);plot(y);title( '滤波前信号波形' );grid;subplot(2, 2, 4);plot( z22 );title( '滤波后的信号波形' );grid;wavwrite(y, 'E:IIR高通.wav');（5）FIR带通滤波%打开含噪语音信号y, fs, bits = wavread( 'E:miao2.wav' );sound(y, fs, bits); n = length (y) ; Y = fft(y, n);%

37、 FIR凯撒窗口带通滤波器Fp1 = 1200; Fc1 = 1000;Fp2 = 3000; Fc2 = 3200;wp1 = 2 * Fp1 / Ft;wc1 = 2 * Fc1 / Ft;wp2 = 2 * Fp2 / Ft;wc2 = 2 * Fc2 / Ft;p = 1 - 10.( -1 / 20); s = 10.( -100 / 20);fpts2 = wc1 wp1 wp2 wc2;mag1 = 0 1 0;dev1 = s p s;n23, wc23, beta3, ftype3 = kaiserord(fpts2, mag1, dev1);b23 = fir1(n23,

38、wc23, ftype3, kaiser(n23+1, beta3);h, w = freqz(b23, 1);%滤波前后的图像比较z23 = fftfilt(b23, y);sound( z23 );m23 = fft( z23 ); %求滤波后的信号subplot(2, 2, 1);plot( abs(Y);title( '滤波前信号频谱' );grid;subplot(2, 2, 2);plot( abs(m23) );title( '滤波后信号频谱' );grid;subplot(2, 2, 3);plot(y);title( '滤波前信号波形' );grid;subplot(2, 2, 4);plot( z23 );title( '滤波后的信号波形' );grid;wavwrite(y, 'E:FIR带通.wav')（6）IIR带通滤波%打开含噪语音信号