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

—— 用MatLab对语音信号进行频谱分析及滤波 器设计

学院：物联网工程学院

班系：

姓名：赵一畅

QQ号：1320846809

指导老师：

第一章 实验任务

录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法和双线性变换设计滤波器，并画出滤波器的频率响应；然后用自己设计的IIR滤波器和FIR滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号。

第二章 实验原理

2.1 采样频率、位数及采样定理

采样频率[2]，也称为采样速度或者采样率，定义了每秒从连续信号中提取并组成离散信号的采样个数，它用赫兹（Hz）来表示。采样频率的倒数是采样周期或者叫作采样时间，它是采样之间的时间间隔。通俗的讲采样频率是指计算机每秒钟采集多少个声音样本，是描述声音文件的音质、音调，衡量声卡、声音文件的质量标准。采样频率越高，即采样的间隔时间越短，则在单位时间内计算机得到的声音样本数据就越多，对声音波形的表示也越精确。

采样位数可以理解为声卡处理声音的解析度。这个数值越大，解析度就越高，录制和回放的声音就越真实。我们首先要知道：电脑中的声音文件是用数字0和1来表示的。所以在电脑上录音的本质就是把模拟声音信号转换成数字信号。反之，在播放时则是把数字信号还原成模拟声音信号输出。

采样定理又称奈奎斯特定理[2]，在进行模拟/数字信号的转换过程中，当采样频率fs不小于信号中最高频率fm的2倍时，采样之后的数字信号完整地保留了原始信号中的信息，一般实际应用中保证采样频率为信号最高频率的5～10倍。

2.2时域信号的FFT分析

信号的频谱分析就是计算机信号的傅里叶变换[5]。连续信号与系统的傅里叶分析显然不便于用计算机进行计算，使其应用受到限制。而FFT是一种时域和频域均离散化的变换，适合数值运算，成为用计算机分析离散信号和系统的有力工具。对连续信号和系统，可以通过时域采样，应用DFT进行近似谱分析。

2.3 IIR数字滤波器设计原理

利用双线性变换设计IIR滤波器（巴特沃斯数字低通滤波器的设计）[6]，首先要设计出满足指标要求的模拟滤波器的传递函数Ha(s),然后由Ha(s)通过双线性变换可得所要设计的IIR滤波器的系统函数H(z)。如果给定的指标为数字滤波器的指标，则首先要转换成模拟滤波器的技术指标，这里主要是边界频率Wp和Ws的转换，对ap和as指标不作变化。边界频率的转换关系为 ∩=2/T tan(w/2)。接着，按照模拟低通滤波器的技术指标根据相应设计公式求出滤波器的阶数N和3dB截止频率 ∩c ；根据阶数N查巴特沃斯归一化低通滤波器参数表，得到归一化传输函数Ha(p)；最后，将p=s/ ∩c 代入Ha(p)去归一，得到实际的模拟滤波器传输函数Ha(s)。之后，通过双线性变换法转换公式s=2/T((1-1/z)/(1+1/z))得到所要设计的IIR滤波器的系统函数H(z)。

2.4 FIR数字滤波器设计原理

基于窗函数的FIR 数字滤波器的设计方法通常也称之为傅立叶级数法，是用一定宽度窗函数截取无限脉冲响应序列，获得有限长的脉冲响应序列，从而得到FIR 滤波器。它是在时域进行的，由理想滤波器的频率响应word/media/image2.gif推导出其单位冲激响应hd(n)，再设计一个FIR数字滤波器的单位冲激响应h(n)去逼近hd(n)，表示word/media/image3.gif=word/media/image4.gifword/media/image5.gif由此得到的离散滤波器的系统传递函数Hd(z) 为word/media/image2.gif＝word/media/image6.gif，该hd(n) 为无限长序列，因此Hd(z)是物理不可实现的。为了使系统变为物理可实现的，且使实际的FIR 滤波器频率响应尽可能逼近理想滤波器的频率响应，采用窗函数将无限脉冲响应hd(n)截取一段h(n)来近似表示hd(n)，可得：h (n) = hd(n)w(n) ,从而有：式中N 表示窗口长度，这样H(z)就是物理可实现的系统。并且从线性相位FIR滤波器的充要条件可知，为了获得线性相位FIR 数字滤波器的冲激响应h(n) ，那么序列h(n) 应有τ= (N −1) / 2的延迟。由于窗函数的选择对结果起着重要的作用，针对不同的信号和不同的处理目的来确定窗函数的选择才能收到良好的效果。

