信號(hào)（互）相關(guān)及其應(yīng)用 2014-02-16 17:35 6768人閱讀 收藏 舉報(bào) 本文章已收錄于： 分類： DSP（35） 作者同類文章X
版權(quán)聲明：本文為博主原創(chuàng)文章，未經(jīng)博主允許不得轉(zhuǎn)載。
目錄(?) [+]
相關(guān)函數(shù)定義 基于定義的相關(guān)函數(shù)計(jì)算 使用FFT計(jì)算相關(guān)函數(shù) 應(yīng)用 參考資料 在信號(hào)處理中，經(jīng)常要研究?jī)蓚€(gè)信號(hào)的相似性，或者一個(gè)信號(hào)經(jīng)過一段時(shí)間延遲后自身的相似性，以便實(shí)現(xiàn)信號(hào)檢測(cè)、識(shí)別與提取等。
可用于研究信號(hào)相似性的方法稱為相關(guān)，該方法的核心概念是 相關(guān)函數(shù)和互相關(guān)函數(shù) 。
1 相關(guān)函數(shù)定義 無限能量信號(hào)，信號(hào)x(n)與y(n)的 互相關(guān)函數(shù) 定義為
等于將x(n)保持不動(dòng)，y(n)左移m個(gè)抽樣點(diǎn)后，兩個(gè)序列逐點(diǎn)對(duì)應(yīng)相乘的結(jié)果。
當(dāng)x(n)與y(n)不是同一信號(hào)時(shí)，rxy中的x、y順序是不能互換等價(jià)的。
當(dāng)x(n)與y(n)為同一信號(hào)時(shí)，記
為信號(hào)x(n)的 自相關(guān)函數(shù) 在m時(shí)刻的值。自相關(guān)函數(shù)反映了x(n)和其自身發(fā)生m個(gè)采樣點(diǎn)平移后的相似程度。
可以想象， 當(dāng)m=0時(shí)，即原信號(hào)不做任何平移，一一對(duì)應(yīng)的疊加時(shí)rx(m)值最大 ，這個(gè)結(jié)論很重要。
對(duì)于有限能量信號(hào)或周期信號(hào)，設(shè)信號(hào)為復(fù)信號(hào)，自相關(guān)函數(shù)和互相關(guān)函數(shù)可表達(dá)為
注意：
（1）m的取值范圍可以從-(N-1)到(N-1)，對(duì)于N點(diǎn)信號(hào)，rx共可計(jì)算得2N-1點(diǎn)相關(guān)函數(shù)結(jié)果值
（2）對(duì)于給定的m，因?yàn)閷?shí)際信號(hào)總是有限長(zhǎng)的N，所以要計(jì)算rx(m)，n+m=N-1，因此實(shí)際寫程序時(shí)注意n的實(shí)際可取長(zhǎng)度為N-1-m
（3）當(dāng)m值越大時(shí)，對(duì)于N點(diǎn)有限長(zhǎng)信號(hào)，可用于計(jì)算的信號(hào)長(zhǎng)度越短，計(jì)算出的rx(n)性能越差，因此實(shí)際應(yīng)用中常令m<<N
（4）Matlab自帶的xcorr函數(shù)可用于計(jì)算互相關(guān)，但計(jì)算結(jié)果是沒有除以N的結(jié)果。
2 基于定義的相關(guān)函數(shù)計(jì)算
[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);}
}
運(yùn)行輸出結(jié)果如下，r[0].real=4.00, r[0].imag=0.00
從0~8依次存儲(chǔ)的是m=-(N-1)到(N-1)的結(jié)果。 為驗(yàn)證正確性，我們不妨用matlab自帶的xcorr計(jì)算
>> y = [1 2 4 5 6]
>> xcorr(x,y)
結(jié)果一致，只是存儲(chǔ)順序相反。
3 使用FFT計(jì)算相關(guān)函數(shù) 采用暴力的按定義計(jì)算信號(hào)相關(guān)的方法的計(jì)算復(fù)雜度約O(N^2)，當(dāng)數(shù)據(jù)點(diǎn)數(shù)N很大時(shí)，尤其在DSP上跑時(shí)耗時(shí)過長(zhǎng)，因此采用FFT和IFFT計(jì)算互相關(guān)函數(shù)顯得尤為重要。
那么，互相關(guān)函數(shù)與FFT之間又是一種什么樣的關(guān)系呢？
設(shè)y(n)是x(n)與h(n)的互相關(guān)函數(shù)，
即
則，
誒，這不對(duì)啊，不是說兩個(gè)信號(hào)時(shí)域的卷積才對(duì)應(yīng)頻域的乘積嗎？難道時(shí)域的互相關(guān)和時(shí)域的卷積等價(jià)了不成？？
這里說明下，通過推倒可以得到，相關(guān)于卷積的關(guān)系滿足：
不管如何，與直接卷積相差一個(gè)負(fù)號(hào)。這時(shí)，看清楚了，相關(guān)函數(shù)在頻域也不完全是乘積，是一個(gè)信號(hào)的共軛再與原信號(hào)乘積，這就是與“時(shí)域卷積頻域相乘不同的地方”。
所以，請(qǐng)記住這個(gè)有用的結(jié)論，
兩個(gè)信號(hào)的互相關(guān)函數(shù)的頻域等于X信號(hào)頻域的共軛乘以Y信號(hào)的頻域 。
我們就有計(jì)算互相關(guān)的新方法了：將信號(hào)x(n)和h(n)都進(jìn)行FFT，將FFT的結(jié)果相乘計(jì)算得互相關(guān)函數(shù)的FFT，在進(jìn)行逆變換IFFT得到互相關(guān)函數(shù)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的程序參考本博客內(nèi)博文FFT算法的完整DSP實(shí)現(xiàn)。
Matlab中自帶的xcorr也是通過FFT實(shí)現(xiàn)的互相關(guān)函數(shù)計(jì)算，這將互相關(guān)函數(shù)計(jì)算的復(fù)雜度降低到。
4 應(yīng)用 互相關(guān)函數(shù)有許多實(shí)際的用途，比如通過不同的傳感器檢測(cè)不同方向到達(dá)的聲音信號(hào)，通過對(duì)不同方位傳感器間的信號(hào)進(jìn)行互相關(guān)可計(jì)算聲音到達(dá)不同傳感器間的時(shí)延。自相關(guān)函數(shù)還可以用來計(jì)算周期信號(hào)的周期。對(duì)此，有時(shí)間將還會(huì)對(duì)此進(jìn)行詳細(xì)研究。
參考資料
[1]? 《數(shù)字信號(hào)處理——理論、算法與實(shí)現(xiàn)》，胡廣書[2] ? 劉永春，陳琳. 基于廣義互相關(guān)時(shí)延估計(jì)算法聲源定位的研究.[3] ? 金中薇,姜明順等. 基于廣義互相關(guān)時(shí)延估計(jì)算法的聲發(fā)射定位技術(shù). 傳感技術(shù)學(xué)報(bào). 第26卷11期，2013年11月.
頂
4 踩
0 上一篇115家電子科技企業(yè)待遇一覽 下一篇關(guān)于Quartus ii無法識(shí)別Modelsim路徑的問題 我的同類文章 DSP（35） http://blog.csdn.net
?Kalman濾波器從原理到實(shí)現(xiàn)2014-09-26閱讀41330 ?白話壓縮感知（含Matlab代碼）2014-08-25閱讀14392 ?DSP-BIOS使用入門2014-07-24閱讀7698 ?對(duì)功率譜的一點(diǎn)理解2014-07-05閱讀4376 ?燒寫Flash后的DSP程序運(yùn)行不正常的情況分析2014-04-12閱讀3964 ?TMS320C6713燒寫Flash的通用方法2014-03-30閱讀5755 ?自適應(yīng)含噪信號(hào)過零率算法2014-09-11閱讀2758 ?DSP/BIOS使用之初窺門徑——滴答時(shí)鐘及燒寫Flash2014-07-25閱讀2270 ?DSP連接不上CCS3.3的問題討論2014-07-06閱讀2610 ?導(dǎo)出CCS3.3數(shù)據(jù)及使用matlab處理的方法2014-04-22閱讀4245 ?在DSP671x上使用Timer統(tǒng)計(jì)信號(hào)處理算法的時(shí)間消耗2014-03-30閱讀1450 更多文章
總結(jié)
以上是生活随笔 為你收集整理的信号互相关及其应用 的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。
如果覺得生活随笔 網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔 推薦給好友。
