matlab的xcorr函数转成c语言（时域相关？还是频域相关？）

一、MATLAB中相关函数xcorr

在matlab当中可以使用xcorr函数来求序列的自相关和互相关。

使用方法：

c = xcorr(x) 为矢量x的自相关估计。

c = xcorr(x,y)  返回矢量长度为2*N-1互相关函数序列，其中x和y的矢量长度均为N，如果x和y的长度不一样，则在短的序列后补零直到两者长度相等。

c = xcorr(x) 为矢量x的自相关估计。

c = xcorr(x,y,'option') 为有正规化选项的互相关计算；其中选项为"biased"为有偏的互相关函数估计；"unbiased"为无偏的互相关函数估计；"coeff"为0延时的正规化序列的自相关计算；"none"为原始的互相关计算。

在Matalb中，求解xcorr的过程事实上是利用Fourier变换中的卷积定理进行的。

参考：1.【总结】matlab求两个序列的相关性_whu_Paprika_新浪博客

           2.【 MATLAB 】xcorr 函数介绍（互相关）简介_李锐博恩的博客-CSDN博客_互相关函数matlab

二、c语言实现相关

1）时域相关：

如果根据定义，两个序列直接滑动相乘，求和，那么需要用两层for循环实现：

// x长度为 iDataN，y长度为 iSyncLength，假设 iDataN > iSyncLength
// corr长度为 iDataN
int xcorr(double *corr, double *x, double *y, int iDataN, int iSyncLength)
{
	
	double r =0.0;
	int i=0, j=0;
	/* 
	for (i = 0; i < iDataN; i++)
	{
		r=0;
		for(j=0; j < iSyncLength; j++)
		{
			if(i+j<iDataN)
			{	
				r+=x[i+j]*y[j];
			}
			//else
			//{
			//	r+=0;
			//}
			
		}
		corr[i]=r;
	}*/

	for (i = 0; i < iDataN- iSyncLength+1; i++)
	{
		r=0;
		for(j=0; j < iSyncLength; j++)
		{
			r+=x[i+j]*y[j];
			
		}
		corr[i]=r;
	}

	for (i = iDataN- iSyncLength+1; i < iDataN; i++)
	{
		r=0;
		for(j=0; j <iDataN- i; j++)
		{		
			r+=x[i+j]*y[j];
			
		}
		corr[i]=r;
	}

	//cout << corr[1] << endl;
	return 0;

}

复杂度为，非常耗时，意味着实时性能会出现问题。

2）频域相关：

我们想到：

两个能量信号的互相关函数的傅里叶变换是他们的互能量谱密度。

两个功率信号的互相关函数的傅里叶变换是他们的互功率谱密度。

假设两序列分别为

他们的傅里叶变换： ， 

令

则

于是求互相关函数可以转换为先分别求两个序列的傅里叶变换 F1(f) 和 F2(f)，再用F1(f)乘以F2(f)的共轭得到S(f)，最后求S(f)的傅里叶反变换得到R12(τ)，代码实现如下：

#include "fft.h"
#include <math.h>
#include <string.h>
#include <stdint.h>
#include <iostream>
using namespace std;

#define pi (double) 3.14159265359


uint32_t fft_getrealsize(uint32_t size) {
	uint32_t m =0;
	while ((size /= 2) != 0)
		m++;

	return 1 << m;
}

// out should have a size of 2*in_size
void real_to_complex(double * out, double * in, int in_size) {
	uint32_t i;

	// set the real ids in answer to the val, the imaginary ones to 0
	for (i = 0; i < in_size; i++) {
		*(out++) = *(in++);
		*(out++) = 0.0f;
	}
}

void complex_to_real(double * data, int samples){
	double * src = data;
	uint32_t i;
	for (i = 0; i < samples; i++) {
		const double I = *(src++);
		const double Q = *(src++);
		*(data++) = sqrtf(I*I+Q*Q);
	}
}

//data奇数位为幅值abs，偶数位为0
void fft_complex_to_absolute_complex(double * data, int samples) {
	uint32_t fft_size2 = samples * 2;

	uint32_t i;
	for (i = 0; i < fft_size2; i+=2) {
		const int i1 = i+1;
		const double I = data[i];
		const double Q = data[i1];
		data[i] = sqrtf(I*I+Q*Q);
		data[i1] = 0;
	}
}

// the array temp needs to be of size at least 2*size
// the answer will be written in answer at entries fro 0 to 2*size
void fft_autocorrelation(double * answer, double * real, uint32_t size) {

	// set the real ids in answer to the val, the imaginary ones to 0
	//实部为real，虚部置0 
	real_to_complex(answer, real, size);

	// do the fft
	uint32_t fft_size = fft_getrealsize(size);//将size变为2的n次幂，并且fft_size<size

	fft_perform(answer, fft_size, 0);

	// get the abs value
	// answer长度为 size*2
	// answer奇数位为abs(I*I+Q*Q)，偶数位为0
	fft_complex_to_absolute_complex(answer, size);

	// get the ifft
	fft_perform(answer, fft_size, 1);
}

