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