陈龙语音信号去噪处理

《数字信号处理》课程设计

题目：基于MATLAB的语音信号去噪处理

学院：电子与信息工程学院

班级： 2014级电信3班

姓名：陈龙

学号： 201407014330

教师：冯地耘

成绩：__________________

2016年12月26日

1、课程设计题目

基于MATLAB的语音信号去噪处理。

2、课程设计目的

综合运用数字信号处理的理论知识进行频谱分析和滤波器设计，通过理论推导得出相应的结论，再利用MATLAB做为编程工具进行计算机实现，从而加深对所学知识的理解，建立概念。

3、设计原理

语音采集的原理：先在自己的电脑录制一小段声音，然后把录制的声音放入matlab的work中，在matlab的主页面输入一段函数，用wavread函数调用自己录制的语音信号。利用Matlab中的wavread命令，读入（采集）本段语音信号，将它赋值给某一向量。

加噪的原理：加入一个随机噪音信号用一个公式：y=x+nosie。函数randn：产生正态分布的随机数或者矩阵的函数。产生的均值为0，方差σ^2=1,标准差σ=1的正态分布的随机数或矩阵的函数。

用法：Y=randn(n)，返回一个n*n的随机项矩阵。如果n不是个数量，则将返回错误信息。

Y=randn(m,n) 或Y=randn([m n]),返回一个m*n的随机项矩阵。

Y=randn(size(A)),返回一个和A有同样维数大小的随机数组。

randn返回一个每次都变化的数量。

去噪的方法

数字信号处理技术经过几十年的发展，在国内外已经取得了很大的成绩。到目前为止，已经比较成熟的去噪方法比较典型的有切比雪夫去噪法、双线性变换去噪法、窗函数去噪法、谱相减去噪法、巴特沃兹去噪法等有名的去噪方法。

FIR滤波器原理

FIR滤波器通常采用窗函数方法来设计。窗设计的基本思想是，首先选择一个适当的理想选频滤波器（它总是具有一个非因果，无限持续时间脉冲响应），然后加窗它的脉冲响应得到线性相位和因果FIR滤波器。我们用Hd（e^jw）表示理想的选频滤波器，它在通带上具有单位增益和线性相位，在阻带上具有零响应。一个带宽wc

为了从hd(n)得到一个FIR滤波器，必须同时在两边截取hd（n）。而要得到一个因果的线性相位滤波器，它的h(n)长度为N，必须有：

这种操作叫做加窗，h(n)可以看做是hd(n)与窗函数w（n）的乘积：

h(n)=hd(n)w(n)

其中

根据w（n）的不同定义，可以得到不同的窗结构。

在频域中，因果FIR滤波器响应H(e^jw)由Hd(e^jw)和窗响应W(e^jw)的周期卷积得到，即

常用的窗函数有矩形窗、巴特利特（BARTLETT）窗、汉宁（HANNING）窗、海明（HAMMING）窗、布莱克曼（BLACKMAN）窗、凯泽（KAISER）窗等。

FIR：有限脉冲响应滤波器。有限说明其脉冲响应是有限的。与IIR相比，它具有线性相位、容易设计的优点。这也就说明，IIR滤波器具有相位不线性，不容易设计的缺点。而另一方面，IIR却拥有FIR所不具有的缺点，那就是设计同样参数的滤波器，FIR比IIR需要更多的参数。这也就说明，要增加DSP的计算量。DSP需要更多的计算时间，对DSP的实时性有影响。FIR滤波器的设计比较简单，就是要设计一个数字滤波器去逼近一个理想的低通滤波器。通常这个理想的低通滤波器在频域上是一个矩形窗。根据傅里叶变换我们可以知道，此函数在时域上是一个采样函数。通常此函数的表达式为：

sa（n）＝sin（n∩）/n（3-2）

但是这个采样序列是无限的，计算机是无法对它进行计算的。故我们需要对此采样函数进行截断处理。也就是加一个窗函数。就是传说中的加窗。也就是把这个时域采样序列去乘一个窗函数，就把这个无限的时域采样序列截成了有限个序列值。但是加窗后对此采样序列的频域也产生了影响：此时的频域便不在是一个理想的矩形窗，而是成了一个有过渡带，阻带有波动的低通滤波器。通常根据所加的窗函数的不同，对采样信号加窗后，在频域所得的低通滤波器的阻带衰减也不同。通常我们就是根据此阻带衰减去选择一个合适的窗函数。如矩形窗、汉宁窗、汉明窗、BLACKMAN窗、凯撒窗等。

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

