数字信号处理课程设计语音信号处理MATLAB

中原工学院

数字信号处理教程课程设计

电信152

杨耀华

绪论

数字信号处理课程是电子信息类与电气类专业本科生继“信号与系统”课之后的一门必修的专业基础课程。设置本课程的目的在于，使学生通过本课程的学习，了解“数字信号处理”这一技术领域的概貌，初步建立起有关“数字信号处理”的基本概念，掌握基本分析方法，为后续课程及从事信息处理等方面有关的研究工作打下基础。本课程是一门结合实际工程应用的基础理论课程。本次课程设计是在MATLAB平台上，对声音信号进行采集，处理，滤波等最终还原出原无噪声的声音信号。运用本课程所学的理论知识对信号进行谱分析，设计滤波器，得出结论。进一步巩固所学的知识。

目录

一．设计目的 (1)

二．设计要求及任务 (1)

三．课程设计平台 (2)

四．设计原理及计算方法 (2)

五．实验论证方案及结果分析 (3)

六. 结论及心得 (15)

七. 附录：程序代码及注释 (15)

八. 参考文献 (28)

一．设计目的

1.学会MATLAB的使用，掌握MATLAB程序设计方法。

2.掌握在windows环境下语音信号采集的方法。

3.掌握数字信号处理的基本概念，基本理论和基本方法。

4.掌握MATLAB设计IIR数字滤波器的方法。

学会用MATLAB对信号进行分析和处理。

二．设计要求及任务

2.1.语音信号的采集；

本设计利用计算机Windows下的录音机录入一句语音信号，然后在Matlab软件平台下，利用函数waveread对语音信号进行采样，记住采样频率和采样点数。

2.2.语音信号的频谱分析；

在Matlab中，可以利用函数FFT对信号进行快速傅里叶变换，得到信号的频谱特性，然后加入一干扰信号，要求画出语音信号干扰前后的时域波形，并对其频谱进行分析。

2.3.设计数字滤波器，给出性能指标（参考指标）；

（1）低通滤波器的性能指标：fp=1000Hz,fs=1200Hz,As=100dB,Ap=1dB

（2）高通滤波器的性能指标：fs=4800Hz,fp=5000Hz,As=100dB,Ap=1dB

（3）带通滤波器的性能指标：fp1=1200Hz,fp2=3000Hz,fs1=1000Hz,

fps2=3200Hz,As=100dB,Ap=1dB；

采用双线性变换法设计上面一种类型的数字滤波器，要求使用切比雪夫II型滤波器。

2.4.用滤波器对信号进行滤波

画出滤波后信号的时域波形及频谱，并对滤波前后的信号进行对比，分析信号的变化。

2.5.回放语音信号，分析滤波前后的语音变化。

在熟悉数字信号处理课程理论的基础上，通过MATLAB仿真实现语音信号的采集与处理，进一步加深对数字信号处理理论和技术的掌握。

三．课程设计平台

计算机、MATLAB6.5以上

四．设计原理及计算方法

4.1.语音信号采集

在新版本的MATLAB中可采用audioread函数将声音信号提取出来，并通过sound 函数在MATLAB中直接播放，也可以通过audiowrite函数将处理后的数据存为声音文件。

4.2.噪声的叠加

①采用单一频率的信号进行叠加，产生噪音，此种方法容易滤除。

②采用rand函数产生全频率的噪声，然后通过高通滤波器产生高频率的噪声。其高通滤波器选择的是FIR等波纹逼近法设计的滤波器。

③采用rand函数产生全频率的噪声，然后通过带通滤波器产生高频率的噪声。其带通滤波器选择的是椭圆带通滤波器。

4.3.滤波器的选择

根据声音频率选择合适的滤波器，根据题目要求要使用双线性变换法。

巴特沃斯和切比雪夫是基础的滤波器，因此我们采用这两种滤波器。

五．实验论证方案及结果分析

5.1.声音的采样及加噪声

我们通过声音剪辑软件截取了一段10s的音乐，通过audioread函数将声音采样，并通过绘图程序将声音的时域波形和频域波形显示了出来。

部分代码如下：（详细见附录）

[y,Fs]=audioread('Z:\伦桑,橙光音乐 - 山有木兮-橙光游戏《人鱼传说之长生烛》主题曲.wav'); //读取声音文件

audiowrite('Z:\sound.wav',y,44100) //保存声音文件;

%sound(y,Fs);

y1=y(:,1)';y2=y(:,2)';

l1=length(y1);

xk1=fft(y1,l1);xk2=fft(y2,l1); //FFT变换求频谱

T=1/Fs;

D=1/(l1*T); //频率分辨率 hz

其运行的结果图：

结果分析

由图可见，由于采用的是高清的声音源，其最高的频率达到了15khz。

如要完全还原声音信号，滤波器的通带要在15k以上。

所遇到的问题及解决方法

