DSP 相關 Correlation

包含交互相關 cross-correlation 以及自相關 autocorrelation 及其應用。

correlation 相關就是兩個隨機變數之間，是否具有關聯性。相關係數可用來測量兩個隨機變數集合的相關性，通常介於 -1 (負相關) ~ +1 (正相關) 之間。

交互相關 cross correlation

給定兩個函數 \(x(t), h(t)\)，連續時間的 cross correlation 定義為：

\(y(t) = \int_{-∞}^{∞}x^*(τ)h(t+τ)dτ\) 或

\(y(t) = x(t) ⊗ h(t)\)

其中 \(x^*(τ)\) 是共軛複數 complex conjugate 。但數位訊號的輸入資料都是實數，所以 cross correlation 跟卷積運算類似，只是把 \(h(t-τ)\) 換成 \(h(t+τ)\)

數位訊號 x[n], h[n] 的交互相關定義為：

\( y[n] = \sum_{k=-∞}^{∞}x^*[k] \cdot h[n+k] \) 或

\( y[n]=x[n] ⊗ h[n] \)

h[n] 可視為第二個數位訊號，用來進行訊號比對

因輸入訊號的數量 N=7，脈衝響應數量為 M = 5，full cross-correlation 結果長度為 \(M+N-1 = 11\)

計算前先將 \(x[n]\) 兩邊補上 \(M-1=4\) 個 0，稱為 Zero-Padding，但不像卷積運算一樣，不將 h(n) 旋轉 180 度

x	h
0 0 0 0 1 2 4 3 2 1 1 0 0 0 0	1 2 3 1 1
1 2 3 1 1
\(y[0] = 0 \cdot 1 + 0 \cdot 2 + 0 \cdot 3 + 0 \cdot 1 + 1 \cdot 1 = 1 \)

x	h
0 0 0 0 1 2 4 3 2 1 1 0 0 0 0	1 2 3 1 1
1 2 3 1 1
\(y[1] = 0 \cdot 1 + 0 \cdot 2 + 0 \cdot 3 + 1 \cdot 1 + 2 \cdot 1 = 3 \)

最後卷積運算結果為

n	0	1	2	3	4	5	6	7	8	9	10
y[n]	1	3	9	15	22	22	18	11	7	3	1
			✓	✓	✓	✓	✓	✓	✓

實際運算時，會擷取卷積運算的部分結果

n	0	1	2	3	4	5	6
x[n]	1	2	4	3	2	1	1
y[n]	9	15	22	22	18	11	7

結果在 n=2, 3 的位置，數值最大，也就是 x[n] 與 h[n] 在這個位置的相關性最高，波形最相似。

cross correlation 經典的應用是訊號比對 Signal Matching，例如：音樂搜尋系統，聲紋比對系統

numpy 有提供 full, same 兩種 cross correlation 運算的結果

import numpy as np

x = np.array( [ 1, 2, 4, 3, 2, 1, 1 ] )
h = np.array( [ 1, 2, 3, 1, 1 ] )
y = np.correlate( x, h, 'full' )
y1 = np.correlate( x, h, 'same' )

print( "x =", x )
print( "h =", h )
print( "Full Correlation y =", y )
print( "Correlation y =", y1 )

執行結果

$ python correlation.py
x = [1 2 4 3 2 1 1]
h = [1 2 3 1 1]
Full Correlation y = [ 1  3  9 15 22 22 18 11  7  3  1]
Correlation y = [ 9 15 22 22 18 11  7]

如果 h[n] 左右對稱且輸入的數位訊號是實數，則 cross correlation 的結果跟卷積運算一樣。另外 cross correlation 跟卷積一樣，滿足交換率、結合率、分配率。

自相關 autocorrelation

給定兩個函數 \(x(t), h(t)\)，連續時間的 autocorrelation 定義為：

\(y(t) = \int_{-∞}^{∞}x^*(τ)x(t+τ)dτ\) 其中 τ 稱為延遲 lag

