something extra

inverse fourier transform iDFT ifft 4 methods in open office + something extra

http://www.mediafire.com/file/hyz4dbski4w00pb/inverse+fourier+transform+iDFT+ifft+4+methods+in+open+office+++something+extra.ods


http://www.mediafire.com/file/59bpnci966ulec9/DFT+FFT+RADIX-2+DIT+algorytm+Transformacja+Fouriera+analitycznie+v3.4.xlsx

universal mixed radix-2-3-4-5-7-11 inverse iFFT standard version

https://github.com/rewertynpl/mixedFFT

other example how it is possible to calculate inverse fourier transform iFFT

//other example how it is possible to calculate inverse fourier transform iFFT only example withoud implementation FFT algorithm
//author marcin matysek (r)ewertyn.PL

 //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/
    //http://www.cmlab.csie.ntu.edu.tw/cml/dsp/training/coding/transform/fft.html
    //http://dsp.stackexchange.com/questions/3481/radix-4-fft-implementation

    //https://community.arm.com/graphics/b/blog/posts/speeding-up-fast-fourier-transform-mixed-radix-on-mobile-arm-mali-gpu-by-means-of-opencl---part-1
    //book: "Cyfrowe przetwarzanie sygnalow" - Tomasz Zielinski it has quick version of radix-2 because it calculates less sin() and cos()




#include <iostream>
#include "conio.h"
#include <stdlib.h>
#include <math.h>


using namespace std;

void fun_inverse_bits(int N,double *tablica2);
void fast_fourier_transform_FFT(int N,double *tablica2);
//not used:
void inverse_fast_fourier_transform_FFT(int N,double *tablica2);

//method nr 2
//for inverse with standard fourier FFT requires 2 times less storing information
double scaling_re_and_im(int N,double *tablica);
//method 2 works only if audio samples are not complex numbers but result of FFT can be comlex numbers
//that are combinet to normal number


int main()
{
    //assumption: if N==signal period in table tab[] then resolution = 1 Hz but N=2^b;
    //when we have signal period 22000 Hz and N=2^15=32768 then
    //the fundamental frequency is 0,671 Hz
    //that means that in F(1) is 1*0,671 Hz =0,671 Hz
    //in F(2) is 2*0,671 Hz =1,342 Hz
    //in F(9) is 9*0,671 Hz =6,039 Hz
    //in F(30) is 30*0,671 Hz =20,13 Hz that means in F(30) we will see what value has our signal in 20,13 Hz


    int N=16;//N=2^b;
     //tab[A] A=2*N;

    // signal period=16  samples= 16
    //double tab[32]={-0.923879533,0.382683432,1.03153013,0.923879533,0.923879533,1.465075634,1.796896994,0.923879533,-0.923879533,
    //-2.230442498,-1.796896994,-0.158512669,0.923879533,0.382683432,-1.03153013,-1.689246397};

   // signal period=12  samples= 12
    double tab[32]={-0.923879533,0.70664666,0.996551596,0.923879533,1.659378744,1.369473808,-0.923879533,
    -2.29335334,-0.735499212,0.923879533,-0.072672064,-1.630526192};

    //tab[0 to N-1]  -re numbers
    //tab[N-1 to 2N-1]  -im numbers

//tab[21]=-4.93;//im number


cout<<"next method only normal numbers"<<endl<<endl<<endl;
system("pause");
/////////////////////////

cout<<"     method 2 "<<endl<<endl;

//It works only if audio samples are not complex numbers but result of FFT can be comlex numbers
//that are combinet to normal number

 cout<<"signal g(x)"<<endl;

    for(int i=0;i<2;i++)
    {
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab[j+i*16]*1000)/1000<<"  ";
    }
    cout<<endl;
    }
    cout<<endl;

 fun_inverse_bits(N,tab);
 fast_fourier_transform_FFT(N,tab);

 cout<<"frequency Hz"<<endl;

    for(int i=0;i<2;i++)
    {
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab[j+i*N]*100)/100<<"  ";
    }
    cout<<endl;
    }
    cout<<endl;

 scaling_re_and_im(N,tab);

 cout<<"scaled"<<endl;

    for(int i=0;i<2;i++)
    {
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab[j+i*N]*100)/100<<"  ";
    }
    cout<<endl;
    }
    cout<<endl;

 fun_inverse_bits(N,tab);
 //not inverse but makes it inverse
 fast_fourier_transform_FFT(N,tab);//not inverse but makes it inverse
 scaling_re_and_im(N,tab);

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

  cout<<"inverse"<<endl;
    for(int i=0;i<2;i++)
    {
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab[j+i*N]*1000)/1000<<"  ";
    }
    cout<<endl;
    }
    cout<<endl;

    system("pause");
    return 0;
}



void fun_inverse_bits(int N,double *tablica2)  // wielkooa N zawsze potega 2: N^2
{
    double (*(tablica1)[2]);
    tablica1[0]=&tablica2[0];
    tablica1[1]=&tablica2[0+N];

  ////////////////////
//////////////////////
}    // koniec void fun_inverse_bits
/////////////////////////////////////////////////////////////

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

void fast_fourier_transform_FFT(int N,double *tablica2)//+void fun_inverse_bits
{

    double (*(tablica1)[2]);
    tablica1[0]=&tablica2[0];
    tablica1[1]=&tablica2[0+N];

 /////////////////
////////////////////
}    //koniec funkcji: void fast_fourier_transform_FFT
/////////////////////////////////////////////////////////////

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

void inverse_fast_fourier_transform_FFT(int N,double *tablica2)//+void fun_inverse_bits????
{


    double (*(tablica1)[2]);
    tablica1[0]=&tablica2[0];
    tablica1[1]=&tablica2[0+N];


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

}    //koniec funkcji: void inverse_fast_fourier_transform_FFT
/////////////////////////////////////////////////////////////










double scaling_re_and_im(int N,double *tablica)
{
    double (*(tablica1)[2]);
    tablica1[0]=&tablica[0];
    tablica1[1]=&tablica[0+N];

     for (int i=0;i<N;i++)
    {
        tablica1[0][i]=tablica1[0][i]+tablica1[1][i];
        tablica1[1][i]=0;
    }

}

