matlab课程设计报告基于matlab有噪声语音信号处理大学论文.doc

Matlab课程设计报告

题目:基于MATLAB有噪声语音信号处理

系 （院）： 计算机与信息工程学院
专 业： 通信工程
班 级： 10623102

指导教师：
学年学期: 2011 ~ 2012 学年 第 2 学期

简介：
我们通信工程专业在实践中经常碰到需要对已接收信号进行处理的情况，而滤波器设计在数字信号处理中占有极其重要的地位。本课题基于MATLAB有噪音语音信号处理的设计与实现，综合运用数字信号处理的理论知识对加噪语音信号进行时域、频域分析和滤波。通过理论推导得出相应结论，再利用MATLAB作为编程工具进行计算机实现。在设计实现的过程中，我们使用双线性变换法设计IIR数字滤波器，对模拟加噪语音信号进行低通滤波、高通滤波及带通滤波，并利用MATLAB作为辅助工具完成设计中的计算与图形的绘制。
1 绪论：
数字信号处理是利用计算机或专用处理设备，以数值计算的方法对信号进行采集、抽样、变换、综合、估值与识别等加工处理，借以达到提取信息和便于应用的目的。数字滤波器, 是数字信号处理中及其重要的一部分。本课题采用IIR滤波器对加噪声音信号进行处理。
IIR滤波器采用递归型结构，即结构上带有反馈环路。IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。同时，IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，如巴特沃斯、契比雪夫和椭圆滤波器等，有现成的设计数据或图表可查，在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。

2.原始语音信号采集与处理
2.1语音信号的采集
由于MATLAB只识别格式为.wav的声音文件，我们利用PC机上的声卡和WINDOWS操作系统进行数字信号的采集。启动录音机进行录音，以文件名“Orisound”保存入原程序所属的文件夹中。可以看到,文件存储器的后缀默认为.wav ,这是WINDOWS操作系统规定的声音文件存的标准。
程序流程图：
开始
用wavread读取语音信号，并进行采样，建立数据文件，并用plot画出数据文件时域波形图和频谱图
利用MATLAB中的随机函数产生噪声加入到语音信号中，模仿语音信号被污染，并对其进行频谱分析。

回放语音信号得出所设计滤波器在语音处理中的优劣并设计一个信号处理系统界面。
结束
运用数字信号处理理论设计IIR数字滤波器， 并对被噪声污染的语音信号进行滤波，分析滤波后信号的时域和频域。特征。

2.2语音信号的时频分析
利用MATLAB中的“wavread”命令来读入（采集）语音信号，将它赋值给某一向量。再对其进行采样，记住采样频率和采样点数。
对语音信号Orisound.wav进行采样其程序如下：
[y,fs,nbits]=wavread (‘OriSound’); %把语音信号加载入MATLAB仿真软件平台中
画出语音信号的时域波形，再对语音信号进行频谱分析。在本次设计中，我们利用fft函数对语音信号进行快速傅里叶变换，就可以得到信号的频谱特性。程序如下：
[y,fs,nbits]=wavread(‘Orisound’); %语音信号的采集
sound(y,fs,nbits); %语音信号的播放
n=length(y) ;
Y=fft(y,n); %快速傅里叶变换
figure;
subplot(2,1,1);
plot(y);
title(‘原始信号波形’,’fontweight’,’bold’);
axis([ 78000 80000 -1 1]);
grid;
subplot(2,1,2);
plot(abs(Y));
title(‘原始信号频谱’,’fontweight’,’bold’);
axis([ 0 150000 0 4000]);
grid;
程序结果如下图：

2.3语音信号加噪与频谱分析
利用MATLAB中的随机函数Randn（m,n）产生噪声加入到语音信号中，模仿语音信号被污染，并对其频谱分析。主要程序如下：
[y,fs,nbits]=wavread('OriSound');
sound(y,fs,nbits);
n = length (y) ;
Noise=0.2*randn(n,2);
s=y+Noise;
sound(s);
figure;
subplot(2,1,1);
plot(s);
title('加噪语音信号的时域波形','fontweight','bold');
axis([ 78000 80000 -1 1]);
grid;
S=fft(s);
subplot(2,1,2);
plot(abs(S));
title('加噪语音信号的频域波形','fontweight','bold');
axis([ 0 150000 0 4000]);
grid;

程序结果如下图：

