基于MATLAB的数字信号滤波处理

数字信号处理课程设计报告

（选题：基于MATLAB的数字信号滤波处理）

系 （院）： 信息学院 专 业： 10通信工程 班 级： 通信一班 学 号： 14102301239 姓 名： 粟清明

完成日期：2013年5月21日

数字信号处理课程设计

目录

1 课程设计内容求 ...............................1 1.1 设计内容.........................................1 1.2 设计要求.........................................1 2 课程设计理论基础 ..............................1 2.1 FIR滤波器........................................1 2.2 窗函数法设计FIR滤波器的原理......................2 2.3 窗函数类型的选择................................. 3 3 课程设计的具体实现及仿真结果分析...............4 3.1 语音信号采集....................................4 3.2 语音信号的时频分析............................... 4 3.3 语音信号加噪与频谱分析........................4 3.4 设计FIR数字滤波器................................5 3.4.1 低通滤波器的设计............................. 5

3.4.2 带通滤波器的设计..............................6

3.4.3 高通滤波器的设计............................. 7 3.5 用滤波器对噪音信号进行滤波........................8 3.6 比较滤波前后语音信号的波形及频谱..................8 3.6.1 低通滤波前后语音信号的波形及频谱..............8 3.6.2 带通滤波前后语音信号的波形及频谱............. 9 3.6.3 高通滤波前后语音信号的波形及频谱.............10

4 设计总结......................................11

2

数字信号处理课程设计

1 课程设计内容及要求

1.1 设计内容

了解wavread()函数的使用方法，利用wavread()函数对语音信号进行采集，将语音信号转换成计算机能够运算的有限长序列，加一定的噪声，并选择适当的滤波器进行滤波观察滤波后的效果。

1.2 设计要求

① 语音信号采集，采用(.wav)格式

② 运用wavread()函数对采集的语音信号进行时频分析

③ 运通rand()或randn()函数给语音信号加一定的噪声信号，并且对加噪后的语音信号进行时频分析

④用窗函数法设计FIR低通滤波器，其设计具体指标如下： 低通滤波器：fp=1000Hz fc=1200Hz Ap=1DB As=100DB 高通滤波器：fp=3500Hz fc=4000Hz Ap=1DB As=100DB 带通滤波器：fp1=1200Hz fc1=1000Hz Ap=1DB As=100DB fp2=3000Hz fc2=3200Hz

⑤运用fftfilt()函数对加噪语音信号进行滤波处理

2课程设计理论基础

2.1 FIR滤波器

FIR滤波器：有限长单位冲激响应滤波器，是数字信号处理系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。

有限长单位冲激响应（FIR）滤波器有以下特点： (1) 系统的单位冲激响应h(n)在有限个n值处不为零；

(2) 系统函数H(z)在|z|>0处收敛，极点全部在z=0处（因果系统）；

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

3

数字信号处理课程设计

FIR滤波器的系统函数用下式表示：。H(n)就是FIR滤波器的

单位脉冲响应。FIR滤波器最重要的优点就是由于不存在系统极点，FIR滤波器是绝对稳定的系统。

2.2窗函数法设计FIR滤波器的原理

设所希望的滤波器的理想的频率响应函数为Hdej，则其对应的单位脉冲响应为

1hdn2Heed （4.1） jjd窗函数设计法的基本原理是用有限长单位脉冲响应序列hn逼近hdn。由于hdn往往是无限长序列，而且是非因果的，所以用窗函数n将hdn截断，并进行加权处理，得到：

hnhdnn

（4.2）

其频率响应hn就作为实际设计的FIR数字滤波器的单位脉冲响应序列，函数Hej为

Hehnejn0N1jn （4.3）

式中，N为所选窗函数n的长度。

我们知道，用窗函数法设计的滤波器性能取决于窗函数n的类型及窗口长度N的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N 。各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见表4.1。

4

数字信号处理课程设计

表1 各种窗函数的基本参数

窗函数 矩形窗 三角形窗 汉宁窗 哈明窗 不莱克曼窗 凯塞窗(α=7.865) 旁瓣峰值幅度/dB 过渡带宽 阻带最小衰减/dB -13 -25 -31 -41 -57 -57 4π/N 8π/N 8π/N 8π/N 12π/N 10π/N -12 -25 -44 -53 -74 -80 这样选定窗函数类型和长度N之后，求出单位脉冲响应