①在声音采集后出现了两组数列，并且在绘制频谱和时域波形是有两个图形。因此推测此声音为双声道，为了便于观察，将其一声道和二声道分离，从而清晰的看到了结果。

②时域波形显示太过于密集，因此我们调整了绘图的x坐标轴，只显示一部分波形便于以后的对比。

5.2.噪声的产生

我们采用了随机函数rand来产生随机噪声

因为rand产生的是全频率的噪声，为了将声音频率分辨快来，必须加上滤波器来产生高频噪声。我们采用了两种滤波器，一种为FIR滤波器，一种为椭圆滤波器，来避免频率和声音频率相同，导致滤除不掉。

部分代码如下：（详细见附录）

fs=14000; fp=15000;Rp=0.1;As=70; % 滤波器指标

fb=[fs,fp];m=[0,1]; % 计算remezord函数所需参数f,m,dev

dev=[10^(-As/20),(10^(Rp/20)-1)/(10^(Rp/20)+1)];

[n,fo,mo,W]=remezord(fb,m,dev,Fs); % 确定remez函数所需参数

hn=remez(n,fo,mo,W); % 调用remez函数进行设计,用于滤除噪声nt中的低频成分

yt=filter(hn,1,5*nt); %滤除随机噪声中低频成分，生成高通噪声yt

fs1=16000; fp1=18000;fp2=20000;fs2=22000;%模拟指标

Rp=0.1;As=70;%滤波器指标

fp=[fp1,fp2];fs=[fs1,fs2];

wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;%数字指标

OmegaP=2*Fs*tan(wp/2);%频率预畸

OmegaS=2*Fs*tan(ws/2);

%椭圆滤波器

[N,OmegaC]=ellipord(OmegaP,OmegaS,Rp,As,'s');%调用函数求滤波器阶数N和wc

[b,a]=ellip(N,Rp,As,OmegaC,'s');%设计滤波器

[bz,az]=bilinear(b,a,Fs);%双线性变化法

yt=filter(bz,az,5*nt);%滤波

绘形显示：

结果分析

我们设计了两种噪声，如图所示。第一种只能通过低通滤波器来滤除噪声，第二种可以使用低通和带阻滤波器来滤除。

所遇到的问题及解决方法

在加噪声时要注意噪声的长度，否则不满足矩阵的运算规则，产生报错。

5.3.滤波器的设计

根据我们所加的噪声和声音频率，很容易通过计算来得出滤波器的参数fs，和fp。为了比较滤波效果，我们分别设计了巴特沃斯低通滤波器，切比雪夫低通滤波器，巴特沃斯带阻滤波器，切比雪夫带阻滤波器。

部分代码如下：（详细见附录）

fp=12000;fs=15000;%根据声音和噪声频率来确定参数

Fs=44100;

Rp=1;As=70;

wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;

OmegaP=2*Fs*tan(wp/2);%预畸

OmegaS=2*Fs*tan(ws/2);

%切比雪夫Ⅱ型低通滤波器

[N1,OmegaC1]=cheb2ord(OmegaP,OmegaS,Rp,As,'s');%求阶数和wc

[b1,a1]=cheby2(N1,As,OmegaS,'s');%调用函数设计

[bz1,az1]=bilinear(b1,a1,Fs);%双线性变换

[H1,w1]=freqz(bz1,az1);%求频响

dBH1=20*log10(abs(H1)/max(abs(H1)));%转换为分贝

%巴特沃斯低通滤波器

[N2,OmegaC2]=buttord(OmegaP,OmegaS,Rp,As,'s');%求阶数和wc

[b2,a2]=butter(N2,OmegaC2,'s');%调用函数设计

[bz2,az2]=bilinear(b2,a2,Fs);%双线性变换

[H2,w2]=freqz(bz2,az2);

fp=[16000;22000];fs=[18000,20000];%根据声音和噪声频率来确定参数Fs=44100;

Rp=1;As=70;

wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;

OmegaP=2*Fs*tan(wp/2);%预畸

OmegaS=2*Fs*tan(ws/2);

%切比雪夫Ⅱ型带阻滤波器

[N1,OmegaC1]=cheb1ord(OmegaP,OmegaS,Rp,As,'s');%求阶数和wc [b1,a1]=cheby2(N1,As,OmegaC1,'stop','s');

[bz1,az1]=bilinear(b1,a1,Fs);

[H1,w1]=freqz(bz1,az1);

%巴特沃斯带阻滤波器

[N2,OmegaC2]=buttord(OmegaP,OmegaS,Rp,As,'s');

[b2,a2]=butter(N2,OmegaC2,'stop','s');

[bz2,az2]=bilinear(b2,a2,Fs);

[H2,w2]=freqz(bz2,az2);

dBH2=20*log10(abs(H2)/max(abs(H2)));

绘图显示：

结果分析：

