自相关函数和互相关函数计算和作图的整理

自相关函数和互相关函数计算和作图的整理

1. 首先说说自相关和互相关的概念。

--[转版友gghhjj]-------------------------------------------------------------------------------------

这个是信号分析里的概念，他们分别表示的是两个时间序列之间和同一个时间序列在任意两个不同时刻的取值之间的相关程度，即互相关函数是描述随机信号 x(t),y(t)在任意两个不同时刻t1，t2的取值之间的相关程度，自相关函数是描述随机信号x(t)在任意两个不同时刻t1，t2的取值之间的相关程度。

----------------------------------------------------------------------------------------------------------------------------------------------

--[转版友hustyoung]-----------------------------------------------------------------------------------

自相关函数是描述随机信号X(t)在任意两个不同时刻t1，t2的取值之间的相关程度；互相关函数给出了在频域内两个信号是否相关的一个判断指标，把两测点之间信号的互谱与各自的自谱联系了起来。它能用来确定输出信号有多大程度来自输入信号，对修正测量中接入噪声源而产生的误差非常有效。

----------------------------------------------------------------------------------------------------------------------------------------------

事实上，在图象处理中，自相关和互相关函数的定义如下：设原函数是f(t)，则自相关函数定义为R(u)=f(t)*f(-t)，其中*

表示卷积；设两个函数分别是f(t)和g(t)，则互相关函数定义为R(u)=f(t)*g(-t)，它反映的是两个函数在不同的相对位置上互相匹配的程度。

那么，如何在matlab中实现这两个相关并用图像显示出来呢？这个问题happy教授给出了完整答案：

-----------[转happy教授]---------------------

dt=.1;

t=[0:dt:100];

x=cos(t);

[a,b]=xcorr(x,'unbiased');

plot(b*dt,a)

-----------------------------------------------------

上面代码是求自相关函数并作图，对于互相关函数，稍微修改一下就可以了，即把[a,b]=xcorr(x,'unbiased');改为

[a,b]=xcorr(x,y,'unbiased');便可。

2. 实现过程：

在Matalb中，求解xcorr的过程事实上是利用Fourier变换中的卷积定理进行的，即R(u)=ifft(fft(f)×fft(g))，其中×表示乘法，注：此公式仅表示形式计算，并非实际计算所用的公式。当然也可以直接采用卷积进行计算，但是结果会与xcorr 的不同。事实上，两者既然有定理保证，那么结果一定是相同的，只是没有用对公式而已。下面是检验两者结果相同的代码：

dt=.1;

t=[0:dt:100];

x=3*sin(t);

y=cos(3*t);

subplot(3,1,1);

plot(t,x);

subplot(3,1,2);

plot(t,y);

[a,b]=xcorr(x,y);

subplot(3,1,3);

plot(b*dt,a);

yy=cos(3*fliplr(t)); % or use: yy=fliplr(y);

z=conv(x,yy);

pause;

subplot(3,1,3);

plot(b*dt,z,'r');

即在xcorr中不使用scaling。

3. 其他相关问题：

1) 相关程度与相关函数的取值有什么联系？

-------------[转版友

gghhjj]-------------------------------------------------------------------------------------

相关系数只是一个比率，不是等单位量度，无什么单位名称，也不是相关的百分数，一般取小数点后两位来表示。

相关系数的正负号只表示相关的方向，绝对值表示相关的程度。因为不是等单位的度量，因而不能说相关系数0.7是0.35两倍，只能说相关系数为0.7的二列变量相关程度比相关系数为0.35的二列变量相关程度更为密切和更高。也不能说相关系数从0.70到0.80与相关系数从0.30到0.40增加的程度一样大。

对于相关系数的大小所表示的意义目前在统计学界尚不一致，但通常按下是这样认为的：

相关系数相关程度

0.00-±0.30微相关

±0.30-±0.50实相关

±0.50-±0.80显著相关

±0.80-±1.00高度相关

----------------------------------------------------------------------------------------------------------------------------------------------

2) 功率，能量，自相关函数的关系：

---[转happy教授]------------------------------------------------------------------------------------------- 参见/jingpinke/xhst/final/XiTongJiaoCai/chap6/chap6_3/chap6_3_3.htm