3设计数字滤波器
3.1 数字滤波器设计的基本思路
数字滤波器的实现有两个关键步骤：一个从数字域到模拟域间的变换，这个变换实现了数字滤波器技术指标到模拟滤波器技术指标的转换，同样也实现了模拟滤波器系统函数到数字滤波器系统函数的转换；另一个是从模拟滤波器技术指标到满足该指标的模拟滤波器的设计。
3.2 模拟滤波器概述
用模拟—数字变换法设计IIR数字滤波器，首先必须设计一个模拟滤波器，它有许多不同的类型，主要有以下两种类型：
（1）、巴特沃思（Botterworth简写BW）滤波器。BW滤波器是根据幅频特性在通带内具有最拼图特性而定义的滤波器，对一个N阶低通滤波器来说，所谓最平坦特性就是模拟函数的前（2N-1）阶导数在处都为零。BW滤波器的另一个特性是在通带和阻带内的幅频特性始终是频率的单调下降函数，且其模拟函数随阶次N 的增大而更接近于理想低通滤波器。
（2）、切比雪夫（Chbyshev简写为CB）滤波器。CB低通滤波器的模拟函数由切比雪夫多项式定义，且在通带内的幅频响应是波动的，在阻带则单调变化。
3.3设计IIR滤波器
目前IIR数字滤波器设计的最通用的方法是借助于模拟滤波器的设计方法。模拟滤波器设计已经有了一套相当成熟的方法，它不但有完整的设计公式，而且还有较为完整的图表供查询，因此，充分利用这些已有的资源将会给数字滤波器的设计带来很大方便。IIR数字滤波器的设计步骤是：
（1）、按一定规则将给出的数字滤波器的技术指标转换为模拟低通滤波器的技术指标；
（2）、根据转换后的技术指标设计模拟低通滤波器；
（3）、再按一定规则将G(s)转换成H(z)。
若设计的数字滤波器是低通的，那么上述设计工作可以结束，若所设计的是高通，带通或带阻滤波器，那么还有步骤（4）：
（4）、将高通、带通、或带阻数字l不去的技术指标先转化为低通模拟滤波器的技术指标，然后按照上述步骤（2）设计出低通，再将转换为所需的H(z)。
利用模拟滤波器设计IIR数字低通滤波器的步骤：
（1）确定数字低通滤波器的技术指标：通带边界频率、通带最大衰减，阻带截止频率、阻带最小衰减。
（2）将数字低通滤波器的技术指标转换成相应的模拟低通滤波器的技术指标。
（3）按照模拟低通滤波器的技术指标设计及过渡模拟低通滤波器。
（4）用双线性变换法，模拟滤波器系统函数转换成数字低通滤波器系统函数。

程序如下：
Ft=8000;
Fp=1000;
Fs=1200;
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’);
[num11,den11]=bilinear(b11,a11,0.5);
[h,w]=freqz(num11,den11);
figure;
plot(w*8000*0.5/pi,abs(h));
legend(‘IIR低通滤波器’,’Location’,’NorthWest’);
grid;
程序结果如下图：

3.2验证所设计的滤波器：
为了验证滤波器的可使用性，我们用常用的sin函数来进行验证。其具体程序及运行结果如下：
t=[0:1/1023:1];
s=sin(2*pi*t);
N=length(s);
y=s+0.5*rand(1,N);
subplot(2,1,1);
plot(y);
title('加噪语音信号的时域波形','fontweight','bold');
S=fft(y);
subplot(2,1,2);
plot(abs(S));
title('加噪语音信号的频域波形','fontweight','bold');
Ft=8000;
Fp=1000;
Fs=1200;
wp=2*pi*Fp/Ft;
ws=2*pi*Fs/Ft;
[n11,wn11]=buttord(wp,ws,1,50,'s'); %求低通滤波器的阶数和截止频率
[b11,a11]=butter(n11,wn11,'s'); %求S域的频率响应的参数
[num11,den11]=bilinear(b11,a11,0.5); %利用双线性变换实现频率响应S域到Z域的变换
z11=filter(num11,den11,s);
sound(z11);
m11=fft(z11); %求滤波后的信号
figure;
subplot(2,2,1);
plot(abs(S),'g');
title('滤波前信号的频谱','fontweight','bold');
grid;
subplot(2,2,2);
plot(abs(m11),'r');
title('滤波后信号的频谱','fontweight','bold');
grid;
subplot(2,2,3);
plot(y);
title('滤波前信号的波形','fontweight','bold');
grid;
subplot(2,2,4);
plot(z11);
title('滤波后的信号波形','fontweight','bold');
grid;
程序结果如下图：

