数字信号处理课程设计.doc

数字信号处理课程设计.doc

《数字信号处理课程设计.doc》由会员分享，可在线阅读，更多相关《数字信号处理课程设计.doc（26页珍藏版）》请在冰豆网上搜索。

语音信号的数字滤波处理

语音信号的数字滤波处理

摘要

数字信号处理的目的是对真实世界的连续模拟信号进行测量或滤波。

数字信号处理的核心算法是离散傅立叶变换（DFT），是DFT使信号在数字域和频域都实现了离散化，从而可以用通用计算机处理离散信号。

而使数字信号处理从理论走向实用的是快速傅立叶变换（FFT），FFT的出现大大减少了DFT的运算量，使实时的数字信号处理成为可能、极大促进了该学科的发展。

数字信号处理课程设计与《数字信号处理》课程配套，是电子、通信等专业的重要实践环节。

数字信号处理是每一个电子信息科学工作者必须掌握的重要知识。

它采用计算机仿真软件，以数值计算的方法对信号进行分析、变换、滤波、检测、估计与识别等加工处理，以达到提取信息便于使用的目的。

本课程设计旨在运用MATLAB软件对线性卷积和循环卷积进行动态演示，以及完成采样定理，采样时产生的混叠效应的演示，最后通过语音采集，完成滤波器设计及其应用。

关键字：

信号；滤波；采样；混叠效应

目录

1设计目的与内容 1

1.1设计目的 1

1.2设计内容 1

2基本原理 2

2.1线性卷积和循环卷积 2

2.1.1线性卷积 2

2.1.2循环卷积 2

2.2采样定理 2

2.3数字滤波器 3

2.3.1巴特沃斯滤波器 3

2.3.2布莱克曼窗 3

2.3.3数字滤波器的指标 3

3程序的运行结果 4

3.1线性卷积与循环卷积 4

3.1.1线性卷积 4

3.1.2循环卷积 5

3.2采样定理 6

3.3数字滤波器 7

3.3.1总体说明 7

3.3.2巴特沃斯滤波器 8

3.3.3布莱克曼窗 10

4心得体会 12

参考文献 13

附录 14

附录A线性卷积的动态演示 14

附录B循环卷积 15

附录C采样定理 16

附录D巴特沃斯低通 17

附录E巴特沃斯高通 19

附录F巴特沃斯带通 20

附录G布莱克曼窗低通 21

附录H布莱克曼窗高通 23

附录I布莱克曼窗带通 24

1设计目的与内容

1.1设计目的

电子信息工程专业的培养目标是具备电子技术的基本理论和应用技术，能从事电子、信息、通信、电信等领域的工作，具有高素质、宽口径、创新晋升的专业人才。

对本专业学生的培养要进行工程素质培养、拓宽专业口径、注重基础和发展潜力。

特别是培养学生的创新能力，以实现技术为主线多进行实验技能的培养。

学生通过数字信号处理课程设计这一重要环节，可以将本专业的主干课程《数字信号处理》从理论学习到实践应用，对数字信号处理技术有较深的了解，进一步增强学生动手能力和适应实际工作的能力。

本次设计旨在综合运用《数字信号处理》课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。

1.2设计内容

本课程设计的内容分为两个部分：

第一个部分为预习题，它包括两个任务，一个是设计卷积运算的演示程序，一个是编写程序演示采样定理（时域采样、频谱周期延拓），同时演示采样频率小于2fc时，产生的混叠效应。

第二个部分是设计题，它是针对滤波器的设计及其应用，包括两个主要内容，一个是用巴特沃斯模拟滤波器（低通、带通、高通）对有噪声信号的语音进行滤波，一个是用布莱克曼窗（低通、带通、高通）对有噪声信号的语音进行滤波。

2基本原理

2.1线性卷积和循环卷积

2.1.1线性卷积

在实际应用中，为了分析时域离散线性非移变系统或者对序列进行滤波处理等，需要计算两个序列的线性卷积。

线性卷积既可以在时域中直接计算，也可以通过变换在频域中计算得到。

假设h（n）和x（n）都是有限长序列，长度分别为N和M，它们的线性卷积可以表示如下：

2.1.2循环卷积

为了提高线性卷积的速度，希望用DFT（FFT）计算线性卷积。

