信号处理实验.docx

1、信号处理实验目 录绪论 11 离散时间信号和系统分析1.1 离散时间信号产生与运算 21.2 离散时间系统的时域分析 91.3 离散时间系统的频域分析 131.4 离散时间系统频响的零极点确定 142 快速傅立叶变换的应用 2.1 FFT的计算 17 2.2 利用FFT进行谱分析 18 2.3利用FFT实现快速卷积 193 数字滤波器的设计 3.1数字滤波器的结构 23 3.2无限冲激响应(IIR)数字滤波器的设计 25 3.3有限冲激响应(FIR)数字滤波器的设计 274 综合应用举例4.1 语音信号处理 324.2 电话拨号音的合成与识别 32 绪 论 数字信号处理主要研究如何对信号进行分

2、析、变换、综合、估计与识别等加工处理的基本理论和方法。随着计算机技术和大规模集成电路技术的发展,数字信号处理以其方便、灵活等特点引起人们越来越多的重视。在40多年的发展过程中,这门学科基本形成了一套完整的理论体系，其中也包括各种快速、优良的算法,而且数字信号处理的理论和技术也在不断、快速地丰富和完善，新理论和新技术也层出不穷。学习这门课程的过程中，容易使人感到数字信号处理的概念抽象难懂，其中的分析方法与基本理论不容易很好地理解与掌握。因此，如何理解与掌握课程中的基本概念、基本原理、基本分析方法以及综合应用所学知识解决实际问题的能力，是本课程学习中所要解决的关键问题。Matlab是一种面向科学和

3、工程的高级语言，现已成为国际上公认的优秀的科技界应用软件，在世界范围内广为流行和使用。在欧美高等院校里,Matlab已成为大专院校学生、教师的必要基本技能，广泛应用于科学研究、工程计算、教学等。上世纪90年代末和本世纪初Matlab在我国也被越来越多地应用于科研和教学工作中。Matlab是一套功能强大的工程计算及数据处理软件，在工业，电子，医疗和建筑等众多领域均被广泛运用。它是一种面向对象的，交互式程序设计语言，其结构完整又具有优良的可移植性。它在矩阵运算，数字信号处理方面有强大的功能。另外，Matlab提供了方便的绘图功能，便于用户直观地输出处理结果。本文通过Matlab系列仿真，旨在掌握基

4、本的数字信号处理的理论和方法，提高综合运用所学知识，提高Matlab计算机编程的能力。进一步加强独立分析问题、解决问题的能力、综合设计及创新能力的培养，同时注意培养实事求是、严肃认真的科学作风和良好的实验习惯。1. 离散时间信号和系统分析1.1 离散时间信号产生与运算 本节的目的是使读者熟悉Matlab中离散时间信号产生和信号运算的基本命令。几种常用的序列如下：（1）单位抽样序列 在MATLAB中可以利用zeros()函数实现：例如，下列程序N = input (Type in length of sequence = ); n=0:N-1; x=zeros(1,N); x(1)=1; ste