由所得结果可知，所设计的滤波器符合要求。
4 滤波
用设计好的IIR低通滤波器对加噪的语音信号进行滤波，程序如下：
[y,fs,nbits]=wavread (‘OriSound’); %IIR低通
n = length (y) ; %求出语音信号的长度
Noise=0.2*randn(n,2); %随机函数产生噪声
s=y+Noise; %语音信号加入噪声
S=fft(s);
Ft=8000;
Fp=1000;
Fs=1200;
wp=2*pi*Fp/Ft;
ws=2*pi*Fs/Ft;
[n11,wn11]=buttord(wp,ws,1,50,’s’);%求低通滤波器的阶数和截止频率
[b11,a11]=butter(n11,wn11,’s’); %求S域的频率响应的参数
[num11,den11]=bilinear(b11,a11,0.5); %利用双线性变换实现频率响应S域到Z域的变换
z11=filter(num11,den11,s);
sound(z11);
m11=fft(z11); %求滤波后的信号
figure;
subplot(2,2,1);
plot(abs(S),’g’);
title(‘滤波前信号的频谱’,’fontweight’,’bold’);
axis([ 0 150000 0 4000]);
grid;
subplot(2,2,2);
plot(abs(m11),’r’);
title(‘滤波后信号的频谱’,’fontweight’,’bold’);
axis([ 0 150000 0 4000]);
grid;
subplot(2,2,3);
plot(s);
title(‘滤波前信号的波形’,’fontweight’,’bold’);
axis([95000 100000 -1 1]);
grid;
subplot(2,2,4);
plot(z11);
title(‘滤波后的信号波形’,’fontweight’,’bold’);
axis([95000 100000 -1 1]);
grid;
程序结果如下图：

经过以上的加噪处理后，可在Matlab中用函数sound对声音进行回放。其调用格式：sound(y,Fs)，sound(y)和sound(y,Fs,bits)。可以察觉滤波前后的声音有明显的变化。

5.创建GUI界面

总结
本次设计通过选择一个语音信号作为分析的对象，对其进行频谱分析；利用MATLAB中的随机函数产生噪声加入到语音信号中，模仿被噪声干扰的语音信号，并对其进行频谱分析；运用数字信号处理理论设计IIR数字滤波器， 并对被噪声污染的语音信号进行滤波，分析滤波后信号的时域和频域特征，回放语音信号。最后，设计一个信号处理系统界面返回。
我们四个人共同完成这个课程设计，分工是很重要的，首先我们确定了大致的方向，我们要做低通，带通以及高通三种滤波器，此外还有界面的设计，程序的整合，以及课程报告的编写等方面的内容，我们的具体分工如下：
声音信号的读入以及ppt的制作：翁淦泉1062310209
低通，带通，高通滤波器的设计及相关程序的编写：冯竹青1062310202，贺静文1062301203
最后整合及界面的设计：王雅青1062310208
课程报告：贺静文

参考文献
1、张志涌 杨祖樱 MATLAB教程R2011a 北京航空航天大学出版社
2、高西全，丁玉美.数字信号处理（第三版）[M].西安：西安电子科技大学出版社，2008.8
3、蒋瑞艳 振动，噪声处理系统[J]. 大连理工大学学院学报(自然科学版) 2002
4、马超  高世伦 基于MATLAB的噪声信号采集与分析系统研究[J]. 华中科技大学能源与动力工程学院院报 2004.5
5、李正周.MATLAB数字信号处理与应用[M].北京:清华大学出版社，2008
6、高萍,祖静.基于MATLAB小波去噪技术浅析[J].科技信息(学术版),2006(11)1-3.
7、刘智 基于MATLAB的机械噪声分析与处理[J]. 华中江师范学院学报(自然科学版) 2010（4）

附录
附录（I） 设计IIR数字滤波器

%=========================IIR低通滤波器=======================
Ft=8000;
Fp=1000;
Fs=1200;
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’);
[num11,den11]=bilinear(b11,a11,0.5);
[h,w]=freqz(num11,den11);
figure;
plot(w*8000*0.5/pi,abs(h));
legend(‘IIR低通滤波器’,’Location’,’NorthWest’);
grid;
程序结果如下图：

%======================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’);
[num,den]=lp2hp(b13,a13,wn13);
[num13,den13]=bilinear(num,den,0.5);
[h,w]=freqz(num13,den13);
figure;
plot(w*21000*0.5/pi,abs(h));
legend(‘IIR高通滤波器’,’Location’,’NorthWest’);
axis([0 11000 0 1.5]);
grid;
程序结果如下图：

%======================IIR带通滤波器==========================
Fp1=1200;
Fp2=3000;
Fs1=1000;
Fs2=3200;
Ft=8000;
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’);
[num2,den2]=lp2bp(b12,a12,sqrt(wp1*wp2),bw);
[num12,den12]=bilinear(num2,den2,0.5);
[h,w]=freqz(num12,den12);
figure;
plot(w*8000*0.5/pi,abs(h));
axis([0 4500 0 1.5]);
legend(‘IIR带通滤波器’,’Location’,’NorthWest’);
grid;

程序结果如下图：

