课程四设计_基于IIR的有噪声语音信号处理

1、.基于IIR的有噪声语音信号处理第1章 绪论1第2章 基于滤波器的语音信号的处理22.1语音信号的采集22.2设计IIR数字滤波器22.3用滤波器对加噪语音信号进行滤波.7第3章 仿真及其结果分析73.1 语音信号的时频分析73.2 加噪后的语音信号及其频谱分析93.3 比较滤波前后语音信号的波形及频谱.103.4 IIR滤波器处理结果的的分析16第二章：基于滤波器的语音信号的处理选择一个语音信号作为分析的对象，或录制一段语音信号，对其进行频谱分析；利用MATLAB中的随机函数产生噪声加入到语音信号中，模仿语音信号被污染，并对其进行频谱分析；设计IIR数字滤波器， 并对被噪声污染的语音信号进行

2、滤波， 分析滤波后信号的时域和频域特征，最后回放语音信号。2.1 语音信号的采集利用WINDOWS 操作系统可以进行数字信号的采集。将话筒输入计算机的语音输入插口上,启动录音机。按下录音按钮,接着对话筒说话“语音信号处理”,说完后停止录音,屏幕左侧将显示所录声音的长度。点击放音按钮,可以实现所录音的重现。以文件名“speech”保存入g : MATLAB work 中。可以看到,文件存储器的后缀默认为. wav ,这是WINDOWS 操作系统规定的声音文件存的标准。 1、原始信号的采集及分析 源程序：y1,fs,bits=wavread('D:MATLABai2.wav');s

