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

 //source:
//https://www.google.ch/patents/US6957241
//http://www.google.ch/patents/US20020083107
//https://www.beechwood.eu/fft-implementation-r2-dit-r4-dif-r8-dif/
//author marcin matysek (r)ewertyn.PL

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

//haven't try it with other function that cos(x)+jsin(x)=sin(x+pi/2)+jsin(x)


//fourier transform iFFT radix-4 for N=256 algorithm c++ source code implementation

#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_256_official(int N,std::complex<double> tab[]);
void fun_inverse_fourier_transform_FFT_radix_4_N_256_official(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=256;
    std::complex<double> tab2[256]={};
    std::complex<double> tab3[256]={};
    for(int i=0;i<256;i++)
    {
        tab2[i].real()=i;
        tab2[i].imag()=256+i;
    }

    double time2;
    double zmienna=0;

    cout<<"signal="<<endl;
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab2[j].real()*1000)/1000<<"  ";
    }
    cout<<endl;
    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_256_official(N,tab2);
    fun_inverse_bits_radix_4(N,tab2);
    ///////////////////////////////////////////////////////////

    time2=diffclock( start, clock() );

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


    system("pause");
    cout<<endl;

    /////////////////////////////////////////////////////////////////
    fun_inverse_fourier_transform_FFT_radix_4_N_256_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;
    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;
    system("pause");
    return 0;
}



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

}



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

    const double pi=3.141592653589793238462;
    std::complex<double> tab2[256]={};    // tab2[]==N
    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}};
    int flag=0;

        //radix-4
//stage 1
        for(int i=0;i<(N/4);i++)//64
        {
          w1[0].real()=cos((0)*2*pi/(N/1));
          w1[0].imag()=sin((0)*2*pi/(N/1));
          w2[0].real()=cos((0)*2*pi/(N/1));
          w2[0].imag()=sin((0)*2*pi/(N/1));
          w3[0].real()=cos((0)*2*pi/(N/1));
          w3[0].imag()=sin((0)*2*pi/(N/1));
          w4[0].real()=cos((0)*2*pi/(N/1));
          w4[0].imag()=sin((0)*2*pi/(N/1));

          tab2[i]     =(w1[0]*tab[i+0]+w2[0]*tab[i+N/4]+w3[0]*tab[i+N/2]+w4[0]*tab[i+3*N/4]);

          w5[0].real()=0;
          w5[0].imag()=1;
          tab2[i+N/4] =(w1[0]*tab[i+0]-w3[0]*tab[i+N/2]+w5[0]*(w2[0]*tab[i+N/4]-w4[0]*tab[i+3*N/4]));

          tab2[i+N/2] =(w1[0]*tab[i+0]-w2[0]*tab[i+N/4]+w3[0]*tab[i+N/2]-w4[0]*tab[i+3*N/4]);

          w5[0].real()=0;
          w5[0].imag()=-1;
          tab2[i+3*N/4]=(w1[0]*tab[i+0]-w3[0]*tab[i+N/2]+w5[0]*(w2[0]*tab[i+N/4]-w4[0]*tab[i+3*N/4]));
        }
