互相關（cross-correlation）及其在Python中的實現

阿新 • • 發佈：2018-12-10

互相關（cross-correlation）及其在Python中的實現

在這裡我想探討一下“互相關”中的一些概念。正如卷積有線性卷積（linear convolution）和迴圈卷積（circular convolution）之分；互相關也有線性互相關（linear cross-correlation）和迴圈互相關（circular cross-correlation）。線性互相關和迴圈互相關的基本公式是一致的，不同之處在於如何處理邊界資料。其本質的不同在於它們對原始資料的看法不同。通過這篇文章，我想整理一下相關概念，並給出示例。

1. 線性相關（Linear Cross-Correlation）的定義和計算

假設我們手裡有兩組資料，分別為 $M$ 個和 $N$ 個，表示為： ${\bf a} = \{a_i,i \in [0,M]\}$ 和 ${\bf v} = \{v_j,j \in [0,N]\}$ ， $\bf a$ 比 $\bf v$ 長，即 $M \geq N$ 。序列 ${\bf v}$ 和 ${\bf a}$ 之間的線性互相關操作表示為 $R_{\text L}({\bf v},{\bf a})$ ，其結果也是一個序列，表示為 ${\bf R}^\text{L}= \{R^\text{L}_k, k=0,1,2,...\}$ 。具體的操作是用這兩個序列進行的一種類似“滑動點積”的操作，如圖1和圖2所示。

圖1. 線性互相關的計算過程示意

圖2. 線性互相關結果序列中單個值計算示意

得到的互相關序列總長度為 $M+N-1$ ，該序列的前 $N-1$ 和後 $N-1$ 個數值是無效的，有效的資料共 $M-N+1$ 個。線性互相關的有效資料第 $k$ 個分量的值為：

$R^\text{L,eff}_{k}=\sum_{j=0}^{N-1}v_ja_{k+j},k \in[0,M-N+1]$

注意，線性互相關並不滿足交換律，即：

$R_{\text L}({\bf v},{\bf a}) \neq R_{\text L}({\bf a},{\bf v})$

一個簡單的應證是，等式兩側操作所得結果的有效資料個數都不一致。

線性相關的實際意義是，向量 ${\bf a}$ 中的各個與向量 ${\bf v}$ 等長的子向量與向量 ${\bf v}$

的相似程度。這樣， ${\bf R}^\text{L,eff}$ 中值最大的索引就是與向量 ${\bf a}$ 中與 ${\bf v}$ 最相似的子向量的起始索引。通常，為了獲得有效的互相關資料，我們總是用較短的資料去滑動點積較長的資料。

用一個實際的應用例子來驗證一下吧。如圖3的第一個子圖表示雷達聲納發射了一個探測訊號。經過一段時間之後，收到了如圖3的第二個子圖所示的回波（帶有一定的噪聲）。此時我們關注的是如何確定回波中從何時開始是對探測訊號的響應，以便計算目標距雷達的距離，這就需要用到線性互相關。在第三個子圖中的‘Valid’曲線即是有效互相關資料，其中清晰地呈現出兩處與探測訊號相似的回波的位置。

圖3. 相關計算的一個例子：雷達回波分析

線性互相關中，還有一些概念值得注意：
一是補零

。由線性相關的計算式不難發現，為了計算出個完整的相關係數序列（包含那些“無效資料”在內的所有結果），需要用到一些“不存在”的點。這就需要人為地對這些值進行補充，線上性相關的計算中，對這些超出原始資料儲存的區域取值為零。
二是末端效應。由圖1可以發現，一頭一尾的個互相關資料並沒有完全“嵌入”兩個原始陣列的全部資訊，它們或多或少地受到了人為補零的影響。因此一般認為這些資料是不可用的。
三是計算模式的選擇。這個問題其實是由問題二衍生而來的，就Python語言中的函式而言，至少有兩個可以直接計算線性相關：

1	`numpy.correlate(a, v, mode)`

和

1	`scipy.signal.correlate(a, v, mode)`

它們的呼叫引數完全相同。在呼叫時有三種模式可供選擇，它們計算的內容是相同的，但是返回值長度各不相同：
mode = ‘valid’：只返回有效的那一部分相關資料，共$M-N+1$個；
mode = ‘same’：只返回與等長的那一部分相關資料，共$N$個；
mode = ‘full’：返回全部相關資料，共$M+N-1$個。
圖3的第三個子圖展示了這三種模式的計算結果，在那個例子中，‘valid’模式是最合適的。

2. 迴圈互相關（Circular Cross-Correlation）的定義和計算

迴圈互相關是表徵兩組等長的週期性資料之間相似性的操作，其與線性互相關的區別也正由“等長”和“週期性”這個兩特點產生。在迴圈互相關中，被處理的原始資料是等長的，即 ${\bf a} = \{a_i,i \in [0,N]\}$ 和 ${\bf v} = \{v_j,j \in [0,N]\}$ 。序列 ${\bf v}$ 和 ${\bf a}$ 之間的線性互相關操作表示為 $R_{\text C}({\bf v},{\bf a})$ ，其結果也是一個序列，表示為 ${\bf R}^\text{C}= \{R^\text{L}_k, k=0,1,2,...\}$ 。其計算式與線性互相關的寫法是一致的：

