通信原理仿真（四）——信号与系统基础（二）

信号运算

信号相加、相乘

由于两个信号的相加和相乘是指对应时间的相加和相乘。数学表示为：
$y(n)=x_1(n)+x_2(n)$
$y(n)=x_1(n)\ast x_2(n)$
对序列向量$f_1(k),f_2(k)$相加或相乘，可以通过补零的方式使得$f_1(k),f_2(k)$成为具有相同维数的序列向量$s_1(k),s_2(k)$，然后对$s_1(k),s_2(k)$逐点相加或相乘。因此，序列向量可以维数不同。
以下的函数可以实现序列的相加或相乘运算。
相加运算函数：

function [f,k] = sigadd(f1,k1,f2,k2)
%sigadd 实现序列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<=max(k1))==1))=f1;
s2(find((k>=min(k2)) & (k<=max(k2))==1))=f2;
f=s1+s2;% 序列相加
% f=s1-s2;% 序列相减
% f=s1.*s2;% 序列相乘
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)]);
end

相减运算函数：

function [f,k] = sigsub(f1,k1,f2,k2)
%sigsub 实现序列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<=max(k1))==1))=f1;
s2(find((k>=min(k2)) & (k<=max(k2))==1))=f2;
% f=s1+s2;% 序列相加
f=s1-s2;% 序列相减
% f=s1.*s2;% 序列相乘
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)]);
end

相乘运算函数：

function [f,k] = sigmult(f1,k1,f2,k2)
%sigmult 实现序列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<=max(k1))==1))=f1;
s2(find((k>=min(k2)) & (k<=max(k2))==1))=f2;
% f=s1+s2;% 序列相加
% f=s1-s2;% 序列相减
f=s1.*s2;% 序列相乘
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)]);
end

例题1、计算两个序列的和与乘积

%% 信号运算
x=[0.5,1.5,1,-0.5];
n0=-1:2;
x1=2*ones(1,6);
n1=0:5;
[y1,n2]=sigmult(x1,n1,x,n0);% 序列相乘
[y2,n3]=sigadd(x1,n1,x,n0);% 序列相加
subplot(411);stem(n0,x);
ylabel('x');axis([-1,6,-2,3]);
subplot(412);stem(n1,x1);
ylabel('x1');axis([-1,6,0,3]);
subplot(413);stem(n2,y1);
ylabel('y1=x*x1');axis([-1,6,-2,3]);
subplot(414);stem(n3,y2);
ylabel('y2=x+x1');xlabel('n');
axis([-1,6,0,5]);

2、信号的移位及周期延拓

（1）信号时移

设一个离散时间信号为$x(n)$，它的时移$y(n)$表示为：
$y(n)=x(n-m)$
序列的平移可以看做是将序列的时间序号向量平移，而对应原时间序列号的序列样值不变。

function [x,n] = sigshift(f,k,k0)
%sigshift 实现序列平移：x(k)=f(k-k0)
%   f,k,k0是参加运算的序列向量及其时间向量
%   x,n作为返回的移位序列及其时间向量
n=k+k0;
x=f;
end

例题2：

%% 信号的移位和周期延拓
k=0:5;
f=(0.5).^k;% 定义序列f(k)
[f1,k1]=sigshift(f,k,2);% 调用平移子函数
subplot(211)
stem(k,f,'fill');
title('f(k)');
xlabel('k');
axis([0,8,0,1]);
subplot(212)
stem(k1,f1,'fill');
title('f(k-2)');
xlabel('k');
axis([0,8,0,1]);

（2）周期移位

序列周期延拓的数学表示为：$y(n)=x((n){)}_M$，其中M表示延拓周期。
MATLAB实现：

ny=nxs:nxf;y=x(mod(ny,M)+1);

其中，mod为求余函数，求余后加一是因为下标只能从1开始。
序列的循环移位实现：

ny=nxs:nxf;y=x(mod(ny-k,M)+1);

