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 } } }
Looking for a fast linear cross correlation implemented in C in C6713
Started by ●November 10, 2010