fourier transform iFFT radix-4 for N=16 algorithm c++ source code

//this is new in that method:
//when you want to have equal results that are in false modificator in normal FFT then change this:
/*
 fun_fourier_transform_FFT_radix_4_N_256_official
{
    for(int j=0;j<N;j++)
    {
      tab[j].real() =tab[j].real()*2/N;
      tab[j].imag() =tab[j].imag()*2/N;
    }
}
//and:

fun_inverse_fourier_transform_FFT_radix_4_N_256_official
{
    for(int j=0;j<N;j++)
    {
      tab[j].real() =tab[j].real()*0.5;
      tab[j].imag() =tab[j].imag()*0.5;
    }
}

//for official modificator that is only in inverse FFT
*/


//fourier transform iFFT radix-4 for N=16 algorithm c++ source code
#include <iostream>
#include "conio.h"
#include <stdlib.h>
#include <math.h>
#include <cmath>
#include <time.h>
#include <complex>
#include <fstream>

using namespace std;

//complex number method:
void fun_inverse_bits_radix_4(int N,std::complex<double> tab[]);
void fun_fourier_transform_FFT_radix_4_N_16_official(int N,std::complex<double> tab[]);
void fun_inverse_fourier_transform_FFT_radix_4_N_16_official(int N,std::complex<double> tab[]);


void fun_fourier_transform_FFT_radix_4_N_16_method2(int N,std::complex<double> tab[]);
void fun_inverse_fourier_transform_FFT_radix_4_N_16_method2(int N,std::complex<double> tab[]);
void fun_fourier_transform_FFT_radix_4_N_16_method3(int N,std::complex<double> tab[]);
void fun_fourier_transform_FFT_radix_4_N_16_method4(int N,std::complex<double> tab[]);