从而引入循环卷积来运用DFT快速计算线性卷积。

循环卷积运用到离散傅立叶变换的循环移位性质，即时域循环移位定理。

假设h（n）和x（n）都是有限长序列，长度分别为N和M，它们的循环卷积可以表示如下：

2.2采样定理

对连续信号进行等间隔采样形成采样信号，采样信号的频谱是原连续信号的频谱以采样频率为周期进行周期性的拓延形成的。

设连续信号xa（t）属带限信号，最高截止频率为Ωc，如果采样角频率Ωs≥2Ωc，那么让采样信号x^a（t）通过一个增益为T，截止频率为Ωs/2=π/T的理想低通滤波器，可以唯一地恢复出院连续信号xa（t）。

否则Ωs＜2Ωc会造成采样信号中频谱混叠现象，不可能无失真地恢复原连续信号。

2.3数字滤波器

2.3.1巴特沃斯滤波器 

　　从幅频特性提出要求，而不考虑相频特性。

巴特沃斯滤波器具有最大平坦幅度特性，其幅频响应表达式为：

2.3.2布莱克曼窗

   增加一个二次谐波余弦分量，可进一步降低旁瓣，但主瓣宽度进一步增加，增加N可减少过渡带。

频谱的幅度函数为：

+0.04

2.3.3数字滤波器的指标

滤波器的频率参数主要有：

①通带截频为通带与过渡带的边界点，在该点信号增益下降到规定的下限；②阻带截频为阻带与过渡带的边界点，在该点信号衰耗下降到规定的下限；③转折频率为信号功率衰减到1/2（约3dB）时的频率，在很多情况下，也常以fc作为通带或阻带截频；④固有频率当电路没有损耗时，就是其谐振频率，复杂电路往往有多个固有频率。

3程序的运行结果

3.1线性卷积与循环卷积

3.1.1线性卷积

输入任意两个序列x1s、h1s，指定x1s为自己的学号，例如：

x1s=[200784250137]，h1s的内容和长度自选。

本次实验h1s=[123456]。

然后用MATLAB语言编写程序，实现线性卷积的动态演示，演示结果如图3.1，3.2所示。

图3.1线性卷积过程中

图3.2线性卷积结果图

3.1.2循环卷积

输入任意两个序列x1、x2，指定x1为自己的学号，例如：

x1=[2,0,0,7,8,4,2,5,0,1,3,7]，x2的内容和长度自选。

本次实验x2=[1,2,3,4]。

然后用MATLAB语言编写程序，实现循环卷积的演示，演示结果如图3.3所示。

图3.3循环卷积结果图

3.2采样定理

采样信号的频谱是原连续信号的频谱以采样频率为周期进行周期性的拓延形成的。

对下面连续信号进行采样：

A=n，a=Ωo=π，n=学号，A为幅度因子，a为衰减因子，为模拟角频率，其中n为学号（例如：

李浩敏同学n=37）。

要求输入采样频率fs（根据程序处理需要指定范围）后，在时域演示信号波形、采样脉冲及采样后信号；在频域演示不同采样频率下对应信号的频谱。

演示结果如图3.4所示。

图3.4采样定理结果图

在用MATLAB实现采样定理时，调用了时域序列绘图函数tstem（），其程序代码如下：

functiontstem（xn,yn）

%时域序列绘图函数

%xn:

信号数据序列，yn:

绘图信号的纵坐标名称（字符串）

n=0:

length（xn）-1;

stem（n,xn,'*-'）;boxon

xlabel（'n'）;ylabel（yn）;

axis（[0,n（end）+2,min（xn）,1.2*max（xn）]）

3.3数字滤波器

3.3.1总体说明

语音采集，采用Windows自带的声音文件（默认为22050Hz，16位），进行语音信号的采集（*.wav），时间控制在几秒左右。

然后在MATLAB软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。

通过wavread函数的使用，要求理解采样频率、采样位数等概念。

wavread函数调用格式：

y=wavread（file），读取file所规定的wav文件，返回采样值放在向量y中。

[y,fs,nbits]=wavread（file），采样值放在向量y中，fs表示采样频率（Hz），nbits表示采样位数。

