序列的卷积运算与相关运算——MATLAB

一、实验目的

1、掌握有限长序列线性卷积的编程计算原理，并能够利用Matlab或C语言编写算法程序进行线性卷积运算的程序实现;
2、学会线性卷积函数和线性相关函数的使用方法，并能利用二者进行有限长序列的线性卷积与线性相关运算的实现;自相关序列的基本性质;
4、验证线性卷积运算和线性相关运算之间的关系,掌握由线性卷积运算实现线性相关运算的方法;
5、学会利用自相关序列进行信号周期性的检测与判别。

二、实验原理

在数字信号处理中，相关(correlation)可以分为互相关(cross correlation)和自相关(auto-correlation). 互相关是两个数字序列之间的运算；自相关是单个数字序列本身的运算，可以看成是两个相同数字序列的互相关运算．互相关用来度量一个数字序列移位后，与另一个数字序列的相似程度．其数学公式如下：

其中，f和g为数字序列，n为移位的位数，f*表示f序列值的复数共轭，即复数的实部不变，虚部取反而卷积（convolution）与互相关运算相似，定义为将其中一个序列反转并移位后，两个序列的乘积的积分（求和），其数学公式如下：

其中，f和g 为数字序列，n 为移位的位数.
在实数范围内，f的复数共轭f*=f．此时，通过比较上面两式可知：序列f与将序列g反转后的序列的卷积为序列f与序列g的互相关
如下所示:
第一步:确定(n) 有限区间为[N1+N2,M1+M2], 这一步的的工作是为了决定y(n) 中的哪些序列值要计算,因为区间外的序列值都是0，无须计算
第二步:把x(n)和h(n)的有限区间都变为从0开始。
第三步:利用公式计算序列值。计算y(n)的n=0,1,2,…M1+M2 -N1- N2，所对应的y(0),y(1),…,y(M1+ M2 - N1+N2)。求每一个序列值时,乘加运算的结束标志是h(n-k)的n-k
第四步:把y(n)的序号由0开始变为由N1+N2开始,其它序号依次变为N1＋M2为止,就是N2＋1，N1＋N2+2,直到M1+真正的输出序列y(n)。
2、线性相关运算及计算方法
设序列x(n)利y(n)为两个能量有限的离散时间信号，则序列之间的线性互相关序列（运算)
Rxy(m)反映了序列x(n)与右移m个单位的序列y(n)之间的相似程度，某一m处Rxy(m)的绝对值越大,则此时二者之间的相似性就越强。
序列的线性相关运算和线性卷积运算之关系:
Rxy(m)= y(m)x(-m)
表明:序列的相关运算可以序列的线性卷积运算为基础进行实现。
(1）若序列x(n)是实序列，则R.(m)为偶函数，即:
Rxy(m)=Rxy(-m)
(2）Rxx(m)在=0时取得最大值，即:
Rxx(0)≥Rxx(m)
(3）若x(n)是能量信号，则有:
lim(x➡∞) Rxx(m)=0
3、利用自相关序列讲行信号周期性的检测
原序列具有相同的周期。这一性质是实际应用中, 利用自相关运算进行信号周期性检测的基础和理论依据。

三、实验步骤、数据记录及处理

本实验主要进行离散时间信号的线性卷积运算和线性相关运算的程序实现，实验内容和具体的实验步骤如下:
其中:
1、用Matlab生成两个有限长序列x(n)、y(n)

2、编写有限长序列线性卷积算法程序计算x(n)与y(n)的线性卷积，结果令为cl(n)，即c1(n)=x(n)*y(n);利用 Matlab 的 conv()函数计算x(n)与y(n)的线性卷积，结果令为 c2(n)，即c2(n)=x(n)*y(n)。将x(n)、y(n)、cl(n)、c2(n)的序列图形绘制在一幅图上，对结果进行对比分析，检验所实现的线性卷积算法仿真程序是否正确。
3、以 conv()函数为基础，编程计算x(n)与vn)的卷积线性互相关序列，结果令为R1(n):利用 matlab的xcorr()函数计算x(n)与y(n)的线性互相关序列结果令为R2(n)。将x(n)、y(n)、R1(n)、R2(n)绘制在一幅图上,对结果进行对比分析，验线性卷积运算和线性相关运算之间的关系。
4、以下为1970年至1869年这100年间每年所记录到的太阳黑子出现的月平均次数:
101，82，66，35，31，7，20，92，154，125，85，68，38，23,
10，24，83，132，131，118，90，67，60，47，41,21,16，6, 8，
42，28，10，8，2，0，1,5，12，4，7,14，34、45，43，416 , 7,
4，2，8，17，36,50,62,14,35，46，41，30，24,37，24, 67，71，
48，28，8，13，57，122,138,103，86，63，96，66，64, 54，39，
21,7，4,23,11，15，40，62，98，124,47，30，16，7，37,74，55，
94, 96，77，59，44,
然后将该组数据去除掉均值计算该组数据的均值，然后将该组数据去除掉均值后做自相关分析,绘制出自相关序列的图形,并以计算结果为依据确定太阳黑子活动的周期。
实验代码：
(1)

clc;clear all;
x=[1,1,1];N1=1;M1=3;nx=N1:M1;
h=[2,2,2];N2=0;M2=2;nh=N2:M2;
%生成做卷积的两个有限长序列
ny=N1+N2:M1+M2;%确定卷积结果序列非零值区间
nx1=0:M1-N1;nh1=0:M2-N2;%把做卷积得两个序列取非零值得区间都变为从0开始
for n=0:M1+M2-N1-N2
    y(n+1)=0;
    for k=0:M1-N1
        m=n-k;
        if(m>=0&m<=M2-N2)
            y(n+1)=y(n+1)+x(k+1)*h(m+1);
        end
    end
end
%二层循环进行乘加运算，完成卷积运算
c1= y;
c2=conv(x,h);%将x与h进行卷积运算
subplot(4,1,1);stem(nx,x,'x');xlabel('n');ylabel('x(n)');grid on;
subplot(4,1,2);stem(nh,h,'o');xlabel('n');ylabel('y(n)');grid on;
subplot(4,1,3);stem(ny,c1,'*');xlabel('n');ylabel('c1(n)');grid on;
subplot(4,1,4);stem(ny,c2,'.');xlabel('n');ylabel('c2(n)');grid on;

clc;clear all;
x=[1,1,1];N1=1;M1=3;nx=N1:M1;
h=[2,2,2];N2=0;M2=2;nh=N2:M2;
 
z=fliplr(h);R1=conv(x,z);len_R1=length(R1);
R2=xcorr(x,h);len_R2=length(R2);
subplot(2,1,1);stem(1:len_R1,R1,'x');axis([0,len_R2+1,0,max(R1)+1]);xlabel('n');ylabel('R1(n)');grid on;
subplot(2,1,2);stem(1:len_R2,R2,'o');axis([0,len_R2+1,0,max(R1)+1]);xlabel('n');ylabel('R2(n)');grid on;

clc;clear all;
BSY=[101,82,66,35,31,7,20,92,154,125,85,68,38,23,10,24,83,132,131,118,90,67,60,47,41,21,16,6,4,7,14,34,45,43,48,42,28,10,8,2,0,1,5,12,14,35,46,41,30,24,16,7,4,2,8,17,36,50,62,67,71,48,28,8,13,57,122,138,103,86,63,37,24,11,15,40,62,98,124,96,66,64,54,39,21,7,4,23,55,94,96,77,59,44,47,30,16,7,37,74]
t1=BSY;%令t1序列为BSY
t=round(mean(BSY));%对BSY取均值并取整
t1(t1==t)=[];%去除t1中得均值t得值
[acorr,span]=xcorr(t1);%求序列得自相关
figure;plot(span,acorr);xlabel('年');ylabel('平均次数');

四、分析

设x(n)x(n)x(n)和y(n)y(n)y(n)是两个功率信号,其互相关函数定义为：

当x(n)=y(n)x(n)=y(n)x(n)=y(n),功率信号的自相关函数定义为：

如果x(n)x(n)x(n)和y(n)y(n)y(n)都是周期为 N NN 的信号，则(10)、(11)中有限区间上的平均值就等于一个周期上的平均值：

当x ( n ) = y (n ) x( n)=y( n)x(n)=y(n)时，功率信号的自相关函数定义为

由周期信号的定义可得：

互相关函数性质
设两个信号x(n)x(n)x(n)和y(n)y(n)y(n) 均为能量信号，其能量分别为

对于线性非时变离散时间系统来说，若序列x(n)是系统的输入，h(n)是系统在单位脉冲作用下的单位脉冲响应，则由于输入序列x(n)可表示为一系列脉冲的线性组合，所以，根据线性系统的叠加性质，系统的输出在系统初始不储能的条件下(零状态响应)可由图求得。

上式在运算过程存在序列的翻转、移位、相乘和相加，所以称为卷积和。x(n)*h(n)表示两个序列相卷积的运算符号，故式①也就是卷积的定义式。为了与离散傅里叶变换的循环卷积以及周期序列的周期卷积相区别，通常所指的卷积又称为线性卷积。卷积运算符合交换率，可写成另一种等效形式。
线性卷积的计算可以用解析法，也可以用图解法。若两 个序列的长度分别为N1和N2，则卷积结果的总长度应为L=N1+N2-1。
同理，对线性非时变连续系统来说，若连续时间信号x(t)是系统的输入，h(t)是系统在单位脉冲作用下的单位冲激响应，则系统在零状态的输出为它们的卷积积分。
线性卷积是数字信号处理中最常见的一种基本运算，不仅用于系统分析还用于系统设计。如果代表滤波器的脉冲响应则卷积运算就是一种线性滤波，y(n)是信号x(n)通过滤波器后的响应。
卷积性质
（1）结合律：三个序列卷和运算，任意两个序列先卷和运算，再与第3个序列作卷和运算，其运算结果等同。即

（2）交换律：离散序列卷和运算满足交换律，即两序列卷和运算与卷和次序无关，即　
（3）分配律：两个序列先行相加运算再与第3个序列做卷和运算，其结果等于这两个序列分别与第3个序列先做卷和运算，然后二者再相加。

离散线性卷积的定义：设长度为N1的序列x(n)和长度为N2的序列h(n)进行线性卷积，得到长度为N1+N2-1的y(n)。离散圆周卷积的定义：圆周卷积是定义在有限长序列之间的。设有限长序列x(n)和h(n)的长度分别为N1和N2，取N>=max(N1,N2)，定义它们的N点圆周卷积6。圆周卷积与线性卷积之间的关系：当有限长序列x(n)和h(n)的长度分别为N1和N2，取N>=max(N1,N2)，当N>=N1+N2-1，则线性卷积与圆周卷积相同。对于线性卷积，一般直接比较麻烦，由上可知当取点数足够多时（点数不够补零），可求解圆周卷积即可，而圆周卷积又可通过FFT实现，从而实现线性卷积通过FFT和IFFT实现。
自相关检测基本原理:
设混有随机噪声的信号为




将f(t)同时输入到相关接收机的2个通道,其中一路经过延迟电路,使它延迟时间r,然后将经延迟后的f(t-T)和未经延迟的f(t)送入相敏检波器(乘法器),再将乘积积分,取平均值,从而得到f(t)的自相关函数R(T)。
由于信号与噪声是不相关的,根据互相关函数的性质知，Rns(T)=Rsn(T)=0,则有R(T）=Rss(t)+Rnn(t)。假设噪声的均值为零,根据自相关函数的性质知,随着t的增大,噪声自相关函数Rnn(T)衰减至零,结果使包含着s(t)信息的e。Rss(T)相对突出从而达到检测信号的目的。

五、思考题

（1）对于无限长序列序列的线性卷积的实现，能否利用直接编程计算的方式，为什么？

（2）试从区间端点及长度两方面总结有限长序列线性卷积的结果序列的非零值区间与做卷积的两个有限长序列的非零值区的关系。

（3）结合实验或仿真结果说明在对周期信号进行自相关序列的计算时为什么总要先去掉直流成分?

六、总结

通过这次实验，让我知道序列的卷积和可以用matlab编程实现，并且实现的方法比较简单，而且matlab本身有丰富的函数库，画图、求卷积和时只要调用相应的函数即可。这次实验的难点主要就在写conv()函数和各序列长度的计算，如果不注意其中的计算关系则很容易出错，而conv()只要写错，后面的题目的图形就不能准确的绘制出来。对于后面的两个小题只要把相应的序列写出来，写入用matlab表示的矩阵中，在调用dconv()画出图形即可。由于matlab这门课的学习距现在隔的的时间较长，所以，在实验之前应该做好预习准备，熟悉matlab的基本操作，熟悉matlab相应函数的用法，这样，在实验中才能比较快速的完成实验的要求，就不会感到无从下手。