數字訊號處理中的自相關和互相關計算和物理意義（二）
在訊號處理中，經常要研究兩個訊號的相似性，或者一個訊號經過一段時間延遲後自身的相似性，以便實現訊號檢測、識別與提取等。
可用於研究訊號相似性的方法稱為相關，該方法的核心概念是相關函式和互相關函式。
1 相關函式定義
無限能量訊號，訊號x(n)與y(n)的互相關函式定義為
等於將x(n)保持不動，y(n)左移m個抽樣點後，兩個序列逐點對應相乘的結果。
當x(n)與y(n)不是同一訊號時，rxy中的x、y順序是不能互換等價的。
當x(n)與y(n)為同一訊號時，記
為訊號x(n)的自相關函式在m時刻的值。自相關函式反映了x(n)和其自身發生m個取樣點平移後的相似程度。
可以想象，當m=0時，即原訊號不做任何平移，一一對應的疊加時rx(m)值最大
對於有限能量訊號或週期訊號，設訊號為覆信號，自相關函式和互相關函式可表達為
注意：
（1）m的取值範圍可以從-(N-1)到(N-1)，對於N點訊號，rx共可計算得2N-1點相關函式結果值
（2）對於給定的m，因為實際訊號總是有限長的N，所以要計算rx(m)，n+m=N-1，因此實際寫程式時注意n的實際可取長度為N-1-m
（3）當m值越大時，對於N點有限長訊號，可用於計算的訊號長度越短，計算出的rx(n)效能越差，因此實際應用中常令m<<N
（4）Matlab自帶的xcorr函式可用於計算互相關，但計算結果是沒有除以N的結果。
2 基於定義的相關函式計算
[cpp]
view plain
copy
print?
- /*
- * FileName : correl.c
- * Author : xiahouzuoxin [email protected]
- * Date : 2014/2/16
- * Version : v1.0
- * Compiler : gcc
- * Brief :
- */
- #include <stdio.h>
- typedefstruct {
- float real;
- float imag;
- } complex;
- staticvoid assert_param(int32_t x)
- {
- }
- /*---------------------------------------------------------------------
- Routine CORRE1:To estimate the biased cross-correlation function
- of complex arrays x and y. If y=x,then it is auto-correlation.
- input parameters:
- x :n dimensioned complex array.
- y :n dimensioned complex array.
- n :the dimension of x and y.
- lag:point numbers of correlation.
- output parameters:
- r :lag dimensioned complex array, the correlation function is
- stored in r(0) to r(lag-1).
- ---------------------------------------------------------------------*/
- void corre1(complex x[],complex y[],complex r[],int n,int lag)
- {
- int m,j,k,kk;
- assert_param(lag >= 2*n-1);
- for (k=n-1; k>0; k--) { /* -(N-1)~0 PART */
- kk = n-1-k;
- r[kk].real = 0.0;
- r[kk].imag = 0.0;
- for (j=k; j<n; j++) {
- r[kk].real += y[j-k].real*x[j].real+y[j-k].imag*x[j].imag;
- r[kk].imag += y[j-k].imag*x[j].real-y[j-k].real*x[j].imag;
- }
- // r[kk].real=r[kk].real/n;
- // r[kk].imag=r[kk].imag/n;
- }
- for (k=0; k<n; k++) { /* 0~(N-1) PART */
- kk = n-1+k;
- m = n-1-k;
- r[kk].real = 0.0;
- r[kk].imag = 0.0;
- for (j=0; j<=m; j++) {
- r[kk].real += y[j+k].real*x[j].real+y[j+k].imag*x[j].imag;
- r[kk].imag += y[j+k].imag*x[j].real-y[j+k].real*x[j].imag;
- }
- // r[kk].real=r[kk].real/n;
- // r[kk].imag=r[kk].imag/n;
- }
- return;
- }
- #define SIG_N 5
- complex x[SIG_N];
- complex y[SIG_N];
- complex r[2*SIG_N-1];
- int main(void)
- {
- int i = 0;
- x[1].real = 1;
- x[2].real = 2;
- x[3].real = 3;
- x[4].real = 4;
- x[0].real = 5;
- x[1].imag = 0;
- x[2].imag = 0;
- x[3].imag = 0;
- x[4].imag = 0;
- x[0].imag = 0;
- y[1].real = 2;
- y[2].real = 4;
- y[3].real = 5;
- y[4].real = 6;
- y[0].real = 1;
- y[1].imag = 0;
- y[2].imag = 0;
- y[3].imag = 0;
- y[4].imag = 0;
- y[0].imag = 0;
- corre1(x,y,r,5,9);
- for (i=0; i<2*SIG_N-1; i++) {
- printf("r[%d].real=%.2f, r[%d].imag=%.2f\n", i, r[i].real, i, r[i].imag);
- }
- }
r[0].real=4.00, r[0].imag=0.00
r[1].real=11.00, r[1].imag=0.00
r[2].real=24.00, r[2].imag=0.00
r[3].real=37.00, r[3].imag=0.00
r[4].real=54.00, r[4].imag=0.00
r[5].real=42.00, r[5].imag=0.00
r[6].real=37.00, r[6].imag=0.00
r[7].real=31.00, r[7].imag=0.00
r[8].real=30.00, r[8].imag=0.00
從0~8依次儲存的是m=-(N-1)到(N-1)的結果。為驗證正確性，我們不妨用matlab自帶的xcorr計算
>> y = [1 2 4 5 6]
y =
1 2 4 5 6
>> x = [5 1 2 3 4]
x =
5 1 2 3 4
>> xcorr(x,y)
ans =
30.0000 31.0000 37.0000 42.0000 54.0000 37.0000 24.0000 11.0000 4.0000
結果一致，只是儲存順序相反。
3 使用FFT計算相關函式
採用暴力的按定義計算訊號相關的方法的計算複雜度約O(N^2)，當資料點數N很大時，尤其在DSP上跑時耗時過長，因此採用FFT和IFFT計算互相關函式顯得尤為重要。
那麼，互相關函式與FFT之間又是一種什麼樣的關係呢？
設y(n)是x(n)與h(n)的互相關函式，
即
則，
誒，這不對啊，不是說兩個訊號時域的卷積才對應頻域的乘積嗎？難道時域的互相關和時域的卷積等價了不成？？
這裡說明下，通過推倒可以得到，相關於卷積的關係滿足：
不管如何，與直接卷積相差一個負號。這時，看清楚了，相關函式在頻域也不完全是乘積，是一個訊號的共軛再與原訊號乘積，這就是與“時域卷積頻域相乘不同的地方”。
所以，請記住這個有用的結論，
兩個訊號的互相關函式的頻域等於X訊號頻域的共軛乘以Y訊號的頻域。
我們就有計算互相關的新方法了：將訊號x(n)和h(n)都進行FFT，將FFT的結果相乘計算得互相關函式的FFT，在進行逆變換IFFT得到互相關函式y(m)。
[cpp] view plain copy print?
- typedef complex TYPE_CORREL;
- /*
- * @brief To estimate the biased cross-correlation function
- * of TYPE_CORREL arrays x and y.
- * the result will store in x, size of x must be >=2*m
- * @input params
- x : n dimensioned TYPE_CORREL array.
- y : n dimensioned TYPE_CORREL array.
- m : the dimension of x and y.
- n : point numbers of correlation.
- icorrel: icorrel=1, cross-correlation; icorrel=0, auto-correlation
- * @retval None
- *
- * ====
- * TEST OK 2013.01.14
- */
- void zx_xcorrel(TYPE_CORREL x[], TYPE_CORREL y[], int m, int n, int icorrel)
- {
- int s,k;
- TYPE_CORREL z;
- assert_param(n >= 2*m);
- /* n must be power of 2 */
- s = n;
- do {
- s = s >> 1;
- k = s - 2;
- } while (k > 0);
- if (k<0) return;
- /* Padding 0 */
- for (k=m; k<n; k++) {
- x[k].real = 0.;
- x[k].imag = 0.0f;
- }
- fft(x, n);
- if (1 == icorrel) {
- /* Padding 0 */
- for (k=m; k<n; k++) {
- y[k].real = 0.;
- y[k].imag = 0.0f;
- }
- fft(y, n);
- /* conjuction */
- for (k=0; k<n; k++) {
- z.real = x[k].real;
- z.imag = x[k].imag;
- x[k].real = (z.real*y[k].real + z.imag*y[k].imag)/(float)m;
- x[k].imag = (z.real*y[k].imag - z.imag*y[k].real)/(float)m;
- }
- } else {
- for (k=0; k<n; k++) {
- x[k].real = (x[k].real*x[k].real+x[k].imag*x[k].imag) / (float)(m);
- x[k].imag = 0.0f;
- }
- }
- ifft(x, n);
- return;
- }
Matlab中自帶的xcorr也是通過FFT實現的互相關函式計算，這將互相關函式計算的複雜度降低到。
4 應用
互相關函式有許多實際的用途，比如通過不同的感測器檢測不同方向到達的聲音訊號，通過對不同方位感測器間的訊號進行互相關可計算聲音到達不同感測器間的時延。自相關函式還可以用來計算週期訊號的週期。對此，有時間將還會對此進行詳細研究。
參考資料
[1] 《數字訊號處理——理論、演算法與實現》，胡廣書
[2] 劉永春，陳琳. 基於廣義互相關時延估計演算法聲源定位的研究.
[3] 金中薇,姜明順等. 基於廣義互相關時延估計演算法的聲發射定位技術. 感測技術學報. 第26卷11期，2013年11月.