互相關（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語言中的函數而言，至少有兩個可以直接計算線性相關：

 
                 numpy.correlate(a, v, mode)

和

 
                 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)函數出來：

 
                 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中的實現

免責聲明！