fourier transform iFFT radix-4 for N=4 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=4 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_4(int N,std::complex<double> tab[]);
void fun_inverse_fourier_transform_FFT_radix_4_N_4(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=4;
    std::complex<double> tab2[4]={{1,5},    {2,11}    ,{3,18},    {4,6}};
    //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;

    clock_t start = clock();


    fun_inverse_bits_radix_4(N,tab2);
    fun_fourier_transform_FFT_radix_4_N_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_bits_radix_4(N,tab2);
    fun_inverse_fourier_transform_FFT_radix_4_N_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_4(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>  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;
*/
          //radix-4
          w[0]={{0,1}};
          w4[0].real()=cos(0*2*pi/(N/4));
          w4[0].imag()=-sin(0*2*pi/(N/4));
          tab2[0]=(tab[0]+tab[1]+tab[2]+tab[3])*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;
          //w3[0].real()=0;
          //w3[0].imag()=-1;
          tab2[1]=(tab[0]-tab[2]+w2[0]*(-tab[1]+tab[3]))*w[0]*w4[0];
          //tab2[1]=(tab[0]-tab[2]+w3[0]*tab[1]+w2[0]*tab[3])*w[0]*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[2]=(tab[0]-tab[1]+tab[2]-tab[3])*w[0]*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;
          //w3[0].real()=0;
          //w3[0].imag()=-1;
          tab2[3]=(tab[0]-tab[2]+w2[0]*(tab[1]-tab[3]))*w[0]*w4[0];
          //tab2[3]=(tab[0]+w2[0]*tab[1]-tab[2]+w3[0]*tab[3])*w[0]*w4[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_4(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>  w2[1]={{1,1}};
    std::complex<double>  w3[1]={{1,1}};
    std::complex<double>  w4[1]={{1,1}};
/*
for(int j=0;j<N;j++)
{
   cout<< tab[j]<<"  ";
}
cout<<endl<<"........."<<endl;
*/
          //radix-4
          w[0]={{0,1}};
          w4[0].real()=cos(0*2*pi/(N/4));
          w4[0].imag()=sin(0*2*pi/(N/4));
          tab2[0]=(tab[0]+tab[1]+tab[2]+tab[3])*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;
          //w3[0].real()=0;
          //w3[0].imag()=1;
          tab2[1]=(tab[0]-tab[2]+w2[0]*(-tab[1]+tab[3]))*w[0]*w4[0];
          //tab2[1]=(tab[0]+w3[0]*tab[1]-tab[2]+w2[0]*tab[3])*w[0]*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[2]=(tab[0]-tab[1]+tab[2]-tab[3])*w[0]*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;
          //w3[0].real()=0;
          //w3[0].imag()=1;
          tab2[3]=(tab[0]-tab[2]+w2[0]*(tab[1]-tab[3]))*w[0]*w4[0];
          //tab2[3]=(tab[0]+w2[0]*tab[1]-tab[2]+w3[0]*tab[3])*w[0]*w4[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[])
{
//code by Sidney Burrus
//http://dsp.stackexchange.com/questions/3481/radix-4-fft-implementation
    std::complex<double> t;
    //N=4^a;
    // Radix-4 bit-reverse
    double T;
    int j = 0;
    int N2 = N>>2;
    int N1=0;
    for (int i=0; i < N-1; i++) {
        if (i < j) {
            t = tab[i];
            tab[i] = tab[j];
            tab[j] = t;
        }
        N1 = N2;
        while ( j >= 3*N1 ) {
            j -= 3*N1;
            N1 >>= 2;
        }
        j += N1;
    }
}




///////////////////////////////////////////