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

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

《课程四设计基于IIR的有噪声语音信号处理.docx》由会员分享，可在线阅读，更多相关《课程四设计基于IIR的有噪声语音信号处理.docx（20页珍藏版）》请在冰豆网上搜索。

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

《基于IIR的有噪声语音信号处理》

第二章：

基于滤波器的语音信号的处理

选择一个语音信号作为分析的对象，或录制一段语音信号，对其进行频谱分析；利用MATLAB中的随机函数产生噪声加入到语音信号中，模仿语音信号被污染，并对其进行频谱分析；设计IIR数字滤波器，并对被噪声污染的语音信号进行滤波，分析滤波后信号的时域和频域特征，最后回放语音信号。

2.1语音信号的采集

利用WINDOWS操作系统可以进行数字信号的采集。

将话筒输入计算机的语音输入插口上,启动录音机。

按下录音按钮,接着对话筒说话“语音信号处理”,说完后停止录音,屏幕左侧将显示所录声音的长度。

点击放音按钮,可以实现所录音的重现。

以文件名“speech”保存入g:

\MATLAB\work中。

可以看到,文件存储器的后缀默认为.wav,这是WINDOWS操作系统规定的声音文件存的标准。

1、原始信号的采集及分析

源程序：

[y1,fs,bits]=wavread（'D:

\MATLAB\ai2.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滤波器。

对于数字高通、带通滤波器的设计，通用方法为双线性变换法。

可以借助于模拟滤波器的频率转换设计一个所需类型的过渡模拟滤波器，再经过双线性变换将其转换策划那个所需的数字滤波器。

具体设计步骤如下：

（1）确定所需类型数字滤波器的技术指标。

（2）将所需类型数字滤波器的边界频率转换成相应的模拟滤波器的边界频率，转换公式为Ω=2/Ttan（0.5ω）

（3）将相应类型的模拟滤波器技术指标转换成模拟低通滤波器技术指标。

（4）设计模拟低通滤波器。

（5）通过频率变换将模拟低通转换成相应类型的过渡模拟滤波器。

（6）采用双线性变换法将相应类型的过渡模拟滤波器转换成所需类型的数字滤波器。

我们知道，脉冲响应不变法的主要缺点是会产生频谱混叠现象，使数字滤波器的频响偏离模拟滤波器的频响特性。

为了克服之一缺点，可以采用双线性变换法。

下面我们总结一下利用模拟滤波器设计IIR数字低通滤波器的步骤：

（1）确定数字低通滤波器的技术指标：

通带边界频率、通带最大衰减，阻带截止频率、阻带最小衰减。

（2）将数字低通滤波器的技术指标转换成相应的模拟低通滤波器的技术指标。

（3）按照模拟低通滤波器的技术指标设计及过渡模拟低通滤波器。

（4）用双线性变换法，模拟滤波器系统函数转换成数字低通滤波器系统函数。

如前所述，IIR滤波器和FIR滤波器的设计方法有很大的区别。

下面我们着重介绍用窗函数法设计FIR滤波器的步骤。

如下：

（1）根据对阻带衰减及过渡带的指标要求，选择串窗数类型（矩形窗、三角窗、汉宁窗、哈明窗、凯塞窗等），并估计窗口长度N。

先按照阻带衰减选择窗函数类型。

原则是在保证阻带衰减满足要求的情况下，尽量选择主瓣的窗函数。

（2）构造希望逼近的频率响应函数。

（3）计算h（n）.。

（4）加窗得到设计结果。

接下来，我们根据语音信号的特点给出有关滤波器的技术指标：

①低通滤波器的性能指标：

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默认为哈明窗。

其中可选的窗函数有RectangularBarlrttHammingHannBlackman窗，其相应的都有实现函数。

MATLAB信号处理工具箱函数buttpbuttorbutter是巴特沃斯滤波器设计函数，其有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设计'）;

图3IIR低通滤波器

%======================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（[0400001.5]）;

legend（'用butter设计'）;

图4IIR带通滤波器

%======================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设计'）;

图5IIR高通滤波器

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表示采样位数。

（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:

\ProgramFiles\MATLAB\R2009a\speech.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'）;

程序结果如下图：

图5原始信号波形及频谱

3.2加噪后的语音信号及其频谱分析

本文中，利用MATLAB中的随机函数（rand或randn）产生噪声加入到语音信号中，模仿语音信号被污染，并对其频谱分析。

Randn函数有两种基本调用格式：

Randn（n）和Randn（m,n）,前者产生n×n服从标准高斯分布的随机数矩阵，后者产生m×n的随机数矩阵。

在这里，我们选用Randn（m,n）函数。

加噪程序如下：

[y,fs,nbits]=wavread（'D:

\MATLAB\ai2.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（'加噪语音信号的频域波形'）;

程序结果如下图：

图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（'滤波后信号的频谱'）;

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）;

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（'滤波后的信号波形'）;

图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（'滤波前信号的波形'）;

subplot（2,2,4）;

plot（z21）;title（'滤波后的信号波形'）;

图12窗函数法低通滤波

%****************************带通滤波器***************************

[y,fs,nbits]=wavread（'speech'）;

n=length（y）;%求出语音信号的长度

noise=0.01*randn（n,2）;%随机函数产生噪声

s=y+noise;%语音信号加入噪声

S=fft（s）;%傅里叶变换

z22=fftfilt（b22,s）;

sound（z22）;

m22=fft（z22）;%求滤波后的信号

subplot（2,2,1）;

plot（abs（S）,'g'）;

title（'滤波前信号的频谱'）;

subplot（2,2,2）;

plot（abs（m22）,'r'）;

title（'滤波后信号的频谱'）;

subplot（2,2,3）;

plot（s）;

title（'滤波前信号的波形'）;

subplot（2,2,4）;

plot（z22）;

title（'滤波后的信号波形'）;

图13窗函数法带通滤波

%*************************高通滤波器*****************************

[y,fs,nbits]=wavread（'speech'）;

n=length（y）;%求出语音信号的长度

noise=0.01*randn（n,2）;%随机函数产生噪声

s=y+noise;%语音信号加入噪声

S=fft（s）;%傅里叶变换

z23=fftfilt（b23,s）;

sound（z23）;

m23=fft（z23）;%求滤波后的信号

subplot（2,2,1）;

plot（abs（S）,'g'）;

subplot（2,2,2）;

plot（abs（m23）,'r'）;

title（'滤波后信号的频谱'）;

subplot（2,2,3）;

plot（s）;

title（'滤波前信号的波形'）;

subplot（2,2,4）;

plot（z23）;

title（'滤波后的信号波形'）;

图14窗函数法高通滤波