///////////////////////////////////

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

        for(int i=0;i<(N/64);i++)//4
        {
        for(int j=0;j<(N/16);j++)//16
        {
          w1[0].real()=cos((0*i+i*j)*2*pi/(N/1));
          w1[0].imag()=sin((0*i+i*j)*2*pi/(N/1));
          w2[0].real()=cos((16*i+i*j)*2*pi/(N/1));
          w2[0].imag()=sin((16*i+i*j)*2*pi/(N/1));
          w3[0].real()=cos((2*16*i+i*j)*2*pi/(N/1));
          w3[0].imag()=sin((2*16*i+i*j)*2*pi/(N/1));
          w4[0].real()=cos((3*16*i+i*j)*2*pi/(N/1));
          w4[0].imag()=sin((3*16*i+i*j)*2*pi/(N/1));

          w5[0].real()=0;
          w5[0].imag()=1;
          tab[0+j+N/(4*1)*i]         =(w1[0]*tab2[0+j+N/(4*1)*i]+w2[0]*tab2[N/(16*1)+j+N/(4*1)*i]+        w3[0]*tab2[N/(8*1) +j+N/(4*1)*i]+w4[0]*tab2[3*N/(16*1)+j+N/(4*1)*i]);
          tab[N/(16*1)+j+N/(4*1)*i]  =(w1[0]*tab2[0+j+N/(4*1)*i]-w3[0]*tab2[N/(8* 1) +j+N/(4*1)*i]+w5[0]*(w2[0]*tab2[N/(16*1)+j+N/(4*1)*i]-w4[0]*tab2[3*N/(16*1)+j+N/(4*1)*i]));
          tab[N/(8 *1)+j+N/(4*1)*i]  =(w1[0]*tab2[0+j+N/(4*1)*i]-w2[0]*tab2[N/(16*1)+j+N/(4*1)*i]+        w3[0]*tab2[N/(8*1) +j+N/(4*1)*i]-w4[0]*tab2[3*N/(16*1)+j+N/(4*1)*i]);
          w5[0].real()=0;
          w5[0].imag()=-1;
          tab[3*N/(16*1)+j+N/(4*1)*i]=(w1[0]*tab2[0+j+N/(4*1)*i]-w3[0]*tab2[N/(8* 1) +j+N/(4*1)*i]+w5[0]*(w2[0]*tab2[N/(16*1)+j+N/(4*1)*i]-w4[0]*tab2[3*N/(16*1)+j+N/(4*1)*i]));
        }
        }

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

//stage 3

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

        for(int i=0;i<(N/16);i++)//16
        {
        for(int j=0;j<4;j++)
        {
          w1[0].real()=cos((i-4*flag)*(16*0+4*j)*2*pi/(N/1));
          w1[0].imag()=sin((i-4*flag)*(16*0+4*j)*2*pi/(N/1));
          w2[0].real()=cos((i-4*flag)*(1*16+4*j)*2*pi/(N/1));
          w2[0].imag()=sin((i-4*flag)*(1*16+4*j)*2*pi/(N/1));
          w3[0].real()=cos((i-4*flag)*(2*16+4*j)*2*pi/(N/1));
          w3[0].imag()=sin((i-4*flag)*(2*16+4*j)*2*pi/(N/1));
          w4[0].real()=cos((i-4*flag)*(3*16+4*j)*2*pi/(N/1));
          w4[0].imag()=sin((i-4*flag)*(3*16+4*j)*2*pi/(N/1));

          w5[0].real()=0;
          w5[0].imag()=1;
          tab2[0+j+N/(4*4)*i]         =(w1[0]*tab[0+j+N/(4*4)*i]+w2[0]*tab[N/(16*4)+j+N/(4*4)*i]+       w3[0]*tab[N/(8*4) +j+N/(4*4)*i]+w4[0]*tab[3*N/(16*4)+j+N/(4*4)*i]);
          tab2[N/(16*4)+j+N/(4*4)*i]  =(w1[0]*tab[0+j+N/(4*4)*i]-w3[0]*tab[N/(8*4) +j+N/(4*4)*i]+w5[0]*(w2[0]*tab[N/(16*4)+j+N/(4*4)*i]-w4[0]*tab[3*N/(16*4)+j+N/(4*4)*i]));
          tab2[N/(8*4)+j+N/(4*4)*i]   =(w1[0]*tab[0+j+N/(4*4)*i]-w2[0]*tab[N/(16*4)+j+N/(4*4)*i]+       w3[0]*tab[N/(8*4) +j+N/(4*4)*i]-w4[0]*tab[3*N/(16*4)+j+N/(4*4)*i]);
          w5[0].real()=0;
          w5[0].imag()=-1;
          tab2[3*N/(16*4)+j+N/(4*4)*i]=(w1[0]*tab[0+j+N/(4*4)*i]-w3[0]*tab[N/(8*4) +j+N/(4*4)*i]+w5[0]*(w2[0]*tab[N/(16*4)+j+N/(4*4)*i]-w4[0]*tab[3*N/(16*4)+j+N/(4*4)*i]));
         }
         if((i+1)%4==0){flag++;}
        }
