sobota, 15 kwietnia 2017

Shared memory by 2 variables instead of pointers, 2 variables in 1 memory




//Shared memory by 2 variables instead of pointers, 2 variables in 1 memory
//współdzielenie pamięci przez 2 zmienne zamiast wskaźników, 2 zmienne w 1 pamięci
 //source:
//http://stackoverflow.com/questions/4629317/what-does-int-mean
// http://www.cplusplus.com/forum/windows/17153/

//author marcin matysek (r)ewertyn.PL
#include <iostream>
#include <conio.h>
#include <stdlib.h>

using namespace std;
void funkcja(int & nb3);

int main()
{
    int nb1=0;
    int &nb2=nb1;

    cout<<"Variable1 value= "<<nb1<<" Its address is: "<<&nb1<<endl;
    cout<<"Variable2 value= "<<nb2<<" Its address is: "<<&nb2<<endl;
    nb2=10;
    cout<<endl;
    cout<<"After modifying the 2 variables have value:"<<endl;
    cout<<"Variable1 value= "<<nb1<<" Its address is: "<<&nb1<<endl;
    cout<<"Variable2 value= "<<nb2<<" Its address is: "<<&nb2<<endl;

    funkcja(nb1);
    cout<<endl;
    cout<<"After leaving the function:"<<endl;
    cout<<"Variable1 value= "<<nb1<<" Its address is: "<<&nb1<<endl;
    cout<<"Variable2 value= "<<nb2<<" Its address is: "<<&nb2<<endl;

    system("pause");
    return 0;
}

void funkcja(int & nb3)
{
    cout<<endl;
    cout<<"in function Variable 3 have value and address of Variable1:"<<endl;
    cout<<"Variable3 value= "<<nb3<<" Its address is: "<<&nb3<<endl;
    cout<<"after modification Variable 3 ,Variable 3 In function has value:"<<endl;
    nb3=20;
    cout<<endl;
    cout<<"Variable3 value= "<<nb3<<" Its address is: "<<&nb3<<endl;
}

inverse fourier transform iFFT radix-4 for N= points algorithm c++ source code implementation v2



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

//author marcin matysek (r)ewertyn.PL

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

#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_4096(int N,std::complex<double> tab[]);
void fun_inverse_fourier_transform_FFT_radix_4_N_4096(int N,std::complex<double> tab[]);
void fun_fourier_transform_DFT_method5_full_complex(int N,std::complex<double> tab[]);