hnhdnn，并按照式（4.3）求出Hej。Hej是否满足要求，要进行演算。一般在hn尾部加零使长度满足2的整数次幂，以便用FFT计算

Hej。如果要观察细节，补零点数增多即可。如果Hej不满足要求，则要重新选择窗函数类型和长度N ，再次验算，直至满足要求。

如果要求线性相位特性，则hn还必须满足

hnhN1n

根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类，例如，要设计线性相位低通特性，可以选择hnhN1n这一类，而不能选择hnhN1n这一类。

2.2窗函数类型的选择

本次课程设计选用凯塞窗（Kaiser），Kaiser窗可以通过调整参数值来折中选择主办宽度和旁办衰减，采用Kaiser窗设计FIR滤波器具有很大的灵活性。

凯塞窗是一种适应性较强且比较灵活的窗函数，它的表达式为：

5

数字信号处理课程设计

这种窗函数是近似于给定旁办电平，使主办具有最大能量意义下的最佳窗函数。其中β是形状参数，用以调节主办宽度和旁办电平，一般β选为4＜β＜9，此时旁办电平约在30dB和-67dB之间。(当β=5.44时，接近汉明；β=8.5，接近布莱克曼；β=0，为矩形。)

3 课程设计的具体实现及仿真结果分析

3.1 语音信号采集

用录音设备录下一段音频以yr.wav格式存储在电脑的中，以便接下的过程中使用。

3.2语音信号的时频分析