$R^\text{C}_{k}=\sum_{j=0}^{N-1}v_ja_{k+j},k \in[0,N]$

只是得到的互相關序列長度也為 $N$ 。迴圈互相關的計算的具體過程如圖4所示，注意到在計算時要用到超出原始資料索引範圍的資料，其資料補充方式並不是“補零”而是“週期延拓”：即 $a_{N+j} = a_j$ 。這意味著對於迴圈互相關，不存在不同的計算模式之分，所有的資料都是有效資料。

圖4. 迴圈互相關的計算過程示意

注意，迴圈互相關也不滿足交換律。

這裡給出了一個關於迴圈相關的算例。兩路原始資料分別由如下函式生成：

$y_1=100sin(2\pi t)$

$y_1=3sin(2\pi t)+2sin(4\pi t)+Random(t)$

如果視 $y_1$ 為某個線性系統的週期輸入訊號，而視 $y_2$ 為這個線性系統的輸出訊號。由於存在外接干擾，因此輸出訊號不完全由輸入訊號決定。此時，迴圈互相關的實際意義是，分辨輸出訊號中的哪一個部分（頻率成分）是由該輸入訊號產生的。

圖5. 時域資料，從上到下： $y_1$ ， $y_2$ 和他們的迴圈互相關

圖6. 頻譜，從上到下： $y_1$ ， $y_2$ 和他們的迴圈互相關

從圖5和圖6可以看出，迴圈互相關的頻譜準確地說明了那些測試訊號的相關性。

遺憾的是，在Python幾大數值計算庫中，並沒有直接可計算迴圈相關的函式。但是可以採用如下程式碼構造出一個可用的（經過歸一化的）cxcorr(a, v)函數出來：

1 2 3 def cxcorr(a,v): nom = np.linalg.norm(a[:])*np.linalg.norm(v[:]) return fftpack.irfft(fftpack.rfft(a)*fftpack.rfft(v[::-1]))/nom

圖4中的資料就是通過這個函式計算出來的。其中用到了傅立葉變換和反變換來計算迴圈互相關，這是可行的。它們之間的關係在第四小節的QA中專門討論。

3. 用線性互相關處理週期性訊號

實際上，線性相關也可以處理週期訊號，前提是將兩組訊號取樣成不長度差異較大的序列。這樣，其有效線性互相關也可以完美地反應資料之間的相關性。

同樣採用第二節中的例子。這時為了保證足夠的有效線性互相關資料，兩組資料的長度故意不一致（但都足夠表徵其特徵），如圖7所示。它們的頻譜如圖8所示，仍然完美地體現了測試資料的相關性。

圖7. 時域資料，從上到下： $y_1$ ， $y_2$ 和他們的線性互相關

圖8. 頻譜，從上到下： $y_1$ ， $y_2$ 和他們的線性互相關

既然線性互相關也能處理週期性資料，為什麼還要專門搞一個基於等長序列和週期延拓的迴圈互相關呢？實際上，正如後文QA中專門討論的，這是為了利用快速傅利葉變換加速計算。

4. 相關問題QA

至此，兩種常用的互相關評價方法及其計算已經總結完畢。然而其中還有一些細節尚待分辨。例如，序列 ${\bf a}$ 和 ${\bf v}$ 之間的互相關的計算式：

$R({\bf v},{\bf a})_{k}=\sum_{j=0}^{N-1}v_ja_{k+j}$

與卷積（convolution）的定義式：

$Cov({\bf v},{\bf a})_{k}=\sum_{j=0}^{N-1}v_ja_{k-j}$

如此類似，如果再聯想起傅立葉變換的卷積定理，那麼，至少會產生如下的問題：

Q.1：它們之間有更深意義上的聯絡嗎？
A.1：文獻[1]的答覆是堅決的：“不要讓求卷積和互相關的數學相似性迷惑你，它們描述了不同的訊號處理過程。卷積是系統輸入訊號、輸出訊號和衝激響應之間的關係。互相關是一種在噪聲背景下檢測已知訊號的方法。二者在數學上的相似僅僅是一種巧合。”實際上，只要注意到卷積操作是滿足交換律的，而互相關操作並不滿足交換律。僅此一點也許就能說明它們有著本質的不同吧。