quicker version of fourier transform radix-4 without using inside object like complex numbers only tables

  //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/
    //http://www.cmlab.csie.ntu.edu.tw/cml/dsp/training/coding/transform/fft.html
    //http://dsp.stackexchange.com/questions/3481/radix-4-fft-implementation

    //https://community.arm.com/graphics/b/blog/posts/speeding-up-fast-fourier-transform-mixed-radix-on-mobile-arm-mali-gpu-by-means-of-opencl---part-1
    //book: "Cyfrowe przetwarzanie sygnalow" - Tomasz Zielinski it has quick version of radix-2 because it calculates less sin() and cos()




    //author marcin matysek (r)ewertyn.PL
    //quicker version of fourier transform radix-4 without using inside object like complex numbers only tables
    //inverse fast fourier transform iFFT radix-4 combinations witch shift for N= points algorithm c++ source code implementation
    //szybka odwrotna transformacja fouriera iFFT radix-4 z przesunieciem dla N= punktów algorytm c++ kod zródlowy implementacja

   void fun_fourier_transform_FFT_radix_4(int N,std::complex<double> tab[],std::complex<double> w10,std::complex<double> w11,std::complex<double> w20,std::complex<double> w21)
{
    //double *(tab)[1]={};
    //tab[0]=&tab2[0].real();
    const double pi=3.141592653589793238462;
    //std::complex<double>  w111[1]={{1,0}};
    //std::complex<double>  w2[1]={{1,0}};
    //std::complex<double>  w3[1]={{1,0}};
    //std::complex<double>  w4[1]={{1,0}};
    //std::complex<double>  w5[1]={{1,0}};

    double w1[4][2]={};
    double tab1[4][2]={};
    double tmp1[4][2]={};
    double tmp2[4][2]={};
    double tmp3[4][2]={};


    //std::complex<double> tmp01,tmp02,tmp03,tmp04,tmp6,tmp7,tmp8,tmp9;
    //std::complex<double> tmp10,tmp20,tmp30,tmp40;
    int nr1=0,nr2=0,nr3=0,nr4=0,nr5=0;
    int nr10=0,nr20=0,nr31=0,nr32=0,nr33=0;
    int nb_stages=0,rx4=0;
    double tmp5,tmp100,tmp200;
    int increment=0;
    int nn=0;

    //w5[0].real()=0;
    //w5[0].imag()=1;//???
    tab1[0][0]=w10.real();
    tab1[0][1]=w10.imag();
    tab1[1][0]=w11.real();
    tab1[1][1]=w11.imag();
    tab1[2][0]=w20.real();
    tab1[2][1]=w20.imag();
    tab1[3][0]=w21.real();
    tab1[3][1]=w21.imag();


   // cout<<tab1[2][0]<<" "<<tab1[2][1]<<endl;
  //  cout<<tab1[3][0]<<" "<<tab1[3][1]<<endl;system("pause");
  /*
    int M=2;
    int **tab9 = new int *[N];//alokacja pamieci
    for ( int i = 0; i < N; ++i )
    {
        tab9[i] = new int [M]; //alokacja pamieci

    }
*/
    double fi2=0;//always 0

    tmp5=2*pi/(N/1);
    rx4=4;//radix-4
    nn=N;
    for(int i=0;i<100;i++)
    {

        nn=(float)nn/(float)rx4;
         if(nn>=1.0)
         {
         nb_stages++;
         }
    }
//stage 1
          w1[0][0]=cos(0*tmp5+fi2);
          w1[0][1]=sin(0*tmp5+fi2);
          w1[1][0]=cos(0*tmp5+fi2);
          w1[1][1]=sin(0*tmp5+fi2);
          w1[2][0]=cos(0*tmp5+fi2);
          w1[2][1]=sin(0*tmp5+fi2);
          w1[3][0]=cos(0*tmp5+fi2);
          w1[3][1]=sin(0*tmp5+fi2);
        for(int i=0;i<(N/rx4);i++)
        {
            tmp1[0][0]=w1[0][0]*tab[i].real()-w1[0][1]*tab[i].imag();
            tmp1[0][1]=w1[0][0]*tab[i].imag()+w1[0][1]*tab[i].real();
            tmp1[1][0]=w1[1][0]*tab[i+N/4].real()-w1[1][1]*tab[i+N/4].imag();
            tmp1[1][1]=w1[1][0]*tab[i+N/4].imag()+w1[1][1]*tab[i+N/4].real();
            tmp1[2][0]=w1[2][0]*tab[i+N/2].real()-w1[2][1]*tab[i+N/2].imag();
            tmp1[2][1]=w1[2][0]*tab[i+N/2].imag()+w1[2][1]*tab[i+N/2].real();
            tmp1[3][0]=w1[3][0]*tab[i+3*N/4].real()-w1[3][1]*tab[i+3*N/4].imag();
            tmp1[3][1]=w1[3][0]*tab[i+3*N/4].imag()+w1[3][1]*tab[i+3*N/4].real();

            tmp2[0][0]=tab1[0][0]*(tmp1[0][0]-tmp1[2][0])-tab1[0][1]*(tmp1[0][1]-tmp1[2][1]);
            tmp2[0][1]=tab1[0][0]*(tmp1[0][1]-tmp1[2][1])+tab1[0][1]*(tmp1[0][0]-tmp1[2][0]);
            tmp2[1][0]=tab1[1][0]*(tmp1[1][0]-tmp1[3][0])-tab1[1][1]*(tmp1[1][1]-tmp1[3][1]);
            tmp2[1][1]=tab1[1][0]*(tmp1[1][1]-tmp1[3][1])+tab1[1][1]*(tmp1[1][0]-tmp1[3][0]);
            tmp2[2][0]=tmp1[0][0]+tmp1[2][0];
            tmp2[2][1]=tmp1[0][1]+tmp1[2][1];
            tmp2[3][0]=tmp1[1][0]+tmp1[3][0];
            tmp2[3][1]=tmp1[1][1]+tmp1[3][1];

            //radix-4
            tmp3[0][0]=tab1[0][0]*(tmp2[2][0]+tmp2[3][0])-tab1[0][1]*(tmp2[2][1]+tmp2[3][1]);
            tmp3[0][1]=tab1[0][0]*(tmp2[2][1]+tmp2[3][1])+tab1[0][1]*(tmp2[2][0]+tmp2[3][0]);
            tmp3[1][0]=tmp2[0][0]+tmp2[1][0];
            tmp3[1][1]=tmp2[0][1]+tmp2[1][1];
            tmp3[2][0]=tab1[0][0]*(tmp2[2][0]-tmp2[3][0])-tab1[0][1]*(tmp2[2][1]-tmp2[3][1]);
            tmp3[2][1]=tab1[0][0]*(tmp2[2][1]-tmp2[3][1])+tab1[0][1]*(tmp2[2][0]-tmp2[3][0]);
            tmp3[3][0]=tmp2[0][0]-tmp2[1][0];
            tmp3[3][1]=tmp2[0][1]-tmp2[1][1];

            tab[i].real()=tmp3[0][0];
            tab[i].imag()=tmp3[0][1];
            tab[i+N/4].real()=tmp3[1][0];
            tab[i+N/4].imag()=tmp3[1][1];
            tab[i+N/2].real()=tmp3[2][0];
            tab[i+N/2].imag()=tmp3[2][1];
            tab[i+3*N/4].real()=tmp3[3][0];
            tab[i+3*N/4].imag()=tmp3[3][1];

            /*
          tmp01=w111[0]*tab[i+0];
          tmp02=w2[0]*tab[i+N/4];
          tmp03=w3[0]*tab[i+N/2];
          tmp04=w4[0]*tab[i+3*N/4];
            tmp6=w10*(tmp01-tmp03);
            tmp7=w11*(tmp02-tmp04);
            tmp8=tmp01+tmp03;
            tmp9=tmp02+tmp04;
            //radix-4
           tmp10=(w10*(tmp8+tmp9));
           tmp20=(tmp6+tmp7);
           tmp30=(w10*(tmp8-tmp9));
           tmp40=(tmp6-tmp7);

          tab[i]       =tmp10;
          tab[i+N/4]   =tmp20;
          tab[i+N/2]   =tmp30;
          tab[i+3*N/4] =tmp40;
          */
        }
///////////////////////////////////////////


    increment=0;
    int flag2=0;
    for(int stg=1;stg<nb_stages;stg++)
    {
          //  cout<<"..if(flag2%2==0).."<<endl;
          //  system("pause");
         nr1=N/pow(rx4,0+increment);
        nr2=N/pow(rx4,1+increment);
        nr3=N/pow(rx4,2+increment);
        nr4=pow(rx4,(nb_stages-2-increment));
        nr5=pow(rx4,0+increment);
        nr31=1*nr3;
        nr32=2*nr3;
        nr33=3*nr3;
        for(int j=0;j<nr4;j++)
        {
            for(int i=0;i<rx4;i++)
            {
                nr20=nr2*i;
                tmp100=nr5*i*(j+nr4);
                tmp200=nr5*i*j;

                w1[0][0]=cos(tmp200*tmp5+fi2);
                w1[0][1]=-sin(tmp200*tmp5+fi2);
                w1[1][0]=cos(tmp100*tmp5+fi2);
                w1[1][1]=-sin(tmp100*tmp5+fi2);
                w1[2][0]=cos((2*tmp100-  tmp200)*tmp5+fi2);
                w1[2][1]=-sin((2*tmp100-  tmp200)*tmp5+fi2);
                w1[3][0]=cos((3*tmp100-2*tmp200)*tmp5+fi2);
                w1[3][1]=-sin((3*tmp100-2*tmp200)*tmp5+fi2);
                for(int m=0;m<nr5;m++)
                {
                    nr10=nr1*m;

                    tmp1[0][0]=w1[0][0]*tab[     nr20+nr10+j].real()-w1[0][1]*tab[     nr20+nr10+j].imag();
                    tmp1[0][1]=w1[0][0]*tab[     nr20+nr10+j].imag()+w1[0][1]*tab[     nr20+nr10+j].real();
                    tmp1[1][0]=w1[1][0]*tab[nr31+nr20+nr10+j].real()-w1[1][1]*tab[nr31+nr20+nr10+j].imag();
                    tmp1[1][1]=w1[1][0]*tab[nr31+nr20+nr10+j].imag()+w1[1][1]*tab[nr31+nr20+nr10+j].real();
                    tmp1[2][0]=w1[2][0]*tab[nr32+nr20+nr10+j].real()-w1[2][1]*tab[nr32+nr20+nr10+j].imag();
                    tmp1[2][1]=w1[2][0]*tab[nr32+nr20+nr10+j].imag()+w1[2][1]*tab[nr32+nr20+nr10+j].real();
                    tmp1[3][0]=w1[3][0]*tab[nr33+nr20+nr10+j].real()-w1[3][1]*tab[nr33+nr20+nr10+j].imag();
                    tmp1[3][1]=w1[3][0]*tab[nr33+nr20+nr10+j].imag()+w1[3][1]*tab[nr33+nr20+nr10+j].real();

                    //tab1[2][0]=1
                    //tab1[2][1]=0
                    //tmp2[0][0]=tab1[2][0]*(tmp1[0][0]-tmp1[2][0])-tab1[2][1]*(tmp1[0][1]-tmp1[2][1]);
                    //tmp2[0][1]=tab1[2][0]*(tmp1[0][1]-tmp1[2][1])+tab1[2][1]*(tmp1[0][0]-tmp1[2][0]);
                    tmp2[0][0]=tmp1[0][0]-tmp1[2][0];
                    tmp2[0][1]=tmp1[0][1]-tmp1[2][1];
                    tmp2[1][0]=tab1[3][0]*(tmp1[1][0]-tmp1[3][0])-tab1[3][1]*(tmp1[1][1]-tmp1[3][1]);
                    tmp2[1][1]=tab1[3][0]*(tmp1[1][1]-tmp1[3][1])+tab1[3][1]*(tmp1[1][0]-tmp1[3][0]);
                    tmp2[2][0]=tmp1[0][0]+tmp1[2][0];
                    tmp2[2][1]=tmp1[0][1]+tmp1[2][1];
                    tmp2[3][0]=tmp1[1][0]+tmp1[3][0];
                    tmp2[3][1]=tmp1[1][1]+tmp1[3][1];

                    //radix-4
                    //tab1[2][0]=1
                    //tab1[2][1]=0
                    //tmp3[0][0]=tab1[2][0]*(tmp2[2][0]+tmp2[3][0])-tab1[2][1]*(tmp2[2][1]+tmp2[3][1]);
                    //tmp3[0][1]=tab1[2][0]*(tmp2[2][1]+tmp2[3][1])+tab1[2][1]*(tmp2[2][0]+tmp2[3][0]);
                    tmp3[0][0]=tmp2[2][0]+tmp2[3][0];
                    tmp3[0][1]=tmp2[2][1]+tmp2[3][1];
                    tmp3[1][0]=tmp2[0][0]+tmp2[1][0];
                    tmp3[1][1]=tmp2[0][1]+tmp2[1][1];
                    //tmp3[2][0]=tab1[2][0]*(tmp2[2][0]-tmp2[3][0])-tab1[2][1]*(tmp2[2][1]-tmp2[3][1]);
                    //tmp3[2][1]=tab1[2][0]*(tmp2[2][1]-tmp2[3][1])+tab1[2][1]*(tmp2[2][0]-tmp2[3][0]);
                    tmp3[2][0]=tmp2[2][0]-tmp2[3][0];
                    tmp3[2][1]=tmp2[2][1]-tmp2[3][1];
                    tmp3[3][0]=tmp2[0][0]-tmp2[1][0];
                    tmp3[3][1]=tmp2[0][1]-tmp2[1][1];

                    tab[     nr20+nr10+j].real()=tmp3[0][0];
                    tab[     nr20+nr10+j].imag()=tmp3[0][1];
                    tab[nr31+nr20+nr10+j].real()=tmp3[1][0];
                    tab[nr31+nr20+nr10+j].imag()=tmp3[1][1];
                    tab[nr32+nr20+nr10+j].real()=tmp3[2][0];
                    tab[nr32+nr20+nr10+j].imag()=tmp3[2][1];
                    tab[nr33+nr20+nr10+j].real()=tmp3[3][0];
                    tab[nr33+nr20+nr10+j].imag()=tmp3[3][1];


                    //complex numbers version:
                   /*
                    tmp01=w111[0]*tab[     nr20+nr10+j];
                    tmp02=w2[0]*tab[nr31+nr20+nr10+j];
                    tmp03=w3[0]*tab[nr32+nr20+nr10+j];
                    tmp04=w4[0]*tab[nr33+nr20+nr10+j];

                    tmp6=w20*(tmp01-tmp03);
                    tmp7=w21*(tmp02-tmp04);
                    tmp8=tmp01+tmp03;
                    tmp9=tmp02+tmp04;
                   //radix-4
                   tmp10=(w20*(tmp8+tmp9));
                   tmp20=(tmp6+tmp7);
                   tmp30=(w20*(tmp8-tmp9));
                   tmp40=(tmp6-tmp7);

                    tab[     nr20+nr10+j]=tmp10;
                    tab[nr31+nr20+nr10+j]=tmp20;
                    tab[nr32+nr20+nr10+j]=tmp30;
                    tab[nr33+nr20+nr10+j]=tmp40;
                    */
                }
            }
        }
        increment++;
    }
    /*
    //zniszcz tablice 2D
    for ( int i(0); i < N; ++i )
    {
        delete [] tab9[i]; //uwolnienie pamieci
    }
    delete [] tab9; //uwolnienie pamieci
    tab9 = NULL;
    */
    ////////////////////////////////////////////////////
    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_radix_4_basis_for_stage1(std::complex<double> &w0,std::complex<double> &w1)
{
    const double pi=3.141592653589793238462;
    //std::complex<double>    w0,w1,w2,w3,w4,w6,w9;
   
    double tmp5=2*pi/4;
   
   
    double fi2=fi;
   
   
    w0.real()=round(cos(0*tmp5+fi2)*10000000)/10000000;
    w0.imag()=round(-sin(0*tmp5+fi2)*10000000)/10000000;
    w1.real()=round(cos(1*tmp5+fi2)*10000000)/10000000;
    w1.imag()=round(-sin(1*tmp5+fi2)*10000000)/10000000;
    w0=w0;
    w1=w1;
    w0.real()=w0.real();
    w0.imag()=w0.imag();
    w1.real()=w1.real();
    w1.imag()=w1.imag();

    //w2.real()=round(cos(2*tmp5)*1000)/1000;//==-w0
    //w2.imag()=round(-sin(2*tmp5)*1000)/1000;//==-w0
    //w3.real()=round(cos(3*tmp5)*1000)/1000;//==-w1
    //w3.imag()=round(-sin(3*tmp5)*1000)/1000;//==-w1
    //w4.real()=round(cos(4*tmp5)*1000)/1000;//==+w0
    //w4.imag()=round(-sin(4*tmp5)*1000)/1000;//==+w0
    //w6.real()=round(cos(6*tmp5)*1000)/1000;//==-w0
    //w6.imag()=round(-sin(6*tmp5)*1000)/1000;//==-w0
    //w9.real()=round(cos(9*tmp5)*1000)/1000;//==+w1
    //w9.imag()=round(-sin(9*tmp5)*1000)/1000;//==+w1

   //cout<<" "<<w0<<" "<<w0<<" "<<w0<<" "<<w0<<endl;
   //cout<<" "<<w0<<" "<<w1<<" "<<w2<<" "<<w3<<endl;
   //cout<<" "<<w0<<" "<<w2<<" "<<w4<<" "<<w6<<endl;
   //cout<<" "<<w0<<" "<<w3<<" "<<w6<<" "<<w9<<endl;
   //system("pause");
}
//////////////////
void fun_radix_4_basis_for_stage_rest(std::complex<double> &w0,std::complex<double> &w1)
{
    const double pi=3.141592653589793238462;
    //std::complex<double>    w0,w1,w2,w3,w4,w6,w9;
   
    double tmp5=2*pi/4;
   
    double fi2=0;//always 0
   
    w0.real()=round(cos(0*tmp5+fi2)*10000000)/10000000;
    w0.imag()=round(-sin(0*tmp5+fi2)*10000000)/10000000;
    w1.real()=round(cos(1*tmp5+fi2)*10000000)/10000000;
    w1.imag()=round(-sin(1*tmp5+fi2)*10000000)/10000000;
    w0=w0;
    w1=w1;
    w0.real()=w0.real();
    w0.imag()=w0.imag();
    w1.real()=w1.real();
    w1.imag()=w1.imag();

    //w2.real()=round(cos(2*tmp5)*1000)/1000;//==-w0
    //w2.imag()=round(-sin(2*tmp5)*1000)/1000;//==-w0
    //w3.real()=round(cos(3*tmp5)*1000)/1000;//==-w1
    //w3.imag()=round(-sin(3*tmp5)*1000)/1000;//==-w1
    //w4.real()=round(cos(4*tmp5)*1000)/1000;//==+w0
    //w4.imag()=round(-sin(4*tmp5)*1000)/1000;//==+w0
    //w6.real()=round(cos(6*tmp5)*1000)/1000;//==-w0
    //w6.imag()=round(-sin(6*tmp5)*1000)/1000;//==-w0
    //w9.real()=round(cos(9*tmp5)*1000)/1000;//==+w1
    //w9.imag()=round(-sin(9*tmp5)*1000)/1000;//==+w1

   //cout<<" "<<w0<<" "<<w0<<" "<<w0<<" "<<w0<<endl;
   //cout<<" "<<w0<<" "<<w1<<" "<<w2<<" "<<w3<<endl;
   //cout<<" "<<w0<<" "<<w2<<" "<<w4<<" "<<w6<<endl;
   //cout<<" "<<w0<<" "<<w3<<" "<<w6<<" "<<w9<<endl;
   //system("pause");
}

inverse fast fourier transform iFFT radix-4 combinations witch shift for N= points 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/
    //http://www.cmlab.csie.ntu.edu.tw/cml/dsp/training/coding/transform/fft.html
    //http://dsp.stackexchange.com/questions/3481/radix-4-fft-implementation

    //https://community.arm.com/graphics/b/blog/posts/speeding-up-fast-fourier-transform-mixed-radix-on-mobile-arm-mali-gpu-by-means-of-opencl---part-1
    //book: "Cyfrowe przetwarzanie sygnalow" - Tomasz Zielinski it has quick version of radix-2 because it calculates less sin() and cos()




    //author marcin matysek (r)ewertyn.PL

    //inverse fast fourier transform iFFT radix-4 combinations witch shift for N= points algorithm c++ source code implementation
    //szybka odwrotna transformacja fouriera iFFT radix-4 z przesunięciem dla N= punktów algorytm c++ kod zródlowy implementacja

    #include <iostream>
    #include "conio.h"
    #include <stdlib.h>
    #include <math.h>
    #include <cmath>
    #include <time.h>
    #include <complex>
    #include <fstream>

    using namespace std;




           //  shift //
        //////////////////////
        double fi=0.23;///////
        //////////////////////






    //complex number method:
    void fun_inverse_bits_radix_4(int N,std::complex<double> tab[]);
    void fun_fourier_transform_FFT_radix_4_N_4096(int N,std::complex<double> tab[],std::complex<double> w10,std::complex<double> w11,std::complex<double> w20,std::complex<double> w21);
    void fun_inverse_fourier_transform_FFT_radix_4_N_4096(int N,std::complex<double> tab[],std::complex<double> w10,std::complex<double> w11,std::complex<double> w20,std::complex<double> w21);
    void fun_fourier_transform_DFT_full_complex(int N,std::complex<double> tab[]);

    void fun_radix_4_basis_for_stage1(std::complex<double> &w0,std::complex<double> &w1);
    void fun_radix_4_basis_for_stage_rest(std::complex<double> &w0,std::complex<double> &w1);
    void fun_inverse_radix_4_basis_for_stage1(std::complex<double> &w0,std::complex<double> &w1);
    void fun_inverse_radix_4_basis_for_stage_rest(std::complex<double> &w0,std::complex<double> &w1);



    int nb100=0,nb101=0,nb102=0,nb103=0;

    static double diffclock(clock_t clock2,clock_t clock1)
    {
        double diffticks=clock1-clock2;
        double diffms=(diffticks)/(CLOCKS_PER_SEC/1000);
        return diffms;
    }
int main()
{
    const double pi=3.141592653589793238462;
    int N,nn=0;
    //if N==period of signal in table tab[] then resolution = 1 Hz
    std::complex<double> w0;
    std::complex<double> w1;
    std::complex<double> w2;
    std::complex<double> w3;
    std::complex<double> w10;
    std::complex<double> w11;
    std::complex<double> w12;
    std::complex<double> w13;
    int maks=0;



    N=256;;//need to change size of tab2 in function: fun_fourier_transform_DFT_full_complex




    std::complex<double> tab2[4096]={};
    std::complex<double> tab3[4096]={};
    int flag=0;



    double time2;
    double zmienna=0;



    for(int i=0;i<N;i++)
    {
        tab2[i].real()=i+1;
        tab2[i].imag()=-N-i-1;
        tab3[i].real()=i+1;
        tab3[i].imag()=-N-i-1;
    }
    for(int i=0;i<N;i++)
    {
       // tab3[i].real()=sin(i*2*pi/N+pi/2);
       // tab3[i].real()=tab3[i].real()+sin(i*2*2*pi/N+pi/2);
       // tab3[i].real()=tab3[i].real()+sin(i*3*2*pi/N+pi/2);
       // tab3[i].imag()=0;
    }
    cout<<"signal 1="<<endl;
    system("pause");
    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<<"signal 2="<<endl;
     system("pause");
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab3[j].real()*1000)/1000<<"  ";
    }
    cout<<endl;
    cout<<endl;
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab3[j].imag()*1000)/1000<<"  ";
    }
    cout<<endl;
    cout<<endl;
    cout<<"calculations"<<endl;
     system("pause");
    clock_t start = clock();
   fun_fourier_transform_DFT_full_complex(N,tab3);
    //////////////////////////////////////////////////////////





    for(int i100=0;i100<2;i100++)
    {
        if(i100==0){nb100=0;}
        else if(i100==1){nb100=1;}

    for(int i101=0;i101<2;i101++)
    {
        if(i101==0){nb101=0;}
        else if(i101==1){nb101=1;}

    for(int i102=0;i102<2;i102++)
    {
        if(i102==0){nb102=0;}
        else if(i102==1){nb102=1;}

    for(int i103=0;i103<2;i103++)
    {
        if(i103==0){nb103=0;}
        else if(i103==1){nb103=1;}



    fun_radix_4_basis_for_stage1(w0,w1);
    fun_radix_4_basis_for_stage_rest(w2,w3);
    for(int i=0;i<N;i++)
    {
        tab2[i].real()=i+1;
        tab2[i].imag()=-N-i-1;
    }
    fun_fourier_transform_FFT_radix_4_N_4096(N,tab2,w0,w1,w2,w3);
    fun_inverse_bits_radix_4(N,tab2);
    ///////////////////////////////////////////////////////////

    fun_inverse_radix_4_basis_for_stage1(w10,w11);
    fun_inverse_radix_4_basis_for_stage_rest(w12,w13);
    /////////////////////////////////////////////////////////////////
   // fun_inverse_fourier_transform_FFT_radix_4_N_4096(N,tab2,w10,w11,w12,w13);
   // fun_inverse_bits_radix_4(N,tab2);
    ////////////////////////////////////////////////////////////////

    time2=diffclock( start, clock() );
/*
    cout<<endl;
    cout<<endl;
    cout<<"frequency Hz radix-4 real()"<<endl;
    system("pause");
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab2[j].real()*1000)/1000<<"  ";
    }
    cout<<endl;
    cout<<endl;
    cout<<"frequency Hz radix-4 imag()"<<endl;
     system("pause");
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab2[j].imag()*1000)/1000<<"  ";
    }
    cout<<endl;
    cout<<endl;

    cout<<"frequency Hz DFT real()"<<endl;
     system("pause");
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab3[j].real()*1000)/1000<<"  ";
    }
    cout<<endl;
    cout<<endl;
    cout<<"frequency Hz DFT imag()"<<endl;
     system("pause");
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab3[j].imag()*1000)/1000<<"  ";
    }
    cout<<endl;
    cout<<endl;
    cout<<"if radix-4 == DFT tab2[j].real(): "<<endl;system("pause");
    */
    flag=0;

       for(int j=0;j<N;j++)
        {
          if(fabs(tab3[j].real()-tab2[j].real())<=0.03&&fabs(tab3[j].imag()-tab2[j].imag()<=0.03))
          {
            flag++;
            //cout.precision(4);
            //cout<<j<<" "<<round(tab2[j].real()*1000)/1000<<" ok "<<flaga;system("pause");
          }
            else {
                //cout<<j<<" "<<-1<<" .not equal.. ";
                //system("pause");
            }
        }



        if(flag>maks){maks=flag;}
        if(flag>=N)
        {
            nn++;
            cout<<endl<<"found good combination"<<endl<<nn<<" "<<flag;
            cout<<" "<<nb100<<" "<<nb101<<" "<<nb102<<" "<<nb103<<" "<<endl;
            system("pause");
        }

    }}}}
    cout<<"maks= "<<maks<<endl;
    system("pause");
    cout<<endl;

    /*
    /////////////////////////////////////////////////////////////////
    fun_inverse_fourier_transform_FFT_radix_4_N_4096(N,tab2,w0,w1,w2,w3);
    fun_inverse_bits_radix_4(N,tab2);
    ////////////////////////////////////////////////////////////////


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

    system("pause");
    */

    return 0;
}
//////////////////////////////////////////////
void fun_inverse_bits_radix_4(int N,std::complex<double> tab[])
{

//code by Sidney Burrus
    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;
    }
}
///////////////////////////////////////////////////
void fun_fourier_transform_DFT_full_complex(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double> tab2[4096*4]={};    // tab2[]==N
    std::complex<double>  w[1]={{1,1}};


double zmienna1=2*pi/(float)N;
double fi2=fi;

for (int i=0;i<N;i++)
{
    for(int j=0;j<N;j++)
    {
          //complex number method:
          w[0].real()=cos(i*j*zmienna1+fi2);
          w[0].imag()=(-sin(i*j*zmienna1+fi2));
          tab2[i]=tab2[i]+tab[j]*w[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_4096(int N,std::complex<double> tab[],std::complex<double> w10,std::complex<double> w11,std::complex<double> w20,std::complex<double> w21)
{
    const double pi=3.141592653589793238462;
    std::complex<double>  w1[1]={{1,0}};
    std::complex<double>  w2[1]={{1,0}};
    std::complex<double>  w3[1]={{1,0}};
    std::complex<double>  w4[1]={{1,0}};
    std::complex<double>  w5[1]={{1,0}};
    std::complex<double> tmp1,tmp2,tmp3,tmp4,tmp6,tmp7;
    std::complex<double> tmp10,tmp20,tmp30,tmp40;
    int nr1=0,nr2=0,nr3=0,nr4=0,nr5=0;
    int nr10=0,nr20=0,nr31=0,nr32=0,nr33=0;
    int nb_stages=0,rx4=0;
    double tmp5,tmp100,tmp200;
    int increment=0;
    int nn=0;

    //w5[0].real()=0;
    //w5[0].imag()=-1;
    double fi2=0;
    if(nb100==0)
    {
      fi2=0;
    }
    else if(nb100==1)
    {
      fi2=fi;
    }


    tmp5=2*pi/(N/1);
    rx4=4;//radix-4
    nn=N;
    for(int i=0;i<100;i++)
    {

        nn=(float)nn/(float)rx4;
         if(nn>=1.0)
         {
         nb_stages++;
         }
    }
//stage 1
          w1[0].real()=cos(0*tmp5+fi2);
          w1[0].imag()=-sin(0*tmp5+fi2);
          w2[0].real()=cos(0*tmp5+fi2);
          w2[0].imag()=-sin(0*tmp5+fi2);
          w3[0].real()=cos(0*tmp5+fi2);
          w3[0].imag()=-sin(0*tmp5+fi2);
          w4[0].real()=cos(0*tmp5+fi2);
          w4[0].imag()=-sin(0*tmp5+fi2);
        for(int i=0;i<(N/rx4);i++)
        {
          tmp1=w1[0]*tab[i+0];
          tmp2=w2[0]*tab[i+N/4];
          tmp3=w3[0]*tab[i+N/2];
          tmp4=w4[0]*tab[i+3*N/4];
                    if(nb101==0)
                    {
                    tmp6=w10*(tmp1-tmp3);
                    tmp7=w11*(tmp2-tmp4);
                    //radix-4
                   tmp10=(w10*(tmp1+tmp2+tmp3+tmp4));
                   tmp20=(tmp6+tmp7);
                   tmp30=(w10*(tmp1-tmp2+tmp3-tmp4));
                   tmp40=(tmp6-tmp7);
                    }
                    else if(nb101==1)
                    {
                    tmp6=w20*(tmp1-tmp3);
                    tmp7=w21*(tmp2-tmp4);
                    //radix-4
                   tmp10=(w20*(tmp1+tmp2+tmp3+tmp4));
                   tmp20=(tmp6+tmp7);
                   tmp30=(w20*(tmp1-tmp2+tmp3-tmp4));
                   tmp40=(tmp6-tmp7);
                    }
          tab[i]       =tmp10;
          tab[i+N/4]   =tmp20;
          tab[i+N/2]   =tmp30;
          tab[i+3*N/4] =tmp40;
        }
///////////////////////////////////////////

    if(nb102==0)
    {
      fi2=0;
    }
    else if(nb102==1)
    {
      fi2=fi;
    }

    increment=0;
    int flag2=0;
    for(int stg=1;stg<nb_stages;stg++)
    {

          //  cout<<"..if(flag2%2==0).."<<endl;
          //  system("pause");
         nr1=N/pow(rx4,0+increment);
        nr2=N/pow(rx4,1+increment);
        nr3=N/pow(rx4,2+increment);
        nr4=pow(rx4,(nb_stages-2-increment));
        nr5=pow(rx4,0+increment);
        nr31=1*nr3;
        nr32=2*nr3;
        nr33=3*nr3;
        for(int j=0;j<nr4;j++)
        {
            for(int i=0;i<rx4;i++)
            {
                nr20=nr2*i;
                tmp100=nr5*i*(j+nr4);
                tmp200=nr5*i*j;
                w1[0].real()= cos(tmp200*tmp5+fi2);
                w1[0].imag()=-sin(tmp200*tmp5+fi2);
                w2[0].real()= cos(tmp100*tmp5+fi2);
                w2[0].imag()=-sin(tmp100*tmp5+fi2);
                w3[0].real()= cos((2*tmp100-  tmp200)*tmp5+fi2);
                w3[0].imag()=-sin((2*tmp100-  tmp200)*tmp5+fi2);
                w4[0].real()= cos((3*tmp100-2*tmp200)*tmp5+fi2);
                w4[0].imag()=-sin((3*tmp100-2*tmp200)*tmp5+fi2);
                for(int m=0;m<nr5;m++)
                {
                    nr10=nr1*m;
                    tmp1=w1[0]*tab[     nr20+nr10+j];
                    tmp2=w2[0]*tab[nr31+nr20+nr10+j];
                    tmp3=w3[0]*tab[nr32+nr20+nr10+j];
                    tmp4=w4[0]*tab[nr33+nr20+nr10+j];
                    if(nb103==0)
                    {
                    tmp6=w10*(tmp1-tmp3);
                    tmp7=w11*(tmp2-tmp4);
                    //radix-4
                   tmp10=(w10*(tmp1+tmp2+tmp3+tmp4));
                   tmp20=(tmp6+tmp7);
                   tmp30=(w10*(tmp1-tmp2+tmp3-tmp4));
                   tmp40=(tmp6-tmp7);
                    }
                    else if(nb103==1)
                    {
                    tmp6=w20*(tmp1-tmp3);
                    tmp7=w21*(tmp2-tmp4);
                    //radix-4
                   tmp10=(w20*(tmp1+tmp2+tmp3+tmp4));
                   tmp20=(tmp6+tmp7);
                   tmp30=(w20*(tmp1-tmp2+tmp3-tmp4));
                   tmp40=(tmp6-tmp7);
                    }

                    tab[     nr20+nr10+j]=tmp10;
                    tab[nr31+nr20+nr10+j]=tmp20;
                    tab[nr32+nr20+nr10+j]=tmp30;
                    tab[nr33+nr20+nr10+j]=tmp40;
                }
            }
        }
        increment++;
    }
    ////////////////////////////////////////////////////
    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_4096(int N,std::complex<double> tab[],std::complex<double> w10,std::complex<double> w11,std::complex<double> w20,std::complex<double> w21)
{
     const double pi=3.141592653589793238462;
    std::complex<double>  w1[1]={{1,0}};
    std::complex<double>  w2[1]={{1,0}};
    std::complex<double>  w3[1]={{1,0}};
    std::complex<double>  w4[1]={{1,0}};
    std::complex<double>  w5[1]={{1,0}};
    std::complex<double> tmp1,tmp2,tmp3,tmp4,tmp6,tmp7;
    std::complex<double> tmp10,tmp20,tmp30,tmp40;
    int nr1=0,nr2=0,nr3=0,nr4=0,nr5=0;
    int nr10=0,nr20=0,nr31=0,nr32=0,nr33=0;
    int nb_stages=0,rx4=0;
    double tmp5,tmp100,tmp200;
    int increment=0;
    int nn=0;

    //w5[0].real()=0;
    //w5[0].imag()=1;//???

    double fi2=0;
    if(nb100==0)
    {
      fi2=0;
    }
    else if(nb100==1)
    {
      fi2=fi;
    }

    tmp5=2*pi/(N/1);
    rx4=4;//radix-4
    nn=N;
    for(int i=0;i<100;i++)
    {

        nn=(float)nn/(float)rx4;
         if(nn>=1.0)
         {
         nb_stages++;
         }
    }
//stage 1
          w1[0].real()=cos(0*tmp5+fi2);
          w1[0].imag()=sin(0*tmp5+fi2);
          w2[0].real()=cos(0*tmp5+fi2);
          w2[0].imag()=sin(0*tmp5+fi2);
          w3[0].real()=cos(0*tmp5+fi2);
          w3[0].imag()=sin(0*tmp5+fi2);
          w4[0].real()=cos(0*tmp5+fi2);
          w4[0].imag()=sin(0*tmp5+fi2);
        for(int i=0;i<(N/rx4);i++)
        {
          tmp1=w1[0]*tab[i+0];
          tmp2=w2[0]*tab[i+N/4];
          tmp3=w3[0]*tab[i+N/2];
          tmp4=w4[0]*tab[i+3*N/4];
                    if(nb101==0)
                    {
                    tmp6=w10*(tmp1-tmp3);
                    tmp7=w11*(tmp2-tmp4);
                    //radix-4
                   tmp10=(w10*(tmp1+tmp2+tmp3+tmp4));
                   tmp20=(tmp6+tmp7);
                   tmp30=(w10*(tmp1-tmp2+tmp3-tmp4));
                   tmp40=(tmp6-tmp7);
                    }
                    else if(nb101==1)
                    {
                    tmp6=w20*(tmp1-tmp3);
                    tmp7=w21*(tmp2-tmp4);
                    //radix-4
                   tmp10=(w20*(tmp1+tmp2+tmp3+tmp4));
                   tmp20=(tmp6+tmp7);
                   tmp30=(w20*(tmp1-tmp2+tmp3-tmp4));
                   tmp40=(tmp6-tmp7);
                    }
          tab[i]       =tmp10;
          tab[i+N/4]   =tmp20;
          tab[i+N/2]   =tmp30;
          tab[i+3*N/4] =tmp40;
        }
///////////////////////////////////////////

    if(nb102==0)
    {
      fi2=0;
    }
    else if(nb102==1)
    {
      fi2=fi;
    }

    increment=0;
    int flag2=0;
    for(int stg=1;stg<nb_stages;stg++)
    {
          //  cout<<"..if(flag2%2==0).."<<endl;
          //  system("pause");
         nr1=N/pow(rx4,0+increment);
        nr2=N/pow(rx4,1+increment);
        nr3=N/pow(rx4,2+increment);
        nr4=pow(rx4,(nb_stages-2-increment));
        nr5=pow(rx4,0+increment);
        nr31=1*nr3;
        nr32=2*nr3;
        nr33=3*nr3;
        for(int j=0;j<nr4;j++)
        {
            for(int i=0;i<rx4;i++)
            {
                nr20=nr2*i;
                tmp100=nr5*i*(j+nr4);
                tmp200=nr5*i*j;
                w1[0].real()= cos(tmp200*tmp5+fi2);
                w1[0].imag()=sin(tmp200*tmp5+fi2);
                w2[0].real()= cos(tmp100*tmp5+fi2);
                w2[0].imag()=sin(tmp100*tmp5+fi2);
                w3[0].real()= cos((2*tmp100-  tmp200)*tmp5+fi2);
                w3[0].imag()=sin((2*tmp100-  tmp200)*tmp5+fi2);
                w4[0].real()= cos((3*tmp100-2*tmp200)*tmp5+fi2);
                w4[0].imag()=sin((3*tmp100-2*tmp200)*tmp5+fi2);
                for(int m=0;m<nr5;m++)
                {
                    nr10=nr1*m;
                    tmp1=w1[0]*tab[     nr20+nr10+j];
                    tmp2=w2[0]*tab[nr31+nr20+nr10+j];
                    tmp3=w3[0]*tab[nr32+nr20+nr10+j];
                    tmp4=w4[0]*tab[nr33+nr20+nr10+j];
                    if(nb103==0)
                    {
                    tmp6=w10*(tmp1-tmp3);
                    tmp7=w11*(tmp2-tmp4);
                    //radix-4
                   tmp10=(w10*(tmp1+tmp2+tmp3+tmp4));
                   tmp20=(tmp6+tmp7);
                   tmp30=(w10*(tmp1-tmp2+tmp3-tmp4));
                   tmp40=(tmp6-tmp7);
                    }
                    else if(nb103==1)
                    {
                    tmp6=w20*(tmp1-tmp3);
                    tmp7=w21*(tmp2-tmp4);
                   //radix-4
                   tmp10=(w20*(tmp1+tmp2+tmp3+tmp4));
                   tmp20=(tmp6+tmp7);
                   tmp30=(w20*(tmp1-tmp2+tmp3-tmp4));
                   tmp40=(tmp6-tmp7);
                    }

                    tab[     nr20+nr10+j]=tmp10;
                    tab[nr31+nr20+nr10+j]=tmp20;
                    tab[nr32+nr20+nr10+j]=tmp30;
                    tab[nr33+nr20+nr10+j]=tmp40;
                }
            }
        }
        increment++;
    }
    ///////////////////////////////
    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_radix_4_basis_for_stage1(std::complex<double> &w0,std::complex<double> &w1)
{
    const double pi=3.141592653589793238462;
    //std::complex<double>    w0,w1,w2,w3,w4,w6,w9;
    double tmp5=2*pi/4;
    double fi2=fi;
    w0.real()=round(cos(0*tmp5+fi2)*100000)/100000;
    w0.imag()=round(-sin(0*tmp5+fi2)*100000)/100000;
    w1.real()=round(cos(1*tmp5+fi2)*100000)/100000;
    w1.imag()=round(-sin(1*tmp5+fi2)*100000)/100000;
    w0=w0;
    w1=w1;
    w0.real()=w0.real();
    w0.imag()=w0.imag();
    w1.real()=w1.real();
    w1.imag()=w1.imag();

    //w2.real()=round(cos(2*tmp5)*1000)/1000;//==-w0
    //w2.imag()=round(-sin(2*tmp5)*1000)/1000;//==-w0
    //w3.real()=round(cos(3*tmp5)*1000)/1000;//==-w1
    //w3.imag()=round(-sin(3*tmp5)*1000)/1000;//==-w1
    //w4.real()=round(cos(4*tmp5)*1000)/1000;//==+w0
    //w4.imag()=round(-sin(4*tmp5)*1000)/1000;//==+w0
    //w6.real()=round(cos(6*tmp5)*1000)/1000;//==-w0
    //w6.imag()=round(-sin(6*tmp5)*1000)/1000;//==-w0
    //w9.real()=round(cos(9*tmp5)*1000)/1000;//==+w1
    //w9.imag()=round(-sin(9*tmp5)*1000)/1000;//==+w1

   //cout<<" "<<w0<<" "<<w0<<" "<<w0<<" "<<w0<<endl;
   //cout<<" "<<w0<<" "<<w1<<" "<<w2<<" "<<w3<<endl;
   //cout<<" "<<w0<<" "<<w2<<" "<<w4<<" "<<w6<<endl;
   //cout<<" "<<w0<<" "<<w3<<" "<<w6<<" "<<w9<<endl;
   //system("pause");
}
//////////////////
void fun_radix_4_basis_for_stage_rest(std::complex<double> &w0,std::complex<double> &w1)
{
    const double pi=3.141592653589793238462;
    //std::complex<double>    w0,w1,w2,w3,w4,w6,w9;
    double tmp5=2*pi/4;
    double fi2=0;//always 0
    w0.real()=round(cos(0*tmp5+fi2)*100000)/100000;
    w0.imag()=round(-sin(0*tmp5+fi2)*100000)/100000;
    w1.real()=round(cos(1*tmp5+fi2)*100000)/100000;
    w1.imag()=round(-sin(1*tmp5+fi2)*100000)/100000;
    w0=w0;
    w1=w1;
    w0.real()=w0.real();
    w0.imag()=w0.imag();
    w1.real()=w1.real();
    w1.imag()=w1.imag();

    //w2.real()=round(cos(2*tmp5)*1000)/1000;//==-w0
    //w2.imag()=round(-sin(2*tmp5)*1000)/1000;//==-w0
    //w3.real()=round(cos(3*tmp5)*1000)/1000;//==-w1
    //w3.imag()=round(-sin(3*tmp5)*1000)/1000;//==-w1
    //w4.real()=round(cos(4*tmp5)*1000)/1000;//==+w0
    //w4.imag()=round(-sin(4*tmp5)*1000)/1000;//==+w0
    //w6.real()=round(cos(6*tmp5)*1000)/1000;//==-w0
    //w6.imag()=round(-sin(6*tmp5)*1000)/1000;//==-w0
    //w9.real()=round(cos(9*tmp5)*1000)/1000;//==+w1
    //w9.imag()=round(-sin(9*tmp5)*1000)/1000;//==+w1

   //cout<<" "<<w0<<" "<<w0<<" "<<w0<<" "<<w0<<endl;
   //cout<<" "<<w0<<" "<<w1<<" "<<w2<<" "<<w3<<endl;
   //cout<<" "<<w0<<" "<<w2<<" "<<w4<<" "<<w6<<endl;
   //cout<<" "<<w0<<" "<<w3<<" "<<w6<<" "<<w9<<endl;
   //system("pause");
}

//////////////////
void fun_inverse_radix_4_basis_for_stage1(std::complex<double> &w0,std::complex<double> &w1)
{
    const double pi=3.141592653589793238462;
    //std::complex<double>    w0,w1,w2,w3,w4,w6,w9;
    double tmp5=2*pi/4;
    double fi2=fi;
    w0.real()=round(cos(0*tmp5+fi2)*100000)/100000;
    w0.imag()=round(sin(0*tmp5+fi2)*100000)/100000;
    w1.real()=round(cos(1*tmp5+fi2)*100000)/100000;
    w1.imag()=round(sin(1*tmp5+fi2)*100000)/100000;
    w0=w0;
    w1=w1;
    w0.real()=w0.real();
    w0.imag()=w0.imag();
    w1.real()=w1.real();
    w1.imag()=w1.imag();

    //w2.real()=round(cos(2*tmp5)*1000)/1000;//==-w0
    //w2.imag()=round(-sin(2*tmp5)*1000)/1000;//==-w0
    //w3.real()=round(cos(3*tmp5)*1000)/1000;//==-w1
    //w3.imag()=round(-sin(3*tmp5)*1000)/1000;//==-w1
    //w4.real()=round(cos(4*tmp5)*1000)/1000;//==+w0
    //w4.imag()=round(-sin(4*tmp5)*1000)/1000;//==+w0
    //w6.real()=round(cos(6*tmp5)*1000)/1000;//==-w0
    //w6.imag()=round(-sin(6*tmp5)*1000)/1000;//==-w0
    //w9.real()=round(cos(9*tmp5)*1000)/1000;//==+w1
    //w9.imag()=round(-sin(9*tmp5)*1000)/1000;//==+w1

   //cout<<" "<<w0<<" "<<w0<<" "<<w0<<" "<<w0<<endl;
   //cout<<" "<<w0<<" "<<w1<<" "<<w2<<" "<<w3<<endl;
   //cout<<" "<<w0<<" "<<w2<<" "<<w4<<" "<<w6<<endl;
   //cout<<" "<<w0<<" "<<w3<<" "<<w6<<" "<<w9<<endl;
   //system("pause");
}
//////////////////
void fun_inverse_radix_4_basis_for_stage_rest(std::complex<double> &w0,std::complex<double> &w1)
{
    const double pi=3.141592653589793238462;
    //std::complex<double>    w0,w1,w2,w3,w4,w6,w9;
    double tmp5=2*pi/4;
    double fi2=0;//always 0
    w0.real()=round(cos(0*tmp5+fi2)*100000)/100000;
    w0.imag()=round(sin(0*tmp5+fi2)*100000)/100000;
    w1.real()=round(cos(1*tmp5+fi2)*100000)/100000;
    w1.imag()=round(sin(1*tmp5+fi2)*100000)/100000;
    w0=w0;
    w1=w1;
    w0.real()=w0.real();
    w0.imag()=w0.imag();
    w1.real()=w1.real();
    w1.imag()=w1.imag();

    //w2.real()=round(cos(2*tmp5)*1000)/1000;//==-w0
    //w2.imag()=round(-sin(2*tmp5)*1000)/1000;//==-w0
    //w3.real()=round(cos(3*tmp5)*1000)/1000;//==-w1
    //w3.imag()=round(-sin(3*tmp5)*1000)/1000;//==-w1
    //w4.real()=round(cos(4*tmp5)*1000)/1000;//==+w0
    //w4.imag()=round(-sin(4*tmp5)*1000)/1000;//==+w0
    //w6.real()=round(cos(6*tmp5)*1000)/1000;//==-w0
    //w6.imag()=round(-sin(6*tmp5)*1000)/1000;//==-w0
    //w9.real()=round(cos(9*tmp5)*1000)/1000;//==+w1
    //w9.imag()=round(-sin(9*tmp5)*1000)/1000;//==+w1

   //cout<<" "<<w0<<" "<<w0<<" "<<w0<<" "<<w0<<endl;
   //cout<<" "<<w0<<" "<<w1<<" "<<w2<<" "<<w3<<endl;
   //cout<<" "<<w0<<" "<<w2<<" "<<w4<<" "<<w6<<endl;
   //cout<<" "<<w0<<" "<<w3<<" "<<w6<<" "<<w9<<endl;
   //system("pause");
}



//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_4096_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_4096_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:

 fun_fourier_transform_FFT_radix_4_N_4096_official
{

}
fun_inverse_fourier_transform_FFT_radix_4_N_4096_official
{
    for(int i=0;i<N;i++)
    {
        tablica1[0][i]=tablica1[0][i]*1/(float)N;
        tablica1[1][i]=tablica1[1][i]*1/(float)N;
    }
}

*/



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

inverse fast fourier transform FFT radix-4 N points universal final

//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/
//http://www.cmlab.csie.ntu.edu.tw/cml/dsp/training/coding/transform/fft.html
//http://dsp.stackexchange.com/questions/3481/radix-4-fft-implementation

//https://community.arm.com/graphics/b/blog/posts/speeding-up-fast-fourier-transform-mixed-radix-on-mobile-arm-mali-gpu-by-means-of-opencl---part-1
//book: "Cyfrowe przetwarzanie sygnalow" - Tomasz Zielinski it has quick version of radix-2 because it calculates less sin() and cos()


//author marcin matysek (r)ewertyn.PL

//inverse fourier transform iFFT radix-4 for N=4096 points algorithm c++ source code implementation
//szybka odwrotna transformacja fouriera iFFT radix-4 dla N=4096 algorytm c++ kod zródlowy implementacja


#include <iostream>
#include "conio.h"
#include <stdlib.h>
#include <math.h>
#include <cmath>
#include <time.h>
#include <complex>
#include <fstream>

using namespace std;
//const double pi2=3.141592653589793238462;
double fi=0;

//complex number method:
void fun_inverse_bits_radix_4(int N,std::complex<double> tab[]);
void fun_fourier_transform_FFT_radix_4_N_points_universal_final(int N,std::complex<double> tab[],std::complex<double> w10,std::complex<double> w11);
void fun_inverse_fourier_transform_FFT_radix_4_N_points_universal_final(int N,std::complex<double> tab[],std::complex<double> w10,std::complex<double> w11);
void fun_fourier_transform_DFT_full_complex(int N,std::complex<double> tab[]);
void fun_radix_4_podstawa1(std::complex<double> &w0,std::complex<double> &w1);
void fun_radix_4_podstawa2(std::complex<double> &w0,std::complex<double> &w1);


std::complex<double> nb1,nb2,nb3,nb4,nb5,nb9,nb12,nb13,nb15,nb16,nb17;
int nb6=0,nb7=0,nb8=0,nb10=0,nb11=0,nb14=0;

static double diffclock(clock_t clock2,clock_t clock1)
{
    double diffticks=clock1-clock2;
    double diffms=(diffticks)/(CLOCKS_PER_SEC/1000);
    return diffms;
}
int main()
{
    const double pi=3.141592653589793238462;
    int N;
    //if N==period of signal in table tab[] then resolution = 1 Hz
    std::complex<double> w0;
    std::complex<double> w1;
    int nn=0,aa=0,bb=0;
    int maks=0;
    double time1,time2;
    double zmienna=0;
    int flag=0;
    std::complex<double> tab2[4096]={};
    std::complex<double> tab3[4096]={};

    N=4096;//need to change size of tab2 in function: fun_fourier_transform_DFT_full_complex



    for(int i=0;i<N;i++)
    {
        tab2[i].real()=i+1;
        tab2[i].imag()=N+i+1;
        tab3[i].real()=i+1;
        tab3[i].imag()=N+i+1;
    }


    cout<<"signal 1="<<endl;
     system("pause");
    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;
    for(int i=0;i<N;i++)
    {
       // tab3[i].real()=sin(i*2*pi/N+pi/2);
        //tab3[i].real()=tab3[i].real()+sin(i*2*2*pi/N+pi/2);
       // tab3[i].real()=tab3[i].real()+sin(i*3*2*pi/N+pi/2);
       // tab3[i].imag()=0;
    }
    cout<<"signal 2="<<endl;
    system("pause");
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab3[j].real()*1000)/1000<<"  ";
    }
    cout<<endl;
    cout<<endl;
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab3[j].imag()*1000)/1000<<"  ";
    }
    cout<<endl;
    cout<<endl;
    cout<<"      start calculations:"<<endl<<endl;
     system("pause");
     cout<<"calculating first DFT"<<endl;
    fun_fourier_transform_DFT_full_complex(N,tab3);
    //////////////////////////////////////////////////////////
     cout<<"calculating FFT"<<endl;
    clock_t start = clock();
    fun_radix_4_podstawa1(w0,w1);
    fun_fourier_transform_FFT_radix_4_N_points_universal_final(N,tab2,w0,w1);
    fun_inverse_bits_radix_4(N,tab2);
    time1=diffclock( start, clock() );
    ///////////////////////////////////////////////////////////


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

    cout<<"frequency Hz DFT real()"<<endl;
     system("pause");
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab3[j].real()*1000)/1000<<"  ";
    }
    cout<<endl;
    cout<<endl;
    cout<<"frequency Hz DFT imag()"<<endl;
     system("pause");
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab3[j].imag()*1000)/1000<<"  ";
    }
    cout<<endl;
    cout<<endl;
    cout<<"if radix-4 == DFT tab2[j].real(): "<<endl;system("pause");

    flag=0;
       for(int j=0;j<N;j++)
        {
          if(((tab3[j].real()-tab2[j].real()>=-0.03)&&(tab3[j].real()-tab2[j].real()<=0.03)))
          {
            flag++;
            //cout.precision(4);
            cout<<" j= "<<j<<" "<<round(tab2[j].real()*1000)/1000<<" ok ";
            //system("pause");
          }
            else {
                cout<<" j= "<<j<<" "<<-1<<" .not equal.. ";
                //system("pause");
            }
        }
        if(flag>maks){maks=flag;}

    cout<<endl<<"max equal= "<<maks<<endl;
    system("pause");
    cout<<endl;

    /////////////////////////////////////////////////////////////////
    start = clock();
    fun_radix_4_podstawa2(w0,w1);
    fun_inverse_fourier_transform_FFT_radix_4_N_points_universal_final(N,tab2,w0,w1);
    fun_inverse_bits_radix_4(N,tab2);
    time2=diffclock( start, clock() );
    ////////////////////////////////////////////////////////////////


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

    cout<<endl;
    cout<<"time for normal  FFT radix-4 = "<<time1<<endl;
    cout<<"time for inverse FFT radix-4 = "<<time2<<endl;
    system("pause");
    return 0;
}
//////////////////////////////////////////////
void fun_inverse_bits_radix_4(int N,std::complex<double> tab[])
{

//code by Sidney Burrus
    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;
    }
}
///////////////////////////////////////////////////
void fun_fourier_transform_DFT_full_complex(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double> tab2[4096*4]={};    // tab2[]==N
    std::complex<double>  w[1]={{1,1}};


double zmienna1=2*pi/(float)N;
double fi2=fi;

for (int i=0;i<N;i++)
{
    for(int j=0;j<N;j++)
    {
          //complex number method:
          w[0].real()=cos(i*j*zmienna1+fi2);
          w[0].imag()=(-sin(i*j*zmienna1+fi2));
          tab2[i]=tab2[i]+tab[j]*w[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_points_universal_final(int N,std::complex<double> tab[],std::complex<double> w10,std::complex<double> w11)
{
    const double pi=3.141592653589793238462;
    std::complex<double>  w1[1]={{1,0}};
    std::complex<double>  w2[1]={{1,0}};
    std::complex<double>  w3[1]={{1,0}};
    std::complex<double>  w4[1]={{1,0}};
    std::complex<double> tmp1,tmp2,tmp3,tmp4,tmp6,tmp7,tmp8;
    std::complex<double> tmp10,tmp20,tmp30,tmp40;
    int nr1=0,nr2=0,nr3=0,nr4=0,nr5=0;
    int nr10=0,nr20=0,nr31=0,nr32=0,nr33=0;
    int nb_stages=0,rx4=0;
    double tmp5,tmp100,tmp200;
    int increment=0;
    int nn=0;


    tmp5=2*pi/(N/1);
    rx4=4;//radix-4
    nn=N;
    for(int i=0;i<100;i++)
    {

        nn=(float)nn/(float)rx4;
         if(nn>=1.0)
         {
         nb_stages++;
         }
    }
//stage 1
          w1[0].real()=cos(0*tmp5);
          w1[0].imag()=-sin(0*tmp5);
          w2[0].real()=cos(0*tmp5);
          w2[0].imag()=-sin(0*tmp5);
          w3[0].real()=cos(0*tmp5);
          w3[0].imag()=-sin(0*tmp5);
          w4[0].real()=cos(0*tmp5);
          w4[0].imag()=-sin(0*tmp5);
        for(int i=0;i<(N/rx4);i++)
        {

          tmp1=w1[0]*tab[i+0];
          tmp2=w2[0]*tab[i+N/4];
          tmp3=w3[0]*tab[i+N/2];
          tmp4=w4[0]*tab[i+3*N/4];
          tmp6=w10*(tmp1-tmp3);
          tmp7=w11*(tmp2-tmp4);
          tmp8=tmp2+tmp4;
         //radix-4
            tmp10=w10*(tmp1+tmp3+tmp8);
            tmp20=tmp6+tmp7;
            tmp30=w10*(tmp1+tmp3-tmp8);
            tmp40=tmp6-tmp7;
          tab[i]       =tmp10;
          tab[i+N/4]   =tmp20;
          tab[i+N/2]   =tmp30;
          tab[i+3*N/4] =tmp40;
        }
///////////////////////////////////////////

    increment=0;

    for(int stg=1;stg<nb_stages;stg++)
    {
        nr1=N/pow(rx4,0+increment);
        nr2=N/pow(rx4,1+increment);
        nr3=N/pow(rx4,2+increment);
        nr4=pow(rx4,(nb_stages-2-increment));
        nr5=pow(rx4,0+increment);
        nr31=1*nr3;
        nr32=2*nr3;
        nr33=3*nr3;
        for(int j=0;j<nr4;j++)
        {
            for(int i=0;i<rx4;i++)
            {
                nr20=nr2*i;
                tmp100=nr5*i*(j+nr4)*tmp5;
                tmp200=nr5*i*j*tmp5;
                w1[0].real()= cos(tmp200);
                w1[0].imag()=-sin(tmp200);
                w2[0].real()= cos(tmp100);
                w2[0].imag()=-sin(tmp100);
                w3[0].real()= cos(2*tmp100-  tmp200);
                w3[0].imag()=-sin(2*tmp100-  tmp200);
                w4[0].real()=cos(3*tmp100-2*tmp200);
                w4[0].imag()=-sin(3*tmp100-2*tmp200);
                 for(int m=0;m<nr5;m++)
                 {
                    nr10=nr1*m;

                    tmp1=w1[0]*tab[     nr20+nr10+j];
                    tmp2=w2[0]*tab[nr31+nr20+nr10+j];
                    tmp3=w3[0]*tab[nr32+nr20+nr10+j];
                    tmp4=w4[0]*tab[nr33+nr20+nr10+j];
                    tmp6=w10*(tmp1-tmp3);
                    tmp7=w11*(tmp2-tmp4);
                    //radix-4
                    tmp10=w10*(tmp1+tmp2+tmp3+tmp4);
                    tmp20=tmp6+tmp7;
                    tmp30=w10*(tmp1-tmp2+tmp3-tmp4);
                    tmp40=tmp6-tmp7;

                    tab[     nr20+nr10+j]=tmp10;
                    tab[nr31+nr20+nr10+j]=tmp20;
                    tab[nr32+nr20+nr10+j]=tmp30;
                    tab[nr33+nr20+nr10+j]=tmp40;
                }
            }
        }
        increment++;
    }
///////////////////////////////////////////
    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_points_universal_final(int N,std::complex<double> tab[],std::complex<double> w10,std::complex<double> w11)
{
    const double pi=3.141592653589793238462;
    std::complex<double>  w1[1]={{1,0}};
    std::complex<double>  w2[1]={{1,0}};
    std::complex<double>  w3[1]={{1,0}};
    std::complex<double>  w4[1]={{1,0}};
    std::complex<double> tmp1,tmp2,tmp3,tmp4,tmp6,tmp7;
    std::complex<double> tmp10,tmp20,tmp30,tmp40;
    int nr1=0,nr2=0,nr3=0,nr4=0,nr5=0;
    int nr10=0,nr20=0,nr31=0,nr32=0,nr33=0;
    int nb_stages=0,rx4=0;
    double tmp5,tmp100,tmp200;
    int increment=0;
    int nn=0;


    tmp5=2*pi/(N/1);
    rx4=4;//radix-4
    nn=N;
    for(int i=0;i<100;i++)
    {

        nn=(float)nn/(float)rx4;
         if(nn>=1.0)
         {
         nb_stages++;
         }
    }
//stage 1
          w1[0].real()=cos(0*tmp5);
          w1[0].imag()= sin(0*tmp5);
          w2[0].real()=cos(0*tmp5);
          w2[0].imag()= sin(0*tmp5);
          w3[0].real()=cos(0*tmp5);
          w3[0].imag()= sin(0*tmp5);
          w4[0].real()=cos(0*tmp5);
          w4[0].imag()= sin(0*tmp5);
        for(int i=0;i<(N/rx4);i++)
        {

          tmp1=w1[0]*tab[i+0];
          tmp2=w2[0]*tab[i+N/4];
          tmp3=w3[0]*tab[i+N/2];
          tmp4=w4[0]*tab[i+3*N/4];
          tmp6=w10*(tmp1-tmp3);
          tmp7=w11*(tmp2-tmp4);
         //radix-4
        tmp10=w10*(tmp1+tmp2+tmp3+tmp4);
        tmp20=tmp6+tmp7;
        tmp30=w10*(tmp1-tmp2+tmp3-tmp4);
        tmp40=tmp6-tmp7;
          tab[i]       =tmp10;
          tab[i+N/4]   =tmp20;
          tab[i+N/2]   =tmp30;
          tab[i+3*N/4] =tmp40;
        }
///////////////////////////////////////////

    increment=0;

    for(int stg=1;stg<nb_stages;stg++)
    {
        nr1=N/pow(rx4,0+increment);
        nr2=N/pow(rx4,1+increment);
        nr3=N/pow(rx4,2+increment);
        nr4=pow(rx4,(nb_stages-2-increment));
        nr5=pow(rx4,0+increment);
        nr31=1*nr3;
        nr32=2*nr3;
        nr33=3*nr3;
        for(int j=0;j<nr4;j++)
        {

            for(int i=0;i<rx4;i++)
            {
                nr20=nr2*i;
                /*
                tmp100=nr5*i*(j+nr4)*tmp5;
                tmp200=nr5*i*j*tmp5;
                w1[0].real()= cos(tmp200);
                w1[0].imag()= sin(tmp200);
                w2[0].real()= cos(tmp100);
                w2[0].imag()= sin(tmp100);
                w3[0].real()= cos(2*tmp100-  tmp200);
                w3[0].imag()= sin(2*tmp100-  tmp200);
                w4[0].real()= cos(3*tmp100-2*tmp200);
                w4[0].imag()= sin(3*tmp100-2*tmp200);
                */

                tmp100=nr5*i*(j+nr4);
                tmp200=nr5*i*j;
                w1[0].real()= cos(tmp200*tmp5);
                w1[0].imag()= sin(tmp200*tmp5);
                w2[0].real()= cos(tmp100*tmp5);
                w2[0].imag()= sin(tmp100*tmp5);
                w3[0].real()= cos((2*tmp100-  tmp200)*tmp5);
                w3[0].imag()= sin((2*tmp100-  tmp200)*tmp5);
                w4[0].real()= cos((3*tmp100-2*tmp200)*tmp5);
                w4[0].imag()= sin((3*tmp100-2*tmp200)*tmp5);


                 for(int m=0;m<nr5;m++)
                 {
                    nr10=nr1*m;

                    tmp1=w1[0]*tab[     nr20+nr10+j];
                    tmp2=w2[0]*tab[nr31+nr20+nr10+j];
                    tmp3=w3[0]*tab[nr32+nr20+nr10+j];
                    tmp4=w4[0]*tab[nr33+nr20+nr10+j];
                    tmp6=w10*(tmp1-tmp3);
                    tmp7=w11*(tmp2-tmp4);
                    //radix-4
                    tmp10=w10*(tmp1+tmp2+tmp3+tmp4);
                    tmp20=tmp6+tmp7;
                    tmp30=w10*(tmp1-tmp2+tmp3-tmp4);
                    tmp40=tmp6-tmp7;

                    tab[     nr20+nr10+j]=tmp10;
                    tab[nr31+nr20+nr10+j]=tmp20;
                    tab[nr32+nr20+nr10+j]=tmp30;
                    tab[nr33+nr20+nr10+j]=tmp40;
                }
            }
        }
        increment++;
    }
///////////////////////////////
    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_radix_4_podstawa1(std::complex<double> &w0,std::complex<double> &w1)
{
    const double pi=3.141592653589793238462;
    //std::complex<double>    w0,w1,w2,w3,w4,w6,w9;
    double tmp5=2*pi/4;
    double fi2=fi;
    w0.real()=round(cos(0*tmp5+fi2)*1000)/1000;
    w0.imag()=round(-sin(0*tmp5+fi2)*1000)/1000;
    w1.real()=round(cos(1*tmp5+fi2)*1000)/1000;
    w1.imag()=round(-sin(1*tmp5+fi2)*1000)/1000;
    w0=w0;
    w1=w1;
    w0.real()=w0.real();
    w0.imag()=w0.imag();
    w1.real()=w1.real();
    w1.imag()=w1.imag();

    //w2.real()=round(cos(2*tmp5)*1000)/1000;//==-w0
    //w2.imag()=round(-sin(2*tmp5)*1000)/1000;//==-w0
    //w3.real()=round(cos(3*tmp5)*1000)/1000;//==-w1
    //w3.imag()=round(-sin(3*tmp5)*1000)/1000;//==-w1
    //w4.real()=round(cos(4*tmp5)*1000)/1000;//==+w0
    //w4.imag()=round(-sin(4*tmp5)*1000)/1000;//==+w0
    //w6.real()=round(cos(6*tmp5)*1000)/1000;//==-w0
    //w6.imag()=round(-sin(6*tmp5)*1000)/1000;//==-w0
    //w9.real()=round(cos(9*tmp5)*1000)/1000;//==+w1
    //w9.imag()=round(-sin(9*tmp5)*1000)/1000;//==+w1

   //cout<<" "<<w0<<" "<<w0<<" "<<w0<<" "<<w0<<endl;
   //cout<<" "<<w0<<" "<<w1<<" "<<w2<<" "<<w3<<endl;
   //cout<<" "<<w0<<" "<<w2<<" "<<w4<<" "<<w6<<endl;
   //cout<<" "<<w0<<" "<<w3<<" "<<w6<<" "<<w9<<endl;
   //system("pause");
}
void fun_radix_4_podstawa2(std::complex<double> &w0,std::complex<double> &w1)
{
    const double pi=3.141592653589793238462;
    //std::complex<double>    w0,w1,w2,w3,w4,w6,w9;
    double tmp5=2*pi/4;
    double fi2=fi;
    w0.real()=round(cos(0*tmp5+fi2)*1000)/1000;
    w0.imag()=round(sin(0*tmp5+fi2)*1000)/1000;
    w1.real()=round(cos(1*tmp5+fi2)*1000)/1000;
    w1.imag()=round(sin(1*tmp5+fi2)*1000)/1000;
    w0=w0;
    w1=w1;
    w0.real()=w0.real();
    w0.imag()=w0.imag();
    w1.real()=w1.real();
    w1.imag()=w1.imag();


   //cout<<" "<<w0<<" "<<w0<<" "<<w0<<" "<<w0<<endl;
   //cout<<" "<<w0<<" "<<w1<<" "<<w2<<" "<<w3<<endl;
   //cout<<" "<<w0<<" "<<w2<<" "<<w4<<" "<<w6<<endl;
   //cout<<" "<<w0<<" "<<w3<<" "<<w6<<" "<<w9<<endl;
   //system("pause");
}







//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_points_universal_final_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_points_universal_final_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:

 fun_fourier_transform_FFT_radix_4_N_points_universal_final_official
{

}
fun_inverse_fourier_transform_FFT_radix_4_N_points_universal_final_official
{
    for(int i=0;i<N;i++)
    {
        tablica1[0][i]=tablica1[0][i]*1/(float)N;
        tablica1[1][i]=tablica1[1][i]*1/(float)N;
    }
}

*/



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