信号处理实验指导

基于MATLAB的DSP软件仿真目 录绪论 11 离散时间信号和系统分析1.1 离散时间信号产生与运算 21.2 离散时间系统的时域分析 91.3 离散时间系统的频域分析 131.4 离散时间系统频响的零极点确定 142 快速傅立叶变换的应用2.1 FFT的计算 172.2 利用FFT进行谱分析 182.3利用FFT实现快速卷积 193 数字滤波器的设计3.1数字滤波器的结构 233.2无限冲激响应(IIR)数字滤波器的设计 253.3有限冲激响应(FIR)数字滤波器的设计 274 综合应用举例4.1 语音信号处理 324.2 电话拨号音的合成与识别 32 绪 论 数字信号处理主要研究如何对信号进行分析、变换、综合、估计与识别等加工处理的基本理论和方法。随着计算机技术和大规模集成电路技术的发展,数字信号处理以其方便、灵活等特点引起人们越来越多的重视。在40多年的发展过程中,这门学科基本形成了一套完整的理论体系，其中也包括各种快速、优良的算法,而且数字信号处理的理论和技术也在不断、快速地丰富和完善，新理论和新技术也层出不穷。学习这门课程的过程中，容易使人感到数字信号处理的概念抽象难懂，其中的分析方法与基本理论不容易很好地理解与掌握。因此，如何理解与掌握课程中的基本概念、基本原理、基本分析方法以及综合应用所学知识解决实际问题的能力，是本课程学习中所要解决的关键问题。Matlab是一种面向科学和工程的高级语言，现已成为国际上公认的优秀的科技界应用软件，在世界范围内广为流行和使用。在欧美高等院校里,Matlab已成为大专院校学生、教师的必要基本技能，广泛应用于科学研究、工程计算、教学等。上世纪90年代末和本世纪初Matlab在我国也被越来越多地应用于科研和教学工作中。Matlab是一套功能强大的工程计算及数据处理软件，在工业，电子，医疗和建筑等众多领域均被广泛运用。它是一种面向对象的，交互式程序设计语言，其结构完整又具有优良的可移植性。它在矩阵运算，数字信号处理方面有强大的功能。另外，Matlab提供了方便的绘图功能，便于用户直观地输出处理结果。本文通过Matlab系列仿真，旨在掌握基本的数字信号处理的理论和方法，提高综合运用所学知识，提高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; stem(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中：例如，下列程序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 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 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); 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函数在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; 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）序列的卷积在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 sequence =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;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); 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*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('振幅');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计算出离散时间系统的频率响应。可用下列程序计算差分方程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.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);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）绘出零、极点分布图；也可以用函数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);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=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)是由两个正弦信号及白噪声的叠加，产生两个正弦加白噪声，对产生的信号进行谱分析，绘出序列和幅频特性曲线。对连续信号进行谱分析也是连续的，应先对连续信号进行时域采样，再利用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; 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点；（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(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('零点是');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常数0.2500下列程序，可实现一个有限冲激响应传输函数num=input('分子系数向量=');den=input('分母系数向量=');z,p,k=tf2zp(num,den);sos=zp2sos(z,p,k);disp('z是');disp(z);disp('p是');disp(p);disp('k是');disp(k);运行后输入：分子系数向量=2 10 23 34 31 16 4，分母系数向量=1 0 0 0 0 0 0，可得z是-2.0000 -0.5000 + 1.3229i-0.5000 - 1.3229i-1.0000 -0.5000 + 0.5000i-0.5000 - 0.5000ip是0k是23.2 无限冲激响应(IIR)数字滤波器的设计IIR数字滤波器的系统函数IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，如巴特沃斯、契比雪夫和椭圆滤波器等，有现成的设计数据或图表可查，其设计工作量比较小，对计算工具的要求不高。在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。 巴特沃斯IIR滤波器的设计： 在MATLAB下，设计巴特沃斯IIR滤波器可使用butter函数。 Butter函数可设计低通、高通、带通和带阻的数字和模拟IIR滤波器，其特性为使通带内的幅度响应最大限度地平坦，但同时损失截止频率处的下降斜度。在期望通带平滑的情况下，可使用butter函数。 butter函数的用法为：b,a=butter(n,Wn,/ftype/) ，其中n代表滤波器阶数，Wn代表滤波器的截止频率，这两个参数可使用buttord函数来确定。buttord函数可在给定滤波器性能的情况下，求出巴特沃斯滤波器的最小阶数n，同时给出对应的截止频率Wn。buttord函数的用法为：n,Wn= buttord(Wp,Ws,Rp,Rs)，其中Wp和Ws分别是通带和阻带的拐角频率（截止频率），其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹系数。 不同类型（高通、低通、带通和带阻）滤波器对应的Wp和Ws值遵循以下规则： (1)高通滤波器：Wp和Ws为一元矢量且Wp>Ws； (2)低通滤波器：Wp和Ws为一元矢量且Wp<Ws； (3)带通滤波器：Wp和Ws为二元矢量且Wp<Ws，如Wp=0.2,0.7,Ws=0.1,0.8; (4)带阻滤波器：Wp和Ws为二元矢量且Wp>Ws，如Wp=0.1,0.8,Ws=0.2,0.7。契比雪夫I型IIR滤波器的设计： 在期望通带下降斜率大的场合，应使用椭圆滤波器或契比雪夫滤波器。在MATLAB下可使用cheby1函数设计出契比雪夫I型IIR滤波器。 cheby1函数可设计低通、高通、带通和带阻契比雪夫I型滤IIR波器，其通带内为等波纹，阻带内为单调。契比雪夫I型的下降斜度比II型大，但其代价是通带内波纹较大。 cheby1函数的用法为：b,a=cheby1(n,Rp,Wn,/ftype/) ，在使用cheby1函数设计IIR滤波器之前，可使用cheblord函数求出滤波器阶数n和截止频率Wn。cheblord函数可在给定滤波器性能的情况下，选择契比雪夫I型滤波器的最小阶和截止频率Wn。 cheblord函数的用法为：n,Wn=cheblord(Wp,Ws,Rp,Rs) ，其中Wp和Ws分别是通带和阻带的拐角频率（截止频率），其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹系数。例如，下列程序可实现一个巴特沃兹低通数字滤波器的设计wp=100*2*pi;ws=200*2*pi;Rp=2; Rs=15;Fs=500; Ts=1/Fs;N,Wn=buttord(wp,ws,Rp,Rs,'s');Z,P,K=buttap(N);Bap,Aap=zp2tf(Z,P,K);b,a=lp2lp(Bap,Aap,Wn);bz,az=bilinear(b,a,Fs);H,W=freqz(bz,az);subplot(211);plot(W*Fs/(2*pi),abs(H); grid xlabel('频率/Hz'); ylabel('频率响应幅度');运行可得到：3.3有限冲激响应(FIR)数字滤波器的设计数字滤波器的设计是数字信号处理中的一个重要内容。数字滤波器设计包括FIR（有限单位脉冲响应）滤波器与IIR（无限单位脉冲响应）滤波器两种。 与IIR滤波器相比，FIR滤波器在保证幅度特性满足技术要求的同时，很容易做到严格的线性相位特性。设FIR滤波器单位脉冲响应h(n)长度为N，其系统函数为H（z），H（z）是z1的N1次多项式，它在z平面上有N1个零点，原点z=0是N1阶重极点，因此H（z）是永远稳定的。稳定和线性相位特性是FIR滤波器突出的优点。 FIR滤波器的设计任务是选择有限长度的h(n)。使传输函数H(z)满足技术要求。FIR滤波器的设计方法有多种，如窗函数法、频率采样法及其它各种优化设计方法，本实验介绍窗函数法的FIR滤波器设计。 窗函数法是使用矩形窗、三角窗、巴特利特窗、汉明窗、汉宁窗和布莱克曼窗等设计出标准响应的高通、低通、带通和带阻FIR滤波器。 （1）firl函数的使用 在MATLAB下设计标准响应FIR滤波器可使用firl函数。firl函数以经典方法实现加窗线性相位FIR滤波器设计，它可以设计出标准的低通、带通、高通和带阻滤波器。firl函数的用法为： b=firl(n,Wn,/ftype/,Window) 各个参数的含义如下： b滤波器系数。对于一个n阶的FIR滤波器，其n+1个滤波器系数可表示为：b(z)=b(1)+b(2)z1+b(n+1)zn。n滤波器阶数。 Wn截止频率，0Wn1，Wn=1对应于采样频率的一半。当设计带通和带阻滤波器时，Wn=W1 W2,W1W2。ftype当指定ftype时，可设计高通和带阻滤波器。Ftype=high时，设计高通FIR滤波器；ftype=stop时设计带阻FIR滤波器。低通和带通FIR滤波器无需输入ftype参数。 Window窗函数。窗函数的长度应等于FIR滤波器系数个数，即阶数n+1。 （2）窗函数的使用 在MATLAB下，这些窗函数分别为： 1矩形窗：w=boxcar(n)，产生一个n点的矩形窗函数。 2三角窗：w=triang(n)，产生一个n点的三角窗函数。 3巴特利特窗：w=Bartlett(n)，产生一个n点的巴特利特窗函数。 4汉明窗：w=hamming(n)，产生一个n点的汉明窗函数。 5汉宁窗：w=hanning(n)，产生一个n点的汉宁窗函数。6布莱克曼窗：w=Blackman(n),产生一个n点的布莱克曼窗函数。 7凯泽窗：w=Kaiser(n,beta)，产生一个n点的凯泽窗数。8契比雪夫窗：w=chebwin(n,r)产生一个n点的契比雪夫窗函数。下列程序可演示常用窗函数及其频谱特性N=31;n=0:1:(N-1);%矩形窗W_box=boxcar(N);Hbox,W=freqz(W_box,1);subplot(4,2,1);stem(n,W_box);xlabel('n');ylabel('矩形窗');subplot(4,2,2);plot(W/pi,20*log10(abs(Hbox)/abs(Hbox(1);ylabel('矩形窗频谱');%三角窗W_tri=triang(N);Htri,W=freqz(W_tri,1);subplot(4,2,3);stem(n,W_tri);xlabel('n');ylabel('三角窗');subplot(4,2,4);plot(W/pi,20*log10(abs(Htri)/abs(Htri(1);ylabel('三角窗频谱');%汉宁窗W_han=hanning(N);Hhan,W=freqz(W_han,1);subplot(4,2,5);stem(n,W_han);xlabel('n');ylabel('汉宁窗');subplot(4,2,6);plot(W/pi,20*log10(abs(Hhan)/abs(Hhan(1);ylabel('汉宁窗频谱');%汉明窗W_ham=hamming(N);Hham,W=freqz(W_ham,1);subplot(4,2,7);stem(n,W_ham);xlabel('n');ylabel('汉明窗');subplot(4,2,8);plot(W/pi,20*log10(abs(Hham)/abs(Hham(1);ylabel('汉明窗频谱');运行得到：下列程序可实现一个有限冲激响应(FIR)数字滤波器的设计% 通带边缘频率:wp1=0.45*pi,wp2=0.65*pi,通带峰值起伏：rp<=1dB% 阻带边缘频率:ws1=0.3*pi,ws2=0.75*pi,阻峰值起伏：rs>=40dBwp1=0.45*pi;wp2=0.65*pi;ws1=0.3*pi;ws2=0.75*pi;rp=1;rs=40;width=min(wp1-ws1),(ws2-wp2);N=ceil(11*pi/width)+1;n=0:1:(N-1); a=(N-1)/2; m=n-a+eps;w1=(wp1+ws1)/2;wh=(ws2+wp2)/2;hd=(sin(wh*m)-sin(w1*m)./(pi*m);w_bla=(blackman(N)'h=hd.*w_bla;H,w=freqz(h,1);subplot(2,1,1); stem(n,h);ylabel('h(n)'); title('脉冲响应');subplot(2,1,2); plot(w/pi,20*log10(abs(H)/max(H);xlabel('频率'); ylabel('幅频响应'); title('滤波器频响特性');运行得到：4 综合应用举例4.1 语音信号处理语音信号处理综合运用了数字信号处理的理论知识，对信号进行计算及频谱分析，设计滤波器，并对含噪信号进行滤波。具体分为以下步骤：（1）语音信号的采集：利用Windows下的录音机，录制一段话音。然后在Matlab软件平台下，利用函数wavread对语音信号进行采样，播放语音信号，并绘制原始语音信号；（2）对原始信号加入噪声：对原始语音信号加入s=sin(2*pi*f*Ts*n)的噪声，采样后可知Fs = 16000，选择f = 2500，播放加入噪声信号的语音信号，并绘制噪声信号和含噪语音信号；（3）频谱分析：分别对原始语音信号，噪声信号和含噪声的语音信号进行频谱分析，并绘出各频谱图；（4）设计滤波器：计算滤波器的性能指标，设计滤波器，绘制滤波器的特性曲线；（5）滤波器滤波：用自己设计的滤波器对采集的信号进行滤波，得出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化，并回放语音信号，感觉滤波前后的声音有变化。4.2 电话拨号音的合成与识别双音多频 DTMF（ Dual Tone Multi-Frequency ）信号，是用两个特定的单音频率信号的组合来代表数字或功能。在 DTMF 电话机中有 16 个按键，其中 10 个数字键 0 9 ， 6 个功能键 * 、 # 、 A 、 B 、 C 、 D 。其中 12 个按键是我们比较熟悉的按键，另外由第 4 列确定的按键作为保留，作为功能键留为今后他用。 根据 CCITT 建议，国际上采用 697Hz 、 770Hz 、 852Hz 、 94lHz 低频群及 1209Hz 、 1336Hz 、 1477Hz 、 1633Hz 高频群。从低频群和高频群任意各抽出一种频率进行组合，共有 16 种组合，代表 16 种不同的数字键或功能，每个按键唯一地由一组行频和列频组成，如下表所示。利用 MATLAB 软件能够利用矩阵不同的基频合成 0 9 不同按键的拨号音，并能够对不同的拨号音加以正确的识别，实现由拨号音解析出电话号码的过程。进一步利用MATLAB 中的图形用户界面GUI 做出简单的图形操作界面。从而实现对电话拨号音系统的简单的实验仿真。具体实现步骤如下：（1）图形电话拨号面板的制作 利用 GUI 图形用户界面设计工具制作电话拨号面板，把 DTMF 信号和电话机的键盘矩阵对应起来。其中选用我们熟悉的 10 个数字键 0 9 ， 2 个功能键“ * ”、“”，另外为了仿真方便，添加信号识别键和复位键。每个按键可用 （ Push Button ）添加。 最终利用 GUI 图形用户界面设计工具生成的图形电话拨号面板用于拨号音的合成产生部分，如下图所示。这里将其保存为Untitle.fig文件。 （2）DTMF 信号的产生合成现在将对上节制作的图形电话拨号面板上的各控件单位的动作和变化进行设置，即对Untitle.m 文件进行编辑。其主要的功能是使对应的按键，按照表中的对应关系产生相应的拨号音，完成对应行频及列频的叠加输出。此外，对于图形界面的需要，还要使按键的号码数字显示在拨号显示窗口中。 鉴于 CCITT 对 DTMF 信号规定的指标，这里每个数字信号取 1000 个采样点模拟按键信号，并且每两个数字之间用 100 个 0 来表示间隔来模拟静音。以便区别连续的两个按键信号。间隔的静音信号也是在按键时产生的。 以按键 0 为例，简单介绍拨号音产生的过程： % 按键 0 的响应函数 function varargout = pushbutton0_Callback(h, eventdata, handles, varargin) n=1:1000; % 每个数字 1000 个采样点表示 d0=sin(0.7217*n)+sin(1.0247*n); % 对应行频列频叠加 n0=strcat(get(handles.edit1,'string'),'0'); % 获取数字号码 set(handles.edit1,'string',n0); % 显示号码 space=zeros(1,100); %100 个 0 模拟静音信号 global NUM phone=NUM,d0; NUM=phone,space; % 存储连续的拨号音信号 wavplay(d0,8192); % 产生拨号音 程序解释： NUM 为定义的全局变量，用于存储连续的拨号音（ DTMF ）信号，包括数字信号音以及静音信号。 d0=sin(0.7217*n)+sin(1.0247*n) 中的行频与列频是由表 1 中 0 键对应的, 计算得出，已知声音取样频率 则取样后,对于保留的两个功能键“ * ”、“”，按照现行键盘式拨号电话的习惯，将“ * ”作为删除键，“”作为确认键。“ * ”删除键的作用是将前面拨错的号码删除退回，表现为将显示窗口已经显示的错误号码退回一位数字，并且将连续拨号音信号的存储单元 NUM 中退回一位拨号音信号和静音信号。删除可以进行连续的操作。“”确认键的作用是将前面拨过的号码进行确认保留，意味着此时连续拨号音信号的存储单元 NUM 中的信号即为最后用于识别的连续拨号音 DTMF 信号，并在显示窗口中显示“”号作为标记。 % 删除键的响应函数 function varargout = pushbuttonback_Callback(h, eventdata, handles, varargin) n=1:1000; num=get(handles.edit1,'string'); l=length(num); n11=strrep(num,num,num(1:l-1); %去掉末尾号码在面板上的显示 d11=sin(0.7217*n)+sin(0.9273*n); set(handles.edit1,'string',n11); global NUM L=length(NUM); NUM=NUM(1:L-1100); %删除末尾号码在拨号音信号中的存储 wavplay(d11,8192);  （3）DTMF 信号的检测识别要实现电话拨号音（ DTMF ）信号的检测识别，可以通过直接计算付里叶变换得到输入信号的组成频率。这里采用 FFT 算法对信号进行解码分析。首先对接收到的数字信号作 FFT 分析，计算出其幅频谱，进而得到功率谱，组成输入信号的频率必定对应功率谱的峰值。对于连续的双音多频（ DTMF ）信号，需要把有效的数字拨号信号从静音间隔信号中分割提取出来，然后再用 FFT 算法对信号进行解码分析。 DTFM信号的最小周期为40ms，取采样周期为8KHz，每个周期至少采样0.04*8000=320点。DTMF解码器计算采样序列的DFT，从给定的8个基本频中找出与结果最接近的频率。DFT的长度N决定了DFT样本的位置和计算该样本所用时间直接的间隔。对8KHz的采样率，人们发现检测八个基本DTMF音调的DFT长度N的最佳值为205，表中给出N205是最接近的音调频率的DFT指标k。基本音调 (Hz)实际k67918770208522294124120931133634147738163342 参考文献1 程佩青数字信号处理教程北京：清华大学出版社，20022 薛年喜MATLAB在数字信号处理中的应用北京：清华大学出版社，20033 丁玉美等数字信号处理西安:西安电子科技大学出版社,20014 陈怀琛等MATLAB及在电子信息课程中的应用北京:电子工业出版社,20035 蔡莲红，黄德智，蔡锐现代语音技术基础与应用北京：清华大学出版社，20036 张明照，刘政波，刘斌等应用MATLAB 实现信号分析和处理北京：科学出版社，20067 王树勋数字信号处理基础及实验北京：机械工业出版社，19928 胡广书数字信号处理基础及实验北京：机械工业出版社，19929 楼顺天，李博菡基于MATLAB的系统分析与设计信号处理西安：西安电子科技大学出版社，199810 侯朝焕等实用FFT信号处理的应用北京：科学出版社，199111 J.G.Proakis and D.G.Manolakis. Digital Signal Pricessing: Principles, Algorithms and Applications Macmillan. New York. NY. third edition, 199612 Programs for Digital Signal Processing. New York. IEEE Press, 197913 Oppenheim A.V.and Schafer R.W.Digital Signal Processing. Prentice Hall. Inc, 197514 Mischa Schwartz and Leonard Shaw. Signal Processing Discrete Spectral Analysis,ditection and Estimation. McGraw HillBook Company, 197515 Bernard Gold and Lawrence R.Rabiner. Theory and Application of Digital Signal Processing. Prentice-Hall. Inc, 197537