static double diffclock(clock_t clock2,clock_t clock1)
{
    double diffticks=clock1-clock2;
    double diffms=(diffticks)/(CLOCKS_PER_SEC/1000);
    return diffms;
}
int main()
{
 int N;
    //if N==period of signal in table tab[] then resolution = 1 Hz
    N=16;
    std::complex<double> tab2[16]={{1,17},    {2,18}    ,{3,19},    {4,20},
    {5,21},    {6,22}    ,{7,23},    {8,24},
    {9,25},    {10,26}    ,{11,27},    {12,28},
    {13,29},{14,30}    ,{15,31},    {16,32}};


    //std::complex<double> tab2[4]={{0.707106781},    {0.292893219}    ,{-0.707106781},    {-0.292893219}};
    //std::complex<double> tab2[4]={{0,1},    {0,2}    ,{0,3},    {0,4}};
    //std::complex<double> tab2[4]={{0,0},    {0,0}    ,{0,0},    {0,3}};
    //std::complex<double> tab2[4]={{7,0},    {0,0}    ,{0,0},    {0,0}};
    //std::complex<double> tab2[4]={{1},    {2}    ,{3},    {4}};

    double time2;
    double zmienna=0;
    /*
    std::fstream plik;
    plik.open("test.txt", std::ios::in | std::ios::out);
    if( plik.good() == false )
    {
        cout<<"nie otwarto pliku"<<endl;
        system("pause");
    }
*/
    cout<<"signal="<<endl;
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab2[j].real()*1000)/1000<<"  ";
    }
    cout<<endl;
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab2[j].imag()*1000)/1000<<"  ";
    }
    cout<<endl;
    cout<<endl;

    clock_t start = clock();

    //fun_fourier_transform_FFT_radix_4_N_16_method2(N,tab2);
    //fun_fourier_transform_FFT_radix_4_N_16_method3(N,tab2);
    //fun_fourier_transform_FFT_radix_4_N_16_method4(N,tab2);

    //////////////////////////////////////////////////////////
    fun_fourier_transform_FFT_radix_4_N_16_official(N,tab2);
    fun_inverse_bits_radix_4(N,tab2);
    ///////////////////////////////////////////////////////////

    time2=diffclock( start, clock() );

    cout<<"frequency Hz"<<endl;
    zmienna=0;
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
   cout<<round(tab2[j].real()*1000)/1000<<"  ";
    zmienna=zmienna+fabs(round(tab2[j].real()*1000)/1000);
    }
    cout<<endl;
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab2[j].imag()*1000)/1000<<"  ";
    zmienna=zmienna+fabs(round(tab2[j].imag()*1000)/1000);
    }
    cout<<endl;


 // plik.close();
    system("pause");
    cout<<endl;
    //fun_inverse_fourier_transform_FFT_radix_4_N_16_method2(N,tab2);

    /////////////////////////////////////////////////////////////////
    fun_inverse_fourier_transform_FFT_radix_4_N_16_official(N,tab2);
    fun_inverse_bits_radix_4(N,tab2);
    ////////////////////////////////////////////////////////////////

    cout<<"inverse/signal="<<endl;
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab2[j].real()*1000)/1000<<"  ";
    }
    cout<<endl;
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab2[j].imag()*1000)/1000<<"  ";
    }
    cout<<endl;
    cout<<endl;
    cout<<endl;

    //std::complex<double> tab3[4]={{0,5},    {1,11}    ,{2,18},    {3,6}};

    //std::complex<double> tab3[16]={{0,5},    {1,11}    ,{2,18},    {3,6},    {4,6},    {5,6},    {6,6},{7,6}
   // ,{8,6},{9,6},{10,6},{11,6},{12,6},{13,6},{14,6},{15,6}};


     std::complex<double> tab3[64]={{0,5},    {1,11}    ,{2,18},    {3,6},    {4,6},    {5,6},    {6,6},{7,6}
    ,{8,6},{9,6},{10,6},{11,6},{12,6},{13,6},{14,6},{15,6}
    ,{16,6},{17,6},{18,6},{19,6},{20,6},{21,6},{22,6},{23,6}
    ,{24,6},{25,6},{26,6},{27,6},{28,6},{29,6},{30,6},{31,6}
    ,{32,6},{33,6},{34,6},{35,6},{36,6},{37,6},{38,6},{39,6}
    ,{40,6},{41,6},{42,6},{43,6},{44,6},{45,6},{46,6},{47,6}
    ,{48,6},{49,6},{50,6},{51,6},{52,6},{53,6},{54,6},{55,6}
    ,{56,6},{57,6},{58,6},{59,6},{60,6},{61,6},{62,6},{63,6}};


    fun_inverse_bits_radix_4(64,tab3);
    cout<<"fun_inverse_bits_radix_4="<<endl;
    for(int j=0;j<64;j++)
    {
    cout.precision(4);
    cout<<round(tab3[j].real()*1000)/1000<<"  ";
    }
    cout<<endl;
    system("pause");
    return 0;
}


void fun_fourier_transform_FFT_radix_4_N_16_method2(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double> tab2[64]={};    // tab2[]==N
    std::complex<double>  w1[1]={{1,1}};
    std::complex<double>  w2[1]={{1,1}};
    std::complex<double>  w3[1]={{1,1}};
/*
for(int j=0;j<N;j++)
{
   cout<< tab[j]<<"  ";
}
cout<<endl<<"........."<<endl;
*/
    //radix-4
    for(int k=0;k<4;k++)
    {
        for(int n=0;n<4;n++)
        {
          w2[0].real()=cos(n*k*2*pi/(N/4));
          w2[0].imag()=-sin(n*k*2*pi/(N/4));

          tab2[4*k+0]=tab2[4*k+0]+(tab[n]+tab[n+4]+tab[n+8]+tab[n+12])*w2[0];

          w1[0].real()=cos(n*2*pi/N);
          w1[0].imag()=-sin(n*2*pi/N);
          w3[0].real()=0;
          w3[0].imag()=-1;
          tab2[4*k+1]=tab2[4*k+1]+(tab[n]-tab[n+8]+w3[0]*(tab[n+4]-tab[n+12]))*w1[0]*w2[0];

          w1[0].real()=cos(2*n*2*pi/N);
          w1[0].imag()=-sin(2*n*2*pi/N);
          tab2[4*k+2]=tab2[4*k+2]+(tab[n]-tab[n+4]+tab[n+8]-tab[n+12])*w1[0]*w2[0];

          w1[0].real()=cos(3*n*2*pi/N);
          w1[0].imag()=-sin(3*n*2*pi/N);
          w3[0].real()=0;
          w3[0].imag()=1;
          tab2[4*k+3]=tab2[4*k+3]+(tab[n]-tab[n+8]+w3[0]*(tab[n+4]-tab[n+12]))*w1[0]*w2[0];
        }
    }

    for(int j=0;j<N;j++)
    {
      tab[j].real() =tab2[j].real()*2/N;
      tab[j].imag() =tab2[j].imag()*2/N;
    }
}

