摘 要

本设计综合利用数字信号处理的理论知识进行语音信号的频谱分析，通过理论推导得出相应结论，再利用MATLAB作为编程工具进行计算机实现，从而加深对所学知识的理解，建立概念。本次课程设计要求利用MATLAB对语音信号进行分析和处理，要求学生采集语音信号后，在MATLAB软件平台进行频谱分析；并对所采集的语音信号加入干扰噪声，对加入噪声的信号进行频谱分析，设计合适的滤波器滤除噪声，恢复原信号。待处理语音信号是一个在20Hz～20kHz频段的低频信号。采用了高效快捷的开发工具——MATLAB，实现了语音信号的采集，对语音信号加噪声及设计滤波器滤除噪声的一系列工作。利用采样原理设计了高通滤波器、低通滤波器、带通滤波器、带阻滤波器。同学通过查阅资料自己获得程序进行滤波器的设计，能过得到很好的锻炼。

关键词：MATLAB 滤波器 数字信号处理

目录

第一章 绪论

1.1 设计的目的及意义

1.2 设计要求

1.3 设计内容

第二章 系统方案论证

2.1 设计方案分析

2.2 实验原理

第三章 信号频谱分析

3.1 原始信号及频谱分析

3.2 加入干扰噪声后的信号及频谱分析

第四章 数字滤波器的设计与实现

4.1 高通滤波器的设计

4.2 低通滤波器的设计

4.3 带通滤波器的设计

4.4 带阻滤波器的设计

第五章 课程设计总结

参考文献

附录Ⅰ

附录Ⅱ

第一章 绪论

1.1 设计的目的及意义

综合利用数字信号处理的理论知识进行语音信号的频谱分析，通过理论推导得出相应结论，再利用MATLAB作为编程工具进行计算机实现，从而加深对所学知识的理解，建立概念。本设计采用了高效快捷的开发工具——MATLAB，实现了语音信号的采集，对语音信号加噪声及设计滤波器滤除噪声的一系列工作。

1.2 设计要求

基本要求：本次课程设计要求利用MATLAB对语音信号进行分析和处理，要求学生采集语音信号后，在MATLAB软件平台进行频谱分析；并对所采集的语音信号加入干扰噪声，对加入噪声的信号进行频谱分析，设计合适的滤波器滤除噪声，恢复原信号。待处理语音信号是一个在20Hz～20kHz频段的低频信号。

1.3 设计内容

选择一个wav文件作为分析的对象，或录制一段语音信号，对其进行频谱分析，分别对加噪前后的语音信号进行频谱分析，再通过不同滤波器根据信号的频谱特点重构语音信号，选出最佳滤波方案。

原理：

(1)采样定理

在进行模拟/数字信号的转换过程中，当采样频率fs.max大于信号中，最高频率fmax的2倍时，即：fs.max>=2fmax,则采样之后的数字信号完整地保留了原始信号中的信息，一般实际应用中保证采样频率为信号最高频率的5～10倍；采样定理又称奈奎斯特定理：

1924年奈奎斯特(Nyquist)就推导出在理想低通信道的最高大码元传输速率的公式:理想低通信道的最高大码元传输速率=2W*log2 N (其中W是理想低通信道的带宽,N是电平强度)

(2)采样频率

采样频率是指计算机每秒钟采集多少个声音样本，是描述声音文件的音质、音调，衡量声卡、声音文件的质量标准。采样频率越高，即采样的间隔时间越短，则在单位时间内计算机得到的声音样本数据就越多，对声音波形的表示也越精确。采样频率与声音频率之间有一定的关系，根据奎斯特理论，只有采样频率高于声音信号最高频率的两倍时，才能把数字信号表示的声音还原成为原来的声音。这就是说采样频率是衡量声卡采集、记录和还原声音文件的质量标准。

(3)采样位数与采样频率

采样位数即采样值或取样值，用来衡量声音波动变化的参数，是指声卡在采集和播放声音文件时所使用数字声音信号的二进制位数。采样频率是指录音设备在一秒钟内对声音信号的采样次数，采样频率越高声音的还原就越真实越自然。

采样位数和采样率对于音频接口来说是最为重要的两个指标，也是选择音频接口的两个重要标准。无论采样频率如何，理论上来说采样的位数决定了音频数据最大的力度范围。每增加一个采样位数相当于力度范围增加了6dB。采样位数越多则捕捉到的信号越精确。对于采样率来说你可以想象它类似于一个照相机，44.1kHz意味着音频流进入计算机时计算机每秒会对其拍照达441000次。显然采样率越高，计算机摄取的图片越多，对于原始音频的还原也越加精确。