例题3：

%% 周期延拓
N=24;M=8;m=3;
n=0:N-1;
x1=0.8.^n;% 生成指数序列
x2=[(n>=0) & (n<=M)];% 生成矩形序列
x=x1.*x2;
xm=zeros(1,N);
for k=m+1:m+M
    xm(k)=x(k-m);% 产生移位序列
end
xc=x(mod(n,M)+1);% 产生周期延拓序列
xcm=x(mod(n-m,M)+1);%产生移位序列x(n-3)周期延拓x(n-3)
subplot(411)
stem(n,x);ylabel('x(n)');
subplot(412);stem(n,xm);ylabel('x(n-3)');
subplot(413);stem(n,xc);ylabel('x(n)_8');
subplot(414);stem(n,xcm);ylabel('x(n-3)_8');

3、信号的翻褶

对一个信号$x[n]$的翻褶运算数学上表示为：
$y[n]=x[-n]$
方法一：
修改绘图函数plot(t,x)和stem(n,x)的时间变量t和n。即用-t和-n替代原来的t和n，这样绘制出来的图形看起来就是原信号经过时域翻褶后的版本。
方法二：
直接利用原信号与其翻褶信号的数学关系式来实现。
方法三：
使用MATLAB内部函数fliplr()实现信号的翻褶运算。其用法为y=fliplr(x)。其中x为原信号x(t)或x[n]。而y则为x的时域翻褶。需要说明的是，函数fliplr()对信号做时域翻褶，仅仅将信号中各个元素的次序作了一个反转，这种反转处理是独立于实践变量t和n的。因此，如果信号与其时间变量能够用一个数学函数来表达的话，那么建议将时间变量t和n的范围指定在一个正负对称的时间区间即可。

%% 信号的翻褶
t=0:0.02:1;
x=3*t;
f=fliplr(x);
t1=-1:0.02:0;
grid on
plot(t,x,t1,f,'r');
title('信号的翻褶');

4、信号的卷积

（1）离散卷积和

两个序列卷积运算的实现：$s=conv(f_1,f_2);$
$s=conv(f_1,f_2)={\sum }_{i=-\inf }^{i=inf}{f}_1(i)f_2(k-i)$
其中，f_1(k),f_2(k)为离散序列，k=-2,-1,0,1,2,…。但是，conv函数只给出纵轴的序列值大小，而不能给出卷积的横轴序号。
对任意输入：设$f_1(k)$的非零区间为$n_1{n}_2$，长度$L_1=n_2-n_1+1$：$f_2(k)$的非零区间为$m_1{m}_2$，长度$L_2=m_2-m_1+1$，则$s(k)=f_1(k)\ast f_2(k)$的非零区间从$n_1+m_1$开始，长度为$L=L_1+L_2-1$，所以$S(k)$的非零区间为$n_1+m_1{n}_1+m_1+L-1$。
例题5：

%% 离散卷积和
Nx=20;Nh=10;%设定序列x与序列h的长度
n=0:Nx-1;
x=0.9.^n;
nh=0:Nh-1;
h=(nh==3);
ny=-Nh+1:Nx-1;
y=conv(x,h);% 计算序列x与序列h的相关结果
subplot(311);
stem(n,x);
xlabel('n');ylabel('x(n)');
subplot(3,1,2);stem(nh,h);
xlabel('n');ylabel('h(n)');
subplot(3,1,3);
stem(ny,y);
xlabel('n');ylabel('y(n)');

（2）连续卷积和离散卷积关系

因为计算机本身不能直接处理连续信号，只能由离散信号来进行近似。

可以用离散卷积和conv()求连续卷积，只需delta足够小以及在卷积和的基础上乘以delta。

（3）连续卷积坐标的确定

