FFT with down-sampling in time domain for C++

This code is a complete example of the implementation and usage of a function which calculates the Fast Fourier Transform based on the downsampling in time algorithm.

It considers the details of the treatment of vector signals in C++ language.

The development of a function capable of inverting the bits of an integer is required to run the code.

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define N 128
#define PI 3.14159265358979

struct COMPLEX
{float real,imag;};

void qfft(struct COMPLEX *X, int n);   //FFT function
int q_log2(int n);
int qflip(int i, int M);     //self developed function like fliplr() in Matlab for bits.

struct COMPLEX w[N/2];       	    //twiddle factors
struct COMPLEX samples[N];          //input signal
float x1[N];                        //output signal
int i;
float s;		           

main(void)
{	
	for (i = 0 ; i<N/2 ; i++)	          // set up twiddle constants in w
    {
      	   w[i].real = cos(2*PI*i/N);  //Re component of twiddle constants
      	   w[i].imag =-sin(2*PI*i/N);  //Im component of twiddle constants
    }
	
/*Generating test signal******************/	
    for (i=0; i<N; i++)
	{
        s=0+i;
		samples[i].real = -1+(s*(2/128));
		samples[i].imag = 0; 		
	}
	printf("%f,%f",samples[50].real,s);
	
/************************************************/	
	system("PAUSE");
	//fft calculation
	qfft(samples,N);
        //magnitude calculation
	x1[i] = sqrt(samples[i].real*samples[i].real
	         + samples[i].imag*samples[i].imag);
	
    printf("El resultado es:\n");
	for (i=0; i<N; i++)
	{
		printf("%f\t%fi\n",samples[i].real,samples[i].imag);
	}
	
    system("PAUSE");
}

void qfft(struct COMPLEX *X, int n) 
{
	int k,l,m,k2,k21,Mk,index;
	struct COMPLEX op1,op2;
	struct COMPLEX X1;
	//nuber of stages
	int M=q_log2(n);
	//to array elements of the input signal
	for (k=0; k<N/2; k++)
	{
        index=qflip(k,M);
        X1=X[k];
        X[k]=X[index];
        X[index]=X1;
        printf("%d\n",index);
        }
	//entering to graph
	//k holds the stage
	for (k=1; k<=M; k++)
	{
		k2=(int)pow(2,k);
		k21=(int)pow(2,k-1);
		Mk=(int)pow(2,M-k);
		//l holds the butterfly
		for (l=1; l <= (n/k2); l++)
		{
			//m holds the upper element of the butterfly
			for (m=1; m<=k21; m++)
			{
				op1=X[(l-1)*k2+m-1];
				op2=X[(l-1)*k2+m+k21-1];
				//Butterfly equations
				//Complex product
				//(a+jb)(c+jd)=ac+j(ad+bc)-bd
				X[(l-1)*k2+m-1].real = op1.real + w[Mk*(m-1)].real*op2.real 
									 - w[Mk*(m-1)].imag*op2.imag ;
				X[(l-1)*k2+m-1].imag = op1.imag + w[Mk*(m-1)].real*op2.imag
									 + w[Mk*(m-1)].imag*op2.real ;
				//su elemento complemento:
				X[(l-1)*k2+m+k21-1].real = op1.real - w[Mk*(m-1)].real*op2.real 
									     + w[Mk*(m-1)].imag*op2.imag ;
				X[(l-1)*k2+m+k21-1].imag = op1.imag - w[Mk*(m-1)].real*op2.imag
									     - w[Mk*(m-1)].imag*op2.real ;
			}//end m
		}//end l
	}//end k
	return;
}//end function

int q_log2(int n)
{
	int res=0;
	while(n>1)
	{
		n=n/2;
		res=res+1;
	}
	return res;
}

/***Implementation of the qflip function*****/

/**************************************************/