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
}
}
}