MATLAB课程设计---MATLAB语音信号的采集与处理.doc

课 程 设 计 报 告设计名称：MATLAB课程设计系别：学生姓名：班级：学号：成绩：指导教师：开课时间: 2010 至 2011 学年第 2 学期 MATLAB语音信号的采集与处理一、实践的目的和要求1.MATLAB软件功能简介 MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。 MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。 MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，使MATLAB成为一个强大的数学软件。MATLAB特点：1) 高效的数值计算及符号计算功能,能使用户从繁杂的数学运算分析中解脱出来; 2) 具有完备的图形处理功能,实现计算结果和编程的可视化; 3) 友好的用户界面及接近数学表达式的自然化语言,使学者易于学习和掌握; 4) 功能丰富的应用工具箱,为用户提供了大量方便实用的处理工具。MATLAB与其他计算机语言相比，它有很多的优点： 友好的工作平台和编程环境 简单易用的程序语言，运算符丰富，语言简洁紧凑，使用方便灵活，库函数极其丰富。 强大的科学计算机数据处理能力，MATLAB既具有结构化的控制语句（如for循环、while循环、break语句和if语句），又有面向对象编程的特性。 应用广泛的模块集合工具箱 语法限制不严格，程序设计自由度大。 程序的可移植性很好，基本上不做修改就可以在各种型号的计算机和操作系统上运行。 出色的图形处理功能 实用的程序接口和发布平台2.本题目的意义（1）有利于基础知识的理解通过MATLAB及在电子信息类（课程中的应用）的学习，我掌握了MATLAB的一些基本知识和用法，对我以后的学习带来方便。（2）有利于逻辑思维的锻炼程序设计是公认的、最能直接有效地训练学生的创新思维，培养分析问题、解决问题能力的学科之一。即使一个简单的程序，从任务分析、确定算法、界面布局、编写代码到调试运行，整个过程学生都需要有条理地构思，这中间有猜测设想、判断推理的抽象思维训练，也有分析问题、解决问题、预测目标等能力的培养。（3）有利于与其他学科的整合在程序设计中，我们可以解决其它学科有关问题，也利用其它课程的有关知识来解决MATLAB中比较抽象很难理解的知识。在MATLAB课中整合其它学科的知识，发挥MATLAB的优势。（4）有利于治学态度的培养。课程设计中，语句的语法和常量变量的定义都有严格的要求，有时输了一个中文标点、打错了一个字母，程序无法正常运行。因此，程序设计初学阶段，我们经常会犯这样的错误，可能要通过几次乃至十多次的反复修改、才能成功，但这种现象会随着学习的深入而慢慢改观。这当中就有一个严谨治学、一丝不苟的科学精神的培养，又有一个不怕失败、百折不挠品格的锻炼。二、实践内容 1.常用信号的实现 通过编写MATLAB软件仿真程序，掌握MATLAB软件中常用信号的表示方法，绘图函数的使用。内容如下：（1）单位阶跃信号f(t)= 和单位阶跃序列的实现（2）单位冲激信号和单位脉冲序列的实现（3）单边衰减指数信号和实指数序列的实现2.信号的时域基本运算 通过调用MATLAB中的相关函数，编写信号的基本运算仿真程序，实现信号的时域基本运算，达到掌握MATLAB相关运算函数的调用格式及功能的目的。 （1）信号的相加，相乘和移位运算 （2）序列的卷积运算3.系统分析 通过调用MATLAB相关函数，实现对LTI系统的时域和频域进行分析，熟练掌握LTI系统单位冲激响应和零状态响应得MATLAB求解方法。 （1）连续LTI系统的时域分析 （2）连续LTI系统的频域分析4.综合设计 在熟练掌握MATLAB基本的信号处理方法的基础上，按照要求设计特定功能的综合系统，掌握MATLAB信号处理系统的仿真程序设计。综合系统题目（1）MATLAB模拟滤波器的实现（2）MATLAB语音信号的采集与处理三、实践仿真程序设计及结果分析1.常用信号的实现 （1）单位阶跃信号f(t)= 和单位阶跃序列的实现单位阶跃信号f(t)= 的实现程序：t0=0;t1=-1;t2=3;dt=0.01;t=t1:dt:-t0;n=length(t);t3=-t0:dt:t2;n3=length(t3);u=zeros(1,n);u3=ones(1,n3);plot(t,u);hold on;plot(t3,u3);plot(-t0,-t0,0,1);hold off;axis(t1,t2,-0.2,1.5); xlabel(时间(t);ylabel(幅值(f);title(单位阶跃信号); 时域波形：单位阶跃序列的实现程序： k0=0;%单位阶跃出现的位置k1=-3;k2=6; k=k1:k0-1; n=length(k); k3=-k0:k2; n3=length(k3); u=zeros(1,n); u3=ones(1,n3); stem(k,u,filled); hold on; stem(k3,u3,filled); hold off; axis(k1,k2,-0.2,1.5); title(单位阶跃序列); xlabel(时间(k);ylabel(幅值f(k);时域波形：（2）单位冲激信号和单位脉冲序列的实现单位冲激信号的实现程序：clear;t0=0;t1=-1;t2=5;dt=0.1;t=t1:dt:t2;n=length(t);x=zeros(1,n);x(1,(t0-t1)/dt+1)=1/dt;stairs(t,x);axis(t1,t2,0,1/dt);xlabel(时间(t);ylabel(幅值(f);title(单位冲激信号);时域波形：单位脉冲序列的实现程序： k1=-3; k2=6; k=k1:k2; n=3;%单位脉冲出现的位置 f=(k-n)=0; stem(k,f,filled);title(单位脉冲序列) xlabel(时间(k);ylabel(幅值f(k);时域波形：（3）单边衰减指数信号和实指数序列的实现单边衰减指数信号的实现程序：clear t1=-1;t2=10;dt=0.1;t=t1: dt: t2;A1=1;a1=0.5;n=A1*exp(-a1*t);plot(t,n);axis(t1,t2,0,1);xlabel(时间(t);ylabel(幅值(f);title(单边衰减指数信号);时域波形：实指数序列的实现程序：clf;k1=-1;k2=10;k=k1:k2;a=-0.6;A=1;f=A*a.k;stem(k,f,filled);title();xlabel(k);ylabel(f(k);时域波形：2.信号的时域基本运算 （1）信号的相加，相乘和移位运算 实现两个连续信号的相加，即。 程序 clear all; t=0:0.0001:3; b=3; t0=1; u=stepfun(t,t0); n=length(t); for i=1:n u(i)=b*u(i)*(t(i)-t0);end %产生一个斜坡信号 y=sin(2*pi*t); %产生一个正弦信号 f=y+u; %信号相加 plot(t,f); xlabel(时间(t);ylabel(幅值f(t);title(连续信号的相加); 时域波形：实现两个连续信号的相乘，即。程序 clear all;t=0:0.0001:5;b=3;t0=1;u=stepfun(t,t0);n=length(t);for i=1:n u(i)=b*u(i)*(t(i)-t0);endy=sin(2*pi*t); f=y.*u;plot(t,f);xlabel(时间(t);ylabel(幅值f(t);title(连续信号的相乘) 时域波形：实现两个连续信号的移位，即，常数程序： clear all; t=0:0.0001:2;y=sin(2*pi*(t); yl=sin(2*pi*(t-0.2); plot(t,y,-,t,yl,-); ylabel(幅值f(t);xlabel(时间（t）);title(信号的位移);时域波形：（2）序列的卷积运算实现两个连续信号的卷积运算， 程序：%连续函数卷积计算a=1000;t1=-5:1/a:5;f1=stepfun(t1,0);f2=stepfun(t1,-1/a)-stepfun(t1,1/a);subplot(231);plot(t1,f1);axis(-5,5,0,1.2); %xlabel(时间(t);ylabel(f1(t);title(单位阶跃函数);subplot(232);plot(t1,f2);ylabel(f2(t);title(单位冲激函数);y=conv(f1,f2);r=2*length(t1)-1;t=-10:1/a:10;subplot(233);plot(t,y);axis(-5,5,0,1.2);title(f1与f2的卷积);ylabel(y(t);f11=conv(f1,f1);f22=conv(f2,f2);subplot(234);plot(t,f11);title(f1与f1的卷积);ylabel(f11(t);axis(-5,5,0,5000);subplot(235);plot(t,f22);title(f2与f2的卷积);ylabel(f22(t);时域波形：3.系统分析 （1）连续LTI系统的时域分析 求系统y(2)(t)+6y(1)(t)+8y(t)=3x(1)(t)+9x(t)的冲激响应和阶跃响应。i 冲激响应MATLAB程序：clear all;a=1,6,8;b=3,9;sys=tf(b,a);t=0:0.1:10;y=impulse(sys,t);plot(t,y);xlabel(时间(t);ylabel(y(t);title(单位冲激响应);系统的冲激响应图：ii 阶跃响应MATLAB程序：clear all;a=1,6,8;b=3,9;sys=tf(b,a);t=0:0.1:10;y=step(sys,t);plot(t,y);xlabel(时间(t);ylabel(y(t);title(单位阶跃响应);系统的阶跃响应图： 求系统y(2)(t)+y(t)=costu(t)，y(0+)=y(1)(0+)=0的全响应。i 系统在正弦激励下的零状态响应MATLAB程序b=1;a=1 0 1;sys=tf(b,a);t=0:0.1:10;x=cos(t);y=lsim(sys,x,t);plot(t,y);xlabel(时间(t);ylabel(y(t);title(零状态响应);系统的零状态响应如图所示：ii求系统的全响应MATLAB程序：b=1;a=1 0 1;A B C D=tf2ss(b,a);sys=ss(A,B,C,D);t=0:0.1:10;x=cos(t);zi=-1 0;y=lsim(sys,x,t,zi);plot(t,y);xlabel(时间(t);ylabel(y(t);title(系统的全响应);系统的全响应如图所示： 已知某LTI系统的激励为f1=sint(t)，单位冲激响应为h(t)=te-2t(t)， 试给出系统零状态响应yf(t)的数学表达式。MATLAB程序：p=0.1; k1=0:p:10; f1=sin(k1); k2=k1; f2=1*exp(-2*k1); f=conv(f1,f2); f=f*p; k0=k1(1)+k2(1); k3=length(f1)+length(f2)-2; k=k0:p:(k3*p+k0); subplot(2,2,1); plot(k1,f1); title(f1(t); xlabel(t),ylabel(f1(t); subplot(2,2,2); plot(k2,f2); title(f2(t); xlabel(t),ylabel(f2(t); subplot(2,2,3); plot(k,f); h=get(gca,position); h(3)=2.5*h(3); set(gca,position,h); title(f(t)=f1(t)*f2(t); xlabel(t),ylabel(f(t);系统的零状态响应如图所示 程序设计实验计算下述系统在指数函数激励下的零状态响应。H(s)= MATLAB程序：a=1 0.996 463 97.8 12131 8.11 0;b=1.65 -0.331 -576 90.6 19080;sys=tf(b,a);t=0:0.1:10;x=exp(-t);y=lsim(sys,x,t);plot(t,y);xlabel(时间(t);ylabel(y(t);title(零状态响应)系统的响应如图所示：（2）连续LTI系统的频域分析 编程实现信号的傅里叶变换和傅里叶逆变换I傅里叶变换。i 已知连续时间信号f(t)=e-2|t|，通过程序完成信号f(t)的傅里叶变换MATLAB程序：syms t x=exp(-2*abs(t);F=fourier(x);subplot(2,1,1);ezplot(x);subplot(2,1,2);ezplot(F);信号f(t)的傅里叶变换如图所示：ii 试画出信号f(t)= 的波形及其幅频特性曲线。MATLAB程序: syms t ;f=2/3*exp(-3*t)*sym(heaviside(t) ;F=fourier(f) ;subplot(2,1,1) ;ezplot(f) ;subplot(2,1,2) ;ezplot(abs(F) ;信号f(t)= 的波形及其幅频特性曲线如图所示II 傅里叶逆变换。已知f(j)= ，求信号F(j)的逆傅里叶变换。 MATLAB程序R=0.02;t=-2:R:2;syms t wifourier(1/(1+w2),t)运行结果：ans =1/2*exp(-t)*heaviside(t)+1/2*exp(t)*heaviside(-t)III 傅里叶变换数值计算。已知门函数f(t)=g2(t)=(t+1)-(t-1)，试采用数值计算方法确定信号的傅里叶变换F(j)。 MATLAB程序：R=0.02 ; %频谱线宽度t=-2 :R :2 ;f=stepfun(t,-1)-stepfun(t,1) ;W1=2*pi*5 ; %频率宽度N=500 ; %采样数为Nk=0 :N ;W=k*W1/N ; %W频率正半周的采样点F=f*exp(-j*t*W)*R ; %求F（jw）F=real(F) ; %求F（jw）的实部W=-fliplr(W),W(2 :501) ; %频谱图频率宽度左右对称F=fliplr(F),F(2 :501) ; %频谱图幅值左右对称subplot(2,1,1) ;plot(t,f) ;xlabel(t) ;ylabel(f(t) ;axis(-2,2,-0.5,2) ;title(f(t)= (t+1)-(t-1) ;subplot(2,1,2) ;plot(W,F) ;title(f(t)傅氏变换) ;信号的傅里叶变换如图所示：IV 傅里叶变换的时移特性分别绘出信号f(t)= 信号f(t-1)的频谱图，并观察信号时移对信号频谱的影响。if(t)= 的频谱。MATLAB程序：r=0.02;t=-5:r:5;N=200;W=2*pi;k=-N:N;w=k*W/N;f1=1/2*exp(-2*t).*stepfun(t,0);F=r*f1*exp(-j*t*w);F1=abs(F);P1=angle(F);subplot(3,1,1);plot(t,f1);gridxlabel(t);ylabel(f(t);title(f(t);subplot(3,1,2);plot(w,F1);xlabel(w);grid;ylabel(F(jw);subplot(3,1,3)plot(w,P1*180/pi);grid;xlabel(w);ylabel(相位（度）);信号f(t)= 傅里叶变换的时频特性如图所示：ii f(t-1)的频谱。MATLAB程序：r=0.02;t=-5:r:5;N=200;W=2*pi;k=-N:N;w=k*W/N;f1=1/2*exp(-2*(t-1).*stepfun(t,1); %f(t)F=r*f1*exp(-j*t*w); %f(t)的傅里叶变换F1=abs(F);P1=angle(F);subplot(3,1,1);plot(t,f1);hold on;gridxlabel(t);ylabel(f(t);title(f(t-1);subplot(3,1,2);plot(w,F1);hold on;xlabel(w);ylabel(F(jw);grid;subplot(3,1,3);plot(w,P1*180/pi);hold on;grid;xlabel(w);ylabel(相位(度); 傅里叶变换的时移特性如图所示：V 傅里叶变换的频移特性信号f(t)=g2(t)为门函数，试绘出信号f1(t)=f(t) 以及信号f2(t)=f(t) 的频谱图，并与原信号频谱图进行比较。MATLAB程序：R=0.02;t=-2:R:2;f=stepfun(t,-1)-stepfun(t,1);f1=f.*exp(-j*10*t);f2=f.*exp(j*10*t);W1=2*pi*5;N=500;k=-N:N;W=k*W1/N;F1=f1*exp(-j*t*W)*R; %f1(t)傅里叶变换F2=f2*exp(-j*t*W)*R; %f2(t)傅里叶变换F1=real(F1);F2=real(F2);subplot(2,1,1);plot(W,F1);hold on;xlabel(w);ylabel(F1(jw);title(频谱F1(jw);subplot(2,1,2);plot(W,F2); hold on;xlabel(w);ylabel(F2(jw);title(频谱F2(jw);傅里叶变换的频移特性如图所示： 程序设计实验I 试确定下列信号的傅里叶变换的数学表达式。 (a) f(t)=(t+1)-(t-1) MATLAB程序： R=0.02;t=-2:R:2;syms t;f=fourier(Heaviside(t+1)-Heaviside(t-1) 运行结果：f =2/w*sin(w) (b) f(t)= MATLAB程序： R=0.02;T=-2:R:2;syms t;f=fourier(exp(-3*t)*Heaviside(t) 运行结果：f =1/(3+i*w)(2) 试画出信号f(t)= ，f(t-4)以及信号f(t) 的频谱图。MATLAB程序：r=0.02 ;t=-5 :r :5 ;f1=exp(-3*t).*stepfun(t,0) ;f2=exp(-3*(t-4).*stepfun(t,4) ;subplot(3,1,1) ;plot(t,f1) ;grid ;subplot(3,1,2) ;plot(t,f2) ;grid f3=f1.*exp(-j*4*t);subplot(3,1,3);plot(t,f3) ;grid ; 频谱图：4.综合设计 （1）MATLAB模拟滤波器的实现 I设计一个butterworth低通滤波器，满足以下性能指标：通带截止频率w1=10000rad/s，通带最大衰减Ap=3dB，阻带截止频率w2=40000rad/s，阻带最大衰减As=35dB.MATLAB程序：w1=10000;w2=40000;Ap=3;As=35;N,Wn=buttord(w1,w2,Ap,As,s);z,p,k=butter(N,Wn,s);b,a=zp2tf(z,p,k);w=linspace(1,4000,1000);H=freqs(b,a,w);magH=abs(H);phaH=unwrap(angle(H);figure(11);plot(w,20*log10(magH),k);xlabel (频率/ Hz);ylabel(幅度/B);滤波器频率响应图： I设计一个butterworth高通滤波器，满足以下性能指标：通带截止频率w1=40000rad/s，通带最大衰减Ap=3dB，阻带截止频率w2=10000rad/s，阻带最大衰减As=35dB.w1=40000;w2=10000;Ap=3;As=35;Fs=10000;N,Wc=buttord(w1,w2,Ap,As,s);Z,P,K=buttap(N);A,B,C,D=zp2ss(Z,P,K);AT,BT,CT,DT=lp2hp(A,B,C,D,Wc);num1,den1=ss2tf(AT,BT,CT,DT);num2,den2=bilinear(num1,den1,Fs);H,W=freqz(num2,den2);plot(W*Fs,abs(H);xlabel(频率（Hz）),ylabel(幅值);滤波器频率响应图：（2）MATLAB语音信号的采集与处理I.用MATLAB对原始语音信号进行分析，画出它的时域波形和频谱。程序：fs=22050;x,fs,Nbits =wavread(C:Documents and SettingsAdministrator桌面07MATLAB课程设计Matlab滤波器程序Windows XP 登录音.wav) ;t=0:1/22050:(length(x)-1)/22050;sound(x,fs) ; %对加载的语音信号进行回放figure(1)plot(x) %做原始语音信号的时域图形title(原始语音信号);xlabel(time n);ylabel(fuzhi n);figure(2)freqz(x) %绘制原始语音信号的频率响应图title(频率响应图)n=length(x); %求出语音信号的长度y1=fft(x,n) ; %傅里叶变换y2=fftshift(y1); %对频谱图进行平移f=0:fs/n:fs*(n-1)/n; %得出频点figure(3) subplot(2,1,1);plot(abs(y2) %做原始语音信号的FFT频谱图title(原始语音信号FFT频谱)subplot(2,1,2);plot(f,abs(y2); %绘制原始语音信号的频谱图title(原始语音信号频谱)%加噪声 k=1:1:n; %定义k 值,噪声与原始语音长度一致 c1=10*sin(2*pi*5000*k); %噪声信号的函数 figure(4) plot(k,c1) title(噪声信号时域图); %绘制噪声信号的时域波形图 figure(5) freqz(c1) %绘制噪声信号的频率响应图 title(噪声信号频率响应图) s=length(c1); %噪声信号的长度 xc=fft(c1,s); %对噪声信号进行频谱变换 xcc=fftshift(xc); f=0:fs/s:fs*(s-1)/s; figure(6) %绘制噪声语音的FFT频谱图 plot(f,abs(xcc) title(噪声信号的频谱图);图 原始语音信号的时域波形图 原始语音信号的频率响应图图 原始语音信号频谱和FFT频谱对比II. 给原始的语音信号加上一个高频余弦噪声，频率为（ Hz。对加噪后的语音进行分析，并画出其信号时域和频谱图。程序：fs=22050; x1=wavread(C:Documents and SettingsAdministrator桌面splash.wav);f=fs*(0:511)/1024;t=0:1/fs:(length(x1)-1)/fs;Au=0.03;%噪声幅值d=Au*sin(2*pi*5500*t);x2=x1+d;sound(x2,22050);y1=fft(x1,1024);y2=fft(x2,1024);figure(1);subplot(2,1,1);plot(t,x1);title(加噪前的信号);xlabel(time n);ylabel(fuzhi n);subplot(2,1,2);plot(t,x2)title(加噪后的信号);xlabel(time n);ylabel(fuzhi n);figure(2)subplot(2,1,1);plot(f,abs(y1(1:512);title(原始语音信号频谱);xlabel(Hz);ylabel(fuzhi);subplot(2,1,2);plot(f,abs(y2(1:512);title(加噪后的信号频谱);xlabel(Hz);ylabel(fuzhi); 图 加噪语音信号的时域波形 图 原始语音信号和加噪后的信号频谱对比与原始信号对比，区别：从原始语音信号和加噪后的信号频谱对比图中我们可以看出：在原始信号的基础上加噪后，在0-4000Hz的波形趋于平缓，而在5500z时的波形发生突变，即加了高频信号的时候，在高频部分出现了幅值变化。另外纵轴上的数值范围也发生了相应的变化。四 心得体会一周办的课程设计结束了，我接触了MATLAB，认识到了MATLAB模拟仿真对于电路设计的重要性。从看书的那一刻，我就很郁闷，MATLAB及在电子信息类（课程中的应用）到底事学什么的，有什么用。刚开始上课时，还在迷茫这门课程是用来做什么的。不过，这些感受都是在对这门课程不了解的情况下产生的。后来慢慢的接触多了，听老师讲的多了，了解多了，渐渐的产生了兴趣。对我们