第二章 系统方案论证

2.1 设计方案分析

信号的傅立叶表示在信号的分析与处理中起着重要的作用。因为对于线性系统来说，可以很方便地确定其对正弦或复指数和的响应，所以傅立叶分析方法能完善地解决许多信号分析和处理问题。另外，傅立叶表示使信号的某些特性变得更明显，因此，它能更深入地说明信号的各项红物理现象。

由于语音信号是随着时间变化的，通常认为，语音是一个受准周期脉冲或随机噪声源激励的线性系统的输出。输出频谱是声道系统频率响应与激励源频谱的乘积。声道系统的频率响应及激励源都是随时间变化的，因此一般标准的傅立叶表示虽然适用于周期及平稳随机信号的表示，但不能直接用于语音信号。由于语音信号可以认为在短时间内，近似不变，因而可以采用短时分析法。

本实验要求掌握傅里叶分析原理，会利用已学的知识，编写程序估计短时谱、倒谱，画出语谱图，并分析实验结果，在此基础上，借助频域分析方法所求得的参数分析语音信号的基音周期或共振峰。

2.2 实验原理

1、短时傅立叶变换

由于语音信号是短时平稳的随机信号，某一语音信号帧的短时傅立叶变换的定义为：

（2-1）

其中w(n-m)是实窗口函数序列，n表示某一语音信号帧。令n-m=k'，则得到

（2-2）

于是可以得到

（2-3）

假定

（2-4）

则可以得到

（2-5）

同样，不同的窗口函数，将得到不同的傅立叶变换式的结果。由上式可见，短时傅立叶变换有两个变量：n和ω，所以它既是时序n的离散函数，又是角频率ω的连续函数。与离散傅立叶变换逼近傅立叶变换一样，如令ω=2πk/N，则得离散的短时傅立叶吧如下：

(2-6)

2、语谱图

水平方向是时间轴，垂直方向是频率轴，图上的灰度条纹代表各个时刻的语音短时谱。语谱图反映了语音信号的动态频率特性，在语音分析中具有重要的实用价值。被成为可视语言。

语谱图的时间分辨率和频率分辨率是由窗函数的特性决定的。时间分辨率高，可以看出时间波形的每个周期及共振峰随时间的变化，但频率分辨率低，不足以分辨由于激励所形成的细微结构，称为宽带语谱图；而窄带语谱图正好与之相反。

宽带语谱图可以获得较高的时间分辨率，反映频谱的快速时变过程；窄带语谱图可以获得较高的频率分辨率，反映频谱的精细结构。两者相结合，可以提供带两与语音特性相关的信息。语谱图上因其不同的灰度，形成不同的纹路，称之为“声纹”。声纹因人而异，因此可以在司法、安全等场合得到应用。

3、复倒谱和倒谱

复倒谱

是x(n)的Z变换取对数后的逆Z变换，其表达式如下:

（2-7）

倒谱c(n)定义为x(n)取Z变换后的幅度对数的逆Z变换，即

（2-8）

在时域上，语音产生模型实际上是一个激励信号与声道冲激响应的卷积。对于浊音，激励信号可以由周期脉冲序列表示；对于清音，激励信号可以由随机噪声序列表示。声道系统相当于参数缓慢变化的零极点线性滤波器。这样经过同态处理后，语音信号的复倒谱，激励信号的复倒谱，声道系统的复倒谱之间满足下面的关系：

（2-9）

由于倒谱对应于复倒谱的偶部，因此倒谱与复倒谱具有同样的特点，很容易知道语音信号的倒谱，激励信号的倒谱以及声道系统的倒谱之间满足下面关系：

（2-10）

浊音信号的倒谱中存在着峰值，它的出现位置等于该语音段的基音周期，而清音的倒谱中则不存在峰值。利用这个特点我们可以进行清浊音的判断，并且可以估计浊音的基音周期。

第三章 信号频谱分析

3.1 原始信号及频谱分析