语音信号的频谱分析，要求首先画出语音信号的时域波形，然后对语音信号进行频谱分析，在MATLAB中，可以利用函数fft对信号进行快速付立叶变换，得到信号的频谱特性，从而加深对频谱特性的理解。

回放语音信号，在Matlab中，函数sound可以对声音进行回放。

其调用格式：

sound（x,fs,bits），可以感觉滤波前后的声音有变化。

3.3.2巴特沃斯滤波器

巴特沃斯低通滤波器的性能指标：

fp1=1000;fs1=1500;wp1=2*pi*fp1;ws1=2*pi*fs1;Rp1=1;As1=60;其结果如图3.5所示。

图3.5巴特沃斯低通滤波器

巴特沃斯高通滤波器的性能指标：

fp1=1500;fs1=1000;wp1=2*pi*fp1;ws1=2*pi*fs1;Rp1=1;As1=60;其结果如图3.6所示。

图3.6巴特沃斯高通滤波器

巴特沃斯带通滤波器的性能指标：

fp1=1500;fp2=2500;fs1=1000;fs2=3000;wp1=2*pi*[fp1,fp2];ws1=2*pi*[fs1,fs2];Rp1=1;As1=40;其结果如图3.7所示。

图3.7巴特沃斯带通滤波器

3.3.3布莱克曼窗

布莱克曼窗低通滤波器的性能指标：

fp1=1000;fs1=1200;wp1=fp1/4000;ws1=fs1/4000;Rp1=1;As1=100;其结果如图3.8所示。

图3.8布莱克曼窗低通滤波器

布莱克曼窗高通滤波器的性能指标：

fp2=500;fs2=300;wp2=fp2/4000;ws2=fs2/4000;Rp2=1;As2=100;其结果如图3.9所示。

图3.9布莱克曼窗高通滤波器

布莱克曼窗带通滤波器的性能指标：

fpl3=500;fpu3=2200;fsl3=350;fsu3=2450;wpl3=fpl3/4000;wpu3=fpu3/4000;

wsl3=fsl3/4000;wsu3=fsu3/4000;Rp3=1;As3=100;其结果如图3.10所示。

图3.10布莱克曼窗带通滤波器

4心得体会

这次的课程设计比较和以往的不一样，不一样就在，设计的工具与课程是分离的，就是一门课程的设计需要另外的一门课程做辅助，现在的辅助工具是MATLAB，由于以前是选修的MATLAB程序设计这个课程，刚开始拿到这个设计，心里没底，好象有点没从下手，但是，就是在设计的过程中让我更好的学习，对它也有了比较深的体会。

做完这个设计，和同学也聊过，有同学认为，这个过程中走了不少的弯路。

但是我认为，每个软件，每个语言对于同一结果的实验都是有不同的方案可以选择，不同的路径可以走，那是毫无疑问的。

所以对于初学者这也是必然的，捷径也都是靠自己摸索出来的，我所说的更深的体会也就体现在这个环节上的。

即便是今天别人告诉了你或者是书上看到了有这个捷径可以走，但，心里还是有一种“虚”，我认为这，不好！

通过本次课程设计，使我们对数字信号处理相关知识有了更深刻的理解，尤其是对各个滤波器的设计，并且也对Matlab这个软件的使用有了进一步的掌握。

综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，利用MATLAB作为工具进行实现，复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，从实践上初步实现对数字信号的处理。

做课程设计，最大的困难是MATLAB代码学习和查找，和不断的调试代码。

不过在这过程中我学到了很多，虽说调试代码总是反反复复很累也很烦，但使我们学到了课本上无法学到的知识，同时也提高了我们学习的兴趣，使学习变的不再枯燥无味。

参考文献

[1]丁玉美，高西全．数字信号处理（第二版）[M]．西安电子科技大学出版社．

[2]王创新，文卉．数字信号处理试验指导书[M]．长沙理工大学印刷．

[3]陈怀琛等译．数字信号处理及其MATLAB实现[M]．电子工业出版社．

[4]陈怀琛等．MATLAB及在电子信息课程中的应用[M]．电子工业出版社．

[5]A.V.奥本海姆，R.W.谢弗著．数字信号处理[M]．北京：

科学出版社．

[6]胡广书编著．数字信号处理——理论、算法与实现（第二版）[M]．北京：

电子工业出版社，2003．

