数字信号处理实验讲义 -采用.doc

实验一 信号（模拟、数字）的输入输出实验（常见离散信号产生和实现）一、实验目的1加深对常用离散信号的理解；2掌握matlab中一些基本函数的建立方法。二、实验原理1.单位抽样序列 在MATLAB中可以利用zeros()函数实现。如果在时间轴上延迟了k个单位，得到即： 2单位阶跃序列 在MATLAB中可以利用ones()函数实现。x=ones(1,N) 3正弦序列在MATLAB中, n=0:N-1; x=A*sin(2*pi*f*n/Fs+fai)4复指数序列在MATLAB中,n=0:N-1; x=r*exp(j*w*n)5指数序列在MATLAB中,n=0:N-1; x=a.n三、实验内容实现和图形生成1五种基本函数的生成程序如下：(1)单位抽样序列方法一：% 单位抽样序列和延时的单位抽样序列n=0:10;x1=1 zeros(1,10);x2=zeros(1,5) 1 zeros(1,5);subplot(1,2,1);stem(n,x1);xlabel (时间序列n);ylabel(振幅);title(单位抽样序列x1);subplot(1,2,2);stem(n,x2); xlabel(时间序列n);ylabel(振幅);title(延时了5的单位抽样序列);方法二：先在matlab中定义单位抽样序列：function x,n=dwxl(n1,n2,n0)n=n1:n2;x=(n=n0);在运行命令：x,n=dwxl(-5,5,3)stem(n,x);xlabel(n);title(单位抽样序列x); （2）单位阶跃序列方法一：n=0:10;u=ones(1,11);stem(n,u);xlabel (时间序列n);ylabel(振幅);title(单位阶跃序列); 所得的图形如下所示：方法二;先在matlab中定义单位阶跃序列：function x,n=jyxl(n1,n2,n0)n=n1:n2;x=(n=n0);在运行命令：x,n=jyxl(-5,5,3)stem(n,x);xlabel(n);title(单位阶跃序列x);（3）正弦函数n=1:30;x=2*sin(pi*n/6+pi/3);stem(n,x); xlabel (时间序列n);ylabel(振幅);title(正弦函数序列x=2*sin(pi*n/6+pi/3);(4)复指数序列n=1:30;x=2*exp(j*3*n);stem(n,x); xlabel (时间序列n);ylabel(振幅);title(复指数序列x=2*exp(j*3*n); 图形如下：（5）指数序列n=1:30;x=1.2.n;stem(n,x); xlabel (时间序列n);ylabel(振幅);title(指数序列x=1.2.n);2绘出信号的频率是多少？周期是多少？产生一个数字频率为0.9的正弦序列，并显示该信号，说明其周期?程序如下：n=0:40;x1=1.5*sin(2*pi*0.1*n);x2=sin(0.9*n);subplot(1,2,1);stem(n,x1); xlabel (时间序列n);ylabel(振幅);title(正弦序列x1=1.5*sin(2*pi*0.1*n);subplot(1,2,2);stem(n,x2); xlabel (时间序列n);ylabel(振幅);title(正弦序列x2=sin(0.9*n);运行结果如下：由上图看出：x1=1.5*sin(2*pi*0.1*n)的周期是10，而x2=sin(0.9*n)是非周期的。理论计算中对第一个，N2*/（0.1*）10，第二个0.9不是pi的倍数，所以不是周期的。因此可以看出，实验结果和理论相符。3x(n)=2,3,1,1,2,-1,0,3，-2n5；h(n)=2,4,1,-2,0,-1，-3n2，手工计算和MATLAB计算卷积y(n)=x(n)*h(n)。参考程序如果、的起点不为0，则采用conv_m计算卷积；编写conv_m函数：functiony,ny=conv_m(x,nx,h,nh)%改进卷积程序nyb=nx(1)+nh(1);nye=nx(length(x)+nh(length(h);ny=nyb:nye;y=conv(x,h); %MATLAB自带的函数在命令窗口输入：x=3,11,7,0,-1,4,2;nx=-3:3;h=2,3,0,-5,2,1;nh=-1:4;y,ny=conv_m(x,nx,h,nh)stem(ny,y,.);xlabel(时间序号n);title(卷积和y(n)=x(n)*h(n);4function x,n=impseq(n1,n2,n0) function x,n=stepseq(n1,n2,n0) n=n1:n2; n=n1:n2; x=(n-n0)=0; x=(n-n0)=0;四、实验报告要求1报告中要给出实验的MATLAB程序；2简述实验目的和原理；3给出用笔算时差分方程解、卷积和conv计算线性卷积对照图；4给出收获和体会。五、问题讨论与总结：1离散正弦序列的性质：离散正弦序列就是一个连续的正弦信号被一系列冲激函数采样后的结果，原连续正弦函数一定是周期的，但采样后的离散序列却不一定是周期的。对于离散序列x=sin(n*w)来说，只有当2*pi/是一个有理数时，也就是说当是pi的倍数时，此离散序列才是周期的。所以在本实验中x1=1.5*sin(2*pi*0.1*n)的周期是10，而x2=sin(0.9*n)是非周期的。因为0.9不是pi的倍数。2离散复指数序列性质：对于离散复指数函数=a*exp(z*)，只有当z是纯虚数，且纯虚数的系数是pi的倍数时，才是周期的。其它情况下均不是。 实验二 FFT频谱分析及应用一、实验目的1通过实验加深对FFT的理解；2熟悉应用FFT对典型信号进行频谱分析的方法。二、实验内容使用MATLAB程序实现信号频域特性的分析。涉及到离散傅立叶变换(DFT)、快速傅立叶变换(FFT)及信号频率分辨率等知识点。三、实验原理与方法和手段 在各种信号序列中，有限长序列占重要地位。对有限长序列可以利用离散傅立叶变换(DFT)进行分析。DFT不但可以很好的反映序列的频谱特性，而且易于用快速算法(FFT)在计算机上进行分析。有限长序列的DFT是其z变换在单位圆上的等距离采样，或者说是序列傅立叶的等距离采样，因此可以用于序列的谱分析。FFT是DFT的一种快速算法，它是对变换式进行一次次分解，使其成为若干小数据点的组合，从而减少运算量。在MATLAB信号处理工具箱中的函数fft（x,N）,可以用来实现序列的N点快速傅立叶变换。经函数fft求得的序列一般是复序列，通常要求出其幅值和相位。MATLAB中提供了求复数的幅值和相位的函数：abs、angle，这些函数一般和fft同时使用。四、实验步骤在“开始-程序”菜单中，找到MATLAB程序，运行启动；进入MATLAB后 ，在Command Window中输入实验程序，并执行；记录运行结果图形，作分析。具体步骤如下：1模拟信号，以进行采样，求：（1）N40点FFT的幅度频谱，从图中能否观察出信号的2个频谱分量?（2）提高采样点数，如N128，再求该信号的幅度频谱，此时幅度频谱发生了什么变化？信号的2个模拟频率和数字频率各为多少？FFT频谱分析结果与理论上是否一致？2有限长序列x(n)=2,1,0,1,3;h(n)=1,3,2,1，试利用FFT实现由DFT计算线性卷积，并与线性卷积直接计算（conv）的结果进行比较。参考程序1程序1：%xh11N=40;n=0:N-1;t=0.01*n;x=2*sin(4*pi*t)+5*cos(8*pi*t);k=0:N/2;w=2*pi/N*k;X=fft(x,N);magX=abs(X(1:N/2+1);subplot(2,1,1);stem(n,x,.);title(signal x(n);subplot(2,1,2);plot(w/pi,magX);title(FFT N=40);xlabel(w (unit :pi);ylabel(|X|);grid2程序2：%xh12% 用FFT实现由DFT计算线性卷积x=1 2 0 1;h=2 2 1 1;L=length(x)+length(h)-1; XE=fft(x,L);HE=fft(h,L); y1=ifft(XE.*HE);%画出由圆周卷积计算线性卷积结果及误差k=0:L-1;subplot(211);stem(k,real(y1);axis(0 6 0 7);title(Result of Linear Convolution);xlabel(Time index k);ylabel(Amplitude);y2=conv(x,h);error=y1-y2;subplot(212);stem(k,abs(error);xlabel(Time index k);ylabel(Amplitude); title(Error Magnitude);五、实验报告要求1报告中要给出实验的MATLAB程序；2简述实验目的和原理；3按实验步骤附上实验信号序列和幅频特性曲线，分析所得到的图形，回答每一步提出的问题。实验三 IIR数字滤波器的设计一、实验目的1掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理和方法；2观察双线性变换法和脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法和脉冲响应不变法的特点和区别。二、实验内容使用MATLAB编写程序，实现IIR数字滤波器的设计。涉及脉冲响应不变法和双线性变换法设计IIR数字滤波器的方法、不同设计方法得到的IIR滤波器频域特性异同等知识点。三、实验原理与方法和手段1脉冲响应不变法所谓脉冲响应不变法就是使数字滤波器的单位脉冲响应序列h(n)等于模拟滤波器的单位冲激响应和ha(t)的采样值，即：，其中，T为采样周期。在脉冲响应不变法中，模拟角频率和数字角频率的变换关系为：，可见，和之间的变换关系为线性的。在MATLAB中，可用函数impinvar实现从模拟滤波器到数字滤波器的脉冲响应不变映射，调用格式为：B,A=impinvar(b,a,fs1)B,A=impinvar(b,a)其中，b、a分别为模拟滤波器的分子和分母多项式系数向量；fs1为采样频率（Hz），缺省值fs=1Hz；B、A分别为数字滤波器分子和分母多项式系数向量。2双线性变换法：由于s平面和z平面的单值双线性映射关系为s，其中T为采样周期。因此，若已知模拟滤波器的传递函数，将上式代入即可得到数字滤波器的系统函数H(z)。在双线性变换中，模拟角频率和数字角频率的变换关系为：，可见，和之间的变换关系为非线性的。在MATLAB中，可用函数bilinear实现从模拟滤波器到数字滤波器的双线性变换映射，调用格式为：B,A=bilinear(b,a,fs1)3数字滤波器设计（1）定技术指标转换为模拟滤波器设计性能指标。（2）估计满足性能指标的模拟相应滤波器性能阶数和截止频率。利用MATLAB中buttord、cheb1ord、cheb2ord、ellipord等函数，调用格式如：其中，p为通带边界频率，rad/s；s为阻带边界频率，rad/s;p为带通波动，dB;s为阻带衰减，dB；s表示为模拟滤波器；函数返回值N为模拟滤波器的最小阶数；c为模拟滤波器的截止频率(-3dB频率)，rad/s。函数适用低通、高通、带通、带阻滤波器。（3）设计模拟滤波器。MATLAB信号处理工具箱提供了模拟滤波器设计的完全工具函数：butter、cheby1，cheby2、ellip、besself。用户只需一次调用就可完成低通、高通、带通、带阻滤波器设计。调用格式如：b,a=butter(N,c,ftype,s)，其中，ftype为滤波器类型：high表示高通滤波器，截止频率为c；stop表示带阻滤波器，c=1 2 (12)；ftype缺省时表示为低通或带通滤波器。（）利用脉冲响应不变法或双线性不变法，实现模拟滤波器到数字滤波器的映射。四、实验条件1具有WINDOWS 98/2000/NT/XP操作系统的计算机一台；2MATLAB编程软件。五、实验步骤在“开始-程序”菜单中，找到MATLAB程序，运行启动；进入MATLAB后 ，在Command Window中输入自己编写的主程序，并执行；记录运行结果图形，作分析对比。具体步骤如下：1、查看帮助文件，了解相关函数的调用格式。2、实验题目 例1：设计模拟巴特沃斯低通滤波器，fp=300Hz，p=1dB，fs=800Hz，s=20dB。%xh13模拟低通滤波器技术指标p=1; s=20; fp=300;fs=800;ap=2*pi*fp;as=2*pi*fs;%设计模拟巴特沃斯低通滤波器N,ac=buttord(ap,as,p,s,s)b,a=butter(N,ac,s) %设计模拟巴特沃斯低通滤波器，ap为通带边界频率，rad/s；as为阻带边界频率，rad/s;p为通带最大衰减，dB；s为阻带最小衰减，dB；s表示为模拟滤波器；函数返回值N为模拟滤波器的最小阶数；ac为模拟滤波器的截止频率(-3dB频率)，rad/s；b、a分别为模拟滤波器的系统函数分子和分母多项式系数向量；H,=freqs(b,a); %求模拟滤波器的频率响应%绘制频响幅度谱plot(/2/pi,20*log10(abs(H); %横轴为频率，单位：HZ;纵轴频响幅度，单位：dBaxis(0,1500,-50,0);xlabel(频率Hz);ylabel( H幅值dB);例2：设计模拟巴特沃斯高通滤波器，fp=800Hz，p=1dB，fs=300Hz，s=20dB。%xh14模拟高通滤波器技术指标p=1;s=20; fp=800;fs=300;p=2*pi*fp;s=2*pi*fs;N,c=buttord(p,s,p,s,s)b,a=butter(N,c,high,s) %设计模拟巴特沃斯高通滤波器，p为通带边界频率，rad/s；s为阻带边界频率，rad/s；p为通带最大衰减，dB；s为阻带最小衰减，dB；s表示为模拟滤波器；函数返回值N为模拟滤波器的最小阶数；c为模拟滤波器的截止频率(-3dB频率)，rad/s；b、a分别为模拟高通滤波器的系统函数分子和分母多项式系数向量；H,=freqs(b,a); plot(/2/pi,20*log10(abs(H); axis(0,1000,-50,0);xlabel(频率Hz);ylabel( H幅值dB);例3：fp=0.kHZ,p=1dB，fs=0.3kHZ，s=25dB，T=1ms；分别用脉冲响应不变法和双线性变换法设计一个Butterworth数字低通滤波器。%xh15采用冲击响应不变法p=1;s=25;fp=100;fs=300;p=2*pi*fp;s=2*pi*fs; %数字滤波器技术指标要求转化成模拟滤波器技术指标要求fs1=1000; %采样频率%设计模拟滤波器N,c=buttord(p,s,p,s,s);b,a=butter(N,c,s); B,A=impinvar(b,a,fs1) %用冲击响应不变法将模拟滤波器变换成数字滤波器。B、A分别为数字滤波器的系统函数分子和分母多项式系数向量；H1,w=freqz(B,A,whole); %求数字滤波器的频率响应%绘制数字滤波器频响幅度谱subplot(211);plot(w*fs1/2/pi,20*log10(abs(H1); axis(0,1000,-100,0);xlabel(频率Hz);ylabel( H1幅值dB);%采用双线性变换法%数字滤波器的技术指标要求p=1;s=25; fp=100;fs=300;fs1=1000; %采样频率wp=2*pi*fp/fs1; %数字滤波器通带边界频率ws=2*pi*fs/fs1; %数字滤波器阻带边界频率%变换为同类型模拟滤波器的技术指标要求 p=2*fs1*tan(wp/2); %同类模拟滤波器通带边界频率s=2*fs1*tan(ws/2); %同类模拟滤波器阻带边界频率%设计同类型模拟滤波器N,c=buttord(p,s,p,s,s);b,a=butter(N,c,s);%用双线性变换法将模拟滤波器变换成数字滤波器B,A=bilinear(b,a,fs1) H2,w=freqz(B,A,whole); %求数字滤波器的频率响应%绘制数字滤波器频响幅度谱subplot(212);plot(w*fs1/2/pi,20*log10(abs(H2);axis(0,1000,-100,0);xlabel(频率Hz);ylabel(H2幅值dB);六、实验报告要求1报告中要给出实验的MATLAB程序，了解每个语句作用；2简述实验目的和原理；3按实验步骤附上所设计滤波器的H(z)及相应的幅频特性曲线定性分析得到的图形，判断设计是否满足要求；4总结双线性变换法和脉冲不变法的特点和区别；5收获与建议。实验四 FIR数字滤波器的设计一、实验目的1掌握用窗函数法设计FIR数字滤波器的原理和方法；2熟悉线性相位FIR滤波器的幅频特性和相频特性；3了解不同窗函数对滤波器性能的影响。二、实验内容使用MATLAB编写程序，实现FIR数字滤波器的设计。涉及窗函数法和频率采样法设计FIR数字滤波器的方法、线性相位FIR滤波器的幅频特性和相频特性的特点、窗函数选择及其对滤波器性能的影响等知识点。三、实验原理与方法和手段窗函数设计法是一种把一个长序列变成有限长的短序列的设计方法，是在时域进行的。用窗函数法设计FIR数字滤波器时，先根据W c和N 求出相应的的理想滤波器的单位脉冲响应hd(n)。因为hd(n)一般是非因果的，且无限长，物理上是不可实现的。为此可选择适当的窗函数w(n)截取有限长的hd(n)， 即h(n)= hd(n)w(n)，只要阶数足够长，截取的方法合理，总能够满足频域的要求。实际中常用的窗函数有矩形（Boxcar）窗、三角（Bartlett）窗、汉宁(Hanning)窗、汉明(Hamming) 窗和布莱克曼（Blackman）窗。这些窗函数各有优缺点,所以要根据实际情况合理选择窗函数类型。1窗函数法设计线性相位FIR滤波器的一般步骤为：（1）确定理想滤波器的特性；（2）由Hd（ejw）求出；（3）根据过渡带宽度和阻带最小衰减，借助窗函数确定窗的形式及N的大小，即选择适当的窗函数，并根据线性相位条件确定窗函数的长度N；在MATLAB中，可由w=boxcar(N)（矩形窗）、w=hanning(N)（汉宁窗）、w=hamming(N)（汉明窗）、w=Blackman(N)（布莱克曼窗）、w=Kaiser(N,beta)（凯塞窗）等函数来实现窗函数设计法中所需的窗函数。（4）由h(n)= .w(n), 0n N-1，得出单位脉冲响应h(n);（5）对h(n)作离散时间傅立叶变换，得到H()。2在MATLAB中，可以用b=fir1(N,Wn,ftype,taper) 等函数辅助设计FIR数字滤波器。N代表滤波器阶数；Wn代表滤波器的截止频率(归一化频率)，当设计带通和带阻滤波器时，Wn为双元素相量；ftype代表滤波器类型，如high高通，stop带阻等；taper为窗函数，默认为海明窗，窗函数实现需要用窗函数blackman, hamming,hanning chebwin, kaiser产生。四、实验条件具有WINDOWS 2000/XP操作系统的计算机一台，安装MATLAB软件五、实验步骤在“开始-程序”菜单中，找到MATLAB程序，运行启动；进入MATLAB后 ，在Command Window中输入自己编写的主程序，并执行；记录运行结果图形，作分析对比。具体步骤如下：1.用窗函数法设计一线性相位FIR低通滤波器，要求通带截止频率wc=，单位冲击响应h（n）的长度N21。绘出h(n)及其幅频响应曲线. 设计分析:理想滤波器单位冲击响应为；为了满足线性相位FIR滤波器条件h(n)= h(N-1-n)，要求a=10。%xh19N=21;wc=pi/3; %理想低通滤波器参数n=0:N-1;a=(N-1)/2;hdn=si