第三章 实验过程

一、原始语音信号采样后的时域、频域分析

1、源程序如下：

时域波形：

z1=wavread(' F:\Study\数字信号处理\实验\Qnylan.wav')

sound(z1,22050); %播放语音信号

figure(1)

plot(z1);

title('时域图');

xlabel('时间（n）');

ylabel('幅值（n）');

%图见3-1

频域波形：

z2=wavread(' F:\Study\数字信号处理\实验\Qnylan.wav');

y1=z2(1:8192);

y1=fft(y1);

n=0:8191;

sound(z1,22050); %播放语音信号

figure(2)

plot(n,y1);

title('频谱图');

xlabel('时间（n）');

ylabel('幅值（n）');

%图见3-2

2、运行结果如图3-1、3-2所示：

图3-1 原始语音采样后时域信号

图3-2 原始语音信号采样后频率响应

二、采样后信号的FFT变换分析

1、源程序如下：

[x1,fs,bits]=wavread(' F:\Study\数字信号处理\实验\Qnylan.wav');

%读取语音信号的数据，赋给变量x1

y1=fft(x1,1024); %对采样后信号做1024点FFT变换

f=fs*(0:511)/1024; %生成频率序列

figure(1)

subplot(2,1,1);

plot(f,abs(y1(1:512)));

title('原始语音信号频谱')

xlabel('频率/Hz');

ylabel('幅值')

subplot(2,1,2);

plot(abs(y1(1:1024))) %采样后信号的FFT频谱图

title('原始语音信号FFT频谱')

xlabel('点数N');

ylabel('幅值');

2、运行结果如图3-3所示：

图3-3 原始语音信号频谱及FFT频谱

3、频谱分析：

从图3-3可以看出，采样点数为1024，离散点数越多，越接近原始信号频谱。

三、双线性变换法设计IIR数字滤波器

1) 低通滤波器性能指标 fp=1000Hz，fst=1200Hz，Ap=1db，As=100db。

（1）源程序如下：

fs=22050;

[x2]=wavread(' F:\Study\数字信号处理\实验\Qnylan.wav');

Ts=1/fs;R1=10;

wp=2*pi*1000/fs; %通带截止频率

ws=2*pi*1200/fs; %阻带截止频率

Rp=1; %通带衰减

Rl=100; %阻带衰减

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

ws1=2/Ts*tan(ws/2);

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

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

[Bap,Aap]=zp2tf(Z,P,K);

[b,a]=lp2lp(Bap,Aap,Wn); %将模拟原型低通滤波器转换为低通滤波器

[bz,az]=bilinear(b,a,fs); %用双线性变换法实现模拟滤波器到数字滤波器的转换

[H,W]=freqz(bz,az); %绘制频率响应曲线

figure(1)

plot(W*fs/(2*pi),abs(H))

grid

xlabel('频率／Hz')

ylabel('频率响应幅度')

title('IIR低通滤波器')

f1=filter(bz,az,x2);

figure(2)

subplot(2,1,1)

plot(x2) %画出滤波前的时域图

title('IIR低通滤波器滤波前的时域波形');

subplot(2,1,2)

plot(f1); %画出滤波后的时域图

title('IIR低通滤波器滤波后的时域波形');

sound(f1,44100); %播放滤波后的信号

F0=fft(f1,1024);

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

figure(3)

y2=fft(x2,1024);

subplot(2,1,1);

plot(f,abs(y2(1:512))); %画出滤波前的频谱图

title('IIR低通滤波器滤波前的频谱')

xlabel('频率/Hz');

ylabel('幅值');

