实验十一 离散信号时域分析的MATLAB 实现
一、实验目的
1. 熟悉MATLAB 编程方法、常用语句和可视化绘图技术； 2. 掌握序列时域运算的MATLAB 编程方法。
二、实验原理
在用MATLAB 表示离散信号并将其可视化时，由于矩阵元素个数是有限的，因此无法表示无限长序列；另外，离散信号无法进行符号运算。在MATLAB 中，绘制离散序列波形图的专用命令为stem( )。其格式有：
（1）stem(k,f)
在图形窗口中，绘制出样值顶部为空心圆的序列f (k)波形图。
（2）stem(k,f,’fill’)
在图形窗口中，绘制出样值顶部为实心圆的序列f (k)波形图。
下面介绍离散序列的MA TLAB 表示、基本运算（相加、相乘、平移、反转、尺度变换）、卷积和的实现及其图形显示方法。
1. 单位序列δ(k ) 单位序列的定义：
δ(k ) =⎨
1, k =0⎧
⎩0, k ≠0
下面为绘制δ(k-k 0) 波形图的子程序：
function impseq(k1,k2,k0) %单位序列δ(k-k0),k0为时移量
k=k1:k2; %k1,k2为序列的起止序列号 n=length(k); x=zeros(1,n);
x(1,k0-k1+1)=1; %在k0时刻信号赋值为1
stem(k,x,'fill') axis([k1,k2,0,1.1]) title('单位序列d(k-k0)')
输入如下命令，则可获得单位序列δ(k-3) 的波形图，如图11-1所示。
图
11-1
impseq(-1,5,3) 2. 单位阶跃序列ε(k ) 单位序列的定义：
ε(k ) =⎨
1, k ≥0⎧
⎩0, k
下面为绘制ε(k-k 0) 波形图的MA TLAB 子程序。
function stepseq(k1,k2,k0) %单位阶跃序列,k0为时移量 k=k1:k0-1; %k1,k2为序列的起止序列号 kk=length(k);
x=zeros(1,kk); %k0前信号赋值为0 stem(k,x,'fill') %绘出k1~k0-1的波形(0值) hold on n=k0:k2;
nn=length(n);
图11-2
x=ones(1,nn); %k0后信号赋值为1 stem(n,x,'fill') %绘出k0~k2的波形(1值) hold off
axis([k1,k2,0,1.1]) title('单位阶跃序列')
运行如下命令，则可获得单位序列 (k-3) 的波形图，如图11-2所示。
stepseq(-1,10,3)
3. 序列的相加（减）、相乘运算
对序列向量f 1(k)、f 2(k)相加或相乘，可以通过补零的方式使f 1(k)、f 2(k)成为具有相同维数的序列向量s 1(k)、s 2(k)，然后对s 1(k)、s 2(k)相加或相乘。因此，序列向量f 1(k)、f 2(k)的维数可以不同。以下函数可实现序列向量f 1(k)、f 2(k)的相加或相乘运算。
function [f,k]=sigadd(f1,k1,f2,k2)
%实现序列f1,f2的相加, 相减, 相乘, 可据实际需要作选择 %f1,k1;f2,k2是参加运算的序列向量及其时间向量 %f,k作为返回的和(差, 积) 序列及其时间向量 %将f1,f2转换成等长序列s1,s2
k=min(min(k1),min(k2)):max(max(k1),max(k2)); s1=zeros(1,length(k)); s2=s1; %初始化序列 s1(find((k>=min(k1))&(k=min(k2))&(k
stem(k,f,'fill')
axis([(min(min(k1),min(k2))-1),(max(max(k1),max(k2))+1),(min(f)-0.5),(max(f)+0.5)])
↓
↓
例11-1. 已知序列f 。编写M 文件求f (k ) +f (k ) 。(k ) ={-3, 2, 3, 1, 2},f (k ) ={1, 1, 1}1212
解：运行如下M 文件，可实现f (k ) +f (k ) ，结果如图11-1所示。 12k1=-1:3; f1=[-3 2 3 1 2]; k2=-1:1; kk=length(k2); f2=ones(1,kk);
subplot(2,2,1);stem(k1,f1,'fill');title('f1(k)'); subplot(2,2,2);stem(k2,f2,'fill');title('f2(k)');
subplot(2,2,3);[f,k]=sigadd(f1,k1,f2,k2);title('f1(k)+f2(k)')
图11-1
若要实现序列f 1(k)、f 2(k)的相乘或相减运算，只需将xlyunxuan(f1,k1,f2,k2)子程序中的相乘或相减语句设置为有效即可。
4. 序列的平移、反转 （1）序列的平移
序列的平移可以看作是将序列的时间序号向量平移，而对应原时间序号的序列样值不变。要将序列左移k 0个单位时，则将时间序号向量都减小k 0个单位；若要右移k 0个单位时，则将时间序号向量都增大k 0个单位。实现序列平移的子函数如下：
function [x,n]=sigshift(f,k,k0) %实现序列平移:x(k)=f(k-k0)
n=k+2;x=f;
subplot(1,2,1); stem(k,f,'fill');title('f(k)');xlabel('k'); subplot(1,2,2); stem(n,x,'fill');title('f(k-k0)');xlabel('k');
例11-2. 已知指数序列f (k ) =(0. 5) ε(k ) ，绘出f (k -2) 的波形图。
k
解：运行如下M 文件，可得如图11-2所示的结果。
k=0:5;
f=(0.5).^k; %定义序列f(k)
[x,k]=sigshift(f,k,2) %调用平移子函数
（2）序列的反转
图11-2
序列的反转可用MA TLAB 中的fliplr 函数实现。以下是实现序列反转及其结果可视化的函数。
function [x,n]=sigfold(f,k) %实现序列反转:x(k)=f(-k) x=fliplr(f); n=-fliplr(k);
subplot(1,2,1); stem(k,f,'fill');title('f(k)');xlabel('k'); subplot(1,2,2); stem(n,x,'fill');title('f(-k)');xlabel('k'); 5. 序列的卷积运算
序列f 1(k)、f 2(k)的卷积和运算f (k)= f1(k)*f 2(k)，可由MA TLAB 的conv( )函数实现，调用格式为：
f=conv(f1,f2)
k
1, -2≤k ≤2⎧2, 0≤k ≤3⎧
如：已知序列：f ，运行如下M 文件(k ) =，f (k ) =⎨⎨12
0, 其他0, 其他⎩⎩
可求其卷积和：
k1=-2:2;
f1=ones(1,length(k1)); k2=0:3; f2=2.^k2; f=conv(f1,f2) 结果为：
f =
1 3 7 15 15 14 12 8
可见，conv( )函数不需要给定f 1(k)、f 2(k)的非零样值的时间序号，也不返回卷积和序列f (k) 的时间序号；此外，conv( )假定f 1(k)、f 2(k)都是从k=0开始，这就限制了它的应用范围。因此，要对从任意k 值开始的序列进行卷积和运算，同时正确标识出函数conv( )的计算结果各量f ，还须构造序列f 1(k)、f 2(k)和f (k)的对应序号向量。下面是求序列f 1(k)、f 2(k)卷积和的实用函数dconv( )，它可实现序号向量的返回。
function [f,k]=dconv(f1,k1,f2,k2) %求卷积和:f(k)=f1(k)*f2(k) f=conv(f1,f2)
k0=k1(1)+k2(1); %计算序列f 非零样值的起点位置k0 k3=length(k1)+length(k2)-2; %计算序列f 非零样值的宽度 k=k0:k0+k3; %确定序列f 非零样值的序号向量 subplot(2,2,1); stem(k1,f1,'fill');title('f1(k)');xlabel('k'); subplot(2,2,2); stem(k2,f2,'fill');title('f2(k)');xlabel('k'); subplot(2,2,3); stem(k,f,'fill');title('f(k)=f1(k)*f2(k)');xlabel('k'); h=get(gca,'position'); h(3)=2.5*h(3); set(gca,'position',h)
k
0. 5k , -2≤k ≤2⎧2, 0≤k ≤3⎧
例11-3. 已知序列：f ，求其卷积(k ) =，f (k ) =⎨⎨12
0, 其他0, 其他⎩⎩
和。
解：运行如下M 文件
k1=-2:2;
f1=0.5.*k1; %定义序列f1(k) k2=0:3;
f2=2.^k2; %定义序列f2(k) [f,k]=dconv(f1,k1,f2,k2); %求卷积和 图形结果则如图11-3所示，文本结果如下：
f =
Columns 1 through 6
图11-3
-1.0000 -2.5000 -5.0000 -9.5000 -2.0000 4.0000 Columns 7 through 8 8.0000 8.0000 6. 离散信号相关函数的计算
在数字信号处理中，广泛用到信号的相关运算。两个序列的互相关用来度量这两个序列的相似程度。给定两个长度相同、能量受限的序列x(k)和y(k)，它们的互相关函数是另一个序列，定义为：
R (n ) =x (k ) y (k -n ) =x (n ) *y (-n ) ∑xy
k =
-∞
∞
其中，n 称为滞后参数。当x(k)=y(k)时，称为x(k)的自相关函数，它提供了序列位置不同基准情况下自相似程度的度量。由于信号的自相关函数与信号的功率谱密度是一对傅里叶变换对，这为信号的功率谱计算提供了另一重要途径。
由互相关的定义式可知，x(k)和y(k)的互相关函数可以由序列x(n)和y(-n)的卷积和求得，实现互相关函数 correlation (x,k1,y,k2)运算的程序如下：
function [Rxy,n]=correlation(x,k1,y,k2) [y0,k0]=sigfold(y,k2); %y(k)反转, 生成y(-k) [Rxy,n]=dconv(x,k1,y0,k0); hold off
subplot(3,1,1);stem(k1,x,'fill');xlabel('k');ylabel('x(k)'); subplot(3,1,2);stem(k2,y,'fill');xlabel('k');ylabel('y(k)'); subplot(3,1,3);stem(n,Rxy,'fill');xlabel('n');ylabel('Rxy(n)');
↓
例11-4. 设x y(k)为x(k)加入噪声干扰并移位后的序列： (k ) =[2, 3, 5, 1, -1, 3, 1]为原型序列，
y(k)=x(k-1)+w(k)
其中，w(k)为具有零均值和单位方差的高斯序列。计算序列y(k)和x(k)的互相关。
解：运行如下M 文件，可得如图11-4的结果。 k1=-3:3;
x=[2 3 5 1 -1 3 1]; %生成x(k) [y2,k2]=sigshift(x,k1,1); %生成x(k-1) w=randn(1,length(y2)); %模拟噪声信号w(k) kw=k2;
[y,k]=sigadd(y2,k2,w,kw); %生成
y(k)=x(k-1)+w(k)
图11-4
[Rxy,n]=correlation(x,k1,y,k) %互相关运算
三. 实验内容与步骤
πj ⎫⎛k 4⎪1. 编写M 文件，绘制复指数序列x 范围内x (k ) 2≤k ≤2(k ) =(1+j ) =2e ⎪在-⎝⎭
k
的实部、虚部、模值和相角的波形图。
2. 已知序列f 1(k ) 如图11-5（A ）所示，编写M 文件，绘出如下序列的波形： （1）f 1(-k +2)，（2）f 1(k -2) ε(k -1) 。
3. 已知序列f 1(k ) 、f 2(k ) 如图11-5所示，编写M 文件，求离散卷积和f 1(k )*f 2(k ) ，并绘出其波形。
4. 编写M 文件，求如图11-5所示序列f 1(k ) 、f 2(k ) 的互相关函数。
图11-5
四. 实验报告要求
列出M 文件和运行结果。
总结运用函数stem 、conv 、fliplr 以及相关子函数等进行序列表示、运算（相加、相乘、平移、反转、卷积和、相关）和可视化的实现方法。
连续信号又称为模拟信号，其信号存在于整个时间范围内，包括单位冲激信号，单位阶跃信号，斜坡信号，实指数信号，正弦信号，指数调制正弦信号，等等。以下我们将这几种常见的连续信号用MATLAB 实现。
1、单位冲激信号
【例1】t=1/A=1/50时，单位脉冲δ(t)的MA TLAB 实现程序如下： clear all
t1=-0.5:0.001:1; K=3;
n1=length(t1); u=zeros(1,n1); t2=-0.5:0.001:0;
实验十一 离散信号时域分析的MATLAB 实现
一、实验目的
1. 熟悉MATLAB 编程方法、常用语句和可视化绘图技术； 2. 掌握序列时域运算的MATLAB 编程方法。
二、实验原理
在用MATLAB 表示离散信号并将其可视化时，由于矩阵元素个数是有限的，因此无法表示无限长序列；另外，离散信号无法进行符号运算。在MATLAB 中，绘制离散序列波形图的专用命令为stem( )。其格式有：
（1）stem(k,f)
在图形窗口中，绘制出样值顶部为空心圆的序列f (k)波形图。
（2）stem(k,f,’fill’)
在图形窗口中，绘制出样值顶部为实心圆的序列f (k)波形图。
下面介绍离散序列的MA TLAB 表示、基本运算（相加、相乘、平移、反转、尺度变换）、卷积和的实现及其图形显示方法。
1. 单位序列δ(k ) 单位序列的定义：
δ(k ) =⎨
1, k =0⎧
⎩0, k ≠0
下面为绘制δ(k-k 0) 波形图的子程序：
function impseq(k1,k2,k0) %单位序列δ(k-k0),k0为时移量
k=k1:k2; %k1,k2为序列的起止序列号 n=length(k); x=zeros(1,n);
x(1,k0-k1+1)=1; %在k0时刻信号赋值为1
stem(k,x,'fill') axis([k1,k2,0,1.1]) title('单位序列d(k-k0)')
输入如下命令，则可获得单位序列δ(k-3) 的波形图，如图11-1所示。
图
11-1
impseq(-1,5,3) 2. 单位阶跃序列ε(k ) 单位序列的定义：
ε(k ) =⎨
1, k ≥0⎧
⎩0, k
下面为绘制ε(k-k 0) 波形图的MA TLAB 子程序。
function stepseq(k1,k2,k0) %单位阶跃序列,k0为时移量 k=k1:k0-1; %k1,k2为序列的起止序列号 kk=length(k);
x=zeros(1,kk); %k0前信号赋值为0 stem(k,x,'fill') %绘出k1~k0-1的波形(0值) hold on n=k0:k2;
nn=length(n);
图11-2
x=ones(1,nn); %k0后信号赋值为1 stem(n,x,'fill') %绘出k0~k2的波形(1值) hold off
axis([k1,k2,0,1.1]) title('单位阶跃序列')
运行如下命令，则可获得单位序列 (k-3) 的波形图，如图11-2所示。
stepseq(-1,10,3)
3. 序列的相加（减）、相乘运算
对序列向量f 1(k)、f 2(k)相加或相乘，可以通过补零的方式使f 1(k)、f 2(k)成为具有相同维数的序列向量s 1(k)、s 2(k)，然后对s 1(k)、s 2(k)相加或相乘。因此，序列向量f 1(k)、f 2(k)的维数可以不同。以下函数可实现序列向量f 1(k)、f 2(k)的相加或相乘运算。
function [f,k]=sigadd(f1,k1,f2,k2)
%实现序列f1,f2的相加, 相减, 相乘, 可据实际需要作选择 %f1,k1;f2,k2是参加运算的序列向量及其时间向量 %f,k作为返回的和(差, 积) 序列及其时间向量 %将f1,f2转换成等长序列s1,s2
k=min(min(k1),min(k2)):max(max(k1),max(k2)); s1=zeros(1,length(k)); s2=s1; %初始化序列 s1(find((k>=min(k1))&(k=min(k2))&(k
stem(k,f,'fill')
axis([(min(min(k1),min(k2))-1),(max(max(k1),max(k2))+1),(min(f)-0.5),(max(f)+0.5)])
↓
↓
例11-1. 已知序列f 。编写M 文件求f (k ) +f (k ) 。(k ) ={-3, 2, 3, 1, 2},f (k ) ={1, 1, 1}1212
解：运行如下M 文件，可实现f (k ) +f (k ) ，结果如图11-1所示。 12k1=-1:3; f1=[-3 2 3 1 2]; k2=-1:1; kk=length(k2); f2=ones(1,kk);
subplot(2,2,1);stem(k1,f1,'fill');title('f1(k)'); subplot(2,2,2);stem(k2,f2,'fill');title('f2(k)');
subplot(2,2,3);[f,k]=sigadd(f1,k1,f2,k2);title('f1(k)+f2(k)')
图11-1
若要实现序列f 1(k)、f 2(k)的相乘或相减运算，只需将xlyunxuan(f1,k1,f2,k2)子程序中的相乘或相减语句设置为有效即可。
4. 序列的平移、反转 （1）序列的平移
序列的平移可以看作是将序列的时间序号向量平移，而对应原时间序号的序列样值不变。要将序列左移k 0个单位时，则将时间序号向量都减小k 0个单位；若要右移k 0个单位时，则将时间序号向量都增大k 0个单位。实现序列平移的子函数如下：
function [x,n]=sigshift(f,k,k0) %实现序列平移:x(k)=f(k-k0)
n=k+2;x=f;
subplot(1,2,1); stem(k,f,'fill');title('f(k)');xlabel('k'); subplot(1,2,2); stem(n,x,'fill');title('f(k-k0)');xlabel('k');
例11-2. 已知指数序列f (k ) =(0. 5) ε(k ) ，绘出f (k -2) 的波形图。
k
解：运行如下M 文件，可得如图11-2所示的结果。
k=0:5;
f=(0.5).^k; %定义序列f(k)
[x,k]=sigshift(f,k,2) %调用平移子函数
（2）序列的反转
图11-2
序列的反转可用MA TLAB 中的fliplr 函数实现。以下是实现序列反转及其结果可视化的函数。
function [x,n]=sigfold(f,k) %实现序列反转:x(k)=f(-k) x=fliplr(f); n=-fliplr(k);
subplot(1,2,1); stem(k,f,'fill');title('f(k)');xlabel('k'); subplot(1,2,2); stem(n,x,'fill');title('f(-k)');xlabel('k'); 5. 序列的卷积运算
序列f 1(k)、f 2(k)的卷积和运算f (k)= f1(k)*f 2(k)，可由MA TLAB 的conv( )函数实现，调用格式为：
f=conv(f1,f2)
k
1, -2≤k ≤2⎧2, 0≤k ≤3⎧
如：已知序列：f ，运行如下M 文件(k ) =，f (k ) =⎨⎨12
0, 其他0, 其他⎩⎩
可求其卷积和：
k1=-2:2;
f1=ones(1,length(k1)); k2=0:3; f2=2.^k2; f=conv(f1,f2) 结果为：
f =
1 3 7 15 15 14 12 8
可见，conv( )函数不需要给定f 1(k)、f 2(k)的非零样值的时间序号，也不返回卷积和序列f (k) 的时间序号；此外，conv( )假定f 1(k)、f 2(k)都是从k=0开始，这就限制了它的应用范围。因此，要对从任意k 值开始的序列进行卷积和运算，同时正确标识出函数conv( )的计算结果各量f ，还须构造序列f 1(k)、f 2(k)和f (k)的对应序号向量。下面是求序列f 1(k)、f 2(k)卷积和的实用函数dconv( )，它可实现序号向量的返回。
function [f,k]=dconv(f1,k1,f2,k2) %求卷积和:f(k)=f1(k)*f2(k) f=conv(f1,f2)
k0=k1(1)+k2(1); %计算序列f 非零样值的起点位置k0 k3=length(k1)+length(k2)-2; %计算序列f 非零样值的宽度 k=k0:k0+k3; %确定序列f 非零样值的序号向量 subplot(2,2,1); stem(k1,f1,'fill');title('f1(k)');xlabel('k'); subplot(2,2,2); stem(k2,f2,'fill');title('f2(k)');xlabel('k'); subplot(2,2,3); stem(k,f,'fill');title('f(k)=f1(k)*f2(k)');xlabel('k'); h=get(gca,'position'); h(3)=2.5*h(3); set(gca,'position',h)
k
0. 5k , -2≤k ≤2⎧2, 0≤k ≤3⎧
例11-3. 已知序列：f ，求其卷积(k ) =，f (k ) =⎨⎨12
0, 其他0, 其他⎩⎩
和。
解：运行如下M 文件
k1=-2:2;
f1=0.5.*k1; %定义序列f1(k) k2=0:3;
f2=2.^k2; %定义序列f2(k) [f,k]=dconv(f1,k1,f2,k2); %求卷积和 图形结果则如图11-3所示，文本结果如下：
f =
Columns 1 through 6
图11-3
-1.0000 -2.5000 -5.0000 -9.5000 -2.0000 4.0000 Columns 7 through 8 8.0000 8.0000 6. 离散信号相关函数的计算
在数字信号处理中，广泛用到信号的相关运算。两个序列的互相关用来度量这两个序列的相似程度。给定两个长度相同、能量受限的序列x(k)和y(k)，它们的互相关函数是另一个序列，定义为：
R (n ) =x (k ) y (k -n ) =x (n ) *y (-n ) ∑xy
k =
-∞
∞
其中，n 称为滞后参数。当x(k)=y(k)时，称为x(k)的自相关函数，它提供了序列位置不同基准情况下自相似程度的度量。由于信号的自相关函数与信号的功率谱密度是一对傅里叶变换对，这为信号的功率谱计算提供了另一重要途径。
由互相关的定义式可知，x(k)和y(k)的互相关函数可以由序列x(n)和y(-n)的卷积和求得，实现互相关函数 correlation (x,k1,y,k2)运算的程序如下：
function [Rxy,n]=correlation(x,k1,y,k2) [y0,k0]=sigfold(y,k2); %y(k)反转, 生成y(-k) [Rxy,n]=dconv(x,k1,y0,k0); hold off
subplot(3,1,1);stem(k1,x,'fill');xlabel('k');ylabel('x(k)'); subplot(3,1,2);stem(k2,y,'fill');xlabel('k');ylabel('y(k)'); subplot(3,1,3);stem(n,Rxy,'fill');xlabel('n');ylabel('Rxy(n)');
↓
例11-4. 设x y(k)为x(k)加入噪声干扰并移位后的序列： (k ) =[2, 3, 5, 1, -1, 3, 1]为原型序列，
y(k)=x(k-1)+w(k)
其中，w(k)为具有零均值和单位方差的高斯序列。计算序列y(k)和x(k)的互相关。
解：运行如下M 文件，可得如图11-4的结果。 k1=-3:3;
x=[2 3 5 1 -1 3 1]; %生成x(k) [y2,k2]=sigshift(x,k1,1); %生成x(k-1) w=randn(1,length(y2)); %模拟噪声信号w(k) kw=k2;
[y,k]=sigadd(y2,k2,w,kw); %生成
y(k)=x(k-1)+w(k)
图11-4
[Rxy,n]=correlation(x,k1,y,k) %互相关运算
三. 实验内容与步骤
πj ⎫⎛k 4⎪1. 编写M 文件，绘制复指数序列x 范围内x (k ) 2≤k ≤2(k ) =(1+j ) =2e ⎪在-⎝⎭
k
的实部、虚部、模值和相角的波形图。
2. 已知序列f 1(k ) 如图11-5（A ）所示，编写M 文件，绘出如下序列的波形： （1）f 1(-k +2)，（2）f 1(k -2) ε(k -1) 。
3. 已知序列f 1(k ) 、f 2(k ) 如图11-5所示，编写M 文件，求离散卷积和f 1(k )*f 2(k ) ，并绘出其波形。
4. 编写M 文件，求如图11-5所示序列f 1(k ) 、f 2(k ) 的互相关函数。
图11-5
四. 实验报告要求
列出M 文件和运行结果。
总结运用函数stem 、conv 、fliplr 以及相关子函数等进行序列表示、运算（相加、相乘、平移、反转、卷积和、相关）和可视化的实现方法。
连续信号又称为模拟信号，其信号存在于整个时间范围内，包括单位冲激信号，单位阶跃信号，斜坡信号，实指数信号，正弦信号，指数调制正弦信号，等等。以下我们将这几种常见的连续信号用MATLAB 实现。
1、单位冲激信号
【例1】t=1/A=1/50时，单位脉冲δ(t)的MA TLAB 实现程序如下： clear all
t1=-0.5:0.001:1; K=3;
n1=length(t1); u=zeros(1,n1); t2=-0.5:0.001:0;