相关性.docx
- 文档编号:13324473
- 上传时间:2023-06-13
- 格式:DOCX
- 页数:23
- 大小:712.06KB
相关性.docx
《相关性.docx》由会员分享,可在线阅读,更多相关《相关性.docx(23页珍藏版)》请在冰点文库上搜索。
相关性
1、Xcorr相关函数
自相关函数是描述随机信号X(t)在任意两个不同时刻t1,t2的取值之间的相关程度.设原函数是f(t),则自相关函数定义为R(u)=f(t)*f(-t),其中*表示卷积.
互相关函数给出了在频域内两个信号是否相关的一个判断指标,把两测点之间信号的互谱与各自的自谱联系了起来。
它能用来确定输出信号有多大程度来自输入信号,对修正测量中接入噪声源而产生的误差非常有效.设两个函数分别是f(t)和g(t),则互相关函数定义为R(u)=f(t)*g(-t),它反映的是两个函数在不同的相对位置上互相匹配的程度。
用过Matlab的人都知道,Matlab的命令总是能一石三鸟,通过改变输入参数的注释项即可实现不同功能,如今xcorr命令的难点就在于其有四个注释项,这些注释项使得计算的结果各有不同,本文将详细介绍对应每个注释项Matlab是如何计算的,当然本文考虑输入的是一个简单一维序列x=[1,2,3],序列中数据对应的序号依次为1,2,3(请读者在阅读下文时,不要把序号和数据值弄混,这里只是个特例),其他情况读者可以轻松扩展得到。
首先确定一下,该序列的均值为2,方差值为0.6667.这个应该不用说怎么算的吧。
然后读者需要了解的是该命令意在计算序列中间隔不同距离的数据之间的关系。
第一:
缺省注释项,[a,b]=xcorr(x),通过该命令计算的结果为:
a=381483;b=-2-1012.
下面介绍一下,该过程计算机是如何计算的,首先讲b的计算,设一维序列的长度为N,则序列中任意两个数据序号相减,最小值为1-N,最大值为N-1,且能取遍两者之间的所有整数,将这些数从小到大排列得到的就是b;然后讲a的计算,在缺省注释项的情况下,a的计算是这样的,a的每一项是对应b的每一项的
1、当b
(1)=-2时,计算a
(1)时只用到一组数据——(3,1)注意顺序,只有这两个数据的序号相减(后面数据的序号减去前面的)满足b=-2,因此a
(1)的计算公式为:
3*1=3
2、当b
(2)=-1时,计算a
(2)时用到两组数据——(2,1)和(3,2),这两组数据的序号相减(后面数据的序号减去前面的)满足b=-1,因此a
(2)的计算公式为:
2*1+3*2=8
3、当b(3)=0时,计算a(3)时用到三组数据——(1,1)、(2,2)、(3,3),这三组数据的序号相减(后面数据的序号减去前面的)满足b=0,因此a(3)的计算公式为:
1*1+2*2+3*3=14
4、当b(4)=1时,计算a(4)时用到两组数据——(1,2)和(2,3),(读者请对比和情况2的区别),这两组数据的序号相减(后面数据的序号减去前面的)满足b=1,因此a(4)的计算公式为:
1*2+2*3=8
5、当b(5)=2时,计算a(4)时用到一组数据——(1,3),(读者请对比和情况1的区别),这两组数据的序号相减(后面数据的序号减去前面的)满足b=2,因此a(4)的计算公式为:
1*3=3
第二:
注释项为‘unbiased’,[a,b]=xcorr(x,‘unbiased’),通过该命令计算的结果为:
a=344.666743;b=-2-1012.
下面介绍计算机如何计算该过程,b的计算在四种注释项的情况下是相同的,就不再讲述了。
a的计算仍是和b的每一项相对应的。
1、当b
(1)=-2时,计算a
(1)时只用到一组数据(记N=1)——(3,1)注意顺序,只有这两个数据的序号相减(后面数据的序号减去前面的)满足b=-2,因此a
(1)的计算公式为:
(3*1)/1=3
2、当b
(2)=-1时,计算a
(2)时用到两组数据(记N=2)——(2,1)和(3,2),这两组数据的序号相减(后面数据的序号减去前面的)满足b=-1,因此a
(2)的计算公式为:
(2*1+3*2)/2=4
3、当b(3)=0时,计算a(3)时用到三组数据(记N=3)——(1,1)、(2,2)、(3,3),这三组数据的序号相减(后面数据的序号减去前面的)满足b=0,因此a(3)的计算公式为:
(1*1+2*2+3*3)/3=4.6667
4、当b(4)=1时,计算a(4)时用到两组数据(记N=2)——(1,2)和(2,3),(读者请对比和情况2的区别),这两组数据的序号相减(后面数据的序号减去前面的)满足b=1,因此a(4)的计算公式为:
(1*2+2*3)/2=4
5、当b(5)=2时,计算a(4)时用到一组数据(记N=1)——(1,3),(读者请对比和情况1的区别),这两组数据的序号相减(后面数据的序号减去前面的)满足b=2,因此a(4)的计算公式为:
(1*3)/1=3
第三:
注释项为‘biased’,[a,b]=xcorr(x,‘biased’),通过该命令计算的结果为:
a=1.00002.66674.66672.66671.0000,b=-2-1012。
下面介绍计算机如何计算该过程,注意到本次计算用到的序列x的长度为3,记为M=3。
1、当b
(1)=-2时,计算a
(1)时只用到一组数据——(3,1)注意顺序,只有这两个数据的序号相减(后面数据的序号减去前面的)满足b=-2,因此a
(1)的计算公式为:
(3*1)/M=1
2、当b
(2)=-1时,计算a
(2)时用到两组数据——(2,1)和(3,2),这两组数据的序号相减(后面数据的序号减去前面的)满足b=-1,因此a
(2)的计算公式为:
(2*1+3*2)/M=2.6667
3、当b(3)=0时,计算a(3)时用到三组数据——(1,1)、(2,2)、(3,3),这三组数据的序号相减(后面数据的序号减去前面的)满足b=0,因此a(3)的计算公式为:
(1*1+2*2+3*3)/M=4.6667
4、当b(4)=1时,计算a(4)时用到两组数据——(1,2)和(2,3),(读者请对比和情况2的区别),这两组数据的序号相减(后面数据的序号减去前面的)满足b=1,因此a(4)的计算公式为:
(1*2+2*3)/M=2.6667
5、当b(5)=2时,计算a(4)时用到一组数据——(1,3),(读者请对比和情况1的区别),这两组数据的序号相减(后面数据的序号减去前面的)满足b=2,因此a(4)的计算公式为:
(1*3)/M=1
第四:
注释项为‘coeff’,[a,b]=xcorr(x,‘coeff’),通过该命令计算的结果为:
a=0.21430.57141.00000.57140.2143,b=-2-1012,
下面介绍计算机如何计算该过程,这种情况实际是将第三种情况下得到的结果进行归一化,使得b=0时对应的值为1,a
(1)=1/4.6667=0.2143;a
(2)=2.6667/4.6667=0.5714,a(3)=4.6667/4.6667=1,a(4)=2.6667/4.6667=0.5714,a(5)=1/4.6667=0.2143
另:
xcorr命令在工程上的应用通常是对时间上的采样数据序列x进行处理,当数据点采完之后交给Matlab处理时,Matlab是不知道你的采样时间间隔的,它仅仅根据上文所述的计算过程对输入的数据序列x进行计算,但我们可以自己定义时间间隔,例如dt=0.01,此时t=dt*b即代表相关性计算中的时间延迟,前半部分是超前,后半部分是滞后,若R=xcorr(x,‘unbiased’),则通过命令:
plot(t,R)即可得到该时域信号的自相关函数曲线。
2、Cov
1。
向量的方差与协方差矩阵
cov(x)
求向量x的方差。
cov(x)为一个数值,数值大小计算公式为S(x)。
cov(x,y)
求向量x与y的协方差矩阵。
cov(x,y)为2*2矩阵,
[S(x)C(x,y);
C(y,x)S(y);]
2。
矩阵协方差矩阵
cov(X)
求矩阵X的协方差矩阵。
diag(cov(X))得到每一个列向量的方差。
sqrt(diag(cov(X)))得到每一个列的标准差。
若X大小为M*N,则cov(X)大小为N*N的矩阵。
cov(X)的第(i,j)个元素等于X的第i列向量与第j列向量的方差,即C(Xi,Xj)。
cov(X,Y)
求矩阵X与Y的协方差矩阵。
若X大小为M*N,Y为K*P,则X,Y的大小必须满足M*N=K*P,即X,Y的元素个数相同。
此时,cov(X,Y)等于cov([X(:
)Y(:
)])和cov(X(:
),Y(:
)),即计算两个向量的协方差矩阵,得到的结果为2*2矩阵。
[S(X(:
))C(X(:
),Y(:
));
C(Y(:
),X(:
))S(Y(:
));]
Matlab中cov函数详细解读-钰-计算机视觉·图像处理可知,S(X)=C(X,X).
3。
关于归一化的问题
在上述的S(X),C(X,Y)计算中,采用的归一化参数是1/(N-1),其中N是向量中元素的个数。
而下面的调用形式采用的归一化参数是1/N。
对应的公式如下图所示。
cov(x,1)
求向量x的方差。
计算方法如cov(x),但归一化参数为1/N。
cov(x,y,1)
求向量x与y的协方差矩阵。
计算方法如cov(x,y),但归一化参数为1/N。
Matlab中cov函数详细解读-钰-计算机视觉·图像处理
4。
PS:
为区别对待,
cov(x)又记作cov(x,0)
cov(x,y)又记作cov(x,y,0)
cov(X)又记作cov(X,0)
cov(X,Y)又记作cov(x,y,0)
对于归一化参数为1/(N-1)的情况,当N=1时,自动将参数调整为1/N。
3、Conv
1.我们假设有两个长度有限的任意序列A(n)和B(n),其中A(n)和B(n)的具体数学表达式可以看下图一。
那么这两个有限长序列的卷积就应该为C(n)=A(n)*B(n),其具体表达式请参。
2.相关函数指令
Matlab中的conv和deconv指令不仅可以用于多项式的乘除运算,还可以用于两个有限长序列的卷积和解积运算。
Matlab提供的函数conv,语法格式:
w=conv(u,v),其中u和v分别是有限长度序列向量,w是u和v的卷积结果序列向量。
如果向量u和v的长度分别为N和M,则向量w的长度为N+M-1.如果向量u和v是两个多项式的系数,则w就是这两个多项式乘积的系数。
下面我们看一下deconv指令。
功能:
求向量反褶积和进行多项式除法运算。
语法格式:
[q,r]=deconv(v,u),参数q和r分别返回多项式v除以多项式u的商多项式和余多项式。
具体实例请看下一步。
3.conv和deconv指令实例
具体实例请看下图,这里我们求多项式(x2+2x+1)与多项式(2x2+x+3)的积,再求积与(x2+2x+1)的商。
需要注意的是向量c代表多项式(2x4+5x3+7x2+7x+3)。
两个有限长序列的卷积实例
1.1具体序列的数学形式
在这一步我们将具体的有限长时间序列按数学方式显示,具体请看下图。
2.2解法一:
循环求合法求卷积
在本例中我们将按照原理方法第一步中图二的方式进行卷积计算,即循环求合法求卷积。
具体的代码及结果请看下图。
图一是是生成有限长度时间序列,图二是根据原理方法第一步中图二的方式即循环求合法求卷积的具体代码,图三是是图二的计算结果。
3.3解法二:
0起点序列法
下面就说一下第二种方法,即“0起点序列法”,所采用的指令就是我们在原理方法中介绍的conv函数指令。
具体代码看下图。
4.4解法三:
非平凡区间序列法
下面就说一下第二种方法,即“0起点序列法”,所采用的指令就是我们在原理方法中介绍的conv函数指令。
具体代码看下图。
图一为计算代码,图二为计算结果。
5.5绘图比较
这一步我们将解法二和解法三的计算结果绘制在一张图片中进行比较,其中第一幅是“0起点法”的计算结果图,第二幅是“非平凡区间法”的计算结果图。
其中画图代码为:
subplot(2,1,1),stem(kc,c),text(20,6,'0起点法')%画解法二的结果
CC=[zeros(1,KC
(1)),C];%补零是为了两子图一致
subplot(2,1,2),stem(kc,CC),text(18,6,'非平凡区间法')%画解法三的结果
xlabel('n')
6.6小结
有以上可以得出如下结论:
1、“解法三”最简洁、通用;
2、“解法二”使用于序列起点时刻N1或(和)M1小于0的情况,比较困难;
3、“解法一”最繁琐,效率低下。
4、Corrcoef
相关系数用来衡量两个数据集合是否在一条线上面。
对于一般的矩阵X,执行A=corrcoef(X)后,A中每个值的所在行a和列b,反应的是原矩阵X中相应的第a个列向量和第b个列向量的相似程度(即相关系数)。
计算公式是:
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表示最大的正相关,-1表示绝对值最大的负相关
命令相关系数
函数corrcoef
格式corrcoef(X,Y)%返回列向量X,Y的相关系数,等同于corrcoef([XY])。
corrcoef(A)%返回矩阵A的列向量的相关系数矩阵
例4-48
>>A=[123;40-1;139]
A=
123
40-1
139
>>C1=corrcoef(A)%求矩阵A的相关系数矩阵
C1=
1.0000-0.9449-0.8030
-0.94491.00000.9538
-0.80300.95381.0000
>>C1=corrcoef(A(:
2),A(:
3))%求A的第2列与第3列列向量的相关系数矩阵
C1=
1.00000.9538
0.95381.0000
5、Filter
autocorr
功能:
计算并描绘时间序列的自相关函数
格式:
autocorr(Series,nLags,M,nSTDs)%计算并绘制单变量随机时间序列的样本ACF及置信区间,如果不想绘制置信区间,则设置nSTDs=0
[ACF,lags,bounds]=autocorr(Series,nLags,M,nSTDs)%计算并返回ACF
说明:
Series--时间序列
nLags--延迟,当nLags=[]或缺省时,计算ACF时在延迟点0、1、2、。
。
。
、T处,T=min([20length(Series-1)])
M--非负整数,表示在多大延迟时理论ACF为0.autocorr假设序列为MA(M),并且使用Bartlett估计方法来计算大于M的延迟的标准误差。
如果M=[]或缺省,则为0,函数假设序列为高斯白噪声。
nSTDs--样本ACF估计误差的标准差。
ACF--样本自相关函数
Lags--与ACF(0,1,2,。
。
。
,nLags)相对应的延迟
Bounds--置信区间的近似上下限,假设序列为MA(M)过程。
如
rng('default') %makeoutputreproducible使结果可复现
x=randn(1000,1); %1000Gaussiandeviates~N(0,1).
y=filter([1-11],1,x); %CreateanMA
(2)process.创建MA
(2)过程
[ACF,lags,bounds]=autocorr(y,[],2) %ComputetheACFwith95percentconfidence.就算95%置信度下的相关系数
autocorr(y,[],2) %绘制结果
可以看到在延迟在2以后不再显著,说面这可能是一个MA
(2)过程。
Note:
此函数在判断MA过程时常用。
crosscorr
功能:
计算并描绘两时间序列的互相关函数
格式:
crosscorr(Series1,Series2,nLags,nSTDs)
[ACF,lags,bounds]=crosscorr(Series1,Series2,nLags,nSTDs)
说明:
参数具体说明同autocorr函数
如
%创建100个高斯分布的元素的时间序列
rng('default') %makeoutputreproducible
x=randn(100,1); %100Gaussiandeviates,N(0,1)
%创建上面序列的延迟副本
y=lagmatrix(x,4); %Delayitby4samples
y(isnan(y))=0;
%计算XCF互相关函数
[XCF,Lags,Bounds]=crosscorr(x,y);
crosscorr(x,y)
可以看到延迟在4处很显著。
Note:
在信号处理中经常要研究两个信号的相似性,或研究一个信号经过一段延迟后与自身的相似性,以实现信号的检测、识别与提取等。
相关函数是描述随机信号的重要统计量,有着广泛的用途,例如噪声中信号的检测、信号中隐含周期性的检测、信号相关性的检测、信号时延长度的测量等。
parcorr
功能:
计算并描绘两时间序列的部分自相关函数
格式:
parcorr(Series,nLags,M,nSTDs)
[ACF,lags,bounds]=parcorr(Series,nLags,M,nSTDs)
说明:
参数具体说明同autocorr函数,差别是parcorr是相对于AR模型的
如
rng('default') %makeoutputreproducible
x=randn(1000,1);
y=filter(1,[1-0.60.08],x);%创建AR
(2)过程的序列
[PartialACF,Lags,Bounds]=parcorr(y,[],2);
parcorr(y,[],2)
可以看到在延迟在2以后不再显著,说面这可能是一个AR
(2)过程。
Note:
此函数在判断AR过程时常用。
ksdensity
功能:
时间序列的概率密度估计,即使用非参数方法估计密度来描述数据样本。
格式:
[f,xi]=ksdensity(x) %计算时间序列x中样本数据的概率密度,f为样本点xi处的概率密度
f=ksdensity(x,xi) %f为样本点xi处的概率密度
如
x=[randn(30,1);5+randn(30,1)];
[f,xi]=ksdensity(x);
plot(xi,f);
normplot
功能:
描绘时间序列的正态概率图。
格式:
normplot(x) %描绘时间序列的正态概率图
h=normplot(x)%返回一个图中直线的句柄
如
x=normrnd(10,1,250,1);
normplot(x)
Note:
正态概率图的目的是从图中看序列是否来自正态分布,如果数据是正态分布的,那么正态概率曲线将是一条直线,而其他类型的分布是曲线。
wblplot
功能:
描绘时间序列的威布尔分布图,判断时间序列是否服从威布尔分布。
格式:
wblplot(x) %生成威布尔概率分布图,要求x的值大于0
h=wblplot(x)
如
r=wblrnd(1.2,1.5,50,1);
wblplot(r)
Note:
威布尔分布图的目的是从图中看序列是否来自威布尔分布,如果数据是威布尔分布的,那么威布尔分布曲线将是一条直线,而其他类型的分布是曲线
dt=0.001;
t=0:
0.001:
pi;
x1=sin(t);
x2=sin(t-(pi/2));
[a1,b1]=xcorr(x1,'unbiased');
[a2,b2]=xcorr(x2,'unbiased');
[a3,b3]=xcorr(x1,x2,'unbiased');
figure
subplot(3,1,1),plot(b1*dt,a1),title('xcorr(sin(t))');
subplot(3,1,2),plot(b2*dt,a2),title('xcorr(sin(t-(pi/2)))');
subplot(3,1,3),plot(b3*dt,a3),title('xcorr(sin(t),sin(t-(pi/2)))');
clf;
N=1000;%长度
Fs=500;%采样频率
n=0:
N-1;
t=n/Fs;%时间序列
lag=200;
%信号1
x=cos(2*pi*10*t);
x=x.*hanning(max(size(x)))';
%信号2
y=cos(2*pi*10*t+60*pi/180);
y=y.*hanning(max(size(y)))';
%互相关函数
[c,lags]=xcorr(x,y)%,lag);
subplot(211);
plot(t,x,'r');
holdon;
plot(t,y,':
');
legend('x信号','y信号');
xlabel('时间/s');ylabel('x(t)y(t)');
title('原始信号');gridon;
holdoff
subplot(212);
plot(lags/Fs,c,'r');
title('相关函数');
xlabel('时间/s');ylabel('Rxy(t)');
gridon
clear;
clc;
f0=10000;%用来模拟模拟信号的数字信号的采样频率fs< f=[1050100];%f是模拟信号的频率表max(f)<250; fs=500;%信号的采样频率 N=500;%数字信号的样点数 %模拟信号的生成 s=signal_generate(f,f0,N); subplot(2,1,1);plot(s);axis([1Nmin(s)max(s)]);gridon t=signal_generate1(f,f0,N); subplot(2,1,2);plot(t);axis([1Nmin(s)max(s)]);gridon [c,lags]=xcorr(s,t); [Am,Lm]=max(c); figure (2) plot(c) Lm ds=Lm-N delay=ds/fs %信号t functions=signal_generate(f,f0,N) f0=10000; num=length(f); s=zeros(1,N); fori=1: num s=s+sin(f(i)*2*pi*(1: N)/f0+pi/2); s=s.*hanning(max(siz
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 相关性
![提示](https://static.bingdoc.com/images/bang_tan.gif)