5、m(n,x); xlabel(n);ylabel(x(n); title(单位抽样序列 N取10);输入Type in length of sequence = 10，可产生（2）单位阶越序列 在MATLAB中可以利用ones()函数实现：例如，下列程序N = input (Type in length of sequence = ); n=0:N-1; x=ones(1,N); stem(n,x); xlabel(n);ylabel(x(n);title(单位阶越序列 N取10);输入Type in length of sequence = 10，可产生（3）正弦序列在MATLAB中：例如，

6、下列程序a = input(Type in a = ); b = input(Type in b = ); A = input(Type in the gain constant = ); N = input (Type in length of sequence = ); n = 0:N; x = A*sin(a*pi*n+pi/b); stem(n,x); title(正弦序列); xlabel(Time index n);ylabel(Amplitude);输入Type in a = 0.1，Type in b = 2，Type in the gain constant = 3，Type

7、 in length of sequence = 40，可产生（4）指数序列在MATLAB中：例如，下列程序a = input(Type in exponent = ); K = input(Type in the gain constant = ); N = input (Type in length of sequence = ); n = 0:N; x = K*a.n; stem(n,x); xlabel(Time index n);ylabel(Amplitude); title( 指数序列 alpha = ,num2str(a);输入Type in exponent = 2，Type

8、 in the gain constant = 1，Type in length of sequence = 20，可产生如下结果（5）复指数序列在MATLAB中：例如，下列程序a = input(Type in real exponent = ); b = input(Type in imaginary exponent = ); c = a + b*i; K = input(Type in the gain constant = ); N = input (Type in length of sequence = ); n = 1:N; x = K*exp(c*n);subplot(211

9、); stem(n,real(x); ylabel(Amplitude); title(复指数序列 Real part); subplot(212); stem(n,imag(x); xlabel(Time index n); ylabel(Amplitude); title(复指数序列 Imaginary part);输入Type in real exponent = 0.2，Type in imaginary exponent = 0.2，Type in the gain constant = 2，Type in length of sequence = 40，可产生如下结果（6）Sinc

10、函数在MATLAB中：例如，下列程序t=-10:0.01:10; x=sinc(t); plot(t,x); xlabel(t);ylabel(x(t); title(Sinc函数);可产生（7）随即序列例如，下列程序clf; R=51; d=0.8*(rand(R,1)-0.5); m=0:R-1; stem(m,d,b); title(随机序列); xlabel(k);ylabel(f(k);可产生序列的基本运算有：（1）序列加法和乘法在MATLAB中：x= c+ b；y= c.* b；例如，下列程序%取a=2，1, 3, 4，b=0，1，2, 3, 1 m=1:4; a=2 1 3 4;

11、 c=2 1 3 4 0; n=1:5; b=0 1 2 3 1; c=a zeros(1); x=c+b; y=c.*b; subplot(4,1,1); stem(m,a);xlabel(m);ylabel(a(m); subplot(4,1,2); stem(n,b);xlabel(n);ylabel(b(n); subplot(4,1,3); stem(n,x);xlabel(n);ylabel(x(n); title(序列的加法); subplot(4,1,4); stem(n,y);xlabel(n);ylabel(y(n) ; title(序列的乘法);可产生 （2）序列的卷积在

12、MATLAB中：c=conv(a,b)；例如，下列程序a=input(Type in the first sequence =); b=input(Type in the second sequence =); c=conv(a,b); M=length(c)-1; n=0:1:M; disp(output sequence =); disp(c)；stem(n,c)； xlabel(Time index n); ylabel(Amplitude);title(序列的卷积);输入Type in the first sequence =1 2 3，Type in the second seque

13、nce =4 5 6，可产生：output sequence = 4 13 28 27 18 1.2 离散时间系统的时域分析对线性离散时间系统，若y1n和y2n分别是输入序列x1n和x2n的响应，则输入xn=ax1n+bx2n的输出响应为yn=ay1n+by2n式中叠加性质对任意常数a和b以及任意输入x1n和x2n都成立。反之，则系统称之为非线性。例如，下列程序% yn-0.4yn-1+0.75yn-2=2.2403xn+2.4908xn-1+2.2403xn-2n=0:40;a=2; b=-3;x1=cos(2*pi*0.1*n);x2=sin(2*pi*0.1*n);x=a*x1+b*x2

14、;num=2.2403 2.4908 2.2403; den=1 -0.4 0.75;ic=0 0; %设置零初始条件y1=filter(num,den,x1,ic); %计算输出y1ny2=filter(num,den,x2,ic); %计算输出y2ny=filter(num,den,x,ic); %计算输出ynyt=a*y1+b*y2;d=y-yt; %计算差值输出dn%画出输出和差信号subplot(3,1,1); stem(n,y); ylabel(振幅);title(加权输入:acdot x_1n+bcdot x_2n的输出);subplot(3,1,2); stem(n,yt);

15、ylabel(振幅);title(加权输出:acdot y_1n+bcdot y_2n);subplot(3,1,3); stem(n,d);xlabel(时间序号 n);ylabel(振幅); title(差信号);可产生对于离散时不变系统，若y1n是x1n的响应，则输入xn=x1n-n0的输出响应为yn=y1n-n0式中n0时任意整数。上面的输入输出关系，对任意输入序列及其相应的输出成立。反之，则系统称之为时变的。例如，下列程序% yn-0.4yn-1+0.75yn-2=2.2403xn+2.4908xn-1+2.2403xn-2clf;n=0:40;D=10;a=3.0;b=-2;x=a

16、*cos(2*pi*0.1*n)+b*sin(2*pi*0.1*n);xd=zeros(1,D) x;num=2.2403 2.4908 2.2403;den=1 -0.4 0.75;ic=0 0; %设置零初始条件y=filter(num,den,x,ic); %计算输出ynyd=filter(num,den,xd,ic); %计算输出ydnd=y-yd(1+D:41+D); %计算差值输出dn%画出输出subplot(3,1,1); stem(n,y); ylabel(振幅);title(输出yn);grid;subplot(3,1,2); stem(n,yt(1:41); ylabel(

17、振幅);title(由于延时输入xn,num2str(D),的输出);grid;subplot(3,1,3); stem(n,d); xlabel(时间序号 n);ylabel(振幅);title(差值信号);grid;可产生结果离散时间系统的仿真：线性和非线性系统、时变和非时变系统的仿真离散系统其输入、输出关系可用以下差分方程描述：输入信号分解为冲激信号记系统单位冲激响应，则系统响应为如下的卷积计算式：当时，hn是有限长度的（n：0，M），称系统为FIR系统；反之，称系统为IIR系统。1.3 离散时间系统的频域分析序列xn 的DTFT定义：在MATLAB中，可用freqz计算出离散时间系统的

18、频率响应。可用下列程序计算差分方程y(n)+0.7y(n-1)-0.45y(n-2)-0.6y(n-3)=0.8x(n)-0.44x(n-1)+0.36x(n-2)+0.02x(n-3) 的单位脉冲响应：% x(n)=zeros(1,N-1),0=n=40N=41;a=0.8 -0.44 0.36 0.22;b=1 0.7 -0.45 -0.6;x=1 zeros(1,N-1);k=0:1:N-1;y=filter(a,b,x);stem(k,y);xlabel(n);ylabel(幅度);可产生可用下列程序计算差分方程y(n)+0.7y(n-1)-0.45y(n-2)-0.6y(n-3)=0

19、.8x(n)-0.44x(n-1)+0.36x(n-2)+0.02x(n-3)对应系统函数的DTFT：k=256;num=0.8 -0.44 0.36 0.02;den=1 0.7 -0.45 -0.6;w=0:pi/k:pi;h=freqz(num,den,w);subplot(2,2,1);plot(w/pi,real(h);grid;title(实部);xlabel(omega/pi);ylabel(幅度);subplot(2,2,2);plot(w/pi,imag(h);grid;title(虚部);xlabel(omega/pi);ylabel(幅值);subplot(2,2,3);

20、plot(w/pi,abs(h);grid;title(幅度谱);xlabel(omega/pi);ylabel(幅值);subplot(2,2,4);plot(w/pi,angle(h);grid;title(相位谱);xlabel(omega/pi);ylabel(弧度);可产生1.4 离散时间系统频响的零极点确定离散系统的时域方程为其变换域分析方法如下：（1）频域 系统的频率响应为 （2）Z域 系统的转移函数为 分解因式 其中和称为零、极点。在MATLAB中，可以用函数z，p，K=tf2zp（num，den）求得有理分式形式的系统转移函数的零、极点，用函数zplane（z，p）绘出零、极

21、点分布图；也可以用函数zplane（num，den）直接绘出有理分式形式的系统转移函数的零、极点分布图。可用下列程序，求解已知离散系统H(z)的零极点图，并求解h(k)和H(ejw):b=1 2 1;a=1 -0.5 -0.005 0.3;subplot(311);zplane(b,a);axis(-3,3,-1,1)num=0 1 2 1;den=1 -0.5 -0.005 0.3;h=impz(num,den);subplot(312);stem(h);xlabel(k);ylabel(h(k);H,w=freqz(num,den);subplot(313);plot(w/pi,abs(H

22、);xlabel(/omege);ylabel(abs(H);可产生系统H(z)的零极点图，以及h(k)和H(ejw)：2 FFT的应用2.1 快速傅立叶变换的计算N点序列的DFT和IDFT变换定义式如下：利用旋转因子具有周期性，可以得到快速算法（FFT）。在MATLAB中，可以用函数U=fft（u，N）和u=ifft（U，N）计算N点序列的DFT正、反变换。例如，可用下列程序求x=cos(5*pi*n/16) 的16点序列的16和32点的DFTN1=16; N2=32;n1=0:N1-1; n2=0:N2-1;a=cos(5*pi*n1/16); b=cos(5*pi*n2/16); x1=

23、fft(a,N1); x2=fft(b,N2);subplot(2,1,1);stem(n1,abs(x1),.);axis(0,20,0,20);xlabel(n1);ylabel(|X1(n1)|);title(16点DFT);subplot(2,1,2);stem(n2,abs(x2),.);axis(0,20,0,20);xlabel(n2);ylabel(|X2(n2)|);title(32点DFT);可得到2.2 利用FFT进行谱分析对信号进行谱分析，就是计算信号的傅里叶变换。用Matlab语言编制信号产生子程序，产生典型信号供谱分析用，其中x(n)是由两个正弦信号及白噪声的叠加，

24、产生两个正弦加白噪声，对产生的信号进行谱分析，绘出序列和幅频特性曲线。对连续信号进行谱分析也是连续的，应先对连续信号进行时域采样，再利用DFT对采样序列进行谱分析。对信号进行谱分析，就是计算信号的傅里叶变换，绘出序列和幅频特性曲线。然而，因为对连续信号进行谱分析也是连续的，应先对连续信号进行时域采样，故而会使谱分析引入误差，所以用DFT对连续信号进行谱分析的结果都是近似的。例如，下列程序是利用FFT对x(n)=a1*sin(w*f1(0:N-1)+sin(w*f2*(0:N-1)+randn(1,N) w=2*pi/fs进行谱分析,x(n)是由两个正弦信号及白噪声的叠加N=256; a1=5;

25、 a2=3;f1=.1; f2=.2; fs=1;w=2*pi/fs;x=a1*sin(w*f1*(0:N-1)+sin(w*f2*(0:N-1)+randn(1,N);% 应用FFT求频谱subplot(2,1,1); plot(x(1:N/4); title(原始信号);y=fft(x);subplot(2,1,2); plot(f,(abs(y); title(频域信号);可得到2.3 利用FFT实现快速卷积运用DFT的快速算法FFT，对序列进行卷积，当N很大时，计算速度会快很多，用FFT计算线形卷积步骤为：（1）求H(k)=FFTh(n),N点；（2）求X(k)=DFTx(n),N点；

26、（3）计算Y(k)=X(k)H(k)；（4）求y(n)=IDFTY(k),N点。下列程序可实现序列xn=0.8.n和hn=ones(1,N2)的快速卷积，n=0:1:50;m=0:1:20;N1=length(n);N2=length(m); xn=0.8.n; %生成x(n)hn=ones(1,N2); %生成h(n)N=N1+N2-1;XK=fft(xn,N);HK=fft(hn,N);YK=XK.*HK;yn=ifft(YK,N);if all(imag(xn)=0)&(all(imag(hn)=0) %实序列的循环卷积仍然为实序列yn=real(yn);endx=0:N-1;stem(

27、x,yn,.);得到卷积序列：3 数字滤波器的设计3.1 数字滤波器的结构无限冲击响应传输函数为有限冲激响应传输函数为在MATLAB中无限冲激响应传输函数可用z,p,k=tf2zp(num,den)实现，有限冲激响应传输函数可用r1,p1,k1=residuez(num,den)和r2,p2,k2=residue(num,den)实现。下列程序，可实现一个无限冲激响应传输函数的并联形式num=input(分子系数向量=);den=input(分母系数向量=);r1,p1,k1=residuez(num,den);r2,p2,k2=residue(num,den);disp(并联I型);disp

28、(零点是);disp(r1);disp(极点是);disp(p1);disp(常数);disp(k1);disp(并联II型);disp(零点是);disp(r2);disp(极点是);disp(p2);disp(常数);disp(k2);运行后输入：分子系数向量=1 2 3，分母系数向量=4 5 6，可得并联I型零点是-0.1250 - 0.0148i -0.1250 + 0.0148i极点是-0.6250 + 1.0533i-0.6250 - 1.0533i常数0.5000并联II型零点是0.0938 - 0.1224i0.0938 + 0.1224i极点是-0.6250 + 1.0533i-0.6250 - 1.0533i常