数字信号处理实验徐俊2016年8月目录实验一离散时间系统及系统响应 (2)实验二离散傅立叶变换及其特性验证 (8)实验三时域采样与频域采样 (17)实验四冲激响应不变法IIR数字滤波器设计 (24)实验一离散时间系统及系统响应一、实验目的1、掌握求解离散时间系统脉冲响应和阶跃响应的方法；2、掌握用线性卷积求解离散时间系统响应的基本方法。

二、实验原理与设计方法1、用impz和dstep函数求解离散系统的单位脉冲响应和阶跃响应【例1-1】已知某因果系统的差分方程为系统为零状态，求系统的脉冲响应和阶跃响应。

解：该系统是一个2阶系统，列出b m和a k系数为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.1*min(gn),1.1*max(gn)]);结果如下图所示：2、用conv函数进行卷积计算求系统响应【例1-2】某离散时间系统的脉冲响应为h b(n)=δ(n)+2.5δ(n-1)+2.5δ(n-2)+δ(n-3)激励信号为设A=444.128，α=50，。

试求该系统在输入信号激励下的响应。

解：MATLAB程序如下：n=1:50; %定义序列的长度是50hb=zeros(1,50); %注意：MA TLAB中数组下标从1开始hb(1)=1;hb(2)=2.5;hb(3)=2.5;hb(4)=1;close all;subplot(3,1,1);stem(hb);title('系统h[n]');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('输入信号x[n]');y=conv(x,hb);subplot(3,1,3);stem(y);title('输出信号y[n]');结果如下图所示：3、用filter 函数求系统响应线性常系数差分方程是描述离散时间LTI 系统的另一个时域模型，即系统的输入信号x[n]输出信号y[n]关系可以用下面的差分方程来表达∑∑==-=-Mk k Nk kk n x p k n y d][][为了求得系统响应信号的显式表达式（Explicit expression ），必须对差分方程求解。

在MATLAB 中，是用差分方程的系数来描述一个离散时间LTI 系统的。

例如，一个LTI 离散时间系统的差分方程为y[n]+y[n-1]-6y[n-2]=x[n]MATLAB 则用两个系数向量num = [1]和den = [1 1 -6]来描述该系统，其中num 和den 分别表示系统差分方程右边和左边的系数，按照差分运算的递减排列。

在用num 和den 定义了一个离散时间LTI 系统之后，可以利用MA TLAB 来求解系统响应。

求解离散时间系统的响应的一个非常有用的MATLAB 函数就是filter()。

它可以用来求解系统的在输入信号x[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 系统，它的差分方程为y[n]+0.1y[n-1]-0.06y[n-2]=x[n]若输入信号为x[n] = 0.5n u[n]，系统的初始条件为ic = [0, 1]，编写程序，绘制输入信号x[n]、系统的零状态响应yzs[n]和系统的完全响应信号y[n]的波形图。

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 x[n]'), axis([0,20,-0.5,1.5])subplot(2,2,2)stem(n,yzs,'.'), title('The zero-state response yzs[n]'), 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 y[n]')xlabel('Time index n')结果如下图所示：三、实验内容1、分别用impz和dstep函数求解下面离散时间系统的脉冲响应和阶跃响应。

（1）系统的差分方程为y(n)=0.8y(n-1)-0.64y(n-2)+0.866x(n)（2）系统的系统函数为2、利用第1题求得的系统的脉冲响应求解系统在激励x(n)=u(n-3)下的响应。

3、利用filter函数求解第1题系统在激励x(n)=u(n-3)下的响应。

四、实验参考MATLAB基础操作1、矩阵输入Matlab有很强的数值矩阵处理能力。

实际上，矩阵和矢量这两个词经常换用。

矩阵是以实数或者复数为元素的长方形矢量。

在输入矩阵时，应遵循下列规则：●整个矩阵用中括号括起来；●同一行的数据用空格或逗号隔开；●不同行用分号隔开。

在Matlab中，列矢量可被当作只有一列的矩阵；行矢量也可被当作是只有一行的矩阵；标量可被当作只有一列的矩阵。

在MATLAB中，用一个列向量来表示一个有限长序列，由于一个列向量并不包含位置信息，因此需要用表示位置的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。

由于内存有限，MA TLAB不能表示一个无限序列。

2、利用函数输入矩阵例如：zeros：生成一个元素全部为0的矩阵；ones：生成一个元素全部为1的矩阵；eye：生成一个单位矩阵。

A=zeros(2,3) ：生成一个2行3列的全0矩阵；A=0 0 00 0 03、利用plot，stem进行波形绘制●绘图命令plot绘制x-y坐标图;●在绘图过程中，经常要把几个图形在同一个图形窗口中表现出来，而不是简单地叠加．这就用到函数subplot．其调用格式如下：subplot(m,n,p)。

subplot函数把一个图形窗口分割成m×n个子区域，用户可以通过参数p调用个各子绘图区域进行操作．子绘图区域的编号为按行从左至右编号．●stem绘制离散序列图，常用格式stem(y)和stem(x,y)分别和相应的plot函数的绘图规则相同，只是用stem命令绘制的是离散序列图。

实验二 离散傅立叶变换及其特性验证一、实验目的1、掌握离散时间傅立叶变换（DTFT ）的计算方法和编程技术。

2、掌握离散傅立叶变换（DFT ）的计算方法和编程技术。

3、掌握用FFT 对信号作频谱分析。

二、实验原理与设计方法1、离散时间傅立叶变换如果序列x (n )满足绝对可和的条件，即∞<∑∞-∞=n n x |)(|，则其离散时间傅立叶变换定义为：∑∞-∞=-==n nj j en x n x F e X ωω)()]([)( （1）如果x (n )是无限长的，则不能直接用MATLAB 由x (n )计算X (e j ω)，但可以用它来估计X (e j ω)表达式在[0,π]频率区间的值并绘制它的幅频和相频（或实部和虚部）曲线。

如果x (n )是有限长的，则可以用MATLAB 对任意频率ω处的X (e j ω)进行数值计算。

如果要在[0,π]间按等间隔频点估计X (e j ω)，则（1）式可以用矩阵－向量相乘的运算来实现。

假设序列x (n )在N n n n ≤≤1(即不一定在[0, N -1])有N 个样本，要估计下列各点上的X (e j ω)：M k k Mk ...,2,1,0==，　πω它们是[0,π]之间的(M +1)个等间隔频点，则（1）式可写成： M k n x ee X Nl l kn Mjj l...,2,1,0)()(1==∑=-，　πω （2）将{x (n l )}和{X (e j ωk )}分别排列成向量x 和X ，则有：X=Wx （3） 其中W 是一个(M +1)×N 维矩阵：⎭⎬⎫⎩⎨⎧=≤≤=-M k n n n e N kn M j ...,2,1,0;1，　πW将{k }和{n }排成列向量，则⎥⎦⎤⎢⎣⎡⎪⎭⎫ ⎝⎛-=n k W T M j πexp 在MA TLAB 中，把序列和下标排成行向量，对（3）式取转置得：⎥⎦⎤⎢⎣⎡⎪⎭⎫ ⎝⎛-=k n x X T T T M j πexp其中n T k 是一个N×(M +1)维矩阵。

用MATLAB 实现如下： k=[0:M]; n=[n1:n2]; X=x*(exp(-j*pi/M)).^(n’*k); 2、离散傅立叶变换一个有限长序列的离散傅立叶变换对定义为：10,)()(10-≤≤=∑-=N k W n x k X N n nk N （4） 10,)(1)(1-≤≤=∑-=-N n Wk X Nn x N k kn N （5）以列向量x 和X 形式排列x (n )和X (k )，则式（4）、（5）可写成：X =W N x X W x *1N N=其中矩阵W N 由下式给出：[]⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=∆---→↓-≤≤２）（１ 1111,011111N N N N N NN n N n k knNN W W W W k W W 可由下面的MA TLAB 函数dft 和idft 实现离散傅立叶变换运算。