选择一个wav文件作为分析的对象，可以利用Windows下的录音机或其他软件，录制一段自己的话音，在MATLAB中，[y,fs,bits]=wavread('Blip',[N1 N2]);用于读取语音，采样值放在向量y中，fs表示采样频率(Hz)，bits表示采样位数。[N1 N2]表示读取的值从N1点到N2点的值。 sound(y); 用于对声音的回放。向量y则就代表了一个信号，也即一个复杂的“函数表达式”，也可以说像处理一个信号的表达式一样处理这个声音信号。 下面是语音信号在MATLAB中的语言程序，它实现了语音的读入与打开，并绘出了语音信号时域波形，然后对语音信号进行频谱分析。在MATLAB中，可以利用函数fft对信号进行快速傅里叶变化，得到信号的频谱特性。 在频谱特性中分析最大值的位置（可能有几个），它代表的频率和时域的采样时间有关，相邻的两点之间的距离为。其中，N是离散傅里叶变换用的点数，是采样的时间，前面在读取 wav文件时得到了采样频率。 既然知道了该声波的频谱，按频率就可以反演它的时域值，利用以上分析的主要峰值来重构声波。由于没有考虑相位和其他的频谱分量，所以波形和原来的波形相差甚大，但大体的频率是没错的。如图3-1所示。

图3-1 原始信号频谱图

3.2 加入干扰噪声后的信号及频谱分析

wav语音信号加噪声 在MATLAB软件平台下，给原始的语音信号叠加上噪声，噪声类型分为如下几种：（1）单频噪色（正弦干扰）；（2）高斯随机噪声。绘出加噪声后的语音信号时域和频谱图，在视觉上与原始语音信号图形对比，也可通过Windows播放软件从听觉上进行对比，分析并 体会含噪语音信号频谱和时域波形的改变。本实验采用正弦干扰。

clc;

clear;

fs=22050; %语音信号采样频率为22050 [x,fs,bits]=wavread('C:\Documents and Settings\Administrator\桌面\语音\音频.wav');

%读取语音信号的数据，赋给变量x y1=fft(x,4096);

%对信号做4096点FFT变换

f=fs*(0:511)/4096;

t=(0:length(x)-1)/22050;

x1=[0.05*sin(2*pi*10000*t)]';

x2=x+x1; sound(x2,fs,bits);

figure(1) subplot(2,1,1) plot(x)

%做原始语音信号的时域图形 title('原语音信号时域图')

subplot(2,1,2) plot(x2)

%做原始语音信号的时域图形 title('加高斯噪声后语音信号时域图') xlabel('time n'); ylabel('fudu');

y2=fft(x2,4096);

figure(2) subplot(2,1,1) plot(abs(y1)) title('原始语音信号频谱');

xlabel('Hz');

ylabel('fudu');

subplot(2,1,2) plot(abs(y2)) title('加噪语音信号频谱');

xlabel('Hz');

ylabel('fudu');

axis([0 4500 0 300]);

wavwrite(x2,fs,'C:\Documents and Settings\Administrator\桌面\语音\加噪.wav');

调试分析过程描述：

I.语音信号的分析及处理方法

a.语音的录入与打开

在MATLAB中，[y,fs,bits]=wavread('Blip',[N1 N2]);用于读取语音，采样值放在向量y中，fs表示采样频率(Hz)，bits表示采样位数。[N1 N2]表示读取从N1点到N2点的值（若只有一个N的点则表示读取前N点的采样值）。

sound(x,fs,bits); 用于对声音的回放。向量y则就代表了一个信号（也即一个复杂的“函数表达式”）也就是说可以像处理一个信号表达式一样处理这个声音信号。

b.时域信号的FFT分析

FFT即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。

在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为正整数。函数执行N点的FFT，若x为向量且长度小于N，则函数将x补零至长度N；若向量x的长度大于N，则函数截短x使之长度为N；若x 为矩阵，按相同方法对x进行处理。

II.数字滤波器设计原理

数字滤波是数字信号分析中最重要的组成部分之一，与模拟滤波相比，它具有精度和稳定性高、系统函数容易改变、灵活性强、便于大规模集成和可实现多维滤波等优点。在信号的过滤、检测和参数的估计等方面，经典数字滤波器是使用最广泛的一种线性系统。

数字滤波器的作用是利用离散时间系统的特性对输入信号波形(或频谱)进行加工处理，或者说利用数字方法按预定的要求对信号进行变换。

数字滤波器的设计步骤：

不论是IIR滤波器还是FIR滤波器的设计都包括三个步骤：

a. 按照实际任务的要求，确定滤波器的性能指标。

b. 用一个因果、稳定的离散线性时不变系统的系统函数去逼近这一性能指标。根据不同的要求可以用IIR系统函数，也可以用FIR系统函数去逼近。

c. 利用有限精度算法实现系统函数，包括结构选择、字长选择等。

III.IIR滤波器与FIR滤波器的性能比较

有限冲击响应

FIR:Finite Impulse response，IIR:Infinite Impulse response，无限冲击响应

