信号互相关及其应用

最新推荐文章于 2024-02-27 11:33:35 发布

liusandian

最新推荐文章于 2024-02-27 11:33:35 发布

阅读量2.6k

点赞数

分类专栏：阵列信号处理

阵列信号处理专栏收录该内容

54 篇文章 49 订阅

订阅专栏

信号（互）相关及其应用

2014-02-16 17:35 6768人阅读评论(5) 收藏举报

本文章已收录于：

分类：

DSP（35）

作者同类文章 X

目录(?)[+]

相关函数定义
基于定义的相关函数计算
使用FFT计算相关函数
应用
参考资料

在信号处理中，经常要研究两个信号的相似性，或者一个信号经过一段时间延迟后自身的相似性，以便实现信号检测、识别与提取等。

可用于研究信号相似性的方法称为相关，该方法的核心概念是相关函数和互相关函数。

1 相关函数定义

无限能量信号，信号x(n)与y(n)的互相关函数定义为

$formdata=\Large+r_{xy}(m)=\sum_{n=-\infty}^{\infty}x(n)y(n+m)$

等于将x(n)保持不动，y(n)左移m个抽样点后，两个序列逐点对应相乘的结果。

$formdata=\Large+r_{yx}(m)=\sum_{n=-\infty}^{\infty}y(n)x(n+m)$

当x(n)与y(n)不是同一信号时，rxy中的x、y顺序是不能互换等价的。

当x(n)与y(n)为同一信号时，记

$formdata=\Large+r_{x}(m)=r_{xx}(m)=\sum_{n=-\infty}^{\infty}x(n)x(n+m)$

为信号x(n)的自相关函数在m时刻的值。自相关函数反映了x(n)和其自身发生m个采样点平移后的相似程度。

可以想象，当m=0时，即原信号不做任何平移，一一对应的叠加时rx(m)值最大，这个结论很重要。

对于有限能量信号或周期信号，设信号为复信号，自相关函数和互相关函数可表达为

$formdata=\Large+r_{x}(m)={\frac1_N}\sum_{n=0}^{N-1}x^*(n)x(n+m)$

$formdata=\Large+r_{xy}(m)={\frac1_N}\sum_{n=0}^{N-1}x^*(n)y(n+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 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 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);
}
}

/* 
 * 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 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)的互相关函数，

即

$formdata=\Large+y(m)={\frac1_N}\sum_{n=0}^{N-1}x(n)h(n+m)$

则，

$formdata=\Large+Y(e^{jw})=X^*(e^{jw})H(e^{jw})$

诶，这不对啊，不是说两个信号时域的卷积才对应频域的乘积吗？难道时域的互相关和时域的卷积等价了不成？？

这里说明下，通过推倒可以得到，相关于卷积的关系满足：

$formdata=\Large+r_{xy}(m)=x(-m)\ast{y(m)}$

不管如何，与直接卷积相差一个负号。这时，看清楚了，相关函数在频域也不完全是乘积，是一个信号的共轭再与原信号乘积，这就是与“时域卷积频域相乘不同的地方”。

所以，请记住这个有用的结论，

两个信号的互相关函数的频域等于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;
}

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实现的互相关函数计算，这将互相关函数计算的复杂度降低到 $formdata=\small+\frac+N2+\log_2{N}$ 。