// a and b need to be complex with size samples * 2
// the final answer can be found in answer_out and it is complex
// you may need to take its absolute value to get the crosscorrelation
void fft_crosscorrelation(double * answer_out, double * answer_temp, uint32_t samples) {
	uint32_t i;

	uint32_t fft_size = fft_getrealsize(samples);
	uint32_t fft_size2 = fft_size * 2;

	// perform the ffts
	fft_perform(answer_out, fft_size, 0);
	fft_perform(answer_temp, fft_size, 0);
	
	//cout<<fft_size<<endl;
	/*for(i=0;i< 50; i++){	
		cout<<answer_temp[i]<<endl;
		
	}*/

	// multiply the cojugate乘以共轭
	for (i = 0; i < fft_size2; i+=2) {
		const int i1 = i+1;
		const double aI = answer_out[i];
		const double aQ = answer_out[i1];
		const double bI = answer_temp[i];
		const double bQ = answer_temp[i1];

		answer_out[i]  = aI*bI + aQ*bQ;
		answer_out[i1] = aQ*bI - aI*bQ;//共轭
	}

	// get the ifft
	fft_perform(answer_out, fft_size, 1);
}

// iq must have a size of size*2
void fft_perform(double * iq, uint32_t size, int inverse)
{
	int64_t i,i1,j,k,i2,l,l1,l2;
	double c1,c2,tx,ty,t1,t2,u1,u2,z;

	int m =0;//size>2^m
	while ((size /= 2) != 0)
		m++;

	uint32_t nn = 1 << m;//nn=2^m
	i2 = nn >> 1;//i2=2^(n-1)
	j = 0;

	for (i=0;i<nn-1;i++) {
		if (i < j) {
			const uint32_t Ii = i << 1;
			const uint32_t Qi = Ii + 1;

			const uint32_t Ij = j << 1;
			const uint32_t Qj = Ij + 1;

			tx = iq[Ii];
			ty = iq[Qi];
			iq[Ii] = iq[Ij];
			iq[Qi] = iq[Qj];
			iq[Ij] = tx;
			iq[Qj] = ty;
		}
		k = i2;
		while (k <= j) {
			j -= k;
			k >>= 1;
		}
		j += k;
	}

	c1 = -1.0;
	c2 = 0.0;
	l2 = 1;
	for (l=0;l<m;l++) {
		l1 = l2;
		l2 <<= 1;
		u1 = 1.0f;
		u2 = 0.0f;
		for (j=0; j<l1; j++) {
			for (i=j; i<nn; i+=l2) {
				const uint32_t Ii = i << 1;
				const uint32_t Qi = Ii + 1;

				i1 = i + l1;

				const uint32_t Ii1 = i1 << 1;
				const uint32_t Qi1 = Ii1 + 1;

				t1 = u1 * iq[Ii1] - u2 * iq[Qi1];
				t2 = u1 * iq[Qi1] + u2 * iq[Ii1];
				iq[Ii1] = iq[Ii] - t1;
				iq[Qi1] = iq[Qi] - t2;
				iq[Ii] += t1;
				iq[Qi] += t2;
			}
			z =  u1 * c1 - u2 * c2;
			u2 = u1 * c2 + u2 * c1;
			u1 = z;
		}
		c2 = sqrt((1.0 - c1) / 2.0);
		if (!inverse)
			c2 = -c2;
		c1 = sqrt((1.0 + c1) / 2.0);
	}

	if (!inverse) {
		for (i=0;i<nn;i++) {
			const uint32_t Ii = i << 1;
			const uint32_t Qi = Ii + 1;

			iq[Ii] /= (double)nn;
			iq[Qi] /= (double)nn;
		}
	}
}

借鉴https://github.com/martinmarinov/TempestSDR中自相关的思想，源码来自其中的一部分

复杂度降低至。

MATLAB与c语言互相关输出结果对比：

close all;
a=load('answer_out.txt');
a=a(1:2:end);
b=load('answer_temp.txt');
b=b(1:2:end);
len_a=length(a);
xc=xcorr(a,b);

figure;
subplot 211
plot(abs(xc(len_a:end)));
title('matlab xcorr输出的互相关结果')

xc2=load('fft_out.txt');
I=xc2(1:2:end);
Q=xc2(2:2:end);
xc2cpl=I+Q*1i;

subplot 212
plot(abs(xc2cpl));
title('c语言输出的互相关结果')

测试数据：matlab的xcorr函数转成c语言测试数据-C++文档类资源-CSDN下载