DSP C PROGRAM POWER SPECTRAL DENSITY - Computer Programming

# Computer Programming

C C++ Java Python Perl Programs Examples with Output -useful for Schools & College Students

## Thursday, November 10, 2011

.   POWER SPECTRAL DENSITY

AIM: -
To design and implement a low pass FIR filter using windowing technique.
APPARATUS: -
1.     TMS320C6713 DSK.
2.     USB Cable.
3.     Power Cord
4 .    5V Adapter.

PROGRAM: -

#include <math.h>
#define PTS 128    //# of points for FFT
#define PI 3.14159265358979
typedef struct {float real,imag;} COMPLEX;
void FFT(COMPLEX *Y, int n);     //FFT prototype
float iobuffer[PTS];                          //as input and output buffer
float x1[PTS],x[PTS];                                   //intermediate buffer
short i;                                   //general purpose index variable
short buffercount = 0;           //number of new samples in iobuffer
short flag = 0;                                    //set to 1 by ISR when iobuffer full
float y[128];
COMPLEX w[PTS];                        //twiddle constants stored in w
COMPLEX samples[PTS];                           //primary working buffer

main()
{

float j,sum=0.0 ;
int n,k,i,a;
for (i = 0 ; i<PTS ; i++)        // set up twiddle constants in w
{
w[i].real = cos(2*PI*i/(PTS*2.0)); //Re component of twiddle constants
w[i].imag =-sin(2*PI*i/(PTS*2.0));  /*Im component of twiddle constants*/
}

/****************Input Signal X(n) ********************************************************/

for(i=0,j=0;i<PTS;i++)
{ x[i] = sin(2*PI*5*i/PTS);      // Signal x(Fs)=sin(2*pi*f*i/Fs);
samples[i].real=0.0;
samples[i].imag=0.0;
}

/********************Auto Correlation of X(n)=R(t) *********************************/

for(n=0;n<PTS;n++)
{
sum=0;
for(k=0;k<PTS-n;k++)
{
sum=sum+(x[k]*x[n+k]);     // Auto Correlation R(t)
}
iobuffer[n] = sum;
}

/********************** FFT of R(t) *****************************/

for (i = 0 ; i < PTS ; i++)    //swap buffers
{
samples[i].real=iobuffer[i]; //buffer with new data
}

for (i = 0 ; i < PTS ; i++)
samples[i].imag = 0.0;                //imag components = 0

FFT(samples,PTS);              //call function FFT.c

/******************** PSD *******************************************/

for (i = 0 ; i < PTS ; i++)    //compute magnitude
{
x1[i] = sqrt(samples[i].real*samples[i].real
+ samples[i].imag*samples[i].imag);
}

}                                                        //end of main

void FFT(COMPLEX *Y, int N)      //input sample array, # of points
{
COMPLEX temp1,temp2;             //temporary storage variables
int i,j,k;                       //loop counter variables
int upper_leg, lower_leg;                 //index of upper/lower butterfly leg
int leg_diff;                    //difference between upper/lower leg
int num_stages = 0;              //number of FFT stages (iterations)
int index, step;                 //index/step through twiddle constant
i = 1;                           //log(base2) of N points= # of stages
do
{
num_stages +=1;
i = i*2;
}while (i!=N);
leg_diff = N/2;                                 //difference between upper&lower legs
step = (PTS*2)/N;                                        //step between values in twiddle.h   // 512
for (i = 0;i < num_stages; i++)  //for N-point FFT
{
index = 0;
for (j = 0; j < leg_diff; j++)
{
for (upper_leg = j; upper_leg < N; upper_leg += (2*leg_diff))
{
lower_leg = upper_leg+leg_diff;
temp1.real = (Y[upper_leg]).real + (Y[lower_leg]).real;
temp1.imag = (Y[upper_leg]).imag + (Y[lower_leg]).imag;
temp2.real = (Y[upper_leg]).real - (Y[lower_leg]).real;
temp2.imag = (Y[upper_leg]).imag - (Y[lower_leg]).imag;
(Y[lower_leg]).real = temp2.real*(w[index]).real
-temp2.imag*(w[index]).imag;
(Y[lower_leg]).imag = temp2.real*(w[index]).imag
+temp2.imag*(w[index]).real;
(Y[upper_leg]).real = temp1.real;
(Y[upper_leg]).imag = temp1.imag;
}
index += step;
}
leg_diff = leg_diff/2;
step *= 2;
}
j = 0;
for (i = 1; i < (N-1); i++)     //bit reversal for resequencing data
{
k = N/2;
while (k <= j)
{
j = j - k;
k = k/2;
}
j = j + k;
if (i<j)
{
temp1.real = (Y[j]).real;
temp1.imag = (Y[j]).imag;
(Y[j]).real = (Y[i]).real;
(Y[j]).imag = (Y[i]).imag;
(Y[i]).real = temp1.real;
(Y[i]).imag = temp1.imag;
}
}
return;
}