附录（II）比较滤波前后语音信号的波形及频谱
% ======================双线性变换法=======================
%*************************低通滤波器************************
[y,fs,nbits]=wavread (‘OriSound’); %IIR低通
n = length (y) ; %求出语音信号的长度
Noise=0.2*randn(n,2); %随机函数产生噪声
s=y+Noise; %语音信号加入噪声
S=fft(s);
Ft=8000;
Fp=1000;
Fs=1200;
wp=2*pi*Fp/Ft;
ws=2*pi*Fs/Ft;
[n11,wn11]=buttord(wp,ws,1,50,’s’);%求低通滤波器的阶数和截止频率
[b11,a11]=butter(n11,wn11,’s’); %求S域的频率响应的参数
[num11,den11]=bilinear(b11,a11,0.5); %利用双线性变换实现频率响应S域到Z域的变换
z11=filter(num11,den11,s);
sound(z11);
m11=fft(z11); %求滤波后的信号
figure;
subplot(2,2,1);
plot(abs(S),’g’);
title(‘滤波前信号的频谱’,’fontweight’,’bold’);
axis([ 0 150000 0 4000]);
grid;
subplot(2,2,2);
plot(abs(m11),’r’);
title(‘滤波后信号的频谱’,’fontweight’,’bold’);
axis([ 0 150000 0 4000]);
grid;
subplot(2,2,3);
plot(s);
title(‘滤波前信号的波形’,’fontweight’,’bold’);
axis([95000 100000 -1 1]);
grid;
subplot(2,2,4);
plot(z11);
title(‘滤波后的信号波形’,’fontweight’,’bold’);
axis([95000 100000 -1 1]);
grid;
程序结果如下图：

附II-1 双线性低通滤波器比较

%**********************高通滤波器*****************************
[y,fs,nbits]=wavread (‘OriSound’); %IIR高通
n = length (y) ; %求出语音信号的长度
Noise=0.2*randn(n,2); %随机函数产生噪声
s=y+Noise; %语音信号加入噪声
S=fft(s); %傅里叶变换
Fp1=1200;
Fs1=1000;
Ft=8000;
wp1=tan(pi*Fp1/Ft);
ws1=tan(pi*Fs1/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域转换
z13=filter(num13,den13,s);
sound(z13);
m13=fft(z13); %求滤波后的信号
figure;
subplot(2,2,1);
plot(abs(S),’g’);
title(‘滤波前信号的频谱’,’fontweight’,’bold’);
axis([0 150000 0 4000]);
grid;
subplot(2,2,2);
plot(abs(m13),’r’);
title(‘滤波后信号的频谱’,’fontweight’,’bold’);
axis([0 150000 0 4000]);
grid;
subplot(2,2,3);
plot(s);
title(‘滤波前信号的波形’,’fontweight’,’bold’);
axis([95000 100000 -1 1]);
grid;
subplot(2,2,4);
plot(z13);
title(‘滤波后的信号波形’,’fontweight’,’bold’);
axis([95000 100000 -1 1]);
grid;
程序结果如下图：

附II-2 双线性高通滤波器比较

%**********************带通滤波器*****************************
[y,fs,nbits]=wavread (‘OriSound’); %IIR带通
n = length (y) ; %求出语音信号的长度
Noise=0.2*randn(n,2); %随机函数产生噪声
s=y+Noise; %语音信号加入噪声
S=fft(s); %傅里叶变换
Ft=8000;
Fp=1000;
Fs=1200;
wp=2*Fp/Ft;
ws=2*Fs/Ft;
rp=1;
rs=50;
p=1-10.^(-rp/20); %通带阻带波纹
q=10.^(-rs/20);
fpts=[wp ws];
mag=[1 0];
dev=[p q];
[n21,wn21,beta,ftype]=kaiserord(fpts,mag,dev);%由kaiserord求滤波器的阶数和截止频率
b21=fir1(n21,wn21,kaiser(n21+1,beta)); %由fir1设计滤波器
z21=fftfilt(b21,s);
sound(z21);
m21=fft(z21); %求滤波后的信号
figure(4);
subplot(2,2,1);
plot(abs(S),’g’);
title(‘滤波前信号的频谱’,’fontweight’,’bold’);
axis([0 150000 0 4000]);
grid;
subplot(2,2,2);
plot(abs(m21),’r’);
title(‘滤波后信号的频谱’,’fontweight’,’bold’);
axis([0 150000 0 4000]);
grid;
subplot(2,2,3);
plot(s);
title(‘滤波前信号的波形’,’fontweight’,’bold’);
axis([95000 100000 -1 1]);
grid;
subplot(2,2,4);
plot(z21);
title(‘滤波后的信号波形’,’fontweight’,’bold’);
axis([95000 100000 -1 1]);
grid;
程序结果如下图：

附II-3 双线性带通滤波器比较