在做訊號處理的時候，經常會對訊號做自相關處理，比如對訊號做功率估計，或者是引數擬合。在機器學習領域，如wule-walker方程也會遇到自相互函式的處理。

1  自相關矩陣的基本概念

     首先給出自相關函式的定義：

                              R(s,t)=E[(x(s)x(t))]

其中s,t是x(n)的不同時刻的訊號值。

    在實際的應用當中，我們把測得的資料一般是隨機訊號，比如說要測每天的10點，14點，20點的天氣溫度，測了N天以後，我們計算要10點x1溫度的自相關函式。那麼該怎麼計算？

    根據定義，我們需要知道x1和x2的概率密度f(x1，x2)，然後通過來求得x1和x2的自相關矩陣。顯然，這個方法是行不通的。因為我們不知道f(x)的概率密度。但是我們又需要求出X1的自相關矩陣。

    第二種方法就是測量無數天的10點溫度和12點溫度，然後帶入R(s,t)做平均。顯然這個方法也是行不通。因為這個方法所需要的資料是無窮的。

    第三種方法就是假設每天10點和12點的溫度是平穩隨機過程。利用平穩隨機過程的各態歷經性，即時間平均收斂於集平均。

給出，其中k=s-t。n從0到N-1。

2  自相關矩陣的matlab實現

    在matlab中，實現自相關的函式是xcorr()。

2.1 xcorr()

    計算時先進行b的計算，用序列x中的元素的序號互相做減法，可以得到的所有值的可能集合，按照從小到大順序排列後就得到了b；然後分別根據序號的“差”的情況計算序列a：
     當b(1)=-2時，只有資料(1, 5)作差可以得到，即序號1和序號3的差，因此計算a(1)=1*5=5；
     當b(2)=-1時，涉及到了序號對應的(3, 1)和序號(5, 3)，所以計算a(2)=3*1+5*3=18;
     當b(3)=0時，涉及到了序號對應的(1, 1), (3, 3)和(5, 5)，因此計算a(3)=1*1+3*3+5*5=35;
     當b(4)=1時，涉及到了序號對應的(3, 1)和(5, 3)，計算a(4)=3*1+5*3=18;

     當b(5)=2時，涉及到了序號對應的(5, 1)（後面的資料的序號減去前面資料的序號正好為2），計算a(5)=5*1=5

2.2 xcorr(x, 'unbiased')
     引數'unbiased'的作用在於基於預設引數時的計算結果，每個組的計算再除上該組的序號組數，比如b(1)時組數為1，記為N=1，則a(1)=1*5/N=5；b(2)時就是a(2)=18/N=18/2=9;類似等等；
 

點選(此處)摺疊或開啟

>> [a, b] = xcorr(x, 'unbiased')
 
a =
 
    5.0000 9.0000 11.6667 9.0000 5.0000
 
 
b =
 
    -2 -1 0 1 2
 
>>
3. xcorr(x, 'biased')
     引數'biased'的作用在於預設引數的基礎上除以序列x的長度，即a(1)=5/3；比如：
 

點選(此處)摺疊或開啟

>> [a, b] = xcorr(x, 'biased')
 
a =
 
    1.6667 6.0000 11.6667 6.0000 1.6667
 
 
b =
 
    -2 -1 0 1 2
 
>>
4. xcorr(x, 'coeff')
     此時用於求序列x的自相關序列，其結果是針對'biased'的情況進行歸一化，使得b=0時即中間的值a(3)=1，因此a(1)=5/11.6667，所有的分組資料在'biased'基礎上都通過11.6667歸一運算：
 

點選(此處)摺疊或開啟

>> [a, b] = xcorr(x, 'coeff')
 
a =
 
    0.1429 0.5143 1.0000 0.5143 0.1429
 
 
b =
 
    -2 -1 0 1 2
 
>>
 
 

3.編寫函式實現xocrr函式

根據1中提到的第三種方法計算自相關函式。u是輸入的資料

clear all;

u=[1 2 3 4];

for k=0:length(u)-1     %U的自相關函式,點數多的時候，計算速度慢等價於xcorr，biased
    tmp = 0;
    for i=1:length(u)-k
        tmp = tmp + u(i)*u(i+k);
    end
    rx(k+1) = tmp/length(u);

end

得到結果為：

用matlab的函式xcorr()

clear all;

u=[1 2 3 4];

rx1=xcorr(u,'biased');

rx1=rx1(length(u):end)

可以得到

rx1=xcorr(u,'biased');

rx1=rx1(length(u):end)

才是我們用第三種方法所求的相關函式。
--------------------- 
作者：黨小板 
來源：CSDN 
原文：https://blog.csdn.net/qq_41912125/article/details/82659653 
版權宣告：本文為博主原創文章，轉載請附上博文連結！