void fun_inverse_fourier_transform_FFT_radix_4_N_16_method2(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double> tab2[64]={};    // tab2[]==N
    std::complex<double>  w1[1]={{1,1}};
    std::complex<double>  w2[1]={{1,1}};
    std::complex<double>  w3[1]={{1,1}};
/*
for(int j=0;j<N;j++)
{
   cout<< tab[j]<<"  ";
}
cout<<endl<<"........."<<endl;
*/
          //radix-4
    for(int k=0;k<4;k++)
    {
        for(int n=0;n<4;n++)
        {
          w2[0].real()=cos(n*k*2*pi/(N/4));
          w2[0].imag()=sin(n*k*2*pi/(N/4));

          tab2[4*k+0]=tab2[4*k+0]+(tab[n]+tab[n+4]+tab[n+8]+tab[n+12])*w2[0];

          w1[0].real()=cos(n*2*pi/N);
          w1[0].imag()=sin(n*2*pi/N);
          w3[0].real()=0;
          w3[0].imag()=1;
          tab2[4*k+1]=tab2[4*k+1]+(tab[n]-tab[n+8]+w3[0]*(tab[n+4]-tab[n+12]))*w1[0]*w2[0];

          w1[0].real()=cos(2*n*2*pi/N);
          w1[0].imag()=sin(2*n*2*pi/N);
          tab2[4*k+2]=tab2[4*k+2]+(tab[n]-tab[n+4]+tab[n+8]-tab[n+12])*w1[0]*w2[0];

          w1[0].real()=cos(3*n*2*pi/N);
          w1[0].imag()=sin(3*n*2*pi/N);
          w3[0].real()=0;
          w3[0].imag()=-1;
          tab2[4*k+3]=tab2[4*k+3]+(tab[n]-tab[n+8]+w3[0]*(tab[n+4]-tab[n+12]))*w1[0]*w2[0];
        }
    }

    for(int j=0;j<N;j++)
    {
      tab[j].real() =tab2[j].real()*0.5;
      tab[j].imag() =tab2[j].imag()*0.5;
    }
}


void fun_inverse_bits_radix_4(int N,std::complex<double> tab[])
{

}

