数字信号处理中的自相关和互相关计算和物理意义（二）

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

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

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 plaincopyprint?

/*
* 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 plaincopyprint?

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}$ 。

4 应用

互相关函数有许多实际的用途，比如通过不同的传感器检测不同方向到达的声音信号，通过对不同方位传感器间的信号进行互相关可计算声音到达不同传感器间的时延。自相关函数还可以用来计算周期信号的周期。对此，有时间将还会对此进行详细研究。

参考资料

[1] 《数字信号处理——理论、算法与实现》，胡广书
[2] 刘永春，陈琳. 基于广义互相关时延估计算法声源定位的研究.
[3] 金中薇,姜明顺等. 基于广义互相关时延估计算法的声发射定位技术. 传感技术学报. 第26卷11期，2013年11月.