Reply by hyd198471 November 11, 20102010-11-11
>How much u r talking about :) > >hyd198471 wrote: > >> Hi everyone: >> I am implementing a fast linear correlation in C in DSK C6713. My cross >> correlation in principle should be implemented as following: >> The problem is I cannot get the correct cross correlation. Any help for
the
>> code will be appreciated. > >My dear inquisitive friend, > >I can answer all of your questions if you can give the right answer for >just one question of mine: "Do you have the money?" > > > >Vladimir Vassilevsky >DSP and Mixed Signal Design Consultant >http://www.abvolt.com > >
Reply by Vladimir Vassilevsky November 10, 20102010-11-10

hyd198471 wrote:

> Hi everyone: > I am implementing a fast linear correlation in C in DSK C6713. My cross > correlation in principle should be implemented as following: > The problem is I cannot get the correct cross correlation. Any help for the > code will be appreciated.
My dear inquisitive friend, I can answer all of your questions if you can give the right answer for just one question of mine: "Do you have the money?" Vladimir Vassilevsky DSP and Mixed Signal Design Consultant http://www.abvolt.com
Reply by hyd198471 November 10, 20102010-11-10
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
board.
- 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
2048
- matlab function to implement the linear correlation is ifft
(fft(x,N_FFT).*conj(fft(x1,N_FFT)))  
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
length*************************/
#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] = AIC23_data.channel[LEFT];
		in_bufR[buffercount] = AIC23_data.channel[RIGHT];
	}
	if (buffercount < N_FFT){ 
// the buffer length cross correlation is N_FFT,N_FFT = 2*N
			
		AIC23_data.channel[LEFT] = (short)CCFRE[buffercount];
		AIC23_data.channel[RIGHT] =(short)CCFRE[buffercount];
		
	}
	
	buffercount++;

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

	return;
}

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
	MCBSP_config(DSK6713_AIC23_DATAHANDLE,&AIC23CfgData);

//start data channel 
	MCBSP_start(DSK6713_AIC23_DATAHANDLE, MCBSP_XMIT_START | MCBSP_RCV_START
|
		MCBSP_SRGR_START | MCBSP_SRGR_FRAMESYNC, 220);

	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
in_buf_temp************/

		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
		cfftr2_dit(x1,W,N_FFT); 
		digitrev_index(iData,N_FFT,RADIX); 
// 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
called
		icfftr2_dif(temp,W,N_FFT);
		for (i =0;i<N_FFT;i++){
		
		CCFRE[i] = temp[i].re; // cross correlation
		}

	}

}