根据各种滤波器的频率响应可以看出，各种性能能够满足要求，同时在设计时我们发现，当完成同一滤波器指标时，切比雪夫滤波器的阶数要比巴特沃斯要小。遇到的问题及解决方法

①滤波器无图形，通过各种分析和调试，我们发现当滤波器参数要求过于苛刻时，其所需的阶数会很大，最终超过了上限，系统报错。因此我们更改了合适的过度带，As，和Rp，再次调试无错误。

②采样频率Fs的设定要保持与声音采样频率一致，否则设计的数字滤波器的参数不能与模拟信号一致，导致无法使用。

5.4.滤波效果

调用函数filter函数，对有噪声的信号进行滤波，通过与原声音时域波形和频谱，声音对比，来分析滤波效果。

部分代码如下：（详细代码见附录）

y11=filter(b1,a1,x1);%滤波

y12=filter(b1,a1,x2);

y21=filter(b2,a2,x1);

y22=filter(b2,a2,x2);

yk11=fft(y11);yk12=fft(y12);yk21=fft(y21);yk22=fft(y22);

T=1/Fs;

D=1/(l1*T);

图形如下：

低通滤波器滤波滤除15k以上噪声结果：

带阻滤波器滤除18k~20k噪声结果：

结果分析

对比原声音信号，不管是时域和频域，声音都无明显的改变，证明滤波器设计成功。

遇到的问题及结果分析

在还原声音时域波形时，产生了增益，声音被放大。经过分析代码，发现出错的部分在双线性变换法时进行了归一化处理，相当于乘以一个系数，因此导致了增益产生，处理过之后，无增益，幅值与原信号相同。

六. 结论及心得

本次课程设计顺利完成，总的来说所用到知识皆是课程中所学。在本次课程设计中，从声音的采集，处理，到还原语音信号，用到了数字信号处理课程中的各章节知识，完成了一个系统的学习，将所学的知识进行巩固，时域到频率的变换关系，频域抽样定理，采样定理，滤波器的设计方法等各种知识在本次课程设计中都进行了验证。完善了我们的知识体系。

七. 附录：程序代码及注释

7.1采集声音并加噪声函数

function[xt1,xt2,Fs,l1]=mp3_1(TYPE)

[y,Fs]=audioread('Z:\伦桑,橙光音乐 - 山有木兮-橙光游戏《人鱼传说之长生烛》主题曲.wav');

audiowrite('Z:\sound.wav',y,44100);

%sound(y,Fs);

y1=y(:,1)';y2=y(:,2)';%分离声道

l1=length(y1);%采样点数

xk1=fft(y1,l1);xk2=fft(y2,l1);%fft变换

T=1/Fs;%采样周期

D=1/(l1*T);%频率分辨率

k=floor(0:l1-1);

figure;

subplot(2,2,1);%绘图程序

plot([0:l1-1]/Fs,y1);title('1声道波形');xlabel('t');ylabel('y1'); axis([0,l1/Fs/100,min(y1)-0.5,max(y1)+0.5]);

subplot(2,2,2);

plot([0:l1-1]/Fs,y2);title('2声道波形');xlabel('t');ylabel('y2'); axis([0,l1/Fs/100,min(y2)-0.5,max(y2)+0.5]);

subplot(2,2,3);

plot(k*D,abs(xk1)/max(abs(xk1)));title('1声道波形频谱

');xlabel('f/hz');ylabel('yk1');

axis([0,l1/2*D,-0.1,1.2]);

subplot(2,2,4);

plot(k*D,abs(xk2)/max(abs(xk2)));title('2声道波形频谱

');xlabel('f/hz');ylabel('yk2');

axis([0,l1/2*D,-0.1,1.2]);

nt=2*rand(1,l1); %产生随机噪声nt

switch TYPE

case 1

%产生15k以上频率的噪声

%等波纹逼近法FIR高通滤波器

%用于滤除噪声nt中的低频成分,生成噪声

fs=14000; fp=15000;Rp=0.1;As=70; % 滤波器指标

fb=[fs,fp];m=[0,1]; % 计算remezord函数所需参数f,m,dev

dev=[10^(-As/20),(10^(Rp/20)-1)/(10^(Rp/20)+1)];

[n,fo,mo,W]=remezord(fb,m,dev,Fs); % 确定remez函数所需参数

hn=remez(n,fo,mo,W); % 调用remez函数进行设计,用于滤除噪声nt中的低频成分

yt=filter(hn,1,5*nt); %滤除随机噪声中低频成分，生成高通噪声yt

%------------------------绘图-----------------------------------------------

figure;

[H,w]=freqz(hn,1);%频率响应

dBH1=20*log10(abs(H)/max(abs(H)));%转换DB显示

plot(w/pi,dBH1);%绘图

xlabel('\omega/\pi');%坐标标签

ylabel('dB');title('等波纹逼近FIR滤波器');