subplot(2,1,2)

F1=plot(f,abs(F0(1:512))); %画出滤波后的频谱图

title('IIR低通滤波器滤波后的频谱')

xlabel('频率/Hz');

ylabel('幅值');

（2）运行结果如图3-4、3-5和3-6所示：

图3-4 IIR低通滤波器

图3-5 IIR低通滤波器滤波前后时域波形

图3-6 IIR低通滤波器滤波前后的频谱

（3）频谱分析：

从图3-6可以看出，经过IIR低通滤波器滤波后将高频部分滤除了。

四、窗函数法设计FIR数字滤波器

1）低通滤波器性能指标 fp=1000Hz，fst=1200Hz，Ap=1db，As=100db。

（1）源程序如下：

fs=10000;

x1=wavread(' F:\Study\数字信号处理\实验\Qnylan.wav');

wp=2*pi*1000/fs;

ws=2*pi*1200/fs;

Rp=1;

Rs=100;

wdelta=ws-wp;

N=ceil(8*pi/wdelta); %取整

wn=(wp+ws)/2;

[b,a]=fir1(N,wn/pi,hamming(N+1)); %选择窗函数，并归一化截止频率

figure(1)

freqz(b,a,512);

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

f2=filter(b,a,x1);

figure(2)

subplot(2,1,1)

plot(x1)

title('FIR低通滤波器滤波前的时域波形');

subplot(2,1,2)

plot(f2);

title('FIR低通滤波器滤波后的时域波形');

sound(f2,44100); %播放滤波后的语音信号

F0=fft(f2,1024);

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

figure(3)

y2=fft(x1,1024);

subplot(2,1,1);

plot(f,abs(y2(1:512)));

title('FIR低通滤波器滤波前的频谱')

xlabel('频率/Hz');

ylabel('幅值');

subplot(2,1,2)

F2=plot(f,abs(F0(1:512)));

title('FIR低通滤波器滤波后的频谱')

xlabel('频率/Hz');

ylabel('幅值');

（2）运行结果如图3-7、3-8和3-9所示：

图3-7 FIR低通滤波器

图3-8 FIR低通滤波器滤波前后时域波形

图3-9 FIR低通滤波器滤波前后频谱

（3）频谱分析：

从图3-9可以看出，经过FIR低通滤波器滤波后，将高频部分滤除了

第四章 心得体会

本次的数字信号处理综合实验的题目是应用Matlab对语音信号进行频谱分析及滤波，首先通过网络和书籍查找有关本次综合实验的资料，编写相关程序，并通过Matlab软件运行得到相关波形频谱图。

实验中利用双线性变换法设计IIR数字滤波器，利用窗函数设计FIR数字滤波器，可以是低通、高通和带通滤波器的设计，此次试验中选用的是低通滤波器，而低通滤波器是基于Butterworth滤波器来设计的，通过设计的滤波器对语音信号进行滤波，再对得出的频谱图进行分析。在实验中遇到一些困难，在设计数字滤波器的时候，通带频率和阻带频率的选取要满足低通的要求，以及通带允许的最大衰减和阻带应达到的最小衰减。

在这次的综合实验中，总的来说并不是很辛苦，实验将上课所学的理论知识运用到实践中。通过这次应用Matlab对语音信号进行频谱分析及滤波的综合实验，让我对Matlab的应用以及数字滤波器的设计有了更深层次的理解，每个程序中的语句表示什么意思也有了很清楚的了解。在实践中增强了我的动手能力，以及增强了我的团队意识，并提高了我的综合能力，使自身得到了很大的锻炼。最后要感谢老师的悉心指导和帮助，

参考文献

[1] 丁玉美《数字信号处理(第三版)》[M]. 西安电子科技大学出版社，2008.8

[2] 程佩青.《数字信号处理教程》[M]. 北京：清华大学出版社,2007.2

[3] 余成波、陶红艳等.《数字信号处理及MATLAB实现（第二版）》[M]. 北京：清华大学出版社，2007

应用MatLab对语音信号进行频谱分析及滤波