信號（互）相關及其應用 2014-02-16 17:35 6768人閱讀 收藏 舉報 本文章已收錄于： 分類： DSP（35） 作者同類文章X
版權聲明：本文為博主原創文章，未經博主允許不得轉載。
目錄(?) [+]
相關函數定義 基于定義的相關函數計算 使用FFT計算相關函數 應用 參考資料 在信號處理中，經常要研究兩個信號的相似性，或者一個信號經過一段時間延遲后自身的相似性，以便實現信號檢測、識別與提取等。
可用于研究信號相似性的方法稱為相關，該方法的核心概念是 相關函數和互相關函數 。
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 plaincopy print?
? ? ? ? ? ? ? ?? #include?<stdio.h> ???? typedef ?struct ?{??????float ?real;?? ????float ?imag;?? }?complex;?? ?? ?? static ?void ?assert_param(int32_t?x)??{?? ?? }?? ?? ? ? ? ? ? ? ? ? ? ? ? ?? 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--)?{???? ????????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;?? ????????}?? ?? ?? ????}?? ????for ?(k=0;?k<n;?k++)?{???? ????????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;?? ????????}?? ?? ?? ????}?? ?? ????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);?? ????}?? }??? /* * FileName : correl.c* Author : xiahouzuoxin xiahouzuoxin@163.com* Date : 2014/2/16* Version : v1.0* Compiler : gcc* Brief : */
#include <stdio.h>typedef struct {float real;float imag;
} complex;static void assert_param(int32_t x)
{}/*---------------------------------------------------------------------Routine CORRE1:To estimate the biased cross-correlation functionof 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 isstored 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
從0~8依次存儲的是m=-(N-1)到(N-1)的結果。 為驗證正確性，我們不妨用matlab自帶的xcorr計算
>> y = [1 2 4 5 6]
>> xcorr(x,y)
結果一致，只是存儲順序相反。
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 plaincopy print?
typedef ?complex?TYPE_CORREL;???? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? 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);?? ?? ?????? ????s?=?n;?? ????do ?{?? ????????s?=?s?>>?1;?? ????????k?=?s?-?2;?? ????}?while ?(k?>?0);?? ????if ?(k<0)?return ;?? ?? ?????? ????for ?(k=m;?k<n;?k++)?{?? ????????x[k].real?=?0.;?? ????????x[k].imag?=?0.0f;?? ????}?? ????fft(x,?n);?? ???????? ????if ?(1?==?icorrel)?{???? ?????????? ????????for ?(k=m;?k<n;?k++)?{?? ????????????y[k].real?=?0.;?? ????????????y[k].imag?=?0.0f;?? ????????}?? ????????fft(y,?n);?? ?? ?????????? ????????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 ;????? }?? 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;
}
FFT的程序參考本博客內博文FFT算法的完整DSP實現。
Matlab中自帶的xcorr也是通過FFT實現的互相關函數計算，這將互相關函數計算的復雜度降低到。
4 應用 互相關函數有許多實際的用途，比如通過不同的傳感器檢測不同方向到達的聲音信號，通過對不同方位傳感器間的信號進行互相關可計算聲音到達不同傳感器間的時延。自相關函數還可以用來計算周期信號的周期。對此，有時間將還會對此進行詳細研究。
參考資料
[1]? 《數字信號處理——理論、算法與實現》，胡廣書[2] ? 劉永春，陳琳. 基于廣義互相關時延估計算法聲源定位的研究.[3] ? 金中薇,姜明順等. 基于廣義互相關時延估計算法的聲發射定位技術. 傳感技術學報. 第26卷11期，2013年11月.
頂
4 踩
0 上一篇115家電子科技企業待遇一覽 下一篇關于Quartus ii無法識別Modelsim路徑的問題 我的同類文章 DSP（35） http://blog.csdn.net
?Kalman濾波器從原理到實現2014-09-26閱讀41330 ?白話壓縮感知（含Matlab代碼）2014-08-25閱讀14392 ?DSP-BIOS使用入門2014-07-24閱讀7698 ?對功率譜的一點理解2014-07-05閱讀4376 ?燒寫Flash后的DSP程序運行不正常的情況分析2014-04-12閱讀3964 ?TMS320C6713燒寫Flash的通用方法2014-03-30閱讀5755 ?自適應含噪信號過零率算法2014-09-11閱讀2758 ?DSP/BIOS使用之初窺門徑——滴答時鐘及燒寫Flash2014-07-25閱讀2270 ?DSP連接不上CCS3.3的問題討論2014-07-06閱讀2610 ?導出CCS3.3數據及使用matlab處理的方法2014-04-22閱讀4245 ?在DSP671x上使用Timer統計信號處理算法的時間消耗2014-03-30閱讀1450 更多文章
總結
以上是生活随笔 為你收集整理的信号互相关及其应用 的全部內容，希望文章能夠幫你解決所遇到的問題。
如果覺得生活随笔 網站內容還不錯，歡迎將生活随笔 推薦給好友。