附录

附录A线性卷积的动态演示

-2-

%线性卷积的动态演示

clear;

x1s=[200784250137];

lx=length（x1s）;

h1s=[123456];

lh=length（h1s）;

lmax=max（lx,lh）;

iflx>lhnx=0;nh=lx-lh;

elseiflx

elsenx=0;nh=0;

end

dt=0.5;

lt=lmax;

x=[zeros（1,lt）,x1s,zeros（1,nx）,zeros（1,lt）];

t1=（-lt+1:

2*lt）*dt;

h=[zeros（1,2*lt）,h1s,zeros（1,nh）];

hf=fliplr（h）;

y=zeros（1,3*lt）;

fork=0:

2*lt

p=[zeros（1,k）,hf（1:

end-k）];

y1=x.*p;

yk=sum（y1）;

y（k+lt+1）=yk;

figure

（1）;

subplot（2,2,1）;stem（t1,x）

axis（[-lt*dt,2*lt*dt,min（x）,max（x）]）,holdon;

ylabel（'x（t）'）;

subplot（2,2,2）;stairs（t1,p）

axis（[-lt*dt,2*lt*dt,min（p）,max（p）]）,

ylabel（'h（k-t）'）;

subplot（2,2,3）;stem（t1,y1）

axis（[-lt*dt,2*lt*dt,min（y1）,max（y1）+eps]）,

ylabel（'s=x*h（k-t）'）;

subplot（2,2,4）;stem（k*dt,yk）

axis（[-lt*dt,2*lt*dt,floor（min（y）+eps）,ceil（max（y）+eps）]）,

holdon,ylabel（'y（k）=sum（s）*dt）'）;

ifk==round（0.8*lt）disp（'pausepressanykeytocountinue'）,pause

elsepause

（1）;

end

end

figure

（2）;%结果比较

x=[200784250137];

h=[123456];

y=conv（x,h）;

stem（y）;

axis（[1,20,0,120]）;

ylabel（'y（n）=x（n）*h（n）'）;

附录B循环卷积

%循环卷积

x1=[2,0,0,7,8,4,2,5,0,1,3,7];

x2=[1,2,3,4];

n1=length（x1）;

n2=length（x2）;

ifn1>n2

x22=zeros（n1）;

x2（1,n1）=0;

fori=0:

n1-1

x22（:

i+1）=（circshift（x1,[-1,i]））';

end

y=x22*x2';

disp（y）;

stem（y）;

elseifn1==n2

x22=zeros（n1）;

fori=0:

n1-1

x22（:

i+1）=（circshift（x1,[-1,i]））';

end

y=x22*x2';

disp（y）;

stem（y）;

elseifn1

x11=zeros（n2）;

x1（1,n2）=0;

fori=0:

n2-1

x11（:

i+1）=（circshift（x2,[-1,i]））';

end

y=x11*x2';

disp（y）;

stem（y）;

end

附录C采样定理

-16-

%时域采样理论验证

Tp=128/1500;

%产生M长采样序列x（n）

%Fs=1500;T=1/Fs;

Fs=1500;T=1/Fs;

M=Tp*Fs;n=0:

M-1;

A=37;alph=pi*37*sqrt

（2）;omega=pi*37*sqrt

（2）;

xnt=A*exp（-alph*n*T）.*sin（omega*n*T）;

Xk=T*fft（xnt,M）;%M点FFT[xnt）]

yn='xa（nT）';

X=fft（xnt）;

k=0:

M-1;fk=k/Tp;

figure

（1）;

subplot（2,1,1）;stem（xnt）;title（'波形'）;

subplot（2,1,2）;plot（fk,abs（X）,'r'）;title（'幅频'）;

figure

（2）;

subplot（3,2,1）;tstem（xnt,yn）;boxon;title（'（a）Fs=1500Hz'）;

k=0:

M-1;fk=k/Tp;

subplot（3,2,2）;plot（fk,abs（Xk）,'r'）;title（'（a）T*FFT[xa（nT）],Fs=1500Hz'）;

xlabel（'f（Hz）'）;ylabel（'幅频'）;axis（[0,Fs,0,1.2*max（abs（Xk））]）

%Fs=300;T=1/Fs;

Fs=300;T=1/Fs;