设$f_1(t)$的非零值坐标范围为$t_1{t}_2$，间隔P：非零值坐标范围为$t{t}_{1}t{t}_2$，$s(t)=f_1(t)\ast f_2(t)$的非零值坐标为$t_1+t{t}_1{t}_2+t{t}_2+1$。
MATLAB实现过程：
（1）将连续信号$f_1(t),f_2(t)$以时间间隔delta进行取样，得离散序列$f_1(k\delta ),f_2(k\delta )$：
（2）构造$f_1(k\delta ),f_2(k\delta )$与相对应的时间向量$k_1$和$k_2$。
（3）调用conv（）函数计算卷积积分$f(t)$的近似向量$f(k\delta )$；
（4）构造$f(k\delta )$对应的时间向量k。
计算出卷积积分的函数：

function [f,k] = sconv(f1,f2,k1,k2,p)
%sconv 计算连续信号卷积积分f(t)=f1(t)*f2(t)
%   f:卷积积分f(t)对应的非零值向量
%   k:f(t)的对应时间向量
%   f1,f2:f1(t),f2(t)的非零样值向量
%   k1,k2:f1(t),f2(t)的对应时间向量
%   p:取样时间间隔
f=conv(f1,f2);%计算序列f1,f2的卷积和f
f=f*p;
k0=k1(1)+k2(1);%计算序列非零样值的起始点位置
k3=length(f1)+length(f2)-2;%计算卷积和f的非零样值的宽度
k=k0:p:k3*p;%计算卷积和f的非零样值的时间向量
subplot(221);
plot(k1,f1);title('f1(t)');xlabel('t');ylabel('f1(t)');
subplot(222)
plot(k2,f2);title('f2(t)');xlabel('t');ylabel('f2(t)');
subplot(223)
plot(k,f);%画出卷积f(t)的波形
title('f(t)=f1(t)*f2(t)');xlabel('t');ylabel('f(t)');
h=get(gca,'position');
h(3)=2.5*h(3);
set(gca,'position',h);%将第三个子图的横坐标扩展为原来的2.5倍
end

例题6：

%% 连续卷积和离散卷积的关系
p=0.01;k1=0:p:2;
f1=exp(-k1);
k2=0:p:3;
f2=ones(1,length(k2));
[f,k]=sconv(f1,f2,k1,k2,p);

5、信号的相关运算

两序列相关运算的数字描述为：
$y(m)={\sum }_n{x}_1(n)x_2(n-m)$
两序列的相关运算为：

y=xcorr(x1,x2);

说明，若两序列具有相同的长度M，则相关序列$y(n)$的长度为2M-1。若两序列的长度不同，则短者将自动填充0。
例题7：

%% 信号的相关运算
Nx=20;Nh=10;
n=0:nx-1;
x=0.9.^n;%产生序列x
nh=0:Nh-1;
h=ones(1,Nh);%产生序列h
y=xcorr(x,h);%计算序列x与序列h的相关结果
subplot(311);stem(n,x);
xlabel('n');ylabel('x(n)');
subplot(312);stem(nh,h);
xlabel('n');ylabel('h(n)');
subplot(313);stem(ny,y);
xlabel('n');ylabel('y(n)');

6、信号的累加

序列累加的数学描述：
$y(n)={\sum }_{i=n_s}^{n}x(i)$
MATLAB中可以用cumsum函数来实现：

y=cumsum(x);

7、信号的能量

序列能量的数学描述为：
$E={\sum }_{n=0}^{N}x[n]x^{\ast }[n]$

E=sum(abs(x).^2);

8、信号的功率

序列功率的数学表述为：
$P=1/N{\sum }_{n=0}^{N-1}|x(n){|}^2$

E=sum(abs(x).^2)/N;

例题8：

%% 信号的平均功率
dt=0.001;
t=0:dt:2*pi;
x=cos(t);
P=sum(abs(x).^2*dt)./(2*pi);

参考文献

赵鸿图,茅艳. 通信原理MATLAB仿真教程. 人民邮电出版社，2010.