How to Perform Discrete Fourier Transform on N-Point Sequence using STM32L476G

Posted

in

by

x(n) = { 0, 1, 0, 1}

Here x(n) is a 4 point sequence

Now to perform Fourier transform on this sequence.

DFT Formula

Here X(k) is the DFT of x(n)

‘k’ is the index representing individual frequency component

‘N’ is the Length of the sample sequence

‘n’ is an index of the element of the sequence in x(n)

For example: if we calculate the value for k =1

This will give us the X(1)

So, if the sequence is sampled at 4Hz and the sampling frequency is equal to two times the max frequency. Then X(1) gives us the 1 Hz frequency component.

To calculate Euler identity

Here is the code:

/*
* Author: ABHAY KANT
* Discrete Fourier Transform using STM32L476vg on STM32L476g-Disco board.
* computes N point DFT
* computes magnitude of N/2 Points of dft
* computes phase of N/2 points of dft
* NOTE: DFT array is named xfft. It DOES NOT Calcultes DFT using FFT algo.
* Limitations: for a 2000 point sequence it takes a lot of time.
               And as the number of points increase the computational time also increases.
a rough estimate of time for a 2000 point dft it takes roughly ~ 60 seconds. 
*/
asm("nop");
//int xn[6] = {1,1,2,2,3,3};

int16_t Len = 1024; 
float xn[Len];

float xfft[Len * 2];

/*
 * calculating phase and mag for only N/2 since the value generated after N/2 are conjugate 
 */
float phase[Len/2], mag[Len/2];
int	k = 1;
int	N = 0;
int	n = 0;
float	real = 0;
float	imag = 0;
float	param = 0;
int	freq = 10;
int i;

/*
 * These two "for loops" are used to initalize the array with zero's
 */
for(n=0;n<=Len;n++){xn[n] = 0;}
for(n=0;n<=Len*2;n++){xfft[n] = 0;}



/*
 * Generating a sine Wave 
 */

	for(i=0; i<Len/2; i++){
		xn[i+512] = sin((2 * 3.1415926 * i * freq)/Len);
	}

N = Len;				
/* Note: if Len/8 then it will only calculate the 1/8 th term from 0 only.*/
/*
 * This loop calculates the x[n].e^(-j2.PI.n.k/N)
 */

	for(k=0 ; k<N ; k++){

		for(n=0 ; n<N; n++){
		param = ((2 * M_PI * n * k)/N);
		real = real + ( xn[n] * cos(param) );
		imag = imag + ((-1) * xn[n] * sin(param) );
		}
		xfft[2 * k] = real;
		xfft[(2 * k) + 1] = imag;
		real = 0;
		imag = 0;
	}



/*
 * Calculate phase angle of each frequency component
 */
float temp;
for(k=0; k<N/2;k++)
{
	temp = xfft[(2 * k ) + 1] / xfft[2 * k];
	phase[k] = atan(temp);
}

/*
 * Calculate Magnitude of each frequency 
 */

for(k=0; k<N/2;k++)
{
	temp = sqrt( (xfft[(2 * k ) + 1] * xfft[(2 * k ) + 1] ) + ( xfft[2 * k] * xfft[2 * k] ) );
	mag[k] = temp;
}
asm("nop");  //added to support for creating breakpoint during debug

Comments

Leave a Reply

Your email address will not be published. Required fields are marked *