///////////////////////////

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

}

        for(int i=0;i<(N/16);i++)
        {
        for(int j=0;j<4;j++)
        {
          w1[0].real()=cos(0*16*j*2*pi/(N/1));
          w1[0].imag()=sin(0*16*j*2*pi/(N/1));
          w2[0].real()=cos(1*16*j*2*pi/(N/1));
          w2[0].imag()=sin(1*16*j*2*pi/(N/1));
          w3[0].real()=cos(2*16*j*2*pi/(N/1));
          w3[0].imag()=sin(2*16*j*2*pi/(N/1));
          w4[0].real()=cos(3*16*j*2*pi/(N/1));
          w4[0].imag()=sin(3*16*j*2*pi/(N/1));
          w5[0].real()=0;
          w5[0].imag()=1;

          tab[N/(4*4)*i+0+4*j]=(w1[0]*tab2[N/(4*4)*i+0+4*j]+w2[0]*tab2[N/(4*4)*i+1+4*j]+       w3[0]*tab2[N/(4*4)*i+2+4*j]+w4[0]*tab2[N/(4*4)*i+3+4*j]);
          tab[N/(4*4)*i+1+4*j]=(w1[0]*tab2[N/(4*4)*i+0+4*j]-w3[0]*tab2[N/(4*4)*i+2+4*j]+w5[0]*(w2[0]*tab2[N/(4*4)*i+1+4*j]-w4[0]*tab2[N/(4*4)*i+3+4*j]));
          tab[N/(4*4)*i+2+4*j]=(w1[0]*tab2[N/(4*4)*i+0+4*j]-w2[0]*tab2[N/(4*4)*i+1+4*j]+       w3[0]*tab2[N/(4*4)*i+2+4*j]-w4[0]*tab2[N/(4*4)*i+3+4*j]);
          w5[0].real()=0;
          w5[0].imag()=-1;
          tab[N/(4*4)*i+3+4*j]=(w1[0]*tab2[N/(4*4)*i+0+4*j]-w3[0]*tab2[N/(4*4)*i+2+4*j]+w5[0]*(w2[0]*tab2[N/(4*4)*i+1+4*j]-w4[0]*tab2[N/(4*4)*i+3+4*j]));
        }
        }

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

}
//////////////////////////////
//////////////////////////////

void fun_fourier_transform_FFT_radix_4_N_256_official(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double> tab2[256]={};    // tab2[]==N
    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}};
    int flag=0;

        //radix-4
//stage 1
        for(int i=0;i<(N/4);i++)//64
        {
          w1[0].real()=cos((0)*2*pi/(N/1));
          w1[0].imag()=-sin((0)*2*pi/(N/1));
          w2[0].real()=cos((0)*2*pi/(N/1));
          w2[0].imag()=-sin((0)*2*pi/(N/1));
          w3[0].real()=cos((0)*2*pi/(N/1));
          w3[0].imag()=-sin((0)*2*pi/(N/1));
          w4[0].real()=cos((0)*2*pi/(N/1));
          w4[0].imag()=-sin((0)*2*pi/(N/1));

          tab2[i]     =(w1[0]*tab[i+0]+w2[0]*tab[i+N/4]+w3[0]*tab[i+N/2]+w4[0]*tab[i+3*N/4]);

          w5[0].real()=0;
          w5[0].imag()=-1;
          tab2[i+N/4] =(w1[0]*tab[i+0]-w3[0]*tab[i+N/2]+w5[0]*(w2[0]*tab[i+N/4]-w4[0]*tab[i+3*N/4]));

          tab2[i+N/2] =(w1[0]*tab[i+0]-w2[0]*tab[i+N/4]+w3[0]*tab[i+N/2]-w4[0]*tab[i+3*N/4]);

          w5[0].real()=0;
          w5[0].imag()=1;
          tab2[i+3*N/4]=(w1[0]*tab[i+0]-w3[0]*tab[i+N/2]+w5[0]*(w2[0]*tab[i+N/4]-w4[0]*tab[i+3*N/4]));
        }
