相关运算及实现

1.定义

这里以长度为N的离散时间序列x(n),y(n)为例，相关运算定义如下：

1) $R_{xy}(n)$

x(n)保持不动，y(n)相对于x(n)向左移动

$R_{xy}(n)=\sum_{m=0}^{N-1}x(m)y(m+n)=\sum_{m=0}^{N-1}x(m-n)y(m)$

最后面的式子是等效情况下，y(n)不动，x(n)相对于y(n)向右移动

2) $R_{yx}(n)$

y(n)保持不动，x(n)相对于y(n)向左移动

$R_{yx}(n)=\sum_{m=0}^{N-1}y(m)x(m+n)=\sum_{m=0}^{N-1}y(m-n)x(m)$

最后面的式子是等效情况下，x(n)不动，y(n)相对于x(n)向右移动

注意：

1) $R_{xy}(n)$ ， $R_{yx}(n)$ 的定义是不一样的， $R_{xy}(n)=R_{yx}(-n)$

2)n的取值范围为[-(N-1),N-1]

3)相关运算结果长度为2*N-1（2个序列长度和减1）

对比卷积公式：

$z(n)=x(n)\ast y(n)=\sum_{m=0}^{N-1}x(m)y(n-m)=\sum_{m=0}^{N-1}y(m)x(n-m)$

注意：

1)n的取值范围为[0,2*N-2]

2)卷积运算结果长度为2*N-1（2个序列长度和减1）

3)后面2个等式成立使用的是卷积交换律

对比卷积公式和相关运算公式，可以知道，无论是 $R_{xy}(n)$ 还是 $R_{yx}(n)$ ，它们与卷积运算差别仅在运算时x(n)或y(n)是否需要折叠（转换成x(-n)或y(-n)）。因此，在做相关运算时可以将输入x(n)或y(n)先折叠，再经过卷积运算即可求出相关运算结果。

这里用Matlab作下验证，在Matlab命令行下输入：

x=[1,2,3];
y=[3,2,1];

z1=xcorr(x,y);
z2=conv(x, flip(y));

其中flip(y)即对y进行折叠（反转序列），得到的结果是z1和z2是相等的。

2.实现

这里基于C语言来实现。参考代码如下：

static void correlate(float *pSrcA, uint32_t srcALen, float *pSrcB, uint32_t srcBLen, float *pDst);


static void correlate(float *pSrcA, uint32_t srcALen, float *pSrcB, uint32_t srcBLen, float *pDst)
{
    float *pIn1 = pSrcA;                       /* inputA pointer */
    float *pIn2 = pSrcB + (srcBLen - 1U);      /* inputB pointer */
    float sum = 0;                             /* Accumulator */
    uint32_t i = 0U;                           /* loop counter */
    uint32_t j = 0U;                           /* loop counter */
    uint32_t inv = 0U;                         /* Reverse order flag */
    uint32_t tot = 0U;                         /* Length */

    if ((pSrcA == NULL) || (srcALen == 0) || (pSrcB == NULL) || (srcBLen == 0) || (pDst == NULL))
    {
        return ;
    }

    /* The algorithm implementation is based on the lengths of the inputs. */
    /* srcB is always made to slide across srcA. */
    /* So srcBLen is always considered as shorter or equal to srcALen */
    /* But CORR(x, y) is reverse of CORR(y, x) */
    /* So, when srcBLen > srcALen, output pointer is made to point to the end of the output buffer */
    /* and a varaible, inv is set to 1 */
    /* If lengths are not equal then zero pad has to be done to  make the two
     * inputs of same length. But to improve the performance, we assume zeroes
     * in the output instead of zero padding either of the the inputs*/
    /* If srcALen > srcBLen, (srcALen - srcBLen) zeroes has to included in the
     * starting of the output buffer */
    /* If srcALen < srcBLen, (srcALen - srcBLen) zeroes has to included in the
     * ending of the output buffer */
    /* Once the zero padding is done the remaining of the output is calcualted
     * using convolution but with the shorter signal time shifted. */

    /* Calculate the length of the remaining sequence */
    tot = ((srcALen + srcBLen) - 2U);

    if (srcALen > srcBLen)
    {
      /* Calculating the number of zeros to be padded to the output */
      j = srcALen - srcBLen;

      /* Initialise the pointer after zero padding */
      pDst += j;
    }
    else if (srcALen < srcBLen)
    {
      /* Initialization to inputB pointer */
      pIn1 = pSrcB;

      /* Initialization to the end of inputA pointer */
      pIn2 = pSrcA + (srcALen - 1U);

      /* Initialisation of the pointer after zero padding */
      pDst = pDst + tot;

      /* Swapping the lengths */
      j = srcALen;
      srcALen = srcBLen;
      srcBLen = j;

      /* Setting the reverse flag */
      inv = 1;
    }

    /* Loop to calculate convolution for output length number of times */
    for (i = 0U; i <= tot; i++)
    {
      /* Initialize sum with zero to carry on MAC operations */
      sum = 0.0f;

      /* Loop to perform MAC operations according to convolution equation */
      for (j = 0U; j <= i; j++)
      {
        /* Check the array limitations */
        if ((((i - j) < srcBLen) && (j < srcALen)))
        {
          /* z[i] += x[i-j] * y[j] */
          sum += pIn1[j] * pIn2[-((int32_t)i - (int32_t)j)];
        }
      }
      /* Store the output in the destination buffer */
      if (inv == 1)
      {
        *pDst-- = sum;
      }
      else
      {
        *pDst++ = sum;
      }
    }
}

main函数：

int main()
{
    float x[] = {1, 2, 3, 4};
    float y[] = {4, 3, 2, 1};
    float z[7] = {0};
    uint32_t i = 0;

    correlate(x, 4, y, 4, z);

    for (i = 0; i < sizeof(z) / sizeof(z[0]); i++)
    {
        printf("%f ", z[i]);
    }

    printf("\r\n");

    return 0;
}

输出结果：

作为对比，在Matlab命令行下输入：

x=[1,2,3,4];
y=[4,3,2,1];
z=xcorr(x,y);

查看结果和用C语言实现的是一致的。

3.应用

通过相关运算，我们得到了相关序列，为了方便后续数据处理，需要对得到的相关序列进行归一化处理，对于2个互相关序列，有如下2种归一化方法：

1)有偏估计

$\hat{R_{xy,biased}(m)}=\frac{1}{N}\hat{R_{xy}(m)}$

2)无偏估计

$\hat{R_{xy,unbiased}(m)}=\frac{1}{N-\left | m \right |}\hat{R_{xy}(m)}$

为了查看这2者之间差异，在Matlab命令行下，输入如下命令：

fm=1000;
fs=100000;
N=fs/fm;
k=0:1:1000;
theta=pi/8;
x=sin(2*k*pi/N + theta);%原始信号
xn=awgn(x,20,0);%原始信号加噪声
subplot(5,1,1);
plot(k,x);
title('x');
subplot(5,1,2);
plot(k,xn);
title('xn');
[r,lags]=xcorr(xn,x,500);%未缩放的互相关
subplot(5,1,3);
plot(k,r);
title('corr none');
[rb,lags]=xcorr(xn,x,500,'biased');%互相关的有偏估计
subplot(5,1,4);
plot(k,rb);
title('corr biased');
[rub,lags]=xcorr(xn,x,500,'unbiased');%互相关的无偏估计
subplot(5,1,5);
plot(k,rub);
title('corr unbiased');

输出结果：