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

基于基于 IIR 的有噪声语音信号处理的有噪声语音信号处理 第第 1 1 章章 绪论绪论 1 1 第第 2 2 章章 基于滤波器的语音信号的处理基于滤波器的语音信号的处理 2 2 2.12.1 语音信号的采集语音信号的采集 2 2 2.22.2 设计设计 IIRIIR 数字滤波器数字滤波器2 2 2.32.3 用滤波器对加噪语音信号进行滤波用滤波器对加噪语音信号进行滤波.7.7 第第 3 3 章章 仿真及其结果分析仿真及其结果分析 7 7 3 3. .1 1 语音信号的时频分析语音信号的时频分析 7 7 3.23.2 加噪后的语音信号及其频谱分析加噪后的语音信号及其频谱分析 9 9 3.33.3 比较滤波前后语音信号的波形及频谱比较滤波前后语音信号的波形及频谱.10.10 3.43.4 I II IR R 滤滤波器波器处处理结果的的分析理结果的的分析1616 第二章：基于滤波器的语音信号的处理第二章：基于滤波器的语音信号的处理 选择一个语音信号作为分析的对象，或录制一段语音信号，对其进行频谱分析；利 用 MATLAB 中的随机函数产生噪声加入到语音信号中，模仿语音信号被污染，并对其进行 频谱分析；设计 IIR 数字滤波器， 并对被噪声污染的语音信号进行滤波， 分析滤波后 信号的时域和频域特征，最后回放语音信号。 2.1 语音信号的采集 利用WINDOWS 操作系统可以进行数字信号的采集。将话筒输入计算机的语音输入插 口上,启动录音机。按下录音按钮,接着对话筒说话“语音信号处理”,说完后停止录音, 屏幕左侧将显示所录声音的长度。点击放音按钮,可以实现所录音的重现。以文件名 1 “speech”保存入g : MATLAB work 中。可以看到,文件存储器的后缀默认为. wav , 这是WINDOWS 操作系统规定的声音文件存的标准。 1、原始信号的采集及分析 源程序： y1,fs,bits=wavread(D:MATLABai2.wav); sound(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 滤波器设计方法有间接法和直接法，间接法是借助于模拟滤波器的设计方法进行 的。其设计步骤是：先设计过渡模拟滤波器得到系统函数 H（s） ，然后将 H（s）按某种 方法转换成数字滤波器的系统函数 H(z)。FIR 滤波器比鞥采用间接法，常用的方法有窗 函数法、频率采样发和切比雪夫等波纹逼近法。对于线性相位滤波器，经常采用 FIR 滤 2 波器。 对于数字高通、带通滤波器的设计，通用方法为双线性变换法。可以借助于模拟滤 波器的频率转换设计一个所需类型的过渡模拟滤波器，再经过双线性变换将其转换策划 那个所需的数字滤波器。具体设计步骤如下： （1）确定所需类型数字滤波器的技术指标。 （2）将所需类型数字滤波器的边界频率转换成相应的模拟滤波器的边界频率，转换 公式为 =2/T tan(0.5) (3)将相应类型的模拟滤波器技术指标转换成模拟低通滤波器技术指标。 （4）设计模拟低通滤波器。 （5）通过频率变换将模拟低通转换成相应类型的过渡模拟滤波器。 （6）采用双线性变换法将相应类型的过渡模拟滤波器转换成所需类型的数字滤波器。 我们知道，脉冲响应不变法的主要缺点是会产生频谱混叠现象，使数字滤波器的频 响偏离模拟滤波器的频响特性。为了克服之一缺点，可以采用双线性变换法。 下面我们总结一下利用模拟滤波器设计 IIR 数字低通滤波器的步骤： （1）确定数字低通滤波器的技术指标：通带边界频率、通带最大衰减，阻带截止频 率、阻带最小衰减。 （2）将数字低通滤波器的技术指标转换成相应的模拟低通滤波器的技术指标。 （3）按照模拟低通滤波器的技术指标设计及过渡模拟低通滤波器。 （4）用双线性变换法，模拟滤波器系统函数转换成数字低通滤波器系统函数。 如前所述，IIR 滤波器和 FIR 滤波器的设计方法有很大的区别。下面我们着重介绍 用窗函数法设计 FIR 滤波器的步骤。如下： （1）根据对阻带衰减及过渡带的指标要求，选择串窗数类型（矩形窗、三角窗、汉 宁窗、哈明窗、凯塞窗等） ，并估计窗口长度 N。先按照阻带衰减选择窗函数类型。原则 是在保证阻带衰减满足要求的情况下，尽量选择主瓣的窗函数。 （2）构造希望逼近的频率响应函数。 （3）计算 h(n).。 （4）加窗得到设计结果。 接下来，我们根据语音信号的特点给出有关滤波器的技术指标： 低通滤波器的性能指标： 3 fp=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 默认为哈明窗。其中可选的窗函数有 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=cheby1(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 域的频率响应的参数 num11,den11=bilinear(b11,a11,0.5); %双线性变换实现 S 域到 Z 域的变换 h,w=freqz(num11,den11); %根据参数求出频率响应 plot(w*8000*0.5/pi,abs(h); legend(用 butter 设计); 4 图 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); n12,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(用 butter 设计); 5 图 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=bilinear(num,den,0.5); %利用双线性变换实现 S 域到 Z 域转换 h,w=freqz(num13,den13); plot(w*21000*0.5/pi,abs(h); title(IIR 高通滤波器); legend(用 cheby1 设计); 6 图 5 IIR 高通滤波器 2.3 用滤波器对加噪语音信号进行滤波： 用自己设计的各滤波器分别对加噪的语音信号进行滤波，在 Matlab 中，IIR 滤波器 利用函数 filter 对信号进行滤波。 函数 fftfilt 用的是重叠相加法实现线性卷积的计算。调用格式为： y=fftfilter(h,x,M)。其中，h 是系统单位冲击响应向量；x 是输入序列向量；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,nbits=wavread(file) 功能说明：采样值放在向量 y 中，fs 表示采样频率（hz） ，nbits 表示采样位数。 7 （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 的函数 fft,其调用格式如下： 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; 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); 程序结果如下图： 8 图5 原始信号波形及频谱 3.2加噪后的语音信号及其频谱分析 本文中，利用 MATLAB 中的随机函数(rand 或 randn)产生噪声加入到语音信号中，模 仿语音信号被污染，并对其频谱分析。Randn 函数有两种基本调用格式：Randn(n)和 Randn（m,n）,前者产生 nn 服从标准高斯分布的随机数矩阵，后者产生 mn 的随机数 矩阵。在这里，我们选用 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); %傅里叶变换 subplot(2,1,2); plot(abs(S); title(加噪语音信号的频域波形); 程序结果如下图： 9 图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(z11); %求滤波后的信号 subplot(2,2,1); plot(abs(S),g); title(滤波前信号的频谱); grid; subplot(2,2,2); plot(abs(m11),r); title(滤波后信号的频谱); 10 grid; subplot(2,2,3); plot(s); title(滤波前信号的波形); grid; subplot(2,2,4); plot(z11); title(滤波后的信号波形); 图 9 双线性法低通滤波 %*带通滤波器* y,fs,nbits=wavread (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(滤波后信号的频谱); subplot(2,2,3); plot(s); title(滤波前信号的波形); subplot(2,2,4); 11 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); m13=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(滤波后的信号波形); 12 图11 双线性法高通滤波 %=窗函数法= %*低通滤波器* y,fs,nbits=wavread (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(滤波前信号的波形);