///////////////////////////////////


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

        for(int i=0;i<(N/64);i++)//4
        {
        for(int j=0;j<(N/16);j++)//16
        {
          w1[0].real()=cos((0*i+i*j)*2*pi/(N/1));
          w1[0].imag()=-sin((0*i+i*j)*2*pi/(N/1));
          w2[0].real()=cos((16*i+i*j)*2*pi/(N/1));
          w2[0].imag()=-sin((16*i+i*j)*2*pi/(N/1));
          w3[0].real()=cos((2*16*i+i*j)*2*pi/(N/1));
          w3[0].imag()=-sin((2*16*i+i*j)*2*pi/(N/1));
          w4[0].real()=cos((3*16*i+i*j)*2*pi/(N/1));
          w4[0].imag()=-sin((3*16*i+i*j)*2*pi/(N/1));

          w5[0].real()=0;
          w5[0].imag()=-1;
          tab[0+j+N/(4*1)*i]         =(w1[0]*tab2[0+j+N/(4*1)*i]+w2[0]*tab2[N/(16*1)+j+N/(4*1)*i]+        w3[0]*tab2[N/(8*1) +j+N/(4*1)*i]+w4[0]*tab2[3*N/(16*1)+j+N/(4*1)*i]);
          tab[N/(16*1)+j+N/(4*1)*i]  =(w1[0]*tab2[0+j+N/(4*1)*i]-w3[0]*tab2[N/(8* 1) +j+N/(4*1)*i]+w5[0]*(w2[0]*tab2[N/(16*1)+j+N/(4*1)*i]-w4[0]*tab2[3*N/(16*1)+j+N/(4*1)*i]));
          tab[N/(8 *1)+j+N/(4*1)*i]  =(w1[0]*tab2[0+j+N/(4*1)*i]-w2[0]*tab2[N/(16*1)+j+N/(4*1)*i]+        w3[0]*tab2[N/(8*1) +j+N/(4*1)*i]-w4[0]*tab2[3*N/(16*1)+j+N/(4*1)*i]);
          w5[0].real()=0;
          w5[0].imag()=1;
          tab[3*N/(16*1)+j+N/(4*1)*i]=(w1[0]*tab2[0+j+N/(4*1)*i]-w3[0]*tab2[N/(8* 1) +j+N/(4*1)*i]+w5[0]*(w2[0]*tab2[N/(16*1)+j+N/(4*1)*i]-w4[0]*tab2[3*N/(16*1)+j+N/(4*1)*i]));
        }
        }

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

//stage 3
int tab4[256]={0};
for(int k=0;k<N;k++)
{
    tab2[k]={0,0};
}

        for(int i=0;i<(N/16);i++)//16
        {
        for(int j=0;j<4;j++)
        {
          w1[0].real()=cos((i-4*flag)*(16*0+4*j)*2*pi/(N/1));
          w1[0].imag()=-sin((i-4*flag)*(16*0+4*j)*2*pi/(N/1));
          w2[0].real()=cos((i-4*flag)*(1*16+4*j)*2*pi/(N/1));
          w2[0].imag()=-sin((i-4*flag)*(1*16+4*j)*2*pi/(N/1));
          w3[0].real()=cos((i-4*flag)*(2*16+4*j)*2*pi/(N/1));
          w3[0].imag()=-sin((i-4*flag)*(2*16+4*j)*2*pi/(N/1));
          w4[0].real()=cos((i-4*flag)*(3*16+4*j)*2*pi/(N/1));
          w4[0].imag()=-sin((i-4*flag)*(3*16+4*j)*2*pi/(N/1));

          w5[0].real()=0;
          w5[0].imag()=-1;
          tab2[0+j+N/(4*4)*i]         =(w1[0]*tab[0+j+N/(4*4)*i]+w2[0]*tab[N/(16*4)+j+N/(4*4)*i]+       w3[0]*tab[N/(8*4) +j+N/(4*4)*i]+w4[0]*tab[3*N/(16*4)+j+N/(4*4)*i]);
          tab2[N/(16*4)+j+N/(4*4)*i]  =(w1[0]*tab[0+j+N/(4*4)*i]-w3[0]*tab[N/(8*4) +j+N/(4*4)*i]+w5[0]*(w2[0]*tab[N/(16*4)+j+N/(4*4)*i]-w4[0]*tab[3*N/(16*4)+j+N/(4*4)*i]));
          tab2[N/(8*4)+j+N/(4*4)*i]   =(w1[0]*tab[0+j+N/(4*4)*i]-w2[0]*tab[N/(16*4)+j+N/(4*4)*i]+       w3[0]*tab[N/(8*4) +j+N/(4*4)*i]-w4[0]*tab[3*N/(16*4)+j+N/(4*4)*i]);
          w5[0].real()=0;
          w5[0].imag()=1;
          tab2[3*N/(16*4)+j+N/(4*4)*i]=(w1[0]*tab[0+j+N/(4*4)*i]-w3[0]*tab[N/(8*4) +j+N/(4*4)*i]+w5[0]*(w2[0]*tab[N/(16*4)+j+N/(4*4)*i]-w4[0]*tab[3*N/(16*4)+j+N/(4*4)*i]));

         tab4[0+j+N/(4*4)*i] =(i-4*flag)*(16*0+4*j);
         tab4[N/(16*4)+j+N/(4*4)*i] =(i-4*flag)*(1*16+4*j);
         tab4[N/(8*4)+j+N/(4*4)*i]=(i-4*flag)*(2*16+4*j);
         tab4[3*N/(16*4)+j+N/(4*4)*i]=(i-4*flag)*(3*16+4*j);
         }

         if((i+1)%4==0){flag++;}
        }

        /*
        for(int i2=0;i2<256;i2++)
        {
            cout<<tab4[i2]<<endl;
          if((i2+1)%16==0){system("pause");}
        }
        */
