1. 程式人生 > >Matlab 自相關檢測 :自相關函式xcorr

Matlab 自相關檢測 :自相關函式xcorr

原文:http://blog.chinaunix.net/uid-26275986-id-4342906.html

最近因為工作的關係需要使用matlab作為資料統計的工具,其中一個關鍵是使用其自相關函式獲得資料的估計。自己只在本科時候馬馬虎虎地學習了一點matlab,這次仗著有C/C++的基礎迅速地過了一遍自己需要的matlab的語法,原來這門語言很像指令碼啊,同Python一樣都是弱型別的,語法也不嚴格。瞭解了語法大概後,立刻在Help?文件中查詢xcorr函式,其介紹如下:

     看上去語法也不難,直接運算不就好了麼?可是運算出來的結果自己卻搞不懂,因為自己沒有多少統計的知識,於是又去巴拉數學的材料,想去搞明白xcorr函式的原理或公式。最後還是去matlab論壇找到了自己想找的答案,這裡就來分析下matlab的互相關函式xcorr。
     matlab中的引數都是以陣列的形式儲存的,標量可以看作是一維陣列。我們採用序列x = [1, 3, 5]作為實驗物件,經過xcorr()函式運算,分析結果:
1. xcorr()

點選(此處)摺疊或開啟

  1. >> x = [1 3 5]
  2. =
  3.      1 3 5
  4. >> [a,b] = xcorr(x)
  5. =
  6.      5 18 35 18 5
  7. =
  8.     --1 0 1 2
  9. >>
      也許你對這個結果感到困惑,不急,待我慢慢道來。計算時先進行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. xcorr(x, 'unbiased')
     引數'unbiased'的作用在於基於預設引數時的計算結果,每個組的計算再除上該組的序號組數,比如b(1)時組數為1,記為N=1,則a(1)=1*5/N=5;b(2)時就是a(2)=18/N=18/2=9;類似等等;

點選(此處

)摺疊或開啟

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

點選(此處)摺疊或開啟

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

點選(此處

)摺疊或開啟

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

      由於xcorr多用於工程上針對時間訊號取樣,但是計算時將採集到的資料一起送給matlab,因此matlab本身並不知道時間間隔,我們可以使用dt=0.1, t=b*dt,plot(t, a)進行作圖,前半部分是超前,後半部分是滯後,如: