数字信号处理.docx
《数字信号处理.docx》由会员分享，可在线阅读，更多相关《数字信号处理.docx（22页珍藏版）》请在冰豆网上搜索。
数字信号处理
*****************
实践教学
*****************
兰州理工大学
计算机与通信学院
2013年春季学期
《信号处理》课程设计
题目：
基于MATLAB的IIR滤波器的语音信号去噪
专业班级：
姓名：
学号：
指导教师：
成绩：
摘要
IIR低通滤波器是数字信号处理中很重要的一部分。
本次课程设计主要是录制一段语音信号对其进行加噪处理，然后利用IIR低通滤波器对加有随机噪声的语音信号进行滤波处理及时频谱分析，画出滤波之后的频谱图与时域波形，并对信号滤波处理前后进行分析比较，分析信号的变化。
关键词：
滤波去噪；MATLAB；滤波器
目录
第一章语音信号去噪的基本原理1
1.1采样定理1
1.1.1时域采样定理1
1.1.2频域采样定理1
1.1.3采样频率1
1.2.2双线性变换法2
1.2.3数字滤波器的设计步骤4
第二章语音信号去噪实现流程图5
第三章语音信号去噪的仿真及结果分析6
3.1原始语音信号的录制6
3.2语音信号的时频域分析7
3.3加随机噪声后的时频域分析7
3.4IIR数字低通滤波器的仿真10
3.5滤波前后的时频域比较10
总结11
参考文献12
附录13
致谢17
第一章语音信号去噪的基本原理
1.1采样定理
1.1.1时域采样定理
频带为F的连续信号f（t）可用一系列离散的采样值f（t1）,f（t1±Δt），f（t1±2Δt），...来表示,只要这些采样点的时间间隔Δt≤1/2F，便可根据各采样值完全恢复原来的信号f（t）。
这是时域采样定理的一种表述方式。
时域采样定理的另一种表述方式是：
当时间信号函数f（t）的最高频率分量为fM时,f（t）的值可由一系列采样间隔小于或等于1/2fM的采样值来确定,即采样点的重复频率f≥2fm。
图为模拟信号和采样样本的示意图。
时域采样定理是采样误差理论、随机变量采样理论和多变量采样理论的基础。
1.1.2频域采样定理
对于时间上受限制的连续信号f（t）（即当│t│>T时,f（t）=0,这里T=T2-T1是信号的持续时间），若其频谱为F（ω）,则可在频域上用一系列离散的采样值来表示,只要这些采样点的频率间隔ω≦π/tm。
1.1.3采样频率
采样频率，也称为采样速度或者采样率，单位为赫兹（Hz）。
定义了每秒从连续信号中提取并组成离散信号的采样个数，常用的表示符号是fs。
采样时间，采样频率的倒数是采样周期，是采样之间的时间间隔。
通俗的讲采样频率是指计算机每秒钟采集多少个声音样本，是描述声音文件的音质、音调，衡量声卡、声音文件的质量标准。
通俗的讲采样频率是指计算机每秒钟采集多少个声音样本，是描述声音文件的音质、音调，衡量声卡、声音文件的质量标准。
采样频率越高，即采样的间隔时间越短，则在单位时间内计算机得到的声音样本数据就越多，对声音波形的表示也越精确。
根据采样定理，只有采样频率高于声音信号最高频率的两倍时，才能把数字信号表示的声音还原成为原来的声音。
这就是说采样频率是衡量声卡采集、记录和还原声音文件的质量标准。
采样位数和采样率对于音频接口来说是最为重要的两个指标，也是选择音频接口的两个重要标准。
无论采样频率如何，理论上来说采样的位数决定了音频数据最大的力度范围。
每增加一个采样位数相当于力度范围增加了6dB。
采样位数越多则捕捉到的信号越精确。
1.2数字滤波器的设计原理
数字滤波器[1]可以通过数值运算实现滤波，所以数字滤波器处理精度高、稳定、体积小、重量轻、灵活不存在阻抗匹配问题，可以实现模拟滤波器无法实现的特殊功能。
数字滤波器种类很多，根据其实现的网络结构或者其冲激响应函数的时域特性，可分为两种，即有限冲激响应（FIR，FiniteImpulseResponse）滤波器和无限冲激响应（IIR，InfiniteImpulseResponse）滤波器。
IIR滤波器采用递归型结构，即结构上带有反馈环路。
IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。
同时，IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，如巴特沃斯、契比雪夫和椭圆滤波器等，有现成的设计数据或图表可查，在设计一个IIR数字滤波器时，根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。
1.2.1数字滤波器设计基本思想
一个数字滤波器可用它的系统函数H（z）来描述，如式（1-1）
（1-1）
或者用一个N阶差分方程来描述，如式（1-2）
（1-2）
本实验采用的IIR滤波器采用递归型结构，即结构上带有反馈环路。
IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、级联型、并联型四种结构形式，都具有反馈回路。
IIR滤波器设计方法有间接法和直接法，间接法是借助于模拟滤波器的设计方法进行的。
IIR数字滤波器设计方法也有很多种,但它们可归纳为两类：
一类是模拟—数字转换法，其设计步骤是：
先设计过渡模拟滤波器得到系统函数H（s），然后将H（s）按某种方法转换成数字滤波器的系统函数H（z）。
这种方法简单易行，方便准确，但它只能用来设计低通、高通、带通、带阻等选频滤波器；另一类是直接设计法或称之为计算机辅助设计法，这是一种最优化设计法，其适合于设计复杂的数字滤波器。
1.2.2双线性变换法
脉冲响应不变法[1]的主要缺点是频谱交叠产生的混淆，这是从S平面到Z平面的标准变换z＝esT的多值对应关系导致的,为了克服这一缺点，设想变换分为两步：
第一步：
将整个S平面压缩到S1平面的一条横带里；
第二步：
通过标准变换关系将此横带变换到整个Z平面上去。
如图1-1所示
由此建立S平面与Z平面一一对应的单值关系，消除多值性，也就消除了混淆现象。
图1-1双线性变换法的复平面的映射
为了将s平面的jΩ轴压缩到s1平面jΩ轴上的-
一段上，可通过以下的正切变换实现：
这里C是待定常数，下面会讲到用不同的方法确定C，可使模拟滤波器的频率特性与数字滤波器的频率特性在不同频率点有对应关系。
经过这样的频率变换，当Ω1由
时Ω由
即映射了整个jΩ轴。
将这一关系解析延拓至整个s平面，则得到s平面
平面的映射关系如式（1-3）：
（1-3）
再将s1平面通过标准变换关系映射到z平面，即令
通常取C=2/T
最后得S平面与Z平面的单值映射关系如（1-4）：
（1-4）
双线性换法的主要优点是S平面与Z平面一一单值对应，S平面的虚轴（整个jΩ）对应于Z平面单位圆的一周，S平面的Ω=0处对应于Z平面的ω=0处，
对应
即数字滤波器的频率响应终止于折迭频率处，所以双线性变换不存在混迭效应。
1.2.3数字滤波器的设计步骤
（1）确定所需类型数字滤波器的技术指标：
通带边界频率Fp、通带最大衰减As，阻带截止频率Fc、阻带最小衰减Ap。
（2）将所需类型数字滤波器的边界频率转换成相应的模拟滤波器的边界频率，转换公式为Ω=2/Ttan（0.5ω）。
（3）将相应类型的数字滤波器技术指标转换成模拟低通滤波器技术指标。
（4）设计模拟滤波器。
（5）通过频率变换将模拟低通转换成相应类型的过渡模拟滤波器。
（6）采用双线性变换法将相应类型的过渡模拟滤波器转换成所需类型的数字滤波器。
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滤波器时都是默认的双线性变换法，所以在设计滤波器时只需要代入相应的实现函数即可。
第二章语音信号去噪实现流程图
本次课程设计主要任务是对语音信号的简单处理。
运用数字信号学基本原理实现语音信号的处理，在matlab7.1环境下综合运用录音，傅里叶变换、滤波等程序来进行语音信号处理。
本课题的研究基本步骤如下：
1.在MATLAB平台下进行语音信号的录制。
2.绘制原始语音信号的时域图和频谱图并回放原始语音信号。
3.利用MATLAB编程加入一段随机噪音，设计滤波器去噪。
4.对比滤波前后的时域频域图
语音信号滤波去噪——使用双线性变换法设计的巴特沃斯滤波器的设计流程如图2-1所示：
图2-1双线性变换法巴特沃斯滤波器对语音信号去噪流程
第三章语音信号去噪的仿真及结果分析
3.1原始语音信号的录制
1.基本步骤
（1）设置采样频率Fs=11025；
（2）利用函数“wavrecord（15*Fs,Fs,'int16'）”命令进行录音；
（3）利用函数wavplay（）重新播放，看看是否满足所要的效果；
（4）利用“wavwrite”进行保存，应记住保存路径，以便后面进行调用时能快速找到。
需要指出的是利用MATLAB录的音已经是通过采样得到的数字信号
2.Wavread函数[3]几种调用格式。
（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中。
3.2语音信号的时频域分析
画出语音信号的时域波形，再对语音信号进行频谱分析。
MATLAB提供了快速傅里叶变换算法FFT计算DFT的函数fft,其调用格式如下：
Xk=fft（xn,N）
参数xn为被变换的时域序列向量，N是DFT变换区间长度，当N大于xn的长度时，fft函数自动在xn后面补零。
，当N小于xn的长度时，fft函数计算xn的前N个元素，忽略其后面的元素。
在本次设计中，我们利用fft对语音信号进行快速傅里叶变换，就可以得到信号的频谱特性。
图3-3原始信号的时域频谱图
3.3加随机噪声后的时频域分析
对本次试验中利用MATLAB中的随机函数（rand或randn）[4]产生随机噪声加入到语音信号中，模仿语音信号被污染，并对其频谱分析。
Matlab中函数randn：
产生正态分布的随机数或矩阵的函数产生均值为0，方差σ^2=1，标准差σ=1的正态分布的随机数或矩阵的函数。
用法：
（1）Y=randn（n）：
返回一个n*n的随机项的矩阵。
如果n不是个数量，将返回错误信息。
（2）Y=randn（m,n）或Y=randn（[mn]）：
返回一个m*n的随机项矩阵。
（3）Y=randn（m,n,p,...）或Y=randn（[mnp...]）：
产生随机数组。
（4）Y=randn（size（A））：
返回一个和A有同样维数大小的随机数组。
（5）randn：
返回一个每次都变化的数量。
仿真图如下：
图3-4原始信号时域频域图
图3-5原始信号相位幅度图
图3-5加随机噪声后语音信号的时域频谱图
仿真结果分析:
通过对两张图片的对比，很明显可以看加噪后的语音信号时域波形比原始语音信号浑浊了许多，在时间轴上可以明显看出0—0.5S的幅值增大了；通过对原始语音信号的频谱图与加噪后的语音信号频谱图的对比，也可以看出在频率5000Hz以后的频率幅值发生了明显的增加。
再通过对原始语音信号的回放效果与加噪后的语音信号回放的效果的对比，人耳可以明显辨别出两种语音信号不一样了，加噪后的语音信号在听觉上比原始语音信号要浑浊很多，而且还有吱吱嘎嘎的混杂音。
3.4IIR数字低通滤波器的仿真
低通滤波器的性能指标如下[5]：
通带截止频率fp=1000Hz，阻带截止频率fc=1200Hz，阻带最小衰减As=50db,通带最大衰减Ap=1dB
由此设计出的数字低通滤波器如下：
图4-2低通滤波器
3.5滤波前后的时频域比较
滤波前后的语音信号频谱和时域图有很大区别，若听语音信号，则发现加噪后的语音信号出现混杂音，频谱图中高频成分增多，但经过低通滤波器滤波后，语音信号和原始语音信号差别不大。
仿真图如图4-6所示
图3-6滤波前后比较
仿真结果分析:
语音信号的频率300Hz到4000Hz，它属于低频信号，采用低通滤波器可有效滤除噪声。
信号的失真是最小的。
由图中滤波前后波形比较可看出，经过滤波后的波形比原波形的振幅有所减小，去除了很多由于噪声所产生的干扰；
从滤波前后的频谱比较可以看出经过滤波后除了原本的声音外，中间由于噪声产生的频谱波形已经滤除；经过MATLAB仿真，听滤波前后的声音，可以听出有明显的滤波效果。
因此利用双线性变换法设计的巴特沃斯滤波器已经达到了设计的要求。
总结
课题“基于MATLAB的IIR滤波器语音信号去噪”是将数字信号处理技术应用于对语音及加噪处理。
作为存储于计算机中的语音信号，其本身就是离散化了的向量，只需将这些离散的量提取出来，就可以对其进行处理了。
这一过程的实现，用到了处理数字信号的强有力工具MATLAB。
通过MATLAB里几个命令函数的调用，很轻易的在实际语音与数字信号的理论之间搭了一座桥。
本次课程设计所做的工作如下：
第一步是语音信号的录制，在MATLAB的环境下录制20s的“.wav”格式的语音信号。
并将它保存在目录下，以备后面信号分析的处理。
第二步是分析原始语音信号的时域和频域，频域分析利用傅里叶变换对于20s的语音信号，它的频域范围在4000Hz以内。
第三步是用双线性变换法设计一个低通滤波器，先设计出模拟低通滤波器，在经过双线性变换法转换为符合要求的数字低通滤波器。
其低通滤波器性能指标，fp=1000Hz，fc=1200Hz，As=50db,Ap=1dB。
第四步是语音信号进行加噪处理，在原来语音信号的基础上加一个随机噪声，对两张图片的对比，很明显可以看加噪后的语音信号时域波形比原始语音信号浑浊了许多，在时间轴上可以明显看出0—0.5S的幅值增大了；通过对原始语音信号的频谱图与加噪后的语音信号频谱图的对比，也可以看出在频率5000Hz以后的频率幅值发生了明显的增加。
再通过对原始语音信号的回放效果与加噪后的语音信号回放的效果的对比，人耳可以明显辨别出两种语音信号不一样了，加噪后的语音信号在听觉上比原始语音信号要浑浊很多，而且还有吱吱嘎嘎的混杂音。
第五步对加噪语音信号进行滤波前后对比，由图中滤波前后波形比较可看出，经过滤波后的波形比原波形的振幅有所减小，去除了很多由于噪声所产生的干扰；从滤波前后的频谱比较可以看出经过滤波后除了原本的声音外，中间由于噪声产生的频谱波形已经滤除；经过MATLAB仿真，听滤波前后的声音，可以听出有明显的滤波效果。
因此利用双线性变换法设计的巴特沃斯滤波器已经达到了设计的要求。
参考文献
[1]高西全，丁玉美．数字信号处理．第3版．北京：
西安电子科技大学出版社，2008151:
187
[2]刘泉，阙大顺．数字信号处理原理与实现．北京：
电子工业出版社，2005123:
126
[3]张磊，毕靖，郭莲英．MATLAB实用教程．北京：
人民邮电出版社，2008245:
237
[4]张威．MATLAB基础与编程入门．西安：
西安电子科技大学出版社，200656:
68
[5]周利清，苏菲．数字信号处理基础．北京：
北京邮电大学出版社，2005198-200
附录
1.语音信号的录制
Fs=11025;%采样频率为11025
y=wavrecord（15*Fs,Fs,'int16'）;%进行15秒-16bit的录音
wavplay（y,Fs）;%重新播放录音
wavwrite（y,Fs,'xj.wav'）;保存该录音，且名字为"xj",格式为wav
wavwrite（y,'xj.wav'）;
2.原始语音信号的时频域分析
[y,fs,nbits]=wavread（'xj.wav'）;
sound（y,fs,nbits）;%回放语音信号
N=length（y）;%求出语音信号的长度
Y=fft（y,N）;%对语音信号进行傅里叶变换
subplot（2,1,1）;
plot（y）;title（'原始信号时域图'）;
subplot（2,1,2）;
plot（abs（Y））;
title（'原始信号频域图'）
3.加随机噪声后的语音信号的时频域分析
fs=11025;%语音信号采样频率为11025
[x,fs,bits]=wavread（'C:
\Users\Administrator\Desktop\function\xj.wav'）;
sound（x,fs,bits）;%播放语音信号
X=fft（x,4096）;
magX=abs（X）;
angX=angle（X）;
y1=fft（x,1024）;%对信号做1024点FFT变换
f=fs*（0:
511）/1024;
figure
（1）
subplot（211）;plot（magX）;title（'原始信号幅值'）;
gridon;
subplot（212）;plot（angX）;title（'原始信号相位'）;
gridon;
figure
（2）
subplot（211）;
plot（x）;%绘制原始语音信号的时域波形图
title（'原始语音信号时域波形图'）;
xlabel（'timen'）;
ylabel（'fuzhin'）;
gridon;
Subplot（212）;%绘制原始语音信号的频率响应图
plot（f,abs（y1（1:
512）））;
title（'原始语音信号频谱图'）
xlabel（'Hz'）;
ylabel（'fudu'）;
gridon;
%添加随机噪声及添加噪声后的时域图和频谱图
noise_mu=0;
noise_var=0.05;
x0=randn（size（x））.*sqrt（noise_var）+noise_mu;
x1=x+x0;
ts=1/fs;%绘制在原始信号上加随机噪声的信号图
ta=（length（x）-1）/fs;
t=0:
ts:
ta;
figure（3）;
subplot（211）;
plot（t,x1）;
title（'加随机噪声后语音信号时域图'）;
xlabel（'t'）;
ylabel（'x1'）;
gridon;
y2=fft（x1,1024）;%对信号做1024点FFT变换
f=fs*（0:
511）/1024;
Subplot（212）;%绘制原始语音信号的频率响应图
plot（f,abs（y2（1:
512）））;
title（'加随机噪声后的语音信号频谱图'）
xlabel（'Hz'）;
ylabel（'fudu'）;
gridon;
4.滤波前后的时频域对比
[y,fs,bits]=wavread（'C:
\Users\Administrator\Desktop\function\xj.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,'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）;
plot（y）;
title（'滤波前信号波形'）;
grid;
subplot（2,2,4）;
plot（z11）;
title（'滤波后的信号波形'）;
grid;
wavwrite（y,'C:
\Users\Administrator\Desktop\function\x.wav'）;
致谢
这次课程设计，给我留下了很深的印象。
虽然只是时间短暂，但在这期间，却让我受益匪浅。
通过这次课程设计，使我对语音信号有了全面的认识，对数字信号处理的知识又有了深刻的理解，在之前数字信号与处理的学习以及完成课后的作业的过程中，已经使用过MATLAB，对其有了一些基础的了解和认识。
本次经过一周的课程设计，让我们更熟悉了数字信号处理的基本知识和MATLAB的语言，让我们把课上的理论知识运用到实际中去，让我们更近一步地巩固了课堂上所学的理论知识，并能很好地理解与掌握数字信号处理中的基本概念、基本原理、基本分析方法。
在此论文撰写过程中，要特别感谢我的指导老师的指导与督促，同时感谢她的谅解与包容。
没有老师的帮助也就没有今天的这篇文章。
求学历程是艰苦的，但又是快乐的。
总的来说，通过这次设计使我了解了MATLAB的使用方法，学会分析滤波器的优劣和性能，提高了分析和动手实践能力。
同时我相信，进一步加强对MATLAB的学习与研究对我今后的学习将会起到很大的帮助！
本文参考了大量的文献资料，在此，向各学术界的前辈们致敬！
升级会员 