（1）根据对阻带衰减及过渡带的指标要求，选择串窗数类型（矩形窗、三角窗、汉宁窗、哈明窗、凯塞窗等），并估计窗口长度N。先按照阻带衰减选择窗函数类型。原则是在保证阻带衰减满足要求的情况下，尽量选择主瓣的窗函数。

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

（3）计算h(n).。

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

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

低通滤波器的性能指标：

fp=1000Hz，fc=1200Hz，As=50db ,Ap=1dB

高通滤波器的性能指标：

fp=3500Hz，fc=4000Hz，As=50dB，Ap=1dB

在Matlab中，可以利用函数fir1设计FIR滤波器，利用Matlab中的函数freqz 画出各步步器的频率响应。

FIR的特点：

有限长单位冲激响应（FIR）滤波器有以下特点：

A位冲激响应h (n)在有限个n值处不为零

B系统函数0处收敛，极点全部在z = 0处

C结构上主要是非递归结构，没有输出到输入的反馈，但有些结构中（例如频率抽样结构）也包含有反馈的递归部分。

设FIR滤波器的单位冲激响应h (n)为一个N点序列，0 ≤ n ≤N —1，则滤波器的系统函数为

H(z)=∑h(n)*z^-n （2-3）

就是说，它有（N—1）阶极点在z = 0处，有（N—1）个零点位于有限z 平面的任何位置。

FIR滤波器基本结构：

FIR滤波器有以下几种基本结构：

横截型

式的系统的差分方程表达式为

