信号处理实验指导

目 录绪论 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=nWs； (2)低通滤波器：Wp和Ws为一元矢量且WpWs； (3)带通滤波器：Wp和Ws为二元矢量且WpWs，如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=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