///////////////////////////

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

}

        for(int i=0;i<(N/16);i++)
        {
        for(int j=0;j<4;j++)
        {
          w1[0].real()=cos(0*16*j*2*pi/(N/1));
          w1[0].imag()=-sin(0*16*j*2*pi/(N/1));
          w2[0].real()=cos(1*16*j*2*pi/(N/1));
          w2[0].imag()=-sin(1*16*j*2*pi/(N/1));
          w3[0].real()=cos(2*16*j*2*pi/(N/1));
          w3[0].imag()=-sin(2*16*j*2*pi/(N/1));
          w4[0].real()=cos(3*16*j*2*pi/(N/1));
          w4[0].imag()=-sin(3*16*j*2*pi/(N/1));
          w5[0].real()=0;
          w5[0].imag()=-1;

          tab[N/(4*4)*i+0+4*j]=(w1[0]*tab2[N/(4*4)*i+0+4*j]+w2[0]*tab2[N/(4*4)*i+1+4*j]+       w3[0]*tab2[N/(4*4)*i+2+4*j]+w4[0]*tab2[N/(4*4)*i+3+4*j]);
          tab[N/(4*4)*i+1+4*j]=(w1[0]*tab2[N/(4*4)*i+0+4*j]-w3[0]*tab2[N/(4*4)*i+2+4*j]+w5[0]*(w2[0]*tab2[N/(4*4)*i+1+4*j]-w4[0]*tab2[N/(4*4)*i+3+4*j]));
          tab[N/(4*4)*i+2+4*j]=(w1[0]*tab2[N/(4*4)*i+0+4*j]-w2[0]*tab2[N/(4*4)*i+1+4*j]+       w3[0]*tab2[N/(4*4)*i+2+4*j]-w4[0]*tab2[N/(4*4)*i+3+4*j]);
          w5[0].real()=0;
          w5[0].imag()=1;
          tab[N/(4*4)*i+3+4*j]=(w1[0]*tab2[N/(4*4)*i+0+4*j]-w3[0]*tab2[N/(4*4)*i+2+4*j]+w5[0]*(w2[0]*tab2[N/(4*4)*i+1+4*j]-w4[0]*tab2[N/(4*4)*i+3+4*j]));

         tab4[N/(4*4)*i+0+4*j] =0*16*j;
         tab4[N/(4*4)*i+1+4*j] =1*16*j;
         tab4[N/(4*4)*i+2+4*j]=2*16*j;
         tab4[N/(4*4)*i+3+4*j]=3*16*j;
         }
        }
        /*
         for(int i2=0;i2<256;i2++)
        {
            cout<<tab4[i2]<<endl;
          if((i2+1)%16==0){system("pause");}
        }
        */
    for(int j=0;j<N;j++)
    {
      tab[j].real() =tab[j].real()*2/N;
      tab[j].imag() =tab[j].imag()*2/N;
    }

}