y(n)=∑h(m)x(n-m)( 2-4）

很明显，这就是线性移不变系统的卷积和公式，也是x (n)的延时链的横向结构，称为横截型结构或卷积型结构，也可称为直接型结构。将转置定理用于可得到的转置直接型结构。

FIR滤波器的横截型结构

级联型

其中[N/2]表示取N/2的整数部分。若N为偶数，则N—1为奇数，故系数B2K中有一个为零，这是因为，这时有奇数个根，其中复数根成共轭对必为偶数，必然有奇数个实根。画出N为奇数时，FIR滤波器的级联结构，其中每一个二阶因子用图4-11的横型结构。

这种结构的每一节控制一对零点，因而再需要控制传输零点时，可以采用它。但是这种结构所需要的系数B2k（I = 0，1，2，k，= 1，2，．．．，[N/2]）比卷积型的系数h (n)要多，因而所需的乘法次数也比卷积型的要多。

4、课程设计的要求

(1)熟悉离散信号和系统的时域特性。

(2)掌握序列快速傅里叶变换FFT方法。

(3)学会MATLAB的使用，掌握MATLAB的程序设计方法。

(4)利用MATLAB对语音信号进行频谱分析。

(5)掌握MATLAB设计各种数字滤波器的方法和对信号进行滤波的方法。5、课程设计步骤

（1）语音信号的获取：

单击自己的电脑开始菜单，选择所有程序，接着选择附件，再选择娱乐，最后选择录音。自己录入语音信号，然后保存在MATLAB的work文件夹里面，命名为“chenlong”。

利用Matlab中的wavread命令，读入（采集）本段语音信号，将它赋值给某一向量。

(2)原始语音信号的频谱分析及仿真

在Matlab中使用wavread函数，可得出一个采样频率为22500HZ，并且声音是单声道的信号。利用Sound函数，我们可以清晰的听到“毕业设计声音文件”的语音。采集数据并画出波形图。

在Matlab的信号处理工具箱中,函数FFT和IFFT分别用于快速傅立叶变换和逆变换。前面已简单介绍过，这里详细介绍下这些函数。

函数fft用于序列快速傅立叶变换，是常用的信号处理函数。

函数的一种调用格式：y=fft(x)。

其中，x是序列，y是序列的快速傅里叶变换FFT，x可以为一向量,也可以为一矩阵。若x为一向量，则y是x的fft，且和x为相同长度。若x为一矩阵，y则是对矩阵的每一列向量进行的fft。

若序列x的长度是2的幂次方，则函数fft执行高速基－2fft算法；否则fft 将执行一种混合基的离散傅立叶变换算法，计算速度会较慢。

函数fft的另一种调用格式为y=fft(x,N)。

该式中，x，y的意义同前，N为正整数。函数执行序列x的N点快速傅里叶变换(fft)。若x为向量,且长度小于N，则函数会将x补零至长度N。若向量x的长度大于N，则函数将截短x,使之长度为N。若x为矩阵，则按相同方法对x进行相应处理。

由于经函数fft所求得的序列y一般是复序列，所以通常要求其幅值和相位。Matlab提供求复数的幅值和相位函数：abs，angle，这些函数一般和fft同时使用。函数abs(x)用于计算复向量x的幅值，而函数angle(x)则用于计算复向量的相角，用弧度表示。

（3）下面的一段程序是语音信号在Matlab中的简单体现。通过运行该段程序,我们实现了语音信号的读入与打开，并绘出了语音信号的波形时域图和频谱图。

原语音初始信号的时域波形、频谱图程序：

[y,fs,bits]=wavread('chenlong');%语音信号的采集

sound(y,fs)%语音信号的播放

n=length(y)%计算语音信号的长度

p=fft(y,n); %快速傅里叶变换

f=fs*(0:n/2-1)/n;%对应点的频率

figure(1)

subplot(2,1,1);

plot(y);

title('原始语音信号采样后的时域波形');

xlabel('时间轴')

ylabel('幅值A')

subplot(2,1,2);

plot(f,abs(p(1:n/2)));%语音信号的频谱图

title('原始语音信号采样后的频谱图');

xlabel('频率Hz');

ylabel('频率幅值');

程序运行时可以听到声音，而后得到原始语音信号采样后的时域波形和频谱图

（4）加噪语音信号频谱分析及仿真

下面将介绍在Matlab平台上如何实现在一段原始语音信号加入一个随机噪声信号。

函数randn：产生正态分布的随机数或者矩阵的函数。产生的均值为0，方差σ^2=1,标准差σ=1的正态分布的随机数或矩阵的函数。

用法：Y=randn(n)，返回一个n*n的随机项矩阵。如果n不是个数量，则将返回错误信息。

Y=randn(m,n) 或Y=randn([m n]),返回一个m*n的随机项矩阵。

Y=randn(size(A)),返回一个和A有同样维数大小的随机数组。

randn返回一个每次都变化的数量。

利用下面的一段程序,将实现用randn函数把一段随机噪音信号加入原始语音信号的信号处理功能：

加噪后语音信号的时域波形、频谱图程序：

L=length(y)%计算音频信号的长度

noise=0.1*randn(L,2);%产生随机噪声信号

z=y+noise;%在原始语音信号后面加入噪音得到噪音信号

sound(z,fs)%播放噪音信号

%对加噪后的语音进行分析

n=length(y);%选取变换的点数

zp=fft(z,n);%对n点进行傅里叶变换到频域

f=fs*(0:n/2-1)/n; %对应点的频率

figure(2)

subplot(2,1,1);%绘出加噪语音信号的时域波形

plot(z);

title('加噪语音信号时域波形');

xlabel('时间轴')

ylabel('幅值A')

subplot(2,1,2);%绘出加噪语音信号的频谱图

plot(f,abs(zp(1:n/2)));

title('加噪语音信号频谱图');

xlabel('频率Hz');

ylabel('频率幅值');

再通过对原始语音信号的回放效果与加噪后的语音信号回放的效果的对比,可以

听出加了噪音后的语音信号有吱吱嘎嘎的声音。

（5）去噪及仿真

FIR滤波器法去噪

通过对上一节中加噪语音信号和原始语音信号频谱图对比可以知道，噪音大部分是Hz大于5000的部分，故设计低通滤波器进行滤波处理。接下来我们要用设计的FIR低通滤波器对上一节中加噪语音信号进行滤波处理。

用自己设计的FIR数字低通滤波器对加噪的语音信号进行滤波时，在Matlab 中，FIR滤波器利用函数fft对信号进行滤波。函数fft用的是重叠相加法实现线性卷积的计算。调用格式为：y=fft(h,x,M)。其中，h是系统单位冲击响应向量；x是输入序列向量；y是系统的输出序列向量，M是有用户选择的输入序列的分段长度，缺省时，默认的输入向量的重长度M=512。

设计一个FIR数字低通滤波器对加噪语音信号的滤波程序：

Ft=8000;

Fp=1000;

Fs=1200;

wp=2*Fp/Ft;%通带截止频率

ws=2*Fs/Ft; %阻带截止频率

rp=1;rs=50;%通带允许波动，%阻带衰减

p=1-10.^(-rp/20);%通带阻带波纹

s=10.^(-rs/20);

fpts=[wp ws];

mag=[1 0];dev=[p s];

[n21,wn21,beta,ftype]=kaiserord(fpts,mag,dev);%由kaiserord求滤波器的

阶数和截至频率

b21=fir1(n21,wn21,Kaiser(n21+1,beta));%由fir1设计滤波器

[h,w]=freqz(b21,1);%得到频率响应

plot(w/pi,abs(h));

title('FIR低通滤波器');

x=fftfilt(b21,z);

X=fft(x,n);

figure(4);

subplot(2,2,1);%绘出滤波前信号的频谱

plot(f,abs(zp(1:n/2)));

title('滤波前信号的频谱');

subplot(2,2,2);%绘出滤波后信号的频谱

plot(f,abs(X(1:n/2)));

title('滤波后信号的频谱');

subplot(2,2,3);%绘出滤波前信号的时域波形

plot(z);

title('滤波前信号的时域波形')

subplot(2,2,4);%绘出滤波后信号的时域波形

plot(x);

title('滤波后信号的时域波形')

sound(x,fs,bits)

得到的图像：

Fir低通滤波器

分析：

通过滤波前的信号波形和频谱图的对比，可以明显看出，滤波后的波形开始变得清晰了，开始有点接近原始信号的波形图了。滤波后信号的频谱图也在5000Hz以后开始逐渐接近原始语音信号的频谱图。

再从对语音信号的回放，人耳可以辨别出加噪后的语音信号比较浑浊，还有很明显嘎吱嘎吱的杂音在里面。滤波后，语音信号较加噪后的信号有了明显的改善，基本可以听清楚了，而且杂音也没有那么强烈，但是声音依然没有原始语音信号那么清晰脆耳。

6、语音去噪总结

通过了这次基于matlab的语音信号去噪处理的课程设计，让我有很大的体会。为实现语音信号的基本处理的功能，我采用一系列语音分析和处理技术，经过测试、运行，程序比较圆满的完成了对语音信号的读取与打开，较好的完成了对语音信号的频谱分析。通过fft变换,得出了语音信号的频谱图。滤波这一个环节，我发现我的加噪音后的频谱图在5000HZ

有一定的变化，所以我决定设计一个fir低通滤波器，来让低频的波通过，高频的波不通过。通过比较加噪前后语音的频谱，我们能看到，高于5000Hz的频率噪音已被滤除；通过声音的回放，也能明显的感觉到加入噪声后回放的声音与原始的语音信号有很大的不同,前者有较尖锐的干扰噪声，而后者噪音有所降低。从含噪语音信号的频谱图中可以看出含噪声的语音信号频谱在整个频域范围内分是布均匀。其实，这正是干扰所造成的。

本设计所使用的函数如wavread，fft，sound等基本都是Matlab中自带的库函数，因而具有操作界面简练，处理速度较快，占用内存空间比较小等特点。在进

行语音信号处理的操作界面中，由于使用的是Matlab，因而菜单按键简单、明了，每个功能只对应一个按键，省去了大型软件的操作复杂步骤。另外，本文中所有的程序基本都是固定结构，对于其他的一些不同的音频信号，只需改变其中的相应参数即可再次调用，因而调用快捷、方便。

参考文献

[1] 张虎堂.数字信号处理器(DSP)的发展趋势与应用研究[J].硅谷.2010(23);5-8

[2] 余道衡.数字信号处理的发展与应用[J].世界科技研究与发展.1999(03);25-28

[3] 苏玉萍,郑琼琼,余冬菊.基于MATLAB的FIR滤波器设计[J].中国科技信息.2008(08):14-15

[4] 沈捷,王莉.基于MATLAB的图形交互式数字信号处理教学实验系统[J].电脑开发与应用.2007(09):68

-69.

[5] 陈巍,相广进.基于MATLAB的数字信号处理教学实验的开发[J].中国科教创新导刊.2008(22):30-32

[6] 付青青,吴爱平.基于Matlab的语音信号自相关基音检测[J].长江大学学报(自科版).2006(10):58-59

[7] 张文,蒋猛,陈智.基于MATLAB的语音信号滤波设计与实现[J].山西电子技术.2008(02);33-36

[8] 徐靖涛,王金根.基于MATLAB的语音信号分析和处理[J].重庆科技学院学报(自然科学版).2008(01):

13-14

[9] 王静.基于MATLAB的语音信号LPC技术分析研究[J].微计算机信息.2007(0

9):36-40

[10] 张勇,柳玉海.MA TLAB在数字语音分析中的应用[J].计算机时代.2010(02):8 1-83

[11] 成利香,张桂新.基于MATLAB的语音信号加噪处理[J].软件导刊.2006(23):34-35

[12] 张恒.基于Matlab的小波消噪技术仿真与实现[J].中国西部科技.2010(29):103-107

[13] 高西全,丁玉美.数字信号处理(M).西安电子科技大学出版社.2003.2

[14] 郑君里.信号与系统.高等教育出版社(M).2006.4

[15] 王艳芬.数字信号处理原理及实现第2版清华大学出版社2008.