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

November 10, 20103 comments Coded in 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);

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
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*****/

/**************************************************/``````