利用MATLAB中的“wavread”命令来读入（采集）语音信号，将它赋值给某一向量。再对其进行采样，记住采样频率和采样点数。接下来，对语音信号yr.wav进行采样,具体程序如下： [y,fs,nbits]=wavread ('yr.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));

6

数字信号处理课程设计

title('原始信号频谱') 程序结果图如下：

3.3语音信号加噪与频谱分析

在本次课程设计中，我们是利用MATLAB中的随机函数(rand或randn)产生噪声加入到语音信号中，模仿语音信号被污染，并对其频谱分析。在这里，选用Randn（m,n）函数。语音信号添加噪声及其频谱分析的主要程序如下： [y,fs,nbits]=wavread ('yr.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('加噪语音信号的频域波形'); 程序结果图如下：

7

数字信号处理课程设计

3.4设计FIR数字滤波器

用窗函数法设计FIR滤波器的步骤。如下：

（1）根据对阻带衰减及过渡带的指标要求，选择串窗数类型（矩形窗、三角窗、

汉宁窗、哈明窗、凯塞窗等），并估计窗口长度N。先按照阻带衰减选择窗函数类型。

原则是在保证阻带衰减满足要求的情况下，尽量选择主瓣的窗函数。 （2）构造希望逼近的频率响应函数。 （3）计算h(n).。 （4）加窗得到设计结果。

3.4.1低通滤波器的设计

设计指标为：fp=1000Hz fc=1200Hz Ap=1DB As=100DB 具体设计程序如下： rp=1;

rs=100; （阻带最小衰减）

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

wp=0.2; （通带截止频率） ws=0.28; （阻带截止频率） fpts=[wp ws];

8

数字信号处理课程设计

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

[n21,wn21,beta,ftype] = kaiserord(fpts,mag,dev); %计算出凯塞窗N，beta的值

b21=fir1(n21,wn21,Kaiser(n21+1,beta)); %由fir1设计滤波器 [h,w]=freqz(b21,1); %得到频率响应 plot(w/pi,abs(h)); title('FIR低通滤波器'); 程序结果图如下：

3.4.2带通滤波器的设计

设计指标：fp1=1200Hz fc1=1000Hz Ap=1DB As=100DB fp2=3000Hz fc2=3200Hz 具体设计程序如下：

[n,wn,bta,ftype]=kaiserord([0.16 0.24 0.6 0.68],[0 1 0],[0.01 0.1087 0.01]);%用kaiserord函数估计出滤波器阶数n和beta参数 b22=fir1(n,wn,ftype,kaiser(n+1,bta),'noscale'); [hh1,w1]=freqz(b22,1,256); figure(1) subplot(2,1,1)

plot(w1/pi,20*log10(abs(hh1))) grid

xlabel('归一化频率w');ylabel('幅度/db');

程序结果图如下：

9

数字信号处理课程设计

3.4.3高通滤波器的设计

设计指标：fp=3500Hz fc=4000Hz Ap=1DB As=100DB 具体设计程序如下： rp=1; rs=100;

p=1-10.^(-rp/20); %通带阻带波纹 s=10.^(-rs/20); wp=0.9; ws=0.7; fpts=[ws wp]; mag=[0 1]; dev=[p s];

[n23,wn23,beta,ftype]=kaiserord(fpts,mag,dev);

b23=fir1(n23,wn23,'high',Kaiser(n23+1,beta));%由fir1设计滤波器 [h,w]=freqz(b23,1); %得到频率响应 plot(w*12000*0.5/pi,abs(h)); title('FIR高通滤波器'); axis([3000 6000 0 1.2]); 程序结果图如下：

10

数字信号处理课程设计

3.5用滤波器对噪音信号进行滤波

用自己设计的各滤波器分别对加噪的语音信号进行滤波，在Matlab中，FIR滤波器利用函数fftfilt对信号进行滤波。 函数fftfilt用的是重叠相加法实现线性卷积的计算。调用格式为：y=fftfilter(h,x,M)。其中，h是系统单位冲击响应向量；x是输入序列向量；y是系统的输出序列向量；M是有用户选择的输入序列的分段长度，缺省时，默认的输入向量的重长度M=512。 函数filter的调用格式：yn=filter(B，A.xn)，它是按照直线型结构实现对xn的滤波。其中xn是输入信号向量，yn输出信号向量。

3.6比较滤波前后语音信号的波形及频谱

3.6.1低通滤波前后语音信号的波形及频谱

具体程序设计如下：

[y,fs,nbits]=wavread ('yr.wav');

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('滤波前信号的频谱');

11

数字信号处理课程设计

subplot(2,2,2);

plot(abs(m21),'r'); title('滤波后信号的频谱'); subplot(2,2,3);

plot(s); title('滤波前信号的波形'); subplot(2,2,4);

plot(z21); title('滤波后的信号波形'); 程序结果图如下：

3.6.2带通滤波前后语音信号的波形及频谱

具体程序如下：

[y,fs,nbits]=wavread ('yr.wav');

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

12

数字信号处理课程设计

plot(abs(m22),'r'); title('滤波后信号的频谱'); subplot(2,2,3); plot(s);

title('滤波前信号的波形'); subplot(2,2,4); plot(z22);

title('滤波后的信号波形'); 程序结果图如下：

3.6.3高通滤波前后语音信号的波形及频谱

具体设计程序如下：

%高通滤波器***************************** [y,fs,nbits]=wavread ('yr.wav');

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

13

数字信号处理课程设计

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

数字信号处理课程设计

4设计总结

本次课程设计完成了基于matlab的语音信号滤波处理，用wavread()函数对语音信号进行采集并且进行时频分析，然后用readn()函数对语音信号进行加噪处理并且进行频谱分析，最后用窗函数法设计FIR低通、高通、带通数字滤波器，分别对加噪后的语音信号进行滤波处理，得到滤波处理的频谱图。本次课程设计基本完成了设计的指标和要求，在课程设计过程中遇到了许多的问题，一开始用格式转换器把音乐文件转换成录音文件，用wavread（）函数进行采集后，发现没有时频分析图，初步判断可能是音乐文件过大，后经剪切音乐文件后，再用wavread()函数采集，就可以得到时频分析图了，在设计FIR滤波器的过程中，部分程序参考了网上的设计程序，但是当用matlab编译的时候，会有很多的错误，通过大量的查阅matlab自带的help文件，以及书本的理论知识，慢慢的把出现的错误改正，并且基本设计出了适合本次课程设计的要求的FIR数字滤波器，通过仿真后也基本实现了相应的滤波功能。通过本次课程设计进一步巩固了之前所学习的数字信号处理的理论知识，对用窗函数法设计FIR滤波器的方法以及技巧有了更深刻的理解，同时也更好的掌握了用matlab进行应用设计，深深体会到了matlab软件功能之强大。在课程设计过程中，一开始感觉不知从何下手，感觉是那么的困难，但通过不断的查阅资料以及认真理解所学的理论知识，一步一步的往下做时，感觉并不是那么的难，最终基本完成了本次课程设计，感觉收获还是挺多的，进一步提高了自己的自主学习能力。

15

数字信号处理课程设计

参考文献

[1] 丁玉美，高西全，数字信号处理. 3 版. 西安：西安电子科技大学出

版社，2008.

[2] 刘敏，魏玲．Matlab通信仿真与应用〔M〕．北京：国防工业出版社，2001． [3] 张贤达 编, 《现代信号处理（第二版）》，清华大学出版社，2002 年.

16