数字信号处理实验指导书.docx

1、数字信号处理实验指导书数字信号处理实验指导书王玉富 编著实验一 抽样定理一、实验题目抽样定理二、实验目的掌握MATLAB相关函数的运用；掌握用MATLAB软件来演示抽样定理；加深对抽样定理的理解；三、具体实验内容及相关结论模拟信号与抽样一般的MATLAB不能用来处理和分析模拟信号（具有符号工具箱的MATLAB除外），但是可以用近似理想抽样（远小于T）的方法处理。已知，求出并绘制傅立叶变换图。利用，可以用一个在范围之间的有限长信号近似。因为，当时，有，则可选这样就可以使用MATLAB解题了。程序Example1_1如下：% 模拟信号Dt=0.00005;t=-0.005:Dt:0.005;xa=

2、exp(-1000*abs(t);% 连续时间傅立叶变换Wmax=2*pi*2000;K=500;k=0:1:K;W=k*Wmax/K;Xa=xa*exp(-j*t*W)*Dt;Xa=real(Xa);W=-fliplr(W),W(2:501); % 频率从-Wmax到WmaxXa=fliplr(Xa),Xa(2:501); %Xa介于-Wmax到Wmax之间subplot(211),plot(t*1000,xa);xlabel(t毫秒);ylabel(xa(t)title(模拟信号)subplot(212),plot(W/(2*pi*1000),Xa*1000);xlabel(频率(kHz)

3、;ylabel(Xa(jW)*1000);title(连续时间傅立叶变换)实验结果如图1-1所示。图1-1 example1-1运行结果为了研究抽样结果，用两种不同的抽样频率进行抽样。抽样频率为，对抽样得到，求出并画出；抽样频率为，对抽样得到，求出并画出；因为的带宽为，两倍频率小于，所以不存在混叠现象。程序Example1_2_1如下：%模拟信号Dt=0.00005;t=-0.005:Dt:0.005;xa=exp(-1000*abs(t);% 连续时间傅立叶变换Ts=0.0002;n=-25:1:25;x=exp(-1000*abs(n*Ts);%离散傅立叶变换K=500;k=0:1:K;w

4、=pi*k/K;X=x*exp(-j*n*w);X=real(X);w=-fliplr(w),w(2:K+1); % 频率从-Wmax到WmaxX=fliplr(X),X(2:K+1); %Xa介于-Wmax到Wmax之间subplot(211),plot(t*1000,xa);xlabel(n毫秒);ylabel(x1(n);title(离散信号);hold onstem(n*Ts*1000,x);gtext(Ts=0.2毫秒);hold offsubplot(212),plot(w/pi,X);xlabel(频率(pi);ylabel(X1(w);title(离散时间傅立叶变换)实验结果如

5、图1-2所示。离散信号叠加在模拟信号图上，可以清楚地看出抽样过程，很显然不存在混叠现象。图1-2 example1_2-1运行结果因为的带宽为，两倍频率大于，所以存在混叠现象。因为抽样频率不同，实现时只需将Example1_2_1中的条件改为Ts=0.001； n=-5:1:5；gtext(Ts=1毫秒)就可以了。运行结果如图1-3所示，其形状不同于图1-2就是因为存在混叠现象。图1-3 混叠现象信号恢复理想低通滤波器的时域输出为：给出的是一无限阶内插系统，在实际中是无法实现的，因此只能用有限阶内插逼近。MATLAB提供了在相邻点内插的几种方法。函数：根据这个公式就可以用矩阵与向量的相乘运算实

6、现内插了。利用函数恢复出，程序Example1_3如下：%由x1恢复xa Ts=0.0002; n=-25:1:25; nTs=n*Ts; x1=exp(-1000*abs(nTs); subplot(221),stem(nTs,x1) title(x1); %恢复模拟信号 Dt=0.00005; t=-0.005:Dt:0.005; Fs=1/Ts; xa=x1*sinc(Fs*(ones(length(n),1)*t-nTs*ones(1,length(t); subplot(222),plot(t,xa) title(由x1恢复xa); %校验 error1=max(abs(xa-exp

7、(-1000*abs(t) %由x2恢复xa Ts=0.001; n=-5:1:5; nTs=n*Ts; x=exp(-1000*abs(nTs); subplot(223),stem(nTs,x) title(x2); %模拟信号恢复 Dt=0.00005; t=-0.005:Dt:0.005; Fs=1/Ts; xa=x*sinc(Fs*(ones(length(n),1)*t-nTs*ones(1,length(t); subplot(224),plot(t,xa) title(由x2恢复xa); %校验 error2=max(abs(xa-exp(-1000*abs(t)校验结果：er

8、ror1 =0.0363；error2 =0.1852，运行结果如图1-4所示。图1-4 sinc(x)函数恢复原始信号由程序结果可见，恢复信号与原始信号存在误差，这是由于采用了有限个样本逼近无限阶内插所造成的，对于恢复的，误差只有0.0363，应该说是效果不错的。而由样本恢复出，结果为error2 =0.1852，显然信号恢复效果比较差，这是因为抽样频率不满足奈奎斯特抽样定理，内插后存在混叠的缘故。图解法：图解法就是给定样本值，而后利用stairs函数画出阶梯形的模拟信号逼近，或者利用plot函数画出内插曲线图。程序Example1_4如下：Ts=0.0002;n=-25:1:25;nTs=

9、n*Ts;x=exp(-1000*abs(nTs);%利用stairs函数恢复模拟信号subplot(211),stairs(nTs*1000,x);xlabel(t毫秒);ylabel(xa(t)title(利用stairs函数恢复模拟信号);hold onstem(n*Ts*1000,x);hold off %画出样本值进行比较%利用plot函数恢复模拟信号subplot(212),plot(nTs*1000,x);xlabel(t毫秒);ylabel(xa(t)title(利用plot函数恢复模拟信号);hold onstem(n*Ts*1000,x);hold off %画出样本值进行

10、比较图1-5 图解法恢复原始信号程序运行结果如图1-5所示，由图可见，利用plot函数绘出的图形比用stairs函数绘出的图形好，只要抽样符合奈奎斯特抽样定理，则用图解法可以得到近似解。spline函数：这种方法是采用一组称为3次样条函数的分段多项式实现，3次样条多项式为， 其中是多项式系数，是由样本值的最小二乘分析获得（严格说来，这不是一个因果性运算），这在MATLAB中是很方便得出的。利用spline函数由恢复出。程序Example1_5如下：s=0.0002;n=-25:1:25;nTs=n*Ts;x=exp(-1000*abs(nTs);%模拟信号恢复Dt=0.00005;t=-0.0

11、05:Dt:0.005;xa=spline(nTs,x,t);plot(t,xa)title(由x1恢复xa)%校验error=max(abs(xa-exp(-1000*abs(t)图1-6 spline函数恢复原始信号程序运行结果：error =0.0317， 由图1-6可以看出，使用spline函数恢复模拟信号仍然存在误差，但其效果要好于其它方法。实验二 序列的傅立叶变换一、实验题目序列的傅立叶变换二、实验目的掌握相关的MATLAB函数；用MATLAB来实现序列的傅立叶变换；利用MATLAB进行序列的对称性分析。三、具体实验内容及相关结论已知，用图解法讨论其傅立叶变换的对称性。程序Exam

12、ple2_1如下：图2-1 程序Example2_1运行结果subplot(1,1,1)n=-5:5;x=(-0.9).n;k=-200:200;w=(pi/100)*k;X=x*(exp(-j*pi/100).(n*k);magX=abs(X);angX=angle(X);subplot(211),plot(w/pi,magX);axis(-2,2,0,15)xlabel(频率（以pi为单位）);ylabel(X的绝对值);title(幅度);subplot(2,1,2);plot(w/pi,angX);axis(-2,2,-4,4)xlabel(频率（以pi为单位）);ylabel(弧度/

13、pi)title(相角)程序运行结果如图2-1所示从中可以看出，对对称，因此，分析时只看到0到之间的模和相位就可以了。验证实序列的对称特性，令程序Example2_1如下：n=-5:10;x=sin(pi*n/2); k=-100:100;w=(pi/100)*k; %-pi到pi之间的频率X=x*(exp(-j*pi/100).(n*k); %x的DTFT%信号分解xe,xo,m=evenodd(x,n); %利用evenodd分解奇偶部XE=xe*(exp(-j*pi/100).(m*k); %xe的DTFTXO=xo*(exp(-j*pi/100).(m*k); %xo的DTFT%校验X

14、R=real(X); %X的实部error1=max(abs(XE-XR) %XE与XR是否相等XI=imag(X); %X的虚部error2=max(abs(XO-j*XI) %XO与XI是否相等%用图形验证subplot(221),plot(w/pi,XR);grid onaxis(-1,1,-2,2)xlabel(频率（以pi为单位）);ylabel(Re(X);title(X的实部);subplot(222),plot(w/pi,XI);grid onaxis(-1,1,-10,10)xlabel(频率（以pi为单位）);ylabel(Im(X);title(X的虚部)subplot(

15、223),plot(w/pi,real(XE);grid onaxis(-1,1,-2,2)xlabel(频率（以pi为单位）);ylabel(XE);title(x偶部的DTFT)subplot(224),plot(w/pi,imag(XO);grid onaxis(-1,1,-10,10)xlabel(频率（以pi为单位）);ylabel(XO);title(x奇部的DTFT)程序运行结果为error1 =6.4666e-014，error2 =6.4661e-014，如图2-2所示，可见的实部等于的DTFT，而的虚部等于的DTFT图2-2 程序Example2_2运行结果evenodd函

16、数function xe,xo,m=evenodd(x,n)%将实信号分解为奇偶两部分%检查序列是否为实序列if any(imag(x)=0) error(x不是实序列)endm=-fliplr(n); %fliplr:Flip matrix in left/right direction.m1=min(m,n);m2=max(m,n);m=m1:m2;nm=n(1)-m(1);n1=1:length(n);x1=zeros(1,length(m);x1(n1+nm)=x;x=x1;xe=0.5*(x+fliplr(x);xo=0.5*(x-fliplr(x);实验三 IIR数字滤波器设计一、

17、实验题目IIR数字滤波器设计二、实验目的利用MATLAB软件设计模拟低通滤波器；利用MATLAB软件实现冲激响应不变法设计IIR数字滤波器；利用MATLAB软件实现双线性变换法设计IIR数字滤波器；利用MATLAB软件实现IIR数字滤波器的频带转换。三、具体实验内容及相关结论模拟低通滤波器设计MATLAB提供了N阶归一化（即）巴特沃斯模拟滤波器设计函数z,p,k=buttap(N)，其返回值为数组z和p（零点和极点）以及增益k。MATLAB扩展函数提供的b,a=u_buttap函数提供了一个直接型滤波器结构，而sdir2cas给除了适合于模拟滤波器的直接型到级联型的转换。巴特沃斯低通模拟滤波器

18、的设计过程可以使用MATLAB扩展函数中的afd_butt实现，用freq_m函数显示滤波器的频域图。程序Example3_1如下：Wp=0.2*pi;Ws=0.3*pi;Rp=7;Ar=16; %设计指标b,a=afd_butt(Wp,Ws,Rp,Ar); %模拟滤波器设计wmax=0.5*pi;db,mag,pha,w=freqs_m(b,a,wmax); %计算频率响应和相位响应ha,x,t=impulse(b,a);C,B,A=sdir2cas(b,a) %直接型到级联型转换%画图figure(1)subplot(331)plot(mag)axis(0 500 0 1.05);titl

19、e(频率响应)subplot(332)plot(pha)axis(0 500 -3.2 3.2);title(相位响应)subplot(333)plot(ha)axis(0 100 -0.025 0.21);title(脉冲响应)滤波器响应如图3-1所示。图3-1 程序Example3_1运行结果常用的模拟滤波器设计方法除了巴特沃斯滤波器外，还有切比雪夫型、切比雪夫型、椭圆滤波器等。MATLAB提供了设计N阶归一化切比雪夫型滤波器的cheb1ap函数，设计N阶归一化切比雪夫型滤波器的cheb2ap函数，设计N阶归一化椭圆滤波器的ellipap函数。MATLAB扩展函数中提供了设计N阶未归一化切

20、比雪夫型滤波器的u_cheb1ap函数和根据给定指标设计切比雪夫型滤波器的afd_cheb1函数；设计N阶未归一化切比雪夫型滤波器的u_cheb2ap函数和根据给定指标设计切比雪夫滤型波器的afd_cheb2ap函数；设计N阶椭圆模拟滤波器的u_elipap函数和根据给定指标设计椭圆模拟滤波器的afd_elip函数。afd_butt函数function b,a=afd_butt(Wp,Ws,Rp,Ar)%按给定的技术指标设计模拟巴特沃思低通滤波器%b直接型滤波器分子系数；a直接型滤波器分母系数%Wp通带截止频率，单位弧度/秒；Ws阻带截止频率，单位弧度/秒；%Rp通带中的振幅波动系数（dB）；

21、Ar阻带衰减，（dB）；if Wp=0 error(通带边缘必须大于0)endif Ws=Wp error(阻带边缘必须大于通带边缘)endif (Rp=0)|(Ar1时计算 A=zeros(K+1,3); for n=1:2:2*K Arow=p(n:1:n+1,:); Arow=poly(Arow); A(fix(n+1)/2),:)=real(Arow); end A(K+1,:)=0 real(poly(p(Na);end%分子的二阶因子部分z=cplxpair(roots(b);K=floor(Nb/2);if Nb=0 %当Nb=0时计算 B=0 0 poly(z);elseif

22、K*2=Nb %当Nb为偶时计算 B=zeros(K,3); for n=1:2:Nb Brow=z(n:1:n+1,:); Brow=poly(Brow); B(fix(n+1)/2),:)=real(Brow); endelseif Nb=1 %当Nb=1时计算 B=0 real(poly(z);else %当Nb为奇且1时计算 B=zeros(K+1,3); for n=1:2:2*K Brow=z(n:1:n+1,:); Brow=poly(Brow); B(fix(n+1)/2),:)=real(Brow); end B(K+1,:)=0 real(poly(z(Nb);end脉冲响

23、应不变法利用MATLAB实现脉冲响应不变法，即已知系统函数，利用residue函数得到系统函数零极点，然后根据将每一个模拟极点映射为数字极点。最后利用residuez函数将转换成滤波器直接形式，这一过程用imp_invr函数实现。滤波器设计指标为，用3种方法设计数字滤波器。程序Example3_2如下：Wp=0.2*pi;Ws=0.3*pi;Rp=1;Ar=15;T=1; %数字滤波器设计指标OmegaP=Wp*T; %模拟原型滤波器通带截止频率OmegaS=Ws*T; %模拟原型滤波器阻带截止频率%cs1,ds1=afd_butt(OmegaP,OmegaS,Rp,Ar); %巴特沃思模拟滤

24、波器设计b1,a1=imp_invr(cs1,ds1,T); %用脉冲响应不变法实现模拟滤波器到数字滤波器的变换db1,mag1,pha1,grd1,w1=freqz_m(b1,a1); %计算频率响应和相位响应C1,B1,A1=dir2par(b1,a1) %直接型到并联型转换%画图figure(1)subplot(341)plot(w1/pi,mag1);axis(0 1 0 1.1);ylabel(a)title(幅度响应);subplot(342);plot(w1/pi,db1);axis(0 1 -40 5);title(幅度响应（dB）);subplot(343);plot(w1/

25、pi,pha1/pi);axis(0 1 -1 1);title(相位响应);subplot(344);plot(w1/pi,grd1);axis(0 1 0 10);title(群延迟);%cs2,ds2=afd_chb1(OmegaP,OmegaS,Rp,Ar); %未归一化切比雪夫型模拟低通滤波器设计b2,a2=imp_invr(cs2,ds2,T); %用脉冲响应不变法实现模拟滤波器到数字滤波器的变换db2,mag2,pha2,grd2,w2=freqz_m(b2,a2); %计算频率响应和相位响应C2,B2,A2=dir2par(b2,a2) %直接型到并联型转换%画图subplot(345)plot(w2/pi,mag2);axis(0 1 0 1.1);ylabel(b)subplot(346);plot(w2/pi,db2);axis(0 1 -40 5);subplot(347);plot(w2/pi,pha2/pi);axis(0 1 -1 1