Looking for a fast linear cross correlation implemented in C in C6713

Started by hyd198471 November 10, 2010
Hi everyone:
I am implementing a fast linear correlation in C in DSK C6713. My cross
correlation in principle should be implemented as following:
- two random digital signals from AIC23 codec are the input signals of DSK
- interrupt intser_McBSP1 is used to input two random signals via two
channels(left and right) and output the real part of cross correlation r
- Because of the linear fast correlation, FFTs are used to transform the
two random signals to frequency domain. 
-Be aware of the length of FFT. My input signals x and x1 have the buffer
length N,N = 1024. Due to linear correlation, zero padding is necessary, so
the FFT length is 2*N = 2048 and the length of cross correlation is also
- matlab function to implement the linear correlation is ifft
The corresponding C code to implement the cross correlation is following:
The problem is I cannot get the correct cross correlation. Any help for the
code will be appreciated.

#include "dsk6713_aic23.h"
#include "dsk6713.h"
#include "C6713dskinit.h"
#include <math.h>

/*****************left,right channel and buffer
#define LEFT 1
#define RIGHT 0
#define PI 3.1415926535897
#define N 1024 //length of input samples 
#define N_FFT 2048 //length of FFT
#define RADIX 2 //radix or base  
#define DELTA (2*PI)/N //argument for sine/cosine 

extern MCBSP_Handle DSK6713_AIC23_DATAHANDLE;
static Uint32 CODECEventId;
Uint32 fs=DSK6713_AIC23_FREQ_8KHZ;     //for sampling frequency

union {
	Uint32 both; 
	short channel[2];
      } AIC23_data;  

short i = 0;								
short iTwid[N_FFT/2];		//index for twiddle constants W 
short iData[N_FFT];		//index for bitrev X 
float Xmag[N_FFT];		//magnitude spectrum of x
float Xmag1[N_FFT]; 		// magnitude spectrum of x1
float CCFRe[N_FFT];		// Real part of cross correlation

typedef struct Complex_tag {float re,im;}Complex;  
Complex W[N_FFT/RADIX];	        //array for twiddle constants
Complex x[N_FFT]; 		//N_FFT complex data values 
Complex x1[N_FFT];		//N_FFT complex data values
Complex temp[N_FFT] ;           // N_FFT complex data values for complex 
//multiplication of x and x1

short buffercount = 0;             
short in_bufL[N];
// input buffer contains the input signal from ADC in channel left
short in_bufR[N];
// input buffer contains the input signal from ADC in channel right
/***************** zero padding variables*****/

short in_bufL_temp[N_FFT];
short in_bufR_temp[N_FFT];

#pragma DATA_ALIGN(W,sizeof(Complex))   //align W on boundary 
#pragma DATA_ALIGN(x,sizeof(Complex))	//align input x on boundary

void cfftr2_dit( Complex* x, Complex* w, short n);
void digitrev_index(short *index, int n, int radix);
void bitrev(Complex *x, short *index, int );
void icfftr2_dit(Complex* x, Complex* w, short n);

interrupt void intser_McBSP1()
	AIC23_data.both = MCBSP_read(DSK6713_AIC23_DATAHANDLE); 
// input data
	if (buffercount < N) 
// 2 channel input store in the in_bufL and in_bufR 
		in_bufL[buffercount] =[LEFT];
		in_bufR[buffercount] =[RIGHT];
	if (buffercount < N_FFT){ 
// the buffer length cross correlation is N_FFT,N_FFT = 2*N[LEFT] = (short)CCFRE[buffercount];[RIGHT] =(short)CCFRE[buffercount];

	MCBSP_write(DSK6713_AIC23_DATAHANDLE, AIC23_data.both);   
//output 32 bit


void main() 
 for( i = 0 ; i < N_FFT/RADIX ; i++ ) 
   W[i].re = cos(DELTA*i);//real component of W
   W[i].im = sin(DELTA*i);//neg imag component

  for(i= 0;i<N_FFT;i++)	// in_buf_temp has the length N_FFT (2048)
	in_bufL_temp[i] = 0;
	in_bufR_temp[i] = 0;
 	digitrev_index(iTwid, N_FFT/RADIX, RADIX);
//produces index for bitrev() W 
 	bitrev(W, iTwid, N_FFT/RADIX);		       	
//bit reverse W 	
	IRQ_globalDisable();     //disable interrupts
	DSK6713_init();          //call BSL to init DSK-EMIF,PLL)
// handle(pointer) to codec
	hAIC23_handle=DSK6713_AIC23_openCodec(0, &config);
	DSK6713_AIC23_setFreq(hAIC23_handle, fs);  //set sample rate
//interface 32 bits toAIC23

//start data channel 

	CODECEventId= MCBSP_getXmtEventId(DSK6713_AIC23_DATAHANDLE);//McBSP1 Xmit
	IRQ_map(CODECEventId, 5);	//map McBSP1 Xmit to INT5
	IRQ_reset(CODECEventId);    	//reset codec INT5
	IRQ_globalEnable();       	//globally enable interrupts
	IRQ_nmiEnable();          	//enable NMI interrupt
	IRQ_enable(CODECEventId);	//enable CODEC eventXmit INT5
	IRQ_set(CODECEventId);        //manually start the first interrupt
	// initialization
	for(i=0; i<N_FFT; i++)
 		Xmag[i] = 0;
		Xmag1[i] = 0;
		CCFRE[i] = 0;
		temp[i].re = 0;
		temp[i].im = 0;
		x[i].re = 0;
		x[i].im = 0;
		x1[i].re = 0;
		x1[i].im = 0;
	while(1)	//infinite loop
		while(buffercount < N_FFT); // wait until buffer is full
		buffercount = 0;
		/*************** zero padding copy N samples from in_buf to

		for (i=0;i<N;i++)
			in_bufL_temp[i] = in_bufL[i];
			in_bufR_temp[i] = in_bufR[i];

		for (i=0;i<N_FFT;i++)
			x[i].re = (float)in_bufL_temp[i];
			x1[i].re= (float)in_bufR_temp[i];
			x[i].im = 0.0;
			x1[i].im = 0.0;
		cfftr2_dit(x,W,N_FFT); // TI-floating-pt complex FFT
// produce the index for bitrev
		bitrev(x,iData,N_FFT); // bit-reverse x
		bitrev(x1,iData,N_FFT); // bit-reverse x1

		//calculate complex conjugate of x1

		for (i =0;i<N_FFT;i++)
		  x1[i].im = -1*x1[i].im;
		//complex multiplication of spectrum of x and spectrum of complex
conjugate of x1
		for (i = 0;i<N_FFT;i++)
		temp[i].re = x[i].re * x1[i].re - x[i].im * x1[i].im;
		temp[i].im = x[i].re * x1[i].im + x[i].im * x1[i].re;

		// IFFT,the asm function iffftr2_dif(float* x, float* w, short n) is
		for (i =0;i<N_FFT;i++){
		CCFRE[i] = temp[i].re; // cross correlation