从性能上来说，IIR滤波器传输函数的极点可位于单位圆内的任何地方，因此可用较低的阶数获得高的选择性，所用的存贮单元少，所以经济而效率高。但是这个高效率是以相位的非线性为代价的。选择性越好，则相位非线性越严重。相反，FIR滤波器却可以得到严格的线性相位，然而由于FIR滤波器传输函数的极点固定在原点，所以只能用较高的阶数达到高的选择性；对于同样的滤波器设计指标，FIR滤波器所要求的阶数可以比IIR滤波器高5~10倍，结果，成本较高，信号延时也较大；如果按相同的选择性和相同的线性要求来说，则IIR滤波器就必须加全通网络进行相位较正，同样要大增加滤波器的节数和复杂性。

整体来看，IIR滤波器达到同样效果阶数少，延迟小，但是有稳定性问题，非线性相位；FIR滤波器没有稳定性问题，线性相位，但阶数多，延迟大。如图3-2所示。

图3-2 加入噪声后的信号频谱图

第四章 数字滤波器的设计与实现

4.1 高通滤波器的设计

（1）模拟低通滤波器原型设计。

在此单元模块主要依据设计参数完成模拟低通滤波器的原型设计。 模拟原型滤波器指的是截止频率为1的滤波器。此处的模拟低通滤波器原型以巴特沃斯（Butterworth）低通滤波器为模型进行逼近。MATLAB信号处理工具箱提供Butterworth模拟低通滤波器原型设计函数buttap,函数调用形式为： [z,p,k]=buttap(N) 式中，N为butterworth滤波器阶数；z,p,k分别为滤波器的零点、极点和增益。又Butterworth低通滤波器的幅度函数只由阶数N控制，假定阶数为1，则语句描述为： [z,p,k]=buttap[1]; 再求解模拟低通滤波器的分子分母系数。最后利用MATLAB工具箱提供的模拟滤波器()aHs的频率响应函数freqs和相关绘图函数求出其传递函数图形。其语句描述为：

[H,w]=freqs(B,A,n);

magH2=(abs(H)).^2;

hold on;

plot(w,magH2);

end xlabel('w/wc');

ylabel('|H(jw)|^2');

title('Butterworth模拟原型滤波器');

至此，可得到模拟低通滤波器的原型，该滤波器的截止频率为为1。

（2）频率转换，得到模拟高通滤波器。

此单元模块可以利用MATLAB工具箱提供的模拟频率变换函数求得。即利用函数lp2hp即可由模拟低通滤波器得到模拟高通滤波器。语句描述为： [Bt,At]=lp2hp(B,A,wc); 其中wc为题设所要求的模拟高通滤波器的阻带截止角频率，即 wc=2*pi*f=2*pi*200=400*pi 至此即得到符合设计参数要求的阻带频率为200Hz的模拟高通滤波器。 最后利用MATLAB工具箱提供的模拟滤波器()aHs的频率响应函数freqs和相关绘图函数求出各响应曲线，以便对滤波器的性能有更为直观的了解和认识。语句描述为：

[h,w]=freqs(Bt,At);

subplot();

plot(w,20*log10(abs(h)));

grid;

xlabel(‘w/pi’);

ylabel(‘频率响应/dB’);

title('模拟高通滤波器');

如图4-1所示。

图4-1 高通滤波器频谱图

4.2 低通滤波器的设计

巴特沃斯低通滤波，对加入高斯随机噪声和正弦噪声的语音信号进行滤波。用双线性变换法设计了巴特沃斯数字低通IIR滤波器对两加噪语音信号进行滤波，并绘制了巴特沃斯低通滤波器的幅度图和两加噪语音信号滤波前后的时域图和频谱图。

clear all;

fb = 1000;

fc = 1200;

fs = 22050;

wp=0.1*pi;

ws=0.4*pi;

Rp=1;

Rs=15;

Fs=22050;

Ts=1/Fs;

wp1=2/Ts*tan(wp/2); %将模拟指标转换成数字指标 ws1=2/Ts*t

an(ws/2);

[N,Wn]=buttord(wp1,ws1,Rp,Rs,'s'); %选择滤波器的最小阶数

[Z,P,K]=buttap(N); %创建butterworth模拟滤波器 [Bap,Aap]=zp2tf(Z,P,K);

[b,a]=lp2lp(Bap,Aap,Wn);

[bd,ad]=bilinear(b,a,Fs);

%用双线性变换法实现模拟滤波器到数字滤波器的转换

[h,w]=freqz(bd,ad);

figure(1) subplot(111);

plot(w*fs/(2*pi),abs(h)) grid;

title('滤波器的性能分析');

pause;

figure(2) [x,fs,bits]=wavread('C:\Documents and Settings\Administrator\桌面\语音\加噪.wav');

