數學原理
在數字信號處理中，相關(correlation)可以分為互相關(cross correlation)和自相關(auto-correlation). 互相關是兩個數字序列之間的運算；自相關是單個數字序列本身的運算，可以看成是兩個相同數字序列的互相關運算．互相關用來度量一個數字序列移位后，與另一個數字序列的相似程度．其數學公式如下：
其中，f 和 g 為數字序列，n 為移位的位數，f* 表示 f 序列值的復數共軛，即復數的實部不變，虛部取反．
而卷積（convolution）與互相關運算相似，定義為將其中一個序列反轉並移位后，兩個序列的乘積的積分（求和），其數學公式如下：
其中，f 和 g 為數字序列，n 為移位的位數.
在實數范圍內，f 的復數共軛 f* = f ．此時，通過比較上面兩式可知：序列 f 與將序列 g 反轉后的序列的卷積為序列 f 與序列 g 的互相關
Python 實現
采用兩種方式實現：自定義互相關函數和直接調用 numpy.correlate 或 numpy.convolve．
在 numpy 中, numpy.correlate 函數實現兩個一維數組的互相關操作；numpy.convolve 實現了兩個一維數組的卷積操作．其中定義了三種模式（'valid', 'same','full'）.
設兩個序列長度分別為 M 和 N，則
- 'valid' 模式：輸出長度為 max(M,N)-min(M,N)+1．只返回兩個序列完全重合部分的點的卷積或相關運算；
- 'same' 模式：輸出長度為兩個序列中的較長者，即 max(M,N)；
- 'full' 模式：輸出長度為 M+N-1, 返回所有包含重疊部分的點.
互相關或卷積，實際上，就是計算兩個序列（一維數組）在不同移位情況下，兩個序列逐位相乘之后，求和的結果．不同模式只是返回互相關或卷積結果的不同部分．
注：如果在超出數組的索引范圍，用 0 填充．
下面代碼，采用自定義函數 correlate_func (只適用於實數值) 實現 numpy.correlate 和 numpy.convolve 的三種模式，並進行測試．
#!//usr/bin/env python # -*- coding: utf8 -*-
""" # Author: klchang # Description: correlation or convolution of one-dimensional array with real numbers. # Date: 2018.11 """
from __future__ import print_function import numpy as np def correlate_func(a, b, mode='valid', conv=True): '''correlation or convolution in 1-d array with real numbers'''
if a is None or b is None: return None if len(a) > len(b):# Ensure the length of a is no longer than that of b.
return correlate_func(b, a, mode) # Convert to np.array type
a, b = list(map(np.array, [a, b])) if conv: a = a[::-1] # if convolution is true, reverse the shorter
res = [] min_len, max_len = len(a), len(b) if mode == 'valid': output_length = max_len - min_len + 1 tmp = b elif mode == 'same': output_length = max_len tmp = np.hstack((np.zeros(min_len-1), b)) elif mode == 'full': output_length = max_len + min_len - 1 tmp = np.hstack((np.zeros(min_len-1), b, np.zeros(min_len-1))) else: raise Exception("No such mode {}!".format(mode)) # For each point, get the total sum of element-wise multiplication
for i in range(output_length): val = np.sum(a * tmp[i:min_len+i]) res.append(val) return np.array(res, dtype=a.dtype) def test(): a = [1, 2, 3] b = [1, 2] names = ['numpy.correlate', 'correlate_func', 'numpy.convolve', 'correlate_func(convolution)'] funcs = [np.correlate, correlate_func, np.convolve, lambda *args: correlate_func(*args, conv=True)] for i, (name, func) in enumerate(zip(names, funcs)[:4]): print ('-----' * 30 if i & 0x01 == 0 else '') print ("{} output result: ".format(name)) print (' valid mode: ', func(a, b, 'valid')) print (' same mode: ', func(a, b, 'same')) print (' full mode: ', func(a, b, 'full')) if __name__ == '__main__': test()
除此之外，在 matplotlib.pyplot 模塊中，實現了用於可視化的自相關函數 matplotlib.pyplot.acorr 和互相關函數 matplotlib.pyplot.xcorr, 官方網址提供的一個示例代碼如下:
import matplotlib.pyplot as plt import numpy as np np.random.seed(0) x, y = np.random.randn(2, 100) fig = plt.figure() ax1 = fig.add_subplot(211) ax1.xcorr(x, y, usevlines=True, maxlags=50, normed=True, lw=2) ax1.grid(True) ax1.axhline(0, color='black', lw=2) ax2 = fig.add_subplot(212, sharex=ax1) ax2.acorr(x, usevlines=True, normed=True, maxlags=50, lw=2) ax2.grid(True) ax2.axhline(0, color='black', lw=2) plt.show()
參考資料
[1] Cross-correlation - Wikipedia. https://en.wikipedia.org/wiki/Cross-correlation
[2] Convolution - Wikipedia. https://en.wikipedia.org/wiki/Convolution
[3] Python: Interpretation on XCORR. https://stackoverflow.com/questions/24396589/python-interpretation-on-xcorr
[4] numpy.correlate - Numpy Reference. https://docs.scipy.org/doc/numpy/reference/generated/numpy.correlate.html
[5] numpy.convolve - Numpy Reference. https://docs.scipy.org/doc/numpy/reference/generated/numpy.convolve.html#numpy.convolve