其中 \(x^*(τ)\) 是共軛複數 complex conjugate ，自相關 autocorrelation 是延遲 τ 的函數。但數位訊號的輸入資料都是實數，所以 autocorrelation 可解釋為：訊號與訊號本身的延遲訊號，進行積分運算的結果。

數位訊號的卷積定義為：

\( R[l] = \sum_{k=-∞}^{∞}x^*[k] \cdot x[n+l] \) 其中 l 是延遲 lag

autocorrelation 是延遲 l 的函數，且 l 為正整數。可用來偵測數位訊號是否具有週期性。通常週期性的數位訊號，自關性較強。 noise 因為每個 sample 在時間軸上都是獨立的，所以不具有 autocorrelation

因輸入訊號的數量 N=5，脈衝響應數量為 M = 5，autocorrelation 結果長度為 5

計算前先將 \(x[n]\) 後面補上 \(N-1=4\) 個 0

x	h
1 2 1 2 1 0 0 0 0	1 2 1 2 1
1 2 1 2 1
\(R[0] = 1 \cdot 1 + 2 \cdot 2 + 1 \cdot 1 + 2 \cdot 2 + 1 \cdot 1 = 11 \)

x	h
1 2 1 2 1 0 0 0 0	1 2 1 2 1
1 2 1 2 1
\(y[1] = 2 \cdot 1 + 1 \cdot 2 + 2 \cdot 1 + 1 \cdot 2 + 0 \cdot 1 = 8 \)

最後 autocorrelation 運算結果為

n	0	1	2	3	4
R[n]	11	8	6	4	1

cross correlation 運算結果的後半部分，就是 autocorrelation 的結果

import numpy as np

def autocorr( x ):
    # cross correlation 運算結果的後半部分，就是 autocorrelation
    R = np.correlate( x, x, 'full' )
    return R[ int( R.size / 2 ) : ]

def main( ):
    x = np.array( [ 1, 2, 1, 2, 1 ] )
    R = autocorr( x )
    print( "x =", x )
    print( "Autocorrelation =", R )

main( )

執行結果

$ python autocorrelation.py
x = [1 2 1 2 1]
Autocorrelation = [11  8  6  4  1]

autocorrelation 的應用

週期性訊號的自相關性較強，雜訊不具有自相關性。自相關適合用來偵測訊號的週期性，即使週期性的訊號混入了雜訊，一樣能找出自相關性。

import numpy as np
import numpy.random as random
import matplotlib.pyplot as plt

def autocorr( x ):
    R = np.correlate( x, x, 'full' )
    return R[ int( R.size / 2 ) : ]

def main( ):
    t = np.linspace( 0, 1,  100, endpoint = False )     # 定義時間陣列
    x = 10 * np.cos( 2 * np.pi * 5 * t )                # 原始訊號
    noise = random.uniform( -2, 2, 100 )                # 雜訊(均勻分佈)  
    y = x + noise

    auto_corr1 = autocorr( x )
    auto_corr2 = autocorr( noise )
    auto_corr3 = autocorr( y )

    plt.figure( 1 )                                     
    plt.subplot( 121 )
    plt.plot( t, x )
    plt.xlabel( 't (second)' )
    plt.ylabel( 'Amplitude' )

    plt.subplot( 122 )                                  
    plt.plot( auto_corr1 )
    plt.xlabel( 'Lag' )
    plt.ylabel( 'Auto Correlation' )

    plt.figure( 2 )
    plt.subplot( 121 )
    plt.plot( t, noise )
    plt.xlabel( 't (second)' )
    plt.ylabel( 'Amplitude' )

    plt.subplot( 122 )
    plt.plot( auto_corr2 )
    plt.xlabel( 'Lag' )
    plt.ylabel( 'Auto Correlation' )

    plt.figure( 3 ) 
    plt.subplot( 121 )
    plt.plot( t, y )
    plt.xlabel( 't (second)' )
    plt.ylabel( 'Amplitude' )

    plt.subplot( 122 )
    plt.plot( auto_corr3 )
    plt.xlabel( 'Lag' )
    plt.ylabel( 'Auto Correlation' )

    plt.show( )
    
main( )

左圖 x[n] 是原始訊號，右邊是 autocorrelation 的結果