void fun_fourier_transform_FFT_radix_4_N_16_method3(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double> tab2[64]={};    // tab2[]==N
    std::complex<double>  w1[1]={{1,1}};
    std::complex<double>  w2[1]={{1,1}};
    std::complex<double>  w3[1]={{1,1}};
/*
for(int j=0;j<N;j++)
{
   cout<< tab[j]<<"  ";
}
cout<<endl<<"........."<<endl;
*/
    //radix-4
    for(int k=0;k<4;k++)
    {
        for(int n=0;n<4;n++)
        {
          w2[0].real()=cos(n*k*2*pi/(N/4));
          w2[0].imag()=-sin(n*k*2*pi/(N/4));

          tab2[4*k+0]=tab2[4*k+0]+(tab[n]+tab[n+4]+tab[n+8]+tab[n+12])*w2[0];

          w1[0].real()=cos(n*2*pi/N);
          w1[0].imag()=-sin(n*2*pi/N);
          w3[0].real()=0;
          w3[0].imag()=-1;
          tab2[4*k+1]=tab2[4*k+1]+(w3[0]*(tab[n]-tab[n+8])+(tab[n+4]-tab[n+12]))*w1[0]*w2[0];

          w1[0].real()=cos(2*n*2*pi/N);
          w1[0].imag()=-sin(2*n*2*pi/N);
          tab2[4*k+2]=tab2[4*k+2]+(tab[n]-tab[n+4]+tab[n+8]-tab[n+12])*w1[0]*w2[0];

          w1[0].real()=cos(3*n*2*pi/N);
          w1[0].imag()=-sin(3*n*2*pi/N);
          w3[0].real()=0;
          w3[0].imag()=1;
          tab2[4*k+3]=tab2[4*k+3]+(w3[0]*(tab[n]-tab[n+8])+(tab[n+4]-tab[n+12]))*w1[0]*w2[0];
        }
    }

    for(int j=0;j<N;j++)
    {
      tab[j].real() =tab2[j].real()*2/N;
      tab[j].imag() =tab2[j].imag()*2/N;
    }
}

 void fun_fourier_transform_FFT_radix_4_N_16_method4(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double> tab2[64]={};    // tab2[]==N
    std::complex<double>  w1[1]={{1,1}};
    std::complex<double>  w2[1]={{1,1}};
    std::complex<double>  w3[1]={{1,1}};
/*
for(int j=0;j<N;j++)
{
   cout<< tab[j]<<"  ";
}
cout<<endl<<"........."<<endl;
*/

        for(int i=0;i<4;i++)
        {
          tab2[4*i]=(tab[0+i]+tab[4+i]+tab[8+i]+tab[12+i]);
        }

        for(int i=0;i<4;i++)
        {
          w1[0].real()=cos(i*2*pi/N);
          w1[0].imag()=-sin(i*2*pi/N);
          w2[0].real()=0;
          w2[0].imag()=-1;
          tab2[4*i+1]=(tab[0+i]-tab[i+8]+w2[0]*(tab[i+4]-tab[i+12]))*w1[0];
        }

        for(int i=0;i<4;i++)
        {
          w1[0].real()=cos(2*i*2*pi/N);
          w1[0].imag()=-sin(2*i*2*pi/N);
          tab2[4*i+2]=(tab[i]-tab[i+4]+tab[i+8]-tab[i+12])*w1[0];
        }

        for(int i=0;i<4;i++)
        {
          w1[0].real()=cos(3*i*2*pi/N);
          w1[0].imag()=-sin(3*i*2*pi/N);
          w3[0].real()=0;
          w3[0].imag()=1;
          tab2[4*i+3]=(tab[i]-tab[i+8]+w3[0]*(tab[i+4]-tab[i+12]))*w1[0];
        }

///////////////////////////////////
for(int k=0;k<N;k++)
{
    tab[k]={0,0};
}

        for(int i=0;i<4;i++)
        {
        for(int j=0;j<13;j=j+4)
        {
          w1[0].real()=cos(i*j/4*2*pi/(N/4));
          w1[0].imag()=-sin(i*j/4*2*pi/(N/4));
          tab[4*i]=tab[4*i]+tab2[j]*w1[0];
          tab[4*i+1]=tab[4*i+1]+tab2[j+1]*w1[0];
          tab[4*i+2]=tab[4*i+2]+tab2[j+2]*w1[0];
          tab[4*i+3]=tab[4*i+3]+tab2[j+3]*w1[0];
        }
        }

    for(int j=0;j<N;j++)
    {
      tab[j].real() =tab[j].real()*2/N;
      tab[j].imag() =tab[j].imag()*2/N;
    }
}