需要指出的是，相关和相关函数的概念原本是为描述随机过程的统计特征而引入的，称之为统计相关函数。按照随机过程的理论，要获得一个实际随机过程的统计相关函数是相当困难的，但对于满足各态历经性(遍历性)或广义平稳的随机过程，它们的统计相关函数等于其一个样本函数的时间相关函数。从确定性信号引出相关的概念，是为后续课程的学习打下一个基础。

两信号互相关函数的傅里叶变换等于其中第一个信号变换与第二个信号变换取共轭二者之乘积，这就是相关定理。对于自相关函数，它的傅里叶变换等于原信号幅度谱的平方。

周期余弦信号和它的自相关函数具有相同的角频率，即周期信号的自相关函数仍然是同周期的周期信号。

在实际应用中，有些信号无法求它的傅里叶变换，但是可以用求自相关函数的方法求得信号的功率谱。

--------------------------------------------------------------------------------------------------------------------------------------

3) 与matlab中corrcoef函数的关系：

以两个不同信号（序列）为例，xcorr函数是通过不反折的卷积来衡量这两个信号在不同位置的相似程度，假设两个序列的长度分别是m和n，则得到的是一个长度为2*max(m,n)-1的序列，也就是说，当m和n不相等的时候，在执行xcorr的时候会先对短的那个序列进行0扩充，使得m与n相等；而 corrcoef函数是通过协方差矩阵来衡量这两个信号在不同局部的相似程度，计算公式是：C(1,2)/SQRT(C(1,1)*C(2,2))，其中 C表示矩阵[f,g]的协方差矩阵，假设f和g都是列向量（这两个序列的长度必须一样才能参与运算），则得到的（我们感兴趣的部分）是一个数。以默认的 A=corrcoef(f,g)为例，输出A是一个二维矩阵（对角元恒为1），我们感兴趣的f和g的相关系数就存放在A(1,2)=A(2,1)上，其值在[-1,1]之间，1表示最大的正相关（例如x=[1;2;3], y=[5;7;9]），-1表示绝对值最大的负相关（例如x=[1;2;3], y=[12;7;2]）。对于一般的矩阵X，执行A=corrcoef(X)后，A中每个值的所在行a和列b，反应的是原矩阵X中相应的第a个列向量和第 b个列向量的相似程度（即相关系数）。

4）互相关函数图像的横坐标问题

以下是我编程的例子，主要求两个信号的相位差，按照某篇参考资料的说法，t_max对应的值就应该是它们的相位差，但是这个程序中做出的互相关函数的横坐标不是-40到+40，而是0到1200，请问这个横坐标表示的是什么意思呢？

n=99;%设定每周期数据采集点数

T=6;%采样周期数

t=0:2*pi/(n-1):2*T*pi;%采样数

y1=4*sin(t);%信号1

y2=8*sin(t+pi/6);%信号2，相位差取pi/6

Cc=xcorr(y1,y2);%求互相关函数

[y_max,t_max]=max(Cc)%找出Cc的最大值及对应的t_max

subplot(311); plot(t,y1); grid;

subplot(312); plot(t,y2); grid;

subplot(313); plot(Cc); grid;

解答：

楼主得到的互相关函数，其横坐标是样点数。因为Cc长为1177，绘图中便按1~1177排列。信号y1和y2分别长589，在计算互相关函数时从-588计算到588，共有1177个互相关系数，中心点是589。

计算出的t_max =597应和中心点求差值，算出差几个样点，再进一步求出相应的相位差。我把程序稍作修改为：

n=99; %设定每周期数据采集点数

T=6; %采样周期数

t=0:2*pi/(n-1):2*T*pi; %采样数

N=length(t);

fs=1/(n-1);

y1=4*sin(t); %信号1

y2=8*sin(t+pi/6); %信号2，相位差取pi/6

Cc=xcorr(y1,y2); %求互相关函数

[y_max,t_max]=max(Cc) %找出Cc的最大值及对应的t_max

sample_delay=t_max-N; %计算与中心点相差的样点数

T_sig=(N-1)/6; %求信号一个周期的样点数,一个周期相对应于2*pi

delay1=pi/6 %pi/6的弧度

delay2=2*pi*sample_delay/T_sig %计算与中心点相差的样点数所对应的弧度值

这样计算出pi/6=0.5236，而从相关函数最大值处求出的相位差是0.5129。

对于单频信号而言，时移等于相移，相关系数从初始值变成最大值的时候，说明相位差也从初始相位差变为零