Q.2：可以利用Python中計算卷積的函式來計算互相關嗎？
A.2：可以，但是隻能用以計算線性互相關。Python中的numpy.convolve()函式就可以計算兩個序列之間的卷積。在卷積的計算過程中也會自動進行補零（而不是週期延拓，這就是為什麼只能計算線性相關的原因），這種卷積有時被稱為線性卷積，同樣涉及末端效應、有效資料長度等考慮。具體地，根據相關和卷積的表示式，如果希望計算序列 ${\bf v}$ 和 ${\bf a}$ 之間的線性互相關序列。等效地，只需要計算序列 ${\bf v}$ 和 ${\bf a[::-1]}$ 之間的卷積。 ${\bf a[::-1]}$ 表示序列 ${\bf a}$ 的“反置”，即將序列[1,2,3]反置為[3,2,1]。

Q.3：可以根據傅立葉變換的性質中有卷積定理，利用傅立葉正/逆變換計算互相關嗎？
A.3：可以，但是隻能用於計算迴圈互相關。傅立葉變換的卷積定理中所涉及的卷積是迴圈卷積。與前述的線性卷積是不同的。實際上不同的並不是卷積本身，它們的計算式是一致的，而是在如何看待參與卷積計算的資料，線性卷積認為參與計算的序列之外都是零，而迴圈卷積認為參與計算的序列是一個無限迴圈的資料的一段——這導致了它們對“越界”資料的補齊方式不一樣。正如線性互相關和迴圈互相關的區別！先將迴圈互相關等效為一個迴圈卷積，再利用快速傅立葉變換計算卷積即可。實際上本文給出的cxcorr(a, v)函式正是利用這一性質來計算迴圈相關的。其對計算速度的提升是相當明顯的。

Q.4：怎樣進行歸一化（normalization），以便於比較互相關資料？
A.4：根據參考[4]，用公式：

$R_\text{norm}=\frac{R({\bf v},{\bf a})}{|\bf v||\bf a|}$

5. 參考資料

[1] Steven W. Smith. Digital Signal Processing: A Practical Guide for Engineering and Scientists [M].
張瑞峰, 詹敏晶等譯. 實用數字訊號處理，從原理到應用[M]. 人民郵電出版社, 北京, 2010.
[2] Mark Owen. Practical Signal Processing [M].
丘天爽, 李麗, 趙林譯. 實用訊號處理 [M]. 電子工業出版社, 北京, 2009.
[3] 關於MATLAB中的xcorr() 的論述
http://www.mathworks.cn/cn/help/signal/ref/xcorr.html
[4] 關於MATLAB中的cxcorr() 的論述
http://www.mathworks.com/matlabcentral/fileexchange/4810-circular-cross-correlation
[5] 網路論壇Stackoverflow關於此問題的討論
http://stackoverflow.com/questions/6991471/computing-cross-correlation-function
http://stackoverflow.com/questions/12323959/fast-cross-correlation-method-in-python
http://stackoverflow.com/questions/9281102/n-fold-fft-convolution-and-circular-overlap
http://stackoverflow.com/questions/6855169/convolution-computations-in-numpy-scipy
http://stackoverflow.com/questions/4688715/find-time-shift-between-two-similar-waveforms
[6] 關於Cross-correlation的定義
http://mathworld.wolfram.com/Cross-Correlation.html
http://paulbourke.net/miscellaneous/correlate/
http://en.wikipedia.org/wiki/Cross-correlation
[7] 關於 Circular Cross-correlation的定義
http://en.wikipedia.org/wiki/Circular_convolution
http://cnx.org/content/m22974/latest/

本文轉載自：https://fanyublog.wordpress.com/2015/11/16/corr_python/

互相關（cross-correlation）及其在Python中的實現

互相關（cross-correlation）及其在Python中的實現

互相關（cross-correlation）及其在Python中的實現

機器學習之決策樹（Decision Tree）及其Python程式碼實現

氣泡排序演算法（起泡排序）及其C語言實現

雙目立體匹配——歸一化互相關（NCC）

【尋優演算法】交叉驗證（Cross Validation）引數尋優的python實現：多引數尋優

【尋優演算法】交叉驗證（Cross Validation）引數尋優的python實現：單一引數尋優

Python 資訊熵條件資訊熵互資訊（資訊增益）的理解以及計算程式碼

雙目立體匹配流程——歸一化互相關（NCC）

sugar與阿龍的互懟（第一季）

6.6-2-數組與數據結構（用數組及其函數實現堆棧等數據結構）

Python學習系列（第一章）：Python 的簡介

Python學習筆記（十五）用Python獲取本地數據

顯卡相關（未整理）

hdu 4135 a到b的範圍中多少數與n互質（容斥）

在python中實現線性回歸（linear regression）

快速排序、程式碼實現（python3版）及其時間空間複雜度分析

Python基礎（4）：python中的特性入門篇（索引，切片，連線，重複，成員操作符）

LeetCode練習：蓄水池問題（內附Java和Python的實現方法）

redis進化五（1）：redis叢集的安裝使用（redis-cluster）及其原理

如何在python中實現整數的二進位制迴圈移位（附程式碼）

互相關（cross-correlation）及其在Python中的實現

互相關（cross-correlation）及其在Python中的實現

相關推薦