void fun_fourier_transform_FFT_radix_4_N_16_official(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double> tab2[64]={};    // tab2[]==N
    std::complex<double>  w[1]={{1,1}};
    std::complex<double>  w1[1]={{1,1}};
    std::complex<double>  w2[1]={{1,1}};
    std::complex<double>  w3[1]={{1,1}};
    std::complex<double>  w4[1]={{1,1}};
    std::complex<double>  w5[1]={{1,1}};
/*
for(int j=0;j<N;j++)
{
   cout<< tab[j]<<"  ";
}
cout<<endl<<"........."<<endl;
*/
        w[0]={{0,1}};
        //radix-4
        for(int i=0;i<(N/4);i++)
        {
          w4[0].real()=cos(0*2*pi/(N/4));
          w4[0].imag()=-sin(0*2*pi/(N/4));
          tab2[i]=(tab[i+0]+tab[i+4]+tab[i+8]+tab[i+12]);

          w[0].real()=cos(0*2*pi/N);
          w[0].imag()=-sin(0*2*pi/N);
          w4[0].real()=cos(0*2*pi/(N/4));
          w4[0].imag()=-sin(0*2*pi/(N/4));
          w2[0].real()=0;
          w2[0].imag()=-1;
          tab2[i+4]=(tab[i+0]-tab[i+8]+w2[0]*(tab[i+4]-tab[i+12]));

          w[0].real()=cos(0*2*pi/N);
          w[0].imag()=-sin(0*2*pi/N);
          w4[0].real()=cos(0*2*pi/(N/4));
          w4[0].imag()=-sin(0*2*pi/(N/4));
          tab2[i+8]=(tab[i+0]-tab[i+4]+tab[i+8]-tab[i+12]);

          w[0].real()=cos(0*2*pi/N);
          w[0].imag()=-sin(0*2*pi/N);
          w4[0].real()=cos(0*2*pi/(N/4));
          w4[0].imag()=-sin(0*2*pi/(N/4));
          w2[0].real()=0;
          w2[0].imag()=1;
          tab2[i+12]=(tab[i+0]-tab[i+8]+w2[0]*(tab[i+4]-tab[i+12]));
        }
///////////////////////////////////
for(int k=0;k<N;k++)
{
    tab[k]={0,0};
}
        for(int i=0;i<4;i++)
        {
          //cout<<i<<endl;
          w1[0].real()=cos(0*2*pi/(N/1));
          w1[0].imag()=-sin(0*2*pi/(N/1));
          w2[0].real()=cos(1*i*2*pi/(N/1));
          w2[0].imag()=-sin(1*i*2*pi/(N/1));
          w3[0].real()=cos(2*i*2*pi/(N/1));
          w3[0].imag()=-sin(2*i*2*pi/(N/1));
          w4[0].real()=cos(3*i*2*pi/(N/1));
          w4[0].imag()=-sin(3*i*2*pi/(N/1));
          w5[0].real()=0;
          w5[0].imag()=-1;
          tab[4*i+0]=(w1[0]*tab2[4*i+0]+w2[0]*tab2[4*i+1]+w3[0]*tab2[4*i+2]+w4[0]*tab2[4*i+3]);
          tab[4*i+1]=(w1[0]*tab2[4*i+0]-w3[0]*tab2[4*i+2]+w5[0]*(w2[0]*tab2[4*i+1]-w4[0]*tab2[4*i+3]));
          tab[4*i+2]=(w1[0]*tab2[4*i+0]-w2[0]*tab2[4*i+1]+w3[0]*tab2[4*i+2]-w4[0]*tab2[4*i+3]);
          w5[0].real()=0;
          w5[0].imag()=1;
          tab[4*i+3]=(w1[0]*tab2[4*i+0]-w3[0]*tab2[4*i+2]+w5[0]*(w2[0]*tab2[4*i+1]-w4[0]*tab2[4*i+3]));
        }


    for(int j=0;j<N;j++)
    {
      tab[j].real() =tab[j].real()*2/N;
      tab[j].imag() =tab[j].imag()*2/N;
    }
}