3、ound(y1,44000); %播放语音信号y=fft(y1); %对信号做FFT变换f=fs*(0:511)/1024;figure(1)subplot(2,1,1);plot(abs(y(1:512) %做原始语音信号的FFT频谱图title('原始语音信号FFT频谱')subplot(2,1,2);%做原始语音信号的时域图形plot(y1)title('原始语音信号波形');xlabel('时间 n');ylabel('幅值 n');2.2设计IIR数字滤波器IIR滤波器设计方法有间接法和直接法，间接法是借助于模拟滤波器的

4、设计方法进行的。其设计步骤是：先设计过渡模拟滤波器得到系统函数H（s），然后将H（s）按某种方法转换成数字滤波器的系统函数H(z)。FIR滤波器比鞥采用间接法，常用的方法有窗函数法、频率采样发和切比雪夫等波纹逼近法。对于线性相位滤波器，经常采用FIR滤波器。 对于数字高通、带通滤波器的设计，通用方法为双线性变换法。可以借助于模拟滤波器的频率转换设计一个所需类型的过渡模拟滤波器，再经过双线性变换将其转换策划那个所需的数字滤波器。具体设计步骤如下：（1）确定所需类型数字滤波器的技术指标。（2）将所需类型数字滤波器的边界频率转换成相应的模拟滤波器的边界频率，转换公式为=2/T tan(0.5)(3)

5、将相应类型的模拟滤波器技术指标转换成模拟低通滤波器技术指标。（4）设计模拟低通滤波器。（5）通过频率变换将模拟低通转换成相应类型的过渡模拟滤波器。（6）采用双线性变换法将相应类型的过渡模拟滤波器转换成所需类型的数字滤波器。我们知道，脉冲响应不变法的主要缺点是会产生频谱混叠现象，使数字滤波器的频响偏离模拟滤波器的频响特性。为了克服之一缺点，可以采用双线性变换法。下面我们总结一下利用模拟滤波器设计IIR数字低通滤波器的步骤：（1）确定数字低通滤波器的技术指标：通带边界频率、通带最大衰减，阻带截止频率、阻带最小衰减。（2）将数字低通滤波器的技术指标转换成相应的模拟低通滤波器的技术指标。（3）按照模拟

6、低通滤波器的技术指标设计及过渡模拟低通滤波器。（4）用双线性变换法，模拟滤波器系统函数转换成数字低通滤波器系统函数。如前所述，IIR滤波器和FIR滤波器的设计方法有很大的区别。下面我们着重介绍用窗函数法设计FIR滤波器的步骤。如下：（1）根据对阻带衰减及过渡带的指标要求，选择串窗数类型（矩形窗、三角窗、汉宁窗、哈明窗、凯塞窗等），并估计窗口长度N。先按照阻带衰减选择窗函数类型。原则是在保证阻带衰减满足要求的情况下，尽量选择主瓣的窗函数。（2）构造希望逼近的频率响应函数。（3）计算h(n).。（4）加窗得到设计结果。接下来，我们根据语音信号的特点给出有关滤波器的技术指标：低通滤波器的性能指标：f

7、p=1000Hz，fc=1200Hz，As=100db ,Ap=1dB高通滤波器的性能指标：fp=3500Hz，fc=4000Hz，As=100dB，Ap=1dB；带通滤波器的性能指标：fp1=1200Hz，fp2=3000hZ，fc1=1000Hz，fc2=3200Hz，As=100dB，Ap=1dB在Matlab中，可以利用函数fir1设计FIR滤波器，利用函数butter,cheby1和ellip设计IIR滤波器，利用Matlab中的函数freqz画出各步步器的频率响应。hn=fir1(M，wc，window)，可以指定窗函数向量window。如果缺省window参数，则fir1默认为哈

8、明窗。其中可选的窗函数有Rectangular Barlrtt Hamming Hann Blackman窗，其相应的都有实现函数。MATLAB信号处理工具箱函数buttp buttor butter是巴特沃斯滤波器设计函数，其有5种调用格式，本课程设计中用到的是N,wc=butter(N,wc,Rp,As,s),该格式用于计算巴特沃斯模拟滤波器的阶数N和3dB截止频率wc。MATLAB信号处理工具箱函数cheblap,cheblord和cheeby1是切比雪夫I型滤波器设计函数。我们用到的是cheeby1函数，其调用格式如下：B,A=cheby1(N,Rp,wpo,ftypr)B,A=che

9、by1(N,Rp,wpo,ftypr,s) 函数butter,cheby1和ellip设计IIR滤波器时都是默认的双线性变换法，所以在设计滤波器时只需要代入相应的实现函数即可。下面我们将给出FIR和IIR数字滤波器的主要程序。%=IIR低通滤波器=wp=2*pi*Fp/Ft;ws=2*pi*Fs/Ft;fp=2*Ft*tan(wp/2);fs=2*Fs*tan(wp/2);n11,wn11=buttord(wp,ws,1,50,'s'); %求低通滤波器的阶数和截止频率b11,a11=butter(n11,wn11,'s'); %求S域的频率响应的参数 num1

10、1,den11=bilinear(b11,a11,0.5); %双线性变换实现S域到Z域的变换 h,w=freqz(num11,den11); %根据参数求出频率响应plot(w*8000*0.5/pi,abs(h);legend('用butter设计');图3 IIR低通滤波器%=IIR带通=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.2)/(bw*w);n

11、12,wn12=buttord(wp,ws,1,50,'s'); %求低通滤波器阶数和截止频率b12,a12=butter(n12,wn12,'s'); %求S域的频率响应参数num2,den2=lp2bp(b12,a12,sqrt(wp1*wp2),bw);%将S域低通参数转为带通的num12,den12=bilinear(num2,den2,0.5);%双线性变换实现S域到Z域的转换h,w=freqz(num12,den12); %根据参数求出频率响应plot(w*8000*0.5/pi,abs(h);axis(0 4000 0 1.5);legend(&#

12、39;用butter设计');图4 IIR带通滤波器%=IIR高通=Ft=8000;Fp=4000;Fs=3500;wp1=tan(pi*Fp/Ft); %高通到低通滤波器参数转换ws1=tan(pi*Fs/Ft);wp=1;ws=wp1*wp/ws1;n13,wn13=cheb1ord(wp,ws,1,50,'s'); %求模拟的低通滤波器阶数和截止频率b13,a13=cheby1(n13,1,wn13,'s'); %求S域的频率响应的参数num,den=lp2hp(b13,a13,wn13); %将S域低通参数转为高通的num13,den13=bil

13、inear(num,den,0.5); %利用双线性变换实现S域到Z域转换h,w=freqz(num13,den13);plot(w*21000*0.5/pi,abs(h);title('IIR高通滤波器');legend('用cheby1设计');图5 IIR 高通滤波器2.3用滤波器对加噪语音信号进行滤波：用自己设计的各滤波器分别对加噪的语音信号进行滤波，在Matlab中，IIR滤波器利用函数filter对信号进行滤波。函数fftfilt用的是重叠相加法实现线性卷积的计算。调用格式为：y=fftfilter(h,x,M)。其中，h是系统单位冲击响应向量；x是

14、输入序列向量；y是系统的输出序列向量；M是有用户选择的输入序列的分段长度，缺省时，默认的输入向量的重长度M=512。函数filter的调用格式：yn=filter(B，A.xn)，它是按照直线型结构实现对xn的滤波。其中xn是输入信号向量，yn输出信号向量。第三章 仿真及其结果分析3.1语音信号的时频分析利用MATLAB中的“wavread”命令来读入（采集）语音信号，将它赋值给某一向量。再对其进行采样，记住采样频率和采样点数。下面介绍Wavread 函数几种调用格式。（1）y=wavread（file）功能说明：读取file所规定的wav文件，返回采样值放在向量y中。（2）y,fs,nbit

15、s=wavread(file) 功能说明：采样值放在向量y中，fs表示采样频率（hz），nbits表示采样位数。（3）y=wavread（file，N）功能说明：读取钱N点的采样值放在向量y中。（4）y=wavread（file，N1,N2）功能说明：读取从N1到N2点的采样值放在向量y中。接下来，对语音信号OriSound.wav进行采样。其程序如下：>> y,fs,nbits=wavered (OriSound); %把语音信号加载入Matlab 仿真软件平台中然后，画出语音信号的时域波形，再对语音信号进行频谱分析。MATLAB提供了快速傅里叶变换算法FFT计算DFT的函数ff

16、t,其调用格式如下：Xk=fft(xn,N)参数xn为被变换的时域序列向量，N是DFT变换区间长度，当N大于xn的长度时，fft函数自动在xn后面补零。，当N小于xn的长度时，fft函数计算xn的前N个元素，忽略其后面的元素。在本次设计中，我们利用fft对语音信号进行快速傅里叶变换，就可以得到信号的频谱特性。原始信号程序如下：y1,fs,bits=wavread('C:Program FilesMATLABR2009aspeech.wav'); sound(y1,44000); %播放语音信号y=fft(y1); %对信号做FFT变换f=fs*(0:511)/1024;figu

17、re(1)subplot(2,1,1);plot(abs(y(1:512) %做原始语音信号的FFT频谱图title('原始语音信号FFT频谱')subplot(2,1,2);%做原始语音信号的时域图形plot(y1)title('原始语音信号波形');xlabel('时间 n');ylabel('幅值 n');程序结果如下图：图5 原始信号波形及频谱3.2加噪后的语音信号及其频谱分析本文中，利用MATLAB中的随机函数(rand或randn)产生噪声加入到语音信号中，模仿语音信号被污染，并对其频谱分析。Randn函数有两种基本调

18、用格式：Randn(n)和Randn（m,n）,前者产生n×n服从标准高斯分布的随机数矩阵，后者产生m×n的随机数矩阵。在这里，我们选用Randn（m,n）函数。加噪程序如下：y,fs,nbits=wavread ('D:MATLABai2.wav');N = length (y) ; %求出语音信号的长度Noise=0.01*randn(N,2); %随机函数产生噪声Si=y+Noise; %语音信号加入噪声 sound(Si);subplot(2,1,1);plot(Si);title('加噪语音信号的时域波形');S=fft(Si);

19、%傅里叶变换subplot(2,1,2);plot(abs(S); title('加噪语音信号的频域波形');程序结果如下图：图6 加噪后的波形及频谱分析3.3 比较滤波前后语音信号的波形及频谱% =双线性变换法=%*低通滤波器*y,fs,nbits=wavread ('speech');n = length (y) ; %求出语音信号的长度noise=0.01*randn(n,2); %随机函数产生噪声s=y+noise; %语音信号加入噪声S=fft(s); %傅里叶变换z11=filter(num11,den11,s);sound(z11);m11=fft

20、(z11); %求滤波后的信号subplot(2,2,1);plot(abs(S),'g');title('滤波前信号的频谱');grid;subplot(2,2,2);plot(abs(m11),'r');title('滤波后信号的频谱');grid;subplot(2,2,3);plot(s);title('滤波前信号的波形');grid;subplot(2,2,4);plot(z11);title('滤波后的信号波形');图9 双线性法低通滤波%*带通滤波器*y,fs,nbits=wavrea

21、d ('speech');n = length (y) ; %求出语音信号的长度noise=0.01*randn(n,2); %随机函数产生噪声s=y+noise; %语音信号加入噪声S=fft(s); %傅里叶变换z12=filter(num12,den12,s);sound(z12);m12=fft(z12); %求滤波后的信号subplot(2,2,1);plot(abs(S),'g');title('滤波前信号的频谱');subplot(2,2,2);plot(abs(m12),'r');title('滤波后信号的

22、频谱');subplot(2,2,3);plot(s);title('滤波前信号的波形');subplot(2,2,4);plot(z12);title('滤波后的信号波形');图10 双线性法带通滤波%*高通滤波器*y,fs,nbits=wavread ('speech');n = length (y) ; %求出语音信号的长度noise=0.01*randn(n,2); %随机函数产生噪声s=y+noise; %语音信号加入噪声S=fft(s); %傅里叶变换z13=filter(num13,den13,s);sound(z13);m

23、13=fft(z13); %求滤波后的信号subplot(2,2,1);plot(abs(S),'g');title('滤波前信号的频谱');subplot(2,2,2);plot(abs(m13),'r');title('滤波后信号的频谱');subplot(2,2,3);plot(s); title('滤波前信号的波形');subplot(2,2,4);plot(z13);title('滤波后的信号波形');图11 双线性法高通滤波%=窗函数法=%*低通滤波器*y,fs,nbits=wavrea

24、d ('speech');n = length (y) ; %求出语音信号的长度noise=0.01*randn(n,2); %随机函数产生噪声s=y+noise; %语音信号加入噪声 S=fft(s); %傅里叶变换z21=fftfilt(b21,s);sound(z21);m21=fft(z21); %求滤波后的信号subplot(2,2,1);plot(abs(S),'g');title('滤波前信号的频谱');subplot(2,2,2);plot(abs(m21),'r'); title('滤波后信号的频谱');subplot(2,2,3);plot(s); title('滤波前信号的波形');subplot(2,2,4);plot(z21); title('滤波后的信号波形'); 图12 窗函数法低通滤波%*带通滤波器*y,fs,nbits=wavread ('speech&