n=length(x);

f=fs*(0:(n/2-1))/n;

X=fft(x);

z=filter(bd,ad,x);

subplot(211);

plot(x);

title('原始信号的波形');

subplot(212);

plot(z);

title('滤波后信号的波形');

pause;

figure(3) sound(z,fs,bits);

subplot(211);

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

title('原始信号的频谱');

xlabel('Hz');

Z=fft(z);

subplot(212);

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

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

xlabel('Hz');

wavwrite(z,fs,'C:\Documents and Settings\Administrator\桌面\语音\巴滤.wav');

如图4-2所示。

图4-2 低通滤波器频谱图

4.3 带通滤波器的设计

带通滤波器设计步骤：

1、根据需求选择合适的低通滤波器原型

2、把带通滤波器带宽作为低通滤波器的截止频率，根据抑制点的频率距离带通滤波器中心频点距离的两倍作为需要抑制的频率，换算抑制频率与截止频率的比值，得出的值，然后根据m值选择低通滤波器的原型参数值。

滤波器的时域特性 ：

任何信号通过滤波器都会产生时延。Bessel filter是特殊的滤波器在于对于通带内的所有频率而言，引入的时延都是恒定的。这就意味着相对于输入，输出信号的相位变化与工作的频率是成比例的。而其他类型的滤波器（如Butterworth, Chebyshev,inverse Chebyshev,andCauser）在输出信号中引入的相位变化与频率不成比例。相位随频率变化的速率称之为群延迟（group delay）。群延迟随滤波器级数的增加而增加。如图4-3所示。

图4-3带通滤波器频谱图

4.4 带阻滤波器的设计

带阻滤波器广泛应用于无线通信系统中，用来抑制高功率发射机的杂散输出以及非线性功放或带通滤波器产生的寄生通带等。将带通滤波器的各种拓扑结构应用于带阻滤波器是可行的，而提取各种拓扑结构滤波器的等效电路参数成为设计滤波器的关键。

function y=bands(x,f1,f3,fsl,fsh,rp,rs,Fs)

%带阻滤波

%使用注意事项：通带或阻带的截止频率与采样率的选取范围是不能超过采样率的一半

%即，f1,f3,fs1,fsh,的值小于 Fs/2

%x:需要带通滤波的序列

% f 1：通带左边界

% f 3：通带右边界

% fs1：衰减截止左边界

% fsh：衰变截止右边界

%rp：边带区衰减DB数设置

%rs：截止区衰减DB数设置

%FS：序列x的采样频率

% f1=300;

f3=500;

%通带截止频率上下限

% fsl=200;fsh=600;

%阻带截止频率上下限

% rp=0.1;

rs=30;

%通带边衰减DB值和阻带边衰减DB值

% Fs=2000;%采样率

%wp1=2*pi*f1/Fs;

wp3=2*pi*f3/Fs;

wsl=2*pi*fsl/Fs;

wsh=2*pi*fsh/Fs;

wp=[wp1 wp3];

ws=[wsl wsh]; % % 设计切比雪夫滤波器；

[n,wn]=cheb1ord(ws/pi,wp/pi,rp,rs);

[bz1,az1]=cheby1(n,rp,wp/pi,'stop'); %查看设计滤波器的曲线

[h,w]=freqz(bz1,az1,256,Fs);

h=20*log10(abs(h));

figure;plot(w,h);

title('所设计滤波器的通带曲线');

grid on;

y=filter(bz1,az1,x); end 使用例子 %带阻滤波器测试

fs=1000;

t=(1:fs)/fs;

y=sin(2*pi*100*t)+sin(2*pi*150*t)+sin(2*pi*200*t);

figure;hua_fft(y,fs,1);

z=bands(y,110,190,140,160,0.1,30,fs);

figure;hua_fft(z,fs,1);

如图4-4所示。

图4-4 带阻滤波器频谱图

参考文献

[1]周辉，董正宏，数字信号处理及MATLAB实现，北京希望出版社，2006

[2]王树勋.数字信号处理处理基础及试验.北京：机械工业出版社，1992

[3]井上伸雄.数字信号处理的应用.北京：科学出版社，1991

[4]郑君里，杨为理.信号与系统（第二版），高等教育出版社，1981

[5]程佩青。数字信号处理教程［M］。北京：清华大学出版社，2002。
[6]刘敏，魏玲。Matlab通信仿真与应用［M］．北京：国防工业出版社，2001。

附录Ⅰ

四种滤波器频谱图

附录Ⅱ

GUI界面：

GUI前面板界面

GUT仿真界面