数字信号处理实验指导书2016-通信

1、数字信号处理实验徐俊2016年8月请浏览后下载，资料供参考，期待您的好评与关注！目 录实验一离散时间系统及系统响应2实验二离散傅立叶变换及其特性验证8实验三时域采样与频域采样17实验四冲激响应不变法IIR数字滤波器设计24 实验一 离散时间系统及系统响应一、实验目的1、掌握求解离散时间系统脉冲响应和阶跃响应的方法；2、掌握用线性卷积求解离散时间系统响应的基本方法。二、实验原理与设计方法1、用impz和dstep函数求解离散系统的单位脉冲响应和阶跃响应【例1-1】已知某因果系统的差分方程为yn+0.5yn-1=xn+2x(n-2)系统为零状态，求系统的脉冲响应和阶跃响应。解：该系统是一个2阶系统

2、，列出bm和ak系数为a0=1，a1=0.5，a2=0，b0=1，b1=0，b2=2MALAB程序如下（取16点作图）：a=1,0.5,0;b=1,0,2;n=16;hn=impz(b,a,n); %脉冲响应gn=dstep(b,a,n); %阶跃响应subplot(1,2,1),stem(hn,k);title(系统的单位脉冲响应);ylabel(h(n);xlabel(n);axis(0,n,1.1*min(hn),1.1*max(hn);subplot(1,2,2),stem(gn,k);title(系统的单位阶跃响应);ylabel(g(n);xlabel(n);axis(0,n,1.

3、1*min(gn),1.1*max(gn);结果如下图所示：2、用conv函数进行卷积计算求系统响应【例1-2】某离散时间系统的脉冲响应为hb(n)=(n)+2.5(n-1)+2.5(n-2)+(n-3)激励信号为xt=Ae-nTsin0nT 0n50设A=444.128，=502，0=502。试求该系统在输入信号激励下的响应。解：MATLAB程序如下：n=1:50; %定义序列的长度是50hb=zeros(1,50); %注意：MATLAB中数组下标从1开始hb(1)=1;hb(2)=2.5;hb(3)=2.5;hb(4)=1;close all;subplot(3,1,1);stem(hb

4、);title(系统hn);m=1:50;T=0.001; %定义序列的长度和采样率A=444.128;a=50*sqrt(2.0)*pi; %设置信号有关的参数w0=50*sqrt(2.0)*pi;x=A*exp(-a*m*T).*sin(w0*m*T);subplot(3,1,2);stem(x);title(输入信号xn);y=conv(x,hb);subplot(3,1,3);stem(y);title(输出信号yn);结果如下图所示：3、用filter函数求系统响应线性常系数差分方程是描述离散时间LTI系统的另一个时域模型，即系统的输入信号xn输出信号yn关系可以用下面的差分方程来表

5、达 为了求得系统响应信号的显式表达式（Explicit expression），必须对差分方程求解。在MATLAB中，是用差分方程的系数来描述一个离散时间LTI系统的。例如，一个LTI离散时间系统的差分方程为 yn+yn-1-6yn-2=xnMATLAB则用两个系数向量num = 1和den = 1 1 -6来描述该系统，其中num和den分别表示系统差分方程右边和左边的系数，按照差分运算的递减排列。在用num和den定义了一个离散时间LTI系统之后，可以利用MATLAB来求解系统响应。求解离散时间系统的响应的一个非常有用的MATLAB函数就是filter()。它可以用来求解系统的在输入信号x

6、n作用下的零状态响应，也可以求解系统在这个输入信号作用下的完全响应。其用法描述如下：y = filter(num, den, x)，求系统在输入x作用下的零状态响应y，x与y的长度相同。y = filter(num, den, x, ic)，求系统在输入x作用下的完全响应y，x与y的长度相同。其中ic为系统的初始条件（Initial Condition），即ic = y-1, y-2, y-3, ., y-N。【例1-3】给定一个离散时间LTI系统，它的差分方程为yn+0.1yn-1-0.06yn-2=xn若输入信号为xn = 0.5nun，系统的初始条件为ic = 0, 1，编写程序，绘制输

7、入信号xn、系统的零状态响应yzsn和系统的完全响应信号yn的波形图。MATLAB程序如下：clear,close all,num = 1; den = 1 .1 -0.06;n = 0:20; x = 0.5.n.*ones(1,21); ic = 0 .9;yzs = filter(num,den,x); y = filter(num,den,x,ic); subplot(2,2,1)stem(n,x,.), title(The input sequence xn), axis(0,20,-0.5,1.5)subplot(2,2,2)stem(n,yzs,.), title(The zer

8、o-state response yzsn), axis(0,20,-0.5,1.5), xlabel(Time index n)subplot(2,2,3)stem(n,y,.), axis(0,20,-0.5,1.5), title(The total solution yn)xlabel(Time index n)结果如下图所示：三、实验内容1、分别用impz和dstep函数求解下面离散时间系统的脉冲响应和阶跃响应。（1）系统的差分方程为y(n)=0.8y(n-1)-0.64y(n-2)+0.866x(n)（2）系统的系统函数为Hz=1-0.5z-11-z-1+z-22、利用第1题求得的

9、系统的脉冲响应求解系统在激励x(n)=u(n-3)下的响应。3、利用filter函数求解第1题系统在激励x(n)=u(n-3)下的响应。四、实验参考MATLAB基础操作1、 矩阵输入Matlab有很强的数值矩阵处理能力。实际上，矩阵和矢量这两个词经常换用。矩阵是以实数或者复数为元素的长方形矢量。在输入矩阵时，应遵循下列规则：l 整个矩阵用中括号括起来；l 同一行的数据用空格或逗号隔开；l 不同行用分号隔开。在Matlab中，列矢量可被当作只有一列的矩阵；行矢量也可被当作是只有一行的矩阵；标量可被当作只有一列的矩阵。在MATLAB中，用一个列向量来表示一个有限长序列，由于一个列向量并不包含位置信

10、息，因此需要用表示位置的n和表示量值的x两个向量来表示任意一个序列，如：【例】n=-3,-2,-1,0,1,2,3,4;x=2,1,-1,0,1,4,3,7;plot(n,x);运行结果：如果不对向量的位置进行定义，则MATLAB默认该序列的起始位置为n=0。由于内存有限，MATLAB不能表示一个无限序列。2、 利用函数输入矩阵例如：zeros：生成一个元素全部为0的矩阵；ones：生成一个元素全部为1的矩阵；eye：生成一个单位矩阵。A=zeros(2,3) ：生成一个2行3列的全0矩阵；A= 0 0 0 0 0 03、 利用plot，stem进行波形绘制l 绘图命令plot绘制x-y坐标图

11、;l 在绘图过程中，经常要把几个图形在同一个图形窗口中表现出来，而不是简单地叠加这就用到函数subplot其调用格式如下：subplot(m,n,p)。subplot函数把一个图形窗口分割成mn个子区域，用户可以通过参数p调用个各子绘图区域进行操作子绘图区域的编号为按行从左至右编号l stem绘制离散序列图，常用格式stem(y)和stem(x,y)分别和相应的plot函数的绘图规则相同，只是用stem命令绘制的是离散序列图。实验二 离散傅立叶变换及其特性验证一、实验目的1、掌握离散时间傅立叶变换（DTFT）的计算方法和编程技术。2、掌握离散傅立叶变换（DFT）的计算方法和编程技术。3、掌握用

12、FFT对信号作频谱分析。二、实验原理与设计方法1、离散时间傅立叶变换如果序列x(n)满足绝对可和的条件，即，则其离散时间傅立叶变换定义为： （1）如果x(n)是无限长的，则不能直接用MATLAB由x(n)计算X(ejw)，但可以用它来估计X(ejw)表达式在0,频率区间的值并绘制它的幅频和相频（或实部和虚部）曲线。如果x(n)是有限长的，则可以用MATLAB对任意频率w处的X(ejw)进行数值计算。如果要在0,间按等间隔频点估计X(ejw)，则（1）式可以用矩阵向量相乘的运算来实现。假设序列x(n)在(即不一定在0, N-1)有N个样本，要估计下列各点上的X(ejw)：它们是0,之间的(M+1

13、)个等间隔频点，则（1）式可写成： （2）将x(nl)和X(ejwk)分别排列成向量x和X，则有： X=Wx （3）其中W是一个(M+1)N维矩阵：将k和n排成列向量，则在MATLAB中，把序列和下标排成行向量，对（3）式取转置得：其中nTk是一个N(M+1)维矩阵。用MATLAB实现如下：k=0:M; n=n1:n2;X=x*(exp(-j*pi/M).(n*k);2、离散傅立叶变换一个有限长序列的离散傅立叶变换对定义为： （4） （5）以列向量x和X形式排列x(n)和X(k)，则式（4）、（5）可写成：X=WNx其中矩阵WN由下式给出：可由下面的MATLAB函数dft和idft实现离散傅立

14、叶变换运算。function Xk = dft(xn,N)% Computes Discrete Fourier Transform% -% Xk = dft(xn,N)% Xk = DFT coeff. array over 0 = k = N-1% xn = N-point finite-duration sequence% N = Length of DFT%n = 0:1:N-1; % row vector for nk = 0:1:N-1; % row vecor for kWN = exp(-j*2*pi/N); % Wn factornk = n*k; % creates a N

15、 by N matrix of nk valuesWNnk = WN . nk; % DFT matrixXk = xn * WNnk; % row vector for DFT coefficientsfunction xn = idft(Xk,N)% Computes Inverse Discrete Transform% -% xn = idft(Xk,N)% xn = N-point sequence over 0 = n = N-1% Xk = DFT coeff. array over 0 = k N，则直接截取R点DFT的前N点，若R2fm、fs=2fm、fs2fm三种情况下采样

16、信号的波形以及采样信号的频谱。（2） 用内插公式实现原始信号的恢复。解：MATLAB程序如下：dt=0.01;f0=1;T0=1/f0;%输入基波频率、周期fm=5*f0;Tm=1/fm;%最高频率为基波频率的5倍t=-2:dt:2;f=sin(2*pi*f0*t)+1/3*sin(6*pi*f0*t);%建立原连续信号subplot(4,1,1),plot(t,f,k);axis(min(t) max(t) 1.1*min(f) 1.1*max(f);title(原连续信号和采样信号);for i=1:3 fs=i*fm;Ts=1/fs;%确定采样频率和周期 n=-2:Ts:2; f=sin

17、(2*pi*f0*n)+1/3*sin(6*pi*f0*n);%生成采样信号 subplot(4,1,i+1);stem(n,f,filled,k); axis(min(t) max(t) 1.1*min(f) 1.1*max(f);endfigureN=length(t);%时间轴上采样点数f=sin(2*pi*f0*t)+1/3*sin(6*pi*f0*t);wm=2*pi*fm;k=0:N-1;w1=k*wm/N;%频率轴上采样点F1=f*exp(-j*t*w1)*dt;%对原始信号进行傅里叶变换subplot(4,1,1),plot(w1/(2*pi),abs(F1),k);axis(

18、0 max(4*fm) 1.1*min(abs(F1) 1.1*max(abs(F1);title(原连续信号和采样信号的频谱);for i=1:3 if i=2 c=0,else c=1,end fs=(i+c)*fm;Ts=1/fs; n=-2:Ts:2; f=sin(2*pi*f0*n)+1/3*sin(6*pi*f0*n); N=length(n); wm=2*pi*fs; k=0:N-1; w=k*wm/N; F=f*exp(-j*n*w)*Ts; subplot(4,1,i+1);plot(w/(2*pi),abs(F),k); axis(0 max(4*fm) 1.1*min(a

19、bs(F) 1.1*max(abs(F);endt=0:dt:3*T0;x= sin(2*pi*f0*t)+1/3*sin(6*pi*f0*t);figure;subplot(4,1,1),plot(t,x,k);axis(min(t) max(t) 1.1*min(x) 1.1*max(x);title(用内插公式重建采样信号);for i=1:3fs=i*fm;Ts=1/fs;n=0:(3*T0)/Ts;t1=0:Ts:3*T0;x1= sin(2*pi*n*f0/fs)+1/3*sin(6*pi*n*f0/fs);%生成采样信号T_N=ones(length(n),1)*t1-n*Ts*

20、ones(1.length(t1);xa=x1*sinc(fs*pi*T_N);%内插subplot(4,1,i+1),plot(t1,xa,k);axis(min(t1) max(t1) 1.1*min(xa) 1.1*max(xa);end运行结果：2、频域采样在单位圆上对任意序列x(n)的z变换X(z)进行N点的等间隔采样可得Xk=X(z)z=ej2Nnk k=0,1,N-1上式实现了序列在频域的采样。那么，由频域采样获得的频谱序列能否和连续信号时域采样一样，在一定的条件下能不失真的恢复出原离散序列呢？频域采样定理给出了答案。频域采样定理由下列公式表述xn=r=-x(n+rN)上式表明，

21、对一个频谱采样后经IDTF生成的周期序列xn是原非周期序列x(n)的额周期延拓序列，其时域周期等于频域采样点数N。假设有限长序列x(n)的长度是M，其频域采样点数为N，则原时域离散信号能够不失真的有频域采样恢复的条件如下。如果x(n)为无限长序列，则必然造成混叠现象，无法无失真恢复。如果x(n)为有限长序列，且频域采样点数N小于序列长度M，则x(n)以N为周期进行延拓也将产生混叠，从xn中不能无失真的恢复出x(n)。如果x(n)为有限长序列，且频域采样点数N大于或等于序列长度M，则x(n)以N为周期进行延拓也将产生混叠，从xn中能够无失真的恢复出x(n)，即xNn=xnRNn=r=-xn+rN

22、RNn=x(n)【例3-2】已知序列x(n)=1,1,1;n=0,1,2，对其频谱进行采样，分别取N=2、3、5，观察频域采样造成的混叠现象。解：MATLAB程序如下：Ts=1;N1=2;D1=2*pi/(Ts*N1);kn1=floor(-(N1-1)/2:-1/2);kp1=floor(0:(N1-1)/2);w1=kp1,kn1*D1;X1=1+1*exp(-j*w1)+exp(-j*2*w1);%2点频谱采样结果n=0:N1-1;x1=abs(ifft(X1,N1);%由2点频谱采样结果恢复的时域信号figuresubplot(1,3,1);stem(n*Ts,x1,filled,k)

23、;N2=3;D2=2*pi/(Ts*N2);kn2=floor(-(N2-1)/2:-1/2);kp2=floor(0:(N2-1)/2);w2=kp2,kn2*D2;X2=1+1*exp(-j*w2)+exp(-j*2*w2);%3点频谱采样结果n=0:N2-1;x2=abs(ifft(X2,N2);%由3点频谱采样结果恢复的时域信号subplot(1,3,2);stem(n*Ts,x2,filled,k);N3=5;D3=2*pi/(Ts*N3);kn3=floor(-(N3-1)/2:-1/2);kp3=floor(0:(N3-1)/2);w3=kp3,kn3*D3;X3=1+1*exp

24、(-j*w3)+exp(-j*2*w3);%5点频谱采样结果n=0:N3-1;x3=abs(ifft(X3,N3);%由5点频谱采样结果恢复的时域信号subplot(1,3,3);stem(n*Ts,x3,filled,k);运行结果：三、 实验内容1、已知连续信号xt=sinc(t)，取最高有限带宽频率fm=1Hz。（1）分别显示原信号波形和fs=fm、fs=2fm、fs=3fm三种情况下采样信号的波形。（2）求解原连续信号和采样信号所对应的频谱。（3）用内插公式重建信号。2、已知一个时间序列x(n)=2,4,6,4,2;n=0,1,2,3,4，先求其频谱，再分别取频域采样点数N=3、5、1

25、0，用IFFT计算并求出其时间序列x(n)，用图形显示各时间序列。观察时域序列是否存在混叠，有何规律？实验四 冲激响应不变法IIR数字滤波器设计一、实验目的1、掌握构成一个频率响应与给定的滤波特性相接近的模拟滤波器的设计原理。2、掌握用冲激响应不变法设计IIR数字滤波器的基本原理和算法。3、了解数字滤波器和模拟滤波器的频率响应特性，掌握相应的计算方法，分析用冲激响应不变法获得的数字滤波器频率响应特性中出现的混叠现象。二、实验原理与设计方法、冲激响应不变法设计IIR数字滤波器的基本原理和算法采用冲激响应不变法设计数字滤波器，就是使其单位样值响应与相应的模拟滤波器的冲激响应在抽样点处的量值相等，即

26、 (1)其中T为抽样周期。因此用冲激响应不变法设计IIR数字滤波器的基本步骤，就是首先根据设计要求确定相应的模拟滤波器的传递函数，经Laplace反变换求出冲激响应后，对它进行抽样得到的等于数字滤波器的单位样值响应，再经z变换所得就是数字滤波器的传递函数。如果模拟滤波器的传递函数的N个极点都是单极点，则可以将写成部分分式展开的形式 (2)那么，经Laplace反变换求出的模拟滤波器的冲激响应为相对应的数字滤波器的单位抽样响应为对上式作z变换，得 (3)由上面的推导可见，只要模拟滤波器的传递函数的N个极点都是单极点，当已经求出各个极点值和部分分式的系数后，则可以从模拟滤波器的传递函数的表达式（2

27、）直接得到数字滤波器的传递函数的表达式（3）。冲激响应不变法设计数字滤波器的步骤：（1） 设数字滤波器的技术指标为wp、ws、Rp、Rs，选择采用周期T，将wp、ws、Rp、Rs转换为模拟滤波器参数p、s、Rp、Rs。其中p=wp/T，s= ws/T。（2） 选择一种模拟低通滤波器原型（巴特沃斯、切比雪夫）设计模拟低通滤波器。（3） 用冲激响应不变法转换为数字滤波器。、Butterworth模拟滤波器的设计方法（1）Butterworh原型MATLAB提供了函数z,p,k=buttap(N)用来设计N阶归一化（即c=1）的Butterworth模拟低通滤波器；其中N为巴特沃斯低通滤波器的阶数，

28、输出值z、p、k分别对应零点数组、极点数组及增益数组。如果c1，则需去归一化处理，极点数组变为cp，零点数组变为cz，分子也要乘以cN（k=cNk）（2）按给定技术指标设计Butterworth模拟低通滤波器N,Wn=buttord(Wp,Ws,Rp,Rs,s)函数实现指定参数的Butterworth模拟低通滤波器的设计。其中输入值Wp、Ws为通带及阻带截止频率，单位为rad/s，Rp,Rs为通带最大衰减及阻带最小衰减，单位为dB，s说明设计的是模拟滤波器。输出值N为滤波器阶数，Wn为滤波器的截止频率。【例4-1】设计巴特沃斯模拟低通滤波器：通带截止频率fp=10Hz，通带最大衰减Rp=2dB

29、，阻带截止频率fs=20Hz，阻带最小衰减Rs=12dB。要求：（1） 求滤波器的阶次N及截止频率Wn；（2） 画出系统的幅频特性。解：fp=10;fs=20;Rp=2;Rs=12;Wp=2*pi*fp;Ws=2*pi*fs;N,Wn=buttord(Wp,Ws,Rp,Rs,s);z0,p0,k0=buttap(N); %计算归一化参数p=Wn*p0; %计算极点z=Wn*z0; %计算零点k=WnN*k0;%计算增益b=real(poly(z);b=k*b;%计算零极点所对应的系统函数的分子的系数a=real(poly(p);%计算零极点所对应的系统函数的分子母的系数H,w=freqs(b,

30、a);plot(w/(2*pi),20*log10(abs(H)/max(abs(H); %画系统的幅频特性图title(20*log10(abs(H)/max(abs(H);xlabel(w/Hz);ylabel(20*log10(abs(H);grid on;axis(0 60 -30 0);运行结果：3、切比雪夫型（Chebyshev）模拟滤波器的设计方法（1）切比雪夫型原型MATLAB提供了函数z,p,k= Cheb1ap (N,Rp)用来设计N阶归一化（即c=1）的切比雪夫型模拟低通滤波器；其中N为切比雪夫型低通滤波器的阶数，Rp为通带波动值，输出值z、p、k分别对应零点数组、极点数

31、组及增益数组。如果c1，则需去归一化处理，极点数组变为cp，零点数组变为cz，分子也要乘以cN（k=cN*k）（2）按给定技术指标设计切比雪夫型模拟低通滤波器N,Wn= Cheb1 ord(Wp,Ws,Rp,Rs,s)函数实现指定参数的切比雪夫型模拟低通滤波器的设计。其中输入值Wp、Ws为通带及阻带截止频率，单位为rad/s，Rp,Rs为通带最大衰减及阻带最小衰减，单位为dB，s说明设计的是模拟滤波器。输出值N为滤波器阶数，Wn为滤波器的截止频率。【例4-2】设计一个设计切比雪夫型模拟低通滤波器：通带截止频率fp=10Hz，通带波纹Rp=2dB，阻带截止频率fs=20Hz，阻带最小衰减Rs=1

32、2dB。要求：（1） 求滤波器的阶次N及截止频率Wn；（2） 画出系统的幅频特性。解：fp=10;fs=20;Rp=2;Rs=12;Wp=2*pi*fp;Ws=2*pi*fs;N,Wn=cheb1ord (Wp,Ws,Rp,Rs,s);z0,p0,k0=cheb1ap(N,Rp); %计算归一化参数p=Wn*p0; %计算极点z=Wn*z0; %计算零点k=WnN*k0;%计算增益b=real(poly(z);b=k*b;%计算零极点所对应的系统函数的分子的系数a=real(poly(p);%计算零极点所对应的系统函数的分子母的系数H,w=freqs(b,a);plot(w/(2*pi),20

33、*log10(abs(H)/max(abs(H); %画系统的幅频特性图title(20*log10(abs(H)/max(abs(H);xlabel(w/Hz);ylabel(20*log10(abs(H);grid on;axis(0 60 -30 0);运行结果：4、模拟滤波器转换为数字滤波器函数bz,az=impinvar(bs,as,fs)用来实现冲激响应不变的映射，bs、as分别为模拟滤波器系统函数Ha(s)的分子、分母的系数；fs为抽样频率；bz、az为数字滤波器的分子、分母向量的系数。【例4-3】用巴特沃斯滤波器原型设计一个低通数字滤波器，满足Wp=0.6，Ws=1.2，Rp=2dB，Rs=18dB滤波器采样频率为0.5Hz。采用冲激响应不变法进行变换。要求：（1） 写出滤波器阶次N及数字频率；（2） 模拟滤波器幅频特性；（3） 数字滤波器幅频特性。解：wp=0.6;ws=1.2;fs=0.5;T=1/fs;omegap=wp*fs;omegas=ws*fs;Rp=2;Rs=18;N,omegan=buttord(omegap,omegas,Rp,Rs,s);wn=omegan*T;z0,p0,k0=buttap