int fi=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()
{
    int N;
    //if N==period of signal in table tab[] then resolution = 1 Hz


    N=4096;//need change in fun_fourier_transform_DFT_method5_full_complex
    std::complex<double> tab2[4096]={};
    std::complex<double> tab3[4096]={};


    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;
    }

    double time2;
    double zmienna=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;
     system("pause");
    clock_t start = clock();
    cout<<"calculations"<<endl;
    fun_fourier_transform_DFT_method5_full_complex(N,tab3);
    //////////////////////////////////////////////////////////
    fun_fourier_transform_FFT_radix_4_N_4096(N,tab2);
    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");
       for(int j=0;j<N;j++)
        {
          if(((tab3[j].real()-tab2[j].real()>=-0.03)&&(tab3[j].real()-tab2[j].real()<=0.03)))
          {

            cout.precision(4);
            cout<<round(tab2[j].real()*1000)/1000<<"  ";//system("pause");
          }
            else {
                cout<<-1<<" .not equal.. ";
                system("pause");
            }
        }


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

    /////////////////////////////////////////////////////////////////
    fun_inverse_fourier_transform_FFT_radix_4_N_4096(N,tab2);
    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;
    cout<<endl;

    cout<<endl;
    system("pause");
    return 0;
}
//////////////////////////////////////////////
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;
    }
}
///////////////////////////////////////////////////
void fun_fourier_transform_DFT_method5_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[])
{
    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>  w6[1]={{1,0}};
    std::complex<double> tmp1,tmp2,tmp3,tmp4;
    std::complex<double> tmp10,tmp20,tmp30,tmp40;
    double nr100=0,nr200=0;
    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;
    int increment=0;
    int nn=0;




    w5[0].real()=0;
    w5[0].imag()=-1;
    w6[0].real()=0;
    w6[0].imag()=1;

    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
        for(int i=0;i<(N/rx4);i++)
        {
          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);

          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];
         //radix-4
        tmp10=tmp1+tmp2+       tmp3+tmp4;
        tmp20=tmp1-tmp3+w5[0]*(tmp2-tmp4);
        tmp30=tmp1-tmp2+       tmp3-tmp4;
        tmp40=tmp1-tmp3+w6[0]*(tmp2-tmp4);
          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 m=0;m<nr5;m++)
        {
          nr10=nr1*m;
            for(int i=0;i<rx4;i++)
            {
                nr20=nr2*i;
                for(int j=0;j<nr4;j++)
                {
                    //w1[0].real()= cos((nr5*0*i*(j+nr4*0)+nr5*1*i*j)*tmp5);
                    //w2[0].real()= cos((nr5*1*i*(j+nr4*1)+nr5*0*i*j)*tmp5);

                    nr100=nr5*i*(j+nr4)*tmp5;
                    nr200=nr5*i*j*tmp5;

                    w1[0].real()= cos(       +  nr200);
                    w1[0].imag()=-sin(       +  nr200);
                    w2[0].real()= cos(  nr100        );
                    w2[0].imag()=-sin(  nr100        );
                    w3[0].real()= cos(2*nr100-  nr200);
                    w3[0].imag()=-sin(2*nr100-  nr200);
                    w4[0].real()= cos(3*nr100-2*nr200);
                    w4[0].imag()=-sin(3*nr100-2*nr200);

                    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];
                    //radix-4
                    tmp10=tmp1+tmp2+       tmp3+tmp4;
                    tmp20=tmp1-tmp3+w5[0]*(tmp2-tmp4);
                    tmp30=tmp1-tmp2+       tmp3-tmp4;
                    tmp40=tmp1-tmp3+w6[0]*(tmp2-tmp4);

                    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[])
{
    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>  w6[1]={{1,0}};
    std::complex<double> tmp1,tmp2,tmp3,tmp4;
    std::complex<double> tmp10,tmp20,tmp30,tmp40;
    double nr100=0,nr200=0;
    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;
    int increment=0;
    int nn=0;


    w5[0].real()=0;
    w5[0].imag()=1;
    w6[0].real()=0;
    w6[0].imag()=-1;

    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
        for(int i=0;i<(N/rx4);i++)
        {
          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);

          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];
         //radix-4
        tmp10=tmp1+tmp2+       tmp3+tmp4;
        tmp20=tmp1-tmp3+w5[0]*(tmp2-tmp4);
        tmp30=tmp1-tmp2+       tmp3-tmp4;
        tmp40=tmp1-tmp3+w6[0]*(tmp2-tmp4);
          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 m=0;m<nr5;m++)
        {
          nr10=nr1*m;
            for(int i=0;i<rx4;i++)
            {
                nr20=nr2*i;
                for(int j=0;j<nr4;j++)
                {
                    //w1[0].real()= cos((nr5*0*i*(j+nr4*0)+nr5*1*i*j)*tmp5);
                    //w2[0].real()= cos((nr5*1*i*(j+nr4*1)+nr5*0*i*j)*tmp5);

                    nr100=nr5*i*(j+nr4)*tmp5;
                    nr200=nr5*i*j*tmp5;

                    w1[0].real()= cos(       +  nr200);
                    w1[0].imag()= sin(       +  nr200);
                    w2[0].real()= cos(  nr100        );
                    w2[0].imag()= sin(  nr100        );
                    w3[0].real()= cos(2*nr100-  nr200);
                    w3[0].imag()= sin(2*nr100-  nr200);
                    w4[0].real()= cos(3*nr100-2*nr200);
                    w4[0].imag()= sin(3*nr100-2*nr200);

                    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];
                    //radix-4
                    tmp10=tmp1+tmp2+       tmp3+tmp4;
                    tmp20=tmp1-tmp3+w5[0]*(tmp2-tmp4);
                    tmp30=tmp1-tmp2+       tmp3-tmp4;
                    tmp40=tmp1-tmp3+w6[0]*(tmp2-tmp4);

                    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;
    }

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





//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 fourier transform iFFT radix-4 for N=4096 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

//author marcin matysek (r)ewertyn.PL

//szybka odwrotna transformacja fouriera iFFT radix-4 dla N=4096 algorytm c++ kod zródlowy implementacja
//inverse fourier transform iFFT radix-4 for N=4096 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_4096(int N,std::complex<double> tab[]);
void fun_inverse_fourier_transform_FFT_radix_4_N_4096(int N,std::complex<double> tab[]);
void fun_fourier_transform_DFT_method5_full_complex(int N,std::complex<double> tab[]);

int fi=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()
{
    int N;
    //if N==period of signal in table tab[] then resolution = 1 Hz


    N=4096;
    std::complex<double> tab2[4096]={};
    std::complex<double> tab3[4096]={};


    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;
    }

    double time2;
    double zmienna=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;
     system("pause");
    clock_t start = clock();

    fun_fourier_transform_DFT_method5_full_complex(N,tab3);
    //////////////////////////////////////////////////////////
    fun_fourier_transform_FFT_radix_4_N_4096(N,tab2);
    fun_inverse_bits_radix_4(N,tab2);
    ///////////////////////////////////////////////////////////

    time2=diffclock( start, clock() );

    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");
       for(int j=0;j<N;j++)
        {
          if(((tab3[j].real()-tab2[j].real()>=-0.03)&&(tab3[j].real()-tab2[j].real()<=0.03)))
          {

            cout.precision(4);
            cout<<round(tab2[j].real()*1000)/1000<<"  ";//system("pause");
          }
            else {
                cout<<-1<<" ........ ";
            }
        }


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

    /////////////////////////////////////////////////////////////////
    fun_inverse_fourier_transform_FFT_radix_4_N_4096(N,tab2);
    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;
    cout<<endl;



    cout<<endl;
    system("pause");
    return 0;
}
//////////////////////////////////////////////
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;
    }
}
///////////////////////////////////////////////////
void fun_fourier_transform_DFT_method5_full_complex(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double> tab2[4096]={};    // 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[])
{
    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>  w6[1]={{1,0}};
    std::complex<double> tmp1,tmp2,tmp3,tmp4;
    std::complex<double> tmp10,tmp20,tmp30,tmp40;
    int nr1=0,nr2=0,nr3=0,nr4=0,nr5=0;;
    int nb_stages=0,rx4=0;
    double tmp5;
    int increment=0;
    int nn=0;




    w5[0].real()=0;
    w5[0].imag()=-1;
    w6[0].real()=0;
    w6[0].imag()=1;

    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
        for(int i=0;i<(N/rx4);i++)
        {
          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);

          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];
         //radix-4
        tmp10=tmp1+tmp2+       tmp3+tmp4;
        tmp20=tmp1-tmp3+w5[0]*(tmp2-tmp4);
        tmp30=tmp1-tmp2+       tmp3-tmp4;
        tmp40=tmp1-tmp3+w6[0]*(tmp2-tmp4);
          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=pow(rx4,0+increment);
    nr2=pow(rx4,1+increment);
    nr3=pow(rx4,2+increment);
    nr4=pow(rx4,(nb_stages-2-increment));
    nr5=pow(rx4,0+increment);
    for(int m=0;m<nr5;m++)//stage 2
    {
        for(int i=0;i<rx4;i++)
        {
            for(int j=0;j<nr4;j++)
            {
                w1[0].real()= cos((nr5*0*i*(j+nr4*0)+nr5*1*i*j)*tmp5);
                w1[0].imag()=-sin((nr5*0*i*(j+nr4*0)+nr5*1*i*j)*tmp5);
                w2[0].real()= cos((nr5*1*i*(j+nr4*1)+nr5*0*i*j)*tmp5);
                w2[0].imag()=-sin((nr5*1*i*(j+nr4*1)+nr5*0*i*j)*tmp5);
                w3[0].real()= cos((nr5*2*i*(j+nr4*1)-nr5*1*i*j)*tmp5);
                w3[0].imag()=-sin((nr5*2*i*(j+nr4*1)-nr5*1*i*j)*tmp5);
                w4[0].real()= cos((nr5*3*i*(j+nr4*1)-nr5*2*i*j)*tmp5);
                w4[0].imag()=-sin((nr5*3*i*(j+nr4*1)-nr5*2*i*j)*tmp5);

                tmp1=w1[0]*tab[0*N/(nr3)+N/(nr2)*i+N/(nr1)*m+j];
                tmp2=w2[0]*tab[1*N/(nr3)+N/(nr2)*i+N/(nr1)*m+j];
                tmp3=w3[0]*tab[2*N/(nr3)+N/(nr2)*i+N/(nr1)*m+j];
                tmp4=w4[0]*tab[3*N/(nr3)+N/(nr2)*i+N/(nr1)*m+j];
                //radix-4
                tmp10=tmp1+tmp2+       tmp3+tmp4;
                tmp20=tmp1-tmp3+w5[0]*(tmp2-tmp4);
                tmp30=tmp1-tmp2+       tmp3-tmp4;
                tmp40=tmp1-tmp3+w6[0]*(tmp2-tmp4);

                tab[0*N/(nr3)+N/(nr2)*i+N/(nr1)*m+j]=tmp10;
                tab[1*N/(nr3)+N/(nr2)*i+N/(nr1)*m+j]=tmp20;
                tab[2*N/(nr3)+N/(nr2)*i+N/(nr1)*m+j]=tmp30;
                tab[3*N/(nr3)+N/(nr2)*i+N/(nr1)*m+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[])
{
    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>  w6[1]={{1,0}};
    std::complex<double> tmp1,tmp2,tmp3,tmp4;
    std::complex<double> tmp10,tmp20,tmp30,tmp40;
    int nr1=0,nr2=0,nr3=0,nr4=0,nr5=0;;
    int nb_stages=0,rx4=0;
    double tmp5;
    int increment=0;
    int nn=0;


    w5[0].real()=0;
    w5[0].imag()=1;
    w6[0].real()=0;
    w6[0].imag()=-1;

    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
        for(int i=0;i<(N/rx4);i++)
        {
          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);

          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];
         //radix-4
        tmp10=tmp1+tmp2+       tmp3+tmp4;
        tmp20=tmp1-tmp3+w5[0]*(tmp2-tmp4);
        tmp30=tmp1-tmp2+       tmp3-tmp4;
        tmp40=tmp1-tmp3+w6[0]*(tmp2-tmp4);
          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=pow(rx4,0+increment);
    nr2=pow(rx4,1+increment);
    nr3=pow(rx4,2+increment);
    nr4=pow(rx4,(nb_stages-2-increment));
    nr5=pow(rx4,0+increment);
    for(int m=0;m<nr5;m++)//stage 2
    {
        for(int i=0;i<rx4;i++)
        {
            for(int j=0;j<nr4;j++)
            {
                w1[0].real()= cos((nr5*0*i*(j+nr4*0)+nr5*1*i*j)*tmp5);
                w1[0].imag()= sin((nr5*0*i*(j+nr4*0)+nr5*1*i*j)*tmp5);
                w2[0].real()= cos((nr5*1*i*(j+nr4*1)+nr5*0*i*j)*tmp5);
                w2[0].imag()= sin((nr5*1*i*(j+nr4*1)+nr5*0*i*j)*tmp5);
                w3[0].real()= cos((nr5*2*i*(j+nr4*1)-nr5*1*i*j)*tmp5);
                w3[0].imag()= sin((nr5*2*i*(j+nr4*1)-nr5*1*i*j)*tmp5);
                w4[0].real()= cos((nr5*3*i*(j+nr4*1)-nr5*2*i*j)*tmp5);
                w4[0].imag()= sin((nr5*3*i*(j+nr4*1)-nr5*2*i*j)*tmp5);

                tmp1=w1[0]*tab[0*N/(nr3)+N/(nr2)*i+N/(nr1)*m+j];
                tmp2=w2[0]*tab[1*N/(nr3)+N/(nr2)*i+N/(nr1)*m+j];
                tmp3=w3[0]*tab[2*N/(nr3)+N/(nr2)*i+N/(nr1)*m+j];
                tmp4=w4[0]*tab[3*N/(nr3)+N/(nr2)*i+N/(nr1)*m+j];
                //radix-4
                tmp10=tmp1+tmp2+       tmp3+tmp4;
                tmp20=tmp1-tmp3+w5[0]*(tmp2-tmp4);
                tmp30=tmp1-tmp2+       tmp3-tmp4;
                tmp40=tmp1-tmp3+w6[0]*(tmp2-tmp4);

                tab[0*N/(nr3)+N/(nr2)*i+N/(nr1)*m+j]=tmp10;
                tab[1*N/(nr3)+N/(nr2)*i+N/(nr1)*m+j]=tmp20;
                tab[2*N/(nr3)+N/(nr2)*i+N/(nr1)*m+j]=tmp30;
                tab[3*N/(nr3)+N/(nr2)*i+N/(nr1)*m+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;
    }

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








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

czwartek, 13 kwietnia 2017

fast inverse fourier transform iFFT radix-4 for N=1024 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

//author marcin matysek (r)ewertyn.PL
//fast inverse fourier transform iFFT radix-4 for N=1024 algorithm c++ source code implementation
//szybka odwrotna transformacja fouriera iFFT radix-4 dla N=1024 algorytm c++ kod źródłowy implementacja

#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_1024(int N,std::complex<double> tab[]);
void fun_inverse_fourier_transform_FFT_radix_4_N_1024(int N,std::complex<double> tab[]);
void fun_fourier_transform_DFT_method5_full_complex(int N,std::complex<double> tab[]);

int fi=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()
{
    int N;
    //if N==period of signal in table tab[] then resolution = 1 Hz

    N=1024;
    std::complex<double> tab2[4096]={};
    std::complex<double> tab3[4096]={};
    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;
    }

    double time2;
    double zmienna=0;

    cout<<"signal 1="<<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<<"signal 2="<<endl;
    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;
     system("pause");
    clock_t start = clock();

    fun_fourier_transform_DFT_method5_full_complex(N,tab3);
    //////////////////////////////////////////////////////////
    fun_fourier_transform_FFT_radix_4_N_1024(N,tab2);
    fun_inverse_bits_radix_4(N,tab2);
    ///////////////////////////////////////////////////////////

    time2=diffclock( start, clock() );

    cout<<"frequency Hz radix-4"<<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<<"frequency Hz DFT"<<endl;
    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<<"if radix-4 == DFT tab2[j].real(): "<<endl;system("pause");
       for(int j=0;j<N;j++)
        {
          if(((tab3[j].real()-tab2[j].real()>=-0.03)&&(tab3[j].real()-tab2[j].real()<=0.03)))
          {

            cout.precision(4);
            cout<<round(tab2[j].real()*1000)/1000<<"  ";//system("pause");
          }
            else {
                cout<<-1<<" ........ ";
            }
        }


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

    /////////////////////////////////////////////////////////////////
    fun_inverse_fourier_transform_FFT_radix_4_N_1024(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[])
{

//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;
    }
}
///////////////////////////////////////////////////
void fun_fourier_transform_DFT_method5_full_complex(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double> tab2[4096]={};    // 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_1024(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double>  tab2[4096]={};    // tab2[]==N
    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>  w6[1]={{1,0}};
    std::complex<double> tmp1,tmp2,tmp3,tmp4;
    double tmp5;


    w5[0].real()=0;
    w5[0].imag()=-1;
    w6[0].real()=0;
    w6[0].imag()=1;

    tmp5=2*pi/(N/1);

//stage 1
        for(int i=0;i<(N/4);i++)//256
        {
          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);

          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];
         //radix-4
          tab2[i]       =tmp1+tmp2+tmp3+tmp4;
          tab2[i+N/4]   =tmp1-tmp3+w5[0]*(tmp2-tmp4);
          tab2[i+N/2]   =tmp1-tmp2+tmp3-tmp4;
          tab2[i+3*N/4] =tmp1-tmp3+w6[0]*(tmp2-tmp4);
        }
/////////////////////////////////////////////////////////////////////////////////
    for(int m=0;m<1;m++)//stage 2
    {
        for(int i=0;i<4;i++)
        {
            for(int j=0;j<4*4*4;j++)
            {
                w1[0].real()= cos((1*0*i*(j+4*4*4*0)+1*1*i*j)*tmp5);
                w1[0].imag()=-sin((1*0*i*(j+4*4*4*0)+1*1*i*j)*tmp5);
                w2[0].real()= cos((1*1*i*(j+4*4*4*1)+1*0*i*j)*tmp5);
                w2[0].imag()=-sin((1*1*i*(j+4*4*4*1)+1*0*i*j)*tmp5);
                w3[0].real()= cos((1*2*i*(j+4*4*4*1)-1*1*i*j)*tmp5);
                w3[0].imag()=-sin((1*2*i*(j+4*4*4*1)-1*1*i*j)*tmp5);
                w4[0].real()= cos((1*3*i*(j+4*4*4*1)-1*2*i*j)*tmp5);
                w4[0].imag()=-sin((1*3*i*(j+4*4*4*1)-1*2*i*j)*tmp5);

                tmp1=w1[0]*tab2[0*N/(4*4)+N/(4)*i+N/(1)*m+j];
                tmp2=w2[0]*tab2[1*N/(4*4)+N/(4)*i+N/(1)*m+j];
                tmp3=w3[0]*tab2[2*N/(4*4)+N/(4)*i+N/(1)*m+j];
                tmp4=w4[0]*tab2[3*N/(4*4)+N/(4)*i+N/(1)*m+j];
                //radix-4
                tab[0*N/(4*4)+N/(4)*i+N/(1)*m+j]=tmp1+tmp2+       tmp3+tmp4;
                tab[1*N/(4*4)+N/(4)*i+N/(1)*m+j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
                tab[2*N/(4*4)+N/(4)*i+N/(1)*m+j]=tmp1-tmp2+       tmp3-tmp4;
                tab[3*N/(4*4)+N/(4)*i+N/(1)*m+j]=tmp1-tmp3+w6[0]*(tmp2-tmp4);
            }
        }
    }
////////////////////////////////////////////////////////////////////////////////
    for(int m=0;m<4;m++)//stage 3
    {
        for(int i=0;i<4;i++)
        {
            for(int j=0;j<4*4;j++)
            {
                w1[0].real()= cos((4*0*i*(j+4*4*0)+4*1*i*j)*tmp5);
                w1[0].imag()=-sin((4*0*i*(j+4*4*0)+4*1*i*j)*tmp5);
                w2[0].real()= cos((4*1*i*(j+4*4*1)+4*0*i*j)*tmp5);
                w2[0].imag()=-sin((4*1*i*(j+4*4*1)+4*0*i*j)*tmp5);
                w3[0].real()= cos((4*2*i*(j+4*4*1)-4*1*i*j)*tmp5);
                w3[0].imag()=-sin((4*2*i*(j+4*4*1)-4*1*i*j)*tmp5);
                w4[0].real()= cos((4*3*i*(j+4*4*1)-4*2*i*j)*tmp5);
                w4[0].imag()=-sin((4*3*i*(j+4*4*1)-4*2*i*j)*tmp5);

                tmp1=w1[0]*tab[0*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j];
                tmp2=w2[0]*tab[1*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j];
                tmp3=w3[0]*tab[2*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j];
                tmp4=w4[0]*tab[3*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j];
                //radix-4
                tab2[0*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=tmp1+tmp2+       tmp3+tmp4;
                tab2[1*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
                tab2[2*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=tmp1-tmp2+       tmp3-tmp4;
                tab2[3*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=tmp1-tmp3+w6[0]*(tmp2-tmp4);
            }
        }
    }
////////////////////////////////////////////////////////////////////////////////
    for(int m=0;m<4*4;m++)//stage 4
    {
        for(int i=0;i<4;i++)
        {
            for(int j=0;j<4;j++)
            {
                w1[0].real()= cos((4*4*0*i*(j+4*0)+4*4*1*i*j)*tmp5);
                w1[0].imag()=-sin((4*4*0*i*(j+4*0)+4*4*1*i*j)*tmp5);
                w2[0].real()= cos((4*4*1*i*(j+4*1)+4*4*0*i*j)*tmp5);
                w2[0].imag()=-sin((4*4*1*i*(j+4*1)+4*4*0*i*j)*tmp5);
                w3[0].real()= cos((4*4*2*i*(j+4*1)-4*4*1*i*j)*tmp5);
                w3[0].imag()=-sin((4*4*2*i*(j+4*1)-4*4*1*i*j)*tmp5);
                w4[0].real()= cos((4*4*3*i*(j+4*1)-4*4*2*i*j)*tmp5);
                w4[0].imag()=-sin((4*4*3*i*(j+4*1)-4*4*2*i*j)*tmp5);

                tmp1=w1[0]*tab2[0*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j];
                tmp2=w2[0]*tab2[1*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j];
                tmp3=w3[0]*tab2[2*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j];
                tmp4=w4[0]*tab2[3*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j];
                //radix-4
                tab[0*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=tmp1+tmp2+       tmp3+tmp4;
                tab[1*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
                tab[2*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=tmp1-tmp2+       tmp3-tmp4;
                tab[3*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=tmp1-tmp3+w6[0]*(tmp2-tmp4);
            }
        }
    }
///////////////////////////////////////////////////////////////////////////////
   for(int m=0;m<4*4*4;m++)//stage 5
    {
        for(int i=0;i<4;i++)
        {
            for(int j=0;j<1;j++)
            {
                w1[0].real()= cos((4*4*4*0*i*(j+1*0)+4*4*4*1*i*j)*tmp5);
                w1[0].imag()=-sin((4*4*4*0*i*(j+1*0)+4*4*4*1*i*j)*tmp5);
                w2[0].real()= cos((4*4*4*1*i*(j+1*1)+4*4*4*0*i*j)*tmp5);
                w2[0].imag()=-sin((4*4*4*1*i*(j+1*1)+4*4*4*0*i*j)*tmp5);
                w3[0].real()= cos((4*4*4*2*i*(j+1*1)-4*4*4*1*i*j)*tmp5);
                w3[0].imag()=-sin((4*4*4*2*i*(j+1*1)-4*4*4*1*i*j)*tmp5);
                w4[0].real()= cos((4*4*4*3*i*(j+1*1)-4*4*4*2*i*j)*tmp5);
                w4[0].imag()=-sin((4*4*4*3*i*(j+1*1)-4*4*4*2*i*j)*tmp5);

                tmp1=w1[0]*tab[0*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j];
                tmp2=w2[0]*tab[1*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j];
                tmp3=w3[0]*tab[2*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j];
                tmp4=w4[0]*tab[3*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j];
                //radix-4
                tab2[0*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j]=tmp1+tmp2+       tmp3+tmp4;
                tab2[1*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
                tab2[2*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j]=tmp1-tmp2+       tmp3-tmp4;
                tab2[3*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j]=tmp1-tmp3+w6[0]*(tmp2-tmp4);
            }
        }
    }

    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_1024(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double>  tab2[4096]={};    // tab2[]==N
    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>  w6[1]={{1,0}};
    std::complex<double> tmp1,tmp2,tmp3,tmp4;
    double tmp5;


    w5[0].real()=0;
    w5[0].imag()=1;
    w6[0].real()=0;
    w6[0].imag()=-1;

    tmp5=2*pi/(N/1);

//stage 1
        for(int i=0;i<(N/4);i++)//256
        {
          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);

          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];
         //radix-4
          tab2[i]       =tmp1+tmp2+tmp3+tmp4;
          tab2[i+N/4]   =tmp1-tmp3+w5[0]*(tmp2-tmp4);
          tab2[i+N/2]   =tmp1-tmp2+tmp3-tmp4;
          tab2[i+3*N/4] =tmp1-tmp3+w6[0]*(tmp2-tmp4);
        }
/////////////////////////////////////////////////////////////////////////////////
    for(int m=0;m<1;m++)//stage 2
    {
        for(int i=0;i<4;i++)
        {
            for(int j=0;j<4*4*4;j++)
            {
                w1[0].real()= cos((1*0*i*(j+4*4*4*0)+1*1*i*j)*tmp5);
                w1[0].imag()=sin((1*0*i*(j+4*4*4*0)+1*1*i*j)*tmp5);
                w2[0].real()= cos((1*1*i*(j+4*4*4*1)+1*0*i*j)*tmp5);
                w2[0].imag()=sin((1*1*i*(j+4*4*4*1)+1*0*i*j)*tmp5);
                w3[0].real()= cos((1*2*i*(j+4*4*4*1)-1*1*i*j)*tmp5);
                w3[0].imag()=sin((1*2*i*(j+4*4*4*1)-1*1*i*j)*tmp5);
                w4[0].real()= cos((1*3*i*(j+4*4*4*1)-1*2*i*j)*tmp5);
                w4[0].imag()=sin((1*3*i*(j+4*4*4*1)-1*2*i*j)*tmp5);

                tmp1=w1[0]*tab2[0*N/(4*4)+N/(4)*i+N/(1)*m+j];
                tmp2=w2[0]*tab2[1*N/(4*4)+N/(4)*i+N/(1)*m+j];
                tmp3=w3[0]*tab2[2*N/(4*4)+N/(4)*i+N/(1)*m+j];
                tmp4=w4[0]*tab2[3*N/(4*4)+N/(4)*i+N/(1)*m+j];
                //radix-4
                tab[0*N/(4*4)+N/(4)*i+N/(1)*m+j]=tmp1+tmp2+       tmp3+tmp4;
                tab[1*N/(4*4)+N/(4)*i+N/(1)*m+j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
                tab[2*N/(4*4)+N/(4)*i+N/(1)*m+j]=tmp1-tmp2+       tmp3-tmp4;
                tab[3*N/(4*4)+N/(4)*i+N/(1)*m+j]=tmp1-tmp3+w6[0]*(tmp2-tmp4);
            }
        }
    }
////////////////////////////////////////////////////////////////////////////////
    for(int m=0;m<4;m++)//stage 3
    {
        for(int i=0;i<4;i++)
        {
            for(int j=0;j<4*4;j++)
            {
                w1[0].real()= cos((4*0*i*(j+4*4*0)+4*1*i*j)*tmp5);
                w1[0].imag()=sin((4*0*i*(j+4*4*0)+4*1*i*j)*tmp5);
                w2[0].real()= cos((4*1*i*(j+4*4*1)+4*0*i*j)*tmp5);
                w2[0].imag()=sin((4*1*i*(j+4*4*1)+4*0*i*j)*tmp5);
                w3[0].real()= cos((4*2*i*(j+4*4*1)-4*1*i*j)*tmp5);
                w3[0].imag()=sin((4*2*i*(j+4*4*1)-4*1*i*j)*tmp5);
                w4[0].real()= cos((4*3*i*(j+4*4*1)-4*2*i*j)*tmp5);
                w4[0].imag()=sin((4*3*i*(j+4*4*1)-4*2*i*j)*tmp5);

                tmp1=w1[0]*tab[0*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j];
                tmp2=w2[0]*tab[1*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j];
                tmp3=w3[0]*tab[2*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j];
                tmp4=w4[0]*tab[3*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j];
                //radix-4
                tab2[0*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=tmp1+tmp2+       tmp3+tmp4;
                tab2[1*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
                tab2[2*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=tmp1-tmp2+       tmp3-tmp4;
                tab2[3*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=tmp1-tmp3+w6[0]*(tmp2-tmp4);
            }
        }
    }
////////////////////////////////////////////////////////////////////////////////
    for(int m=0;m<4*4;m++)//stage 4
    {
        for(int i=0;i<4;i++)
        {
            for(int j=0;j<4;j++)
            {
                w1[0].real()= cos((4*4*0*i*(j+4*0)+4*4*1*i*j)*tmp5);
                w1[0].imag()=sin((4*4*0*i*(j+4*0)+4*4*1*i*j)*tmp5);
                w2[0].real()= cos((4*4*1*i*(j+4*1)+4*4*0*i*j)*tmp5);
                w2[0].imag()=sin((4*4*1*i*(j+4*1)+4*4*0*i*j)*tmp5);
                w3[0].real()=cos((4*4*2*i*(j+4*1)-4*4*1*i*j)*tmp5);
                w3[0].imag()=sin((4*4*2*i*(j+4*1)-4*4*1*i*j)*tmp5);
                w4[0].real()= cos((4*4*3*i*(j+4*1)-4*4*2*i*j)*tmp5);
                w4[0].imag()=sin((4*4*3*i*(j+4*1)-4*4*2*i*j)*tmp5);

                tmp1=w1[0]*tab2[0*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j];
                tmp2=w2[0]*tab2[1*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j];
                tmp3=w3[0]*tab2[2*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j];
                tmp4=w4[0]*tab2[3*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j];
                //radix-4
                tab[0*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=tmp1+tmp2+       tmp3+tmp4;
                tab[1*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
                tab[2*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=tmp1-tmp2+       tmp3-tmp4;
                tab[3*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=tmp1-tmp3+w6[0]*(tmp2-tmp4);
            }
        }
    }
///////////////////////////////////////////////////////////////////////////////
   for(int m=0;m<4*4*4;m++)//stage 5
    {
        for(int i=0;i<4;i++)
        {
            for(int j=0;j<1;j++)
            {
                w1[0].real()= cos((4*4*4*0*i*(j+1*0)+4*4*4*1*i*j)*tmp5);
                w1[0].imag()=sin((4*4*4*0*i*(j+1*0)+4*4*4*1*i*j)*tmp5);
                w2[0].real()= cos((4*4*4*1*i*(j+1*1)+4*4*4*0*i*j)*tmp5);
                w2[0].imag()=sin((4*4*4*1*i*(j+1*1)+4*4*4*0*i*j)*tmp5);
                w3[0].real()= cos((4*4*4*2*i*(j+1*1)-4*4*4*1*i*j)*tmp5);
                w3[0].imag()=sin((4*4*4*2*i*(j+1*1)-4*4*4*1*i*j)*tmp5);
                w4[0].real()= cos((4*4*4*3*i*(j+1*1)-4*4*4*2*i*j)*tmp5);
                w4[0].imag()=sin((4*4*4*3*i*(j+1*1)-4*4*4*2*i*j)*tmp5);

                tmp1=w1[0]*tab[0*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j];
                tmp2=w2[0]*tab[1*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j];
                tmp3=w3[0]*tab[2*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j];
                tmp4=w4[0]*tab[3*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j];
                //radix-4
                tab2[0*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j]=tmp1+tmp2+       tmp3+tmp4;
                tab2[1*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
                tab2[2*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j]=tmp1-tmp2+       tmp3-tmp4;
                tab2[3*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j]=tmp1-tmp3+w6[0]*(tmp2-tmp4);
            }
        }
    }

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

}



//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:
 fun_fourier_transform_FFT_radix_4_N_256_official
{

}
fun_inverse_fourier_transform_FFT_radix_4_N_256_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)