void fun_inverse_fourier_transform_FFT_radix_4_N_16_official(int N,std::complex<double> tab[])
{
     const double pi=3.141592653589793238462;
    std::complex<double> tab2[64]={};    // tab2[]==N
    std::complex<double>  w[1]={{1,1}};
    std::complex<double>  w1[1]={{1,1}};
    std::complex<double>  w2[1]={{1,1}};
    std::complex<double>  w3[1]={{1,1}};
    std::complex<double>  w4[1]={{1,1}};
    std::complex<double>  w5[1]={{1,1}};
/*
for(int j=0;j<N;j++)
{
   cout<< tab[j]<<"  ";
}
cout<<endl<<"........."<<endl;
*/
        w[0]={{0,1}};
        //radix-4
        for(int i=0;i<(N/4);i++)
        {
          w4[0].real()=cos(0*2*pi/(N/4));
          w4[0].imag()=sin(0*2*pi/(N/4));
          tab2[i]=(tab[i+0]+tab[i+4]+tab[i+8]+tab[i+12]);

          w[0].real()=cos(0*2*pi/N);
          w[0].imag()=sin(0*2*pi/N);
          w4[0].real()=cos(0*2*pi/(N/4));
          w4[0].imag()=sin(0*2*pi/(N/4));
          w2[0].real()=0;
          w2[0].imag()=1;
          tab2[i+4]=(tab[i+0]-tab[i+8]+w2[0]*(tab[i+4]-tab[i+12]))*w4[0];

          w[0].real()=cos(0*2*pi/N);
          w[0].imag()=sin(0*2*pi/N);
          w4[0].real()=cos(0*2*pi/(N/4));
          w4[0].imag()=sin(0*2*pi/(N/4));
          tab2[i+8]=(tab[i+0]-tab[i+4]+tab[i+8]-tab[i+12])*w4[0];

          w[0].real()=cos(0*2*pi/N);
          w[0].imag()=sin(0*2*pi/N);
          w4[0].real()=cos(0*2*pi/(N/4));
          w4[0].imag()=sin(0*2*pi/(N/4));
          w2[0].real()=0;
          w2[0].imag()=-1;
          tab2[i+12]=(tab[i+0]-tab[i+8]+w2[0]*(tab[i+4]-tab[i+12]))*w4[0];
        }
///////////////////////////////////
for(int k=0;k<N;k++)
{
    tab[k]={0,0};
}
        for(int i=0;i<4;i++)
        {
          //cout<<i<<endl;
          w1[0].real()=cos(0*2*pi/(N/1));
          w1[0].imag()=sin(0*2*pi/(N/1));
          w2[0].real()=cos(1*i*2*pi/(N/1));
          w2[0].imag()=sin(1*i*2*pi/(N/1));
          w3[0].real()=cos(2*i*2*pi/(N/1));
          w3[0].imag()=sin(2*i*2*pi/(N/1));
          w4[0].real()=cos(3*i*2*pi/(N/1));
          w4[0].imag()=sin(3*i*2*pi/(N/1));
          w5[0].real()=0;
          w5[0].imag()=1;
          tab[4*i+0]=(w1[0]*tab2[4*i+0]+w2[0]*tab2[4*i+1]+w3[0]*tab2[4*i+2]+w4[0]*tab2[4*i+3]);
          tab[4*i+1]=(w1[0]*tab2[4*i+0]-w3[0]*tab2[4*i+2]+w5[0]*(w2[0]*tab2[4*i+1]-w4[0]*tab2[4*i+3]));
          tab[4*i+2]=(w1[0]*tab2[4*i+0]-w2[0]*tab2[4*i+1]+w3[0]*tab2[4*i+2]-w4[0]*tab2[4*i+3]);
          w5[0].real()=0;
          w5[0].imag()=-1;
          tab[4*i+3]=(w1[0]*tab2[4*i+0]-w3[0]*tab2[4*i+2]+w5[0]*(w2[0]*tab2[4*i+1]-w4[0]*tab2[4*i+3]));
        }


    for(int j=0;j<N;j++)
    {
      tab[j].real() =tab[j].real()*0.5;
      tab[j].imag() =tab[j].imag()*0.5;
    }
}