三、从小白到入门之功率谱分析（下）（基于Aires大神的功率谱程序浅析自相关函数）

liuyg123

二、总结

对于xcorr、autocorr函数的理解

在matlab中，xcorr与autocorr均可以求自相关函数，最基础、最灵活的方法是根据公式自己编写算法。对于给定的一组时间序列，通过某种变化，xcorr与autocorr所求的自相关可以达成一致，参考这篇博文（http://blog.sina.com.cn/s/blog_50cfd0fc0102v1ij.html）。autocorr是更高级一点的自相关函数，但缺乏灵活性，所以我们先来了解一下xcorr具体是如何计算的。
xcorr函数可以输入三种参数，分别是unbiased、biased、、coeff，以及缺省状态。我会给出两种思路来思考xcorr是如何计算的，这两种方法对于帮助我们理解xcorr有很大的帮助。首先，我们先给出一组时间序列，假设x=[11 2 4]，这是一个时间序列长度为3的时序，我们用大写数字来代表数据对应的序号，即11对应一，2对应二，4对应三。

第一、缺省状态 xcorr的用法为[r,r_n]=xcorr(x),通过此命令求得的结果为r=44 30 141 30 44;r_n=-2 -1 0 1 2,r_n的计算为x的时序长度为3，那么最小值为1-3，最大值为1+3，在最小值与最大值之间取所有整数，结果就是r_n。再说r的计算
1、当b(1)=-2时，我们需要找到x这个时间序列数据序号相减等于-2的时候对应的数据（后面的序号减去前面的），即一减三等于-2，所用到数据为(4,11)所以a(1)=4*11=44（注意顺序，（4,11），4在前）。
2、当b(2)=-1时，即一减二、二减三均等于-1，所以用到了两组数据，(2,11)、(4,2)，所以a(2)=2*11+4*2=30。
3、当b(3)=0时，即一减一、二减二、三减三均要等于0，所以用到了三组数据，(11,11)、(2,2)、(4,4)，所以a(3)=11^2+2^2+4^2=141。
4、当b(4)=1时，即二减一、三减二均等于1，所有用到两组数据，(11,2)、(2,4)，所以a(4)=11*2+2*4=30。
5、当b(5)=2时，即三减1等于2，所用到数据为（11,4），所以a(5)=11*4=44。
以上就是xcorr在缺省状态下的计算思路。
第二、参数为'unbiased'。[r,r_n]=xcorr(x,'unbiased')，通过此命令求得的结果为r=44 15 47 15 44。b的计算是相同的，不再多说。而对于r的计算，由缺省状态可以看出，a(1)、a(5)只用到了一组数据来计算，所以除以1，故是没有变化的。而a(2)、a(4)用到了两组数据，故除以2。a(3)用到了三组数据，所以除以3。
第三、参数为'biased'。[r,r_n]=xcorr(x,'biased')，通过此命令求得的结果为r=14.6667 10 47 10 14.6667。r的计算为所有数据除以时序长度。
第四。参数为'coeff'。[r,r_n]=xcorr(x,'coeff')，通过此命令求得的结果为r=0.3121 0.2182 1 0.2182 0.3121。r的计算为在三的基础上做归一化。

以上给出了xcorr所有参数的计算方式，而autocorr则是将缺省状态下的xcorr先归一化，然后取中间数以后的数据。这就是这两个函数的区别。

下面再给出根据公式编写算法的思路。首先我们要确定一个时滞j，j的取值范围为0到时序的长度减1。在此即为j=0:2
1 当j=0时，即为[11 2 4].*[11 2 4]
2 当j=1时，即为[11 2].*[2 4]
3 当j=2时，[11].*[4]
程序为
for j=0:2
   r1(j+1)=sum((x(1:n-j).*x(j+1:n)))/length(x(1:n-j));
end
通过循环得到一个r1，
再写一个循环
for j=1:2
   r2(j)=sum((x(1:n-j).*x(j+1:n)))/length(x(1:n-j));
end
得到一个r2，
再将r2左右调换，最后得到相关系数r=[r2,r1]。
以上程序是根据参数'unbiased'编写的，缺省状态下要把分母删去，'biased'状态下分母改为length(x)。从上面可以看出，根据数据的不同，在使用matlab自带函数求自相关时，需要根据数据的不同区别不同参数。

通过自相关间接求出功率谱的完整程序如下：

alpha1=0.05;
alpha2=0.05;
m=48;
A=load('D:/gdy.txt');
x=A;

n=length(x);
x=zscore(x,1);%标准化

%求自相关
for j=0:143
   r1(j+1)=sum((x(1:n-j).*x(j+1:n)))/length(x(1:n-j));
end

for j=1:143
   r2(j)=sum((x(1:n-j).*x(j+1:n)))/length(x(1:n-j));
end

r2=fliplr(r2);
r=[r2,r1];

real_r=(length(r)+1)/2;
r=r(real_r:real_r+m);

l=0:m;
k=1:m-1;

a=r(2:end-1).*(1+cos(pi*k/m));
for i=2:m
    S(i)=(r(1)+sum(a.*cos(l(i)*pi*k/m)))/m;
end
S(1)=(r(1)+sum(a.*cos(l(1)*pi*k/m)))/(2*m);
S(m+1)=(r(1)+sum(a.*cos(l(end)*pi*k/m)))/(2*m);

plot(l,S,'-b','linewidth',2);
最后还有一点疑问，在Aires大神的程序中，求自相关的循环只截至到n-7，我一直没想明白为什么，所以我写程序是循环到length(n)-1，还望大神指教。谢谢。

饼脸大人

感谢楼主赐教！

风叶ele

alpha1=0.05;
alpha2=0.05;
m=48;


请问这个m（最大滞后时间长度）是如何确定为48的？
我看书上说是m一般取n/3 ~ n/10为宜。

RippleZ

风叶ele 发表于 2019-3-1 15:55
alpha1=0.05;
alpha2=0.05;
m=48;

自己算一下，143/3≈48

13top

楼主关于自相关函数的介绍很棒了

2daidai

感谢楼主的分享，顺便想问下红噪声检验线怎么弄？