M=ceil（Tp*Fs）;n=0:

M-1;

xnt=A*exp（-alph*n*T）.*sin（omega*n*T）;

Xk=T*fft（xnt,M）;%M点FFT[xnt）]

yn='xa（nT）';

subplot（3,2,3）;tstem（xnt,yn）;boxon;title（'（b）Fs=300Hz'）;

k=0:

M-1;fk=k/Tp;

subplot（3,2,4）;plot（fk,abs（Xk）,'r'）;title（'（b）T*FFT[xa（nT）,Fs=300Hz'）;

xlabel（'f（Hz）'）;ylabel（'幅频'）;axis（[0,Fs,0,1.2*max（abs（Xk））]）

%Fs=200;T=1/Fs;

Fs=200;T=1/Fs;

M=ceil（Tp*Fs）;n=0:

M-1;

xnt=A*exp（-alph*n*T）.*sin（omega*n*T）;

Xk=T*fft（xnt,M）;%M点FFT[xnt）]

yn='xa（nT）';subplot（3,2,5）;tstem（xnt,yn）;boxon;title（'（c）Fs=200Hz'）;

k=0:

M-1;fk=k/Tp;

subplot（3,2,6）;plot（fk,abs（Xk）,'r'）;title（'（c）T*FFT[xa（nT）,Fs=200Hz'）;

xlabel（'f（Hz）'）;ylabel（'幅频'）;axis（[0,Fs,0,1.2*max（abs（Xk））]）;

附录D巴特沃斯低通

%专业：

电子信息工程班级：

电信0701班学号：

200784250137姓名：

李浩敏

%采用Windows自带的声音文件（默认为22050Hz，16位），进行语音信号的采集（*.wav）

clf;closeall;

clearall;

[x,Fs,bits]=wavread（'E:

\LHM\lihaomin.wav'）;%此声音文件的时间为1s

x=x（:

1）;

m=length（x）;

sound（x,Fs,bits）;

y=fft（x,m）;%对声音信号进行傅立叶变换

f=（Fs/m）*[1:

m];

t=[1:

m]/Fs;

subplot（4,2,1）;plot（t,x）;title（'原始信号波形李浩敏'）;xlabel（'time（s）'）;

axis（[0,1.8,-0.5,0.5]）;

subplot（4,2,2）;plot（f,abs（y））;title（'原始信号频谱李浩敏'）;xlabel（'frequency（hz）'）;

axis（[0,11000,0,300]）;

%噪声信号

noise1=0.01*sin（2*pi*3000*t）;%高频噪声

%巴特沃斯低通滤波器

fp1=1000;fs1=1500;

wp1=2*pi*fp1;ws1=2*pi*fs1;Rp1=1;As1=60;%设置滤波器参数

X1=x+noise1';%加了高频噪声的信号

sound（X1,Fs,bits）;%回放噪声信号

Y11=fft（X1,m）;%对加了高频噪声的信号进行傅立叶变换

subplot（4,2,3）;plot（t,X1）;title（'加了高频噪声的信号时域图'）;xlabel（'time（s）'）;

axis（[0,1.8,-0.5,0.5]）;

subplot（4,2,4）;plot（f,abs（Y11））;title（'加了高频噪声的信号频谱图'）;xlabel（'frequency（hz）'）;

axis（[0,11000,0,300]）;

[N,wc]=buttord（wp1,ws1,Rp1,As1,'s'）;%计算滤波器阶数N和3dB截止频率

[B,A]=butter（N,wc,'s'）;%计算滤波器系统函数分子分母多项式系数

[Bz,Az]=bilinear（B,A,Fs）;%impinvar（B,A,Fs）;

k=0:

511;fk=0:

8000/512:

8000;wk=2*pi*fk;

Hk1=freqs（B,A,wk）;

subplot（4,2,5）;plot（fk/1000,20*log10（abs（Hk1）））;title（'低通滤波器幅频'）;xlabel（'Hz'）;

subplot（4,2,6）;plot（fk/1000,unwrap（angle（Hk1）））;title（'低通滤波器相频'）;xlabel（'Hz'）;

x1get=filter（Bz,Az,X1）;

Y1get=fft（x1get,m）;

sound（x1get,Fs,bits）;

subplot（4,2,7）;plot（t,