mixed radix fast fourier transform FFT radix-5 * radix-3 * radix-2 * radix-4 * radix-3 for N=360 points example c++ v2

//mixed radix fast fourier transform FFT radix-5 * radix-3 * radix-2 * radix-4 * radix-3 for N=360 points example c++ v2



void fun_fourier_transform_FFT_mixed_radix_5_3_2_4_3_N_360_points(int N,std::complex<double> tab[])
{
  //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()

//Rabiner L.R., Gold B. Theory and application of digital signal processing p 378



    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>  w34[1]={{1,0}};
    std::complex<double>  w40[1]={{1,0}};
    std::complex<double>  w5[1]={{1,0}};
    std::complex<double>  w6[1]={{1,0}};
    std::complex<double>  w50[1]={{1,0}};
    std::complex<double>  w60[1]={{1,0}};
    std::complex<double>  w7[1]={{1,0}};
    std::complex<double>  w11[1]={{1,0}};
    std::complex<double>  w12[1]={{1,0}};
    std::complex<double>  w13[1]={{1,0}};
    std::complex<double>  w14[1]={{1,0}};
    std::complex<double>  w15[1]={{1,0}};
    std::complex<double>  w16[1]={{1,0}};
    std::complex<double>  w17[1]={{1,0}};
    std::complex<double>  w18[1]={{1,0}};
    std::complex<double>  w19[1]={{1,0}};
    std::complex<double>  w20[1]={{1,0}};
    std::complex<double> tmp1,tmp2,tmp3,tmp4,tmp11;
    double tmp5;
    double tmp6;
    double tmp7;
    double tmp8;

    int rx5=5,rx4=4,rx3=3,rx2=2;
    int stg[100]={};

    int nb1,nb2,nb3,nb4;

    stg[0]=1;
    stg[1]=rx5;
    stg[2]=rx3;
    stg[3]=rx2;
    stg[4]=rx4;
    stg[5]=rx3;



    tmp5=2*pi/(N/1);
    tmp6=2*pi/(N/1);
    tmp7=2*pi/(3/1);
    tmp8=2*pi/(5/1);

//radix 3 fundament
          w34[0].real()=cos(0*tmp7);
          w34[0].imag()=-sin(0*tmp7);
          w5[0].real()=cos(1*tmp7);
          w5[0].imag()=-sin(1*tmp7);
          w6[0].real()=cos(2*tmp7);
          w6[0].imag()=-sin(2*tmp7);
          w7[0].real()=cos(4*tmp7);
          w7[0].imag()=-sin(4*tmp7);


//radix 5 fundament
          w11[0].real()=cos(0*tmp8);
          w11[0].imag()=-sin(0*tmp8);
          w12[0].real()=cos(1*tmp8);
          w12[0].imag()=-sin(1*tmp8);
          w13[0].real()=cos(2*tmp8);
          w13[0].imag()=-sin(2*tmp8);
          w14[0].real()=cos(3*tmp8);
          w14[0].imag()=-sin(3*tmp8);
          w15[0].real()=cos(4*tmp8);
          w15[0].imag()=-sin(4*tmp8);
          w16[0].real()=cos(6*tmp8);
          w16[0].imag()=-sin(6*tmp8);
          w17[0].real()=cos(8*tmp8);
          w17[0].imag()=-sin(8*tmp8);
          w18[0].real()=cos(9*tmp8);
          w18[0].imag()=-sin(9*tmp8);
          w19[0].real()=cos(12*tmp8);
          w19[0].imag()=-sin(12*tmp8);
          w20[0].real()=cos(16*tmp8);
          w20[0].imag()=-sin(16*tmp8);

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

//stage 1  radix-5
          w1[0].real()=cos(0);
          w1[0].imag()=-sin(0);
          w2[0].real()=cos(0);
          w2[0].imag()=-sin(0);
          w3[0].real()=cos(0);
          w3[0].imag()=-sin(0);
          w4[0].real()=cos(0);
          w4[0].imag()=-sin(0);
          w40[0].real()=cos(0);
          w40[0].imag()=-sin(0);
        for(int i=0;i<3*4*2*3;i++)
        {
          tmp1=w1[0]*tab[i+0];
          tmp2=w2[0]*tab[i+1*N/(stg[1])];
          tmp3=w3[0]*tab[i+2*N/(stg[1])];
          tmp4=w4[0]*tab[i+3*N/(stg[1])];
          tmp11=w40[0]*tab[i+4*N/(stg[1])];
         //radix-5
          tab2[i]             =w11[0]*tmp1+w11[0]*tmp2+w11[0]*tmp3+w11[0]*tmp4+w11[0]*tmp11;
          tab2[i+1*N/(stg[1])]   =w11[0]*tmp1+w12[0]*tmp2+w13[0]*tmp3+w14[0]*tmp4+w15[0]*tmp11;
          tab2[i+2*N/(stg[1])]   =w11[0]*tmp1+w13[0]*tmp2+w15[0]*tmp3+w16[0]*tmp4+w17[0]*tmp11;
          tab2[i+3*N/(stg[1])]   =w11[0]*tmp1+w14[0]*tmp2+w16[0]*tmp3+w18[0]*tmp4+w19[0]*tmp11;
          tab2[i+4*N/(stg[1])]   =w11[0]*tmp1+w15[0]*tmp2+w17[0]*tmp3+w19[0]*tmp4+w20[0]*tmp11;
        }

//stage 2 radix-3
    int b=0;
    int i=0;
    nb1=N;
    nb4=1;
    nb1=nb1/stg[0+i];
    nb2=nb1/stg[1+i];
    nb3=nb2/stg[2+i];
    nb4=nb4*stg[0+i];

    for(int b=0;b<stg[1+i];b=b+1)
    {
        for(int j=0;j<nb3;j=j+1)
        {
          w1[0].real()=cos(nb4*b*(0*nb3+j)*tmp6);
          w1[0].imag()=-sin(nb4*b*(0*nb3+j)*tmp6);
          w2[0].real()=cos(nb4*b*(1*nb3+j)*tmp6);
          w2[0].imag()=-sin(nb4*b*(1*nb3+j)*tmp6);
          w3[0].real()=cos(nb4*b*(2*nb3+j)*tmp6);
          w3[0].imag()=-sin(nb4*b*(2*nb3+j)*tmp6);

            for(int i=0;i<nb4;i=i+1)
            {
              tmp1=w1[0]*tab2[i*nb1+b*nb2+0*nb3+j];
              tmp2=w2[0]*tab2[i*nb1+b*nb2+1*nb3+j];
              tmp3=w3[0]*tab2[i*nb1+b*nb2+2*nb3+j];
             //radix-3
              tab[i*nb1+b*nb2+0*nb3+j]  =w34[0]*tmp1+w34[0]*tmp2+w34[0]*tmp3;
              tab[i*nb1+b*nb2+1*nb3+j]  =w34[0]*tmp1+w5[0]*tmp2+w6[0]*tmp3;
              tab[i*nb1+b*nb2+2*nb3+j]  =w34[0]*tmp1+w6[0]*tmp2+w7[0]*tmp3;
            }
        }
    }
//stage 3 radix-2
b=0;
    i=1;
    nb1=nb1/stg[0+i];
    nb2=nb1/stg[1+i];
    nb3=nb2/stg[2+i];
    nb4=nb4*stg[0+i];

    for(int b=0;b<stg[1+i];b=b+1)
    {
        for(int j=0;j<nb3;j=j+1)
        {
            w1[0].real()= cos(nb4*b*(0*nb3+j)*tmp5);
            w1[0].imag()=-sin(nb4*b*(0*nb3+j)*tmp5);
            w2[0].real()= cos(nb4*b*(1*nb3+j)*tmp5);
            w2[0].imag()=-sin(nb4*b*(1*nb3+j)*tmp5);

            for(int i=0;i<nb4;i=i+1)
            {
                tmp1=w1[0]*tab[i*nb1+b*nb2+0*nb3+j];
                tmp2=w2[0]*tab[i*nb1+b*nb2+1*nb3+j];
                tab2[i*nb1+b*nb2+0*nb3+j]=tmp1+tmp2;
                tab2[i*nb1+b*nb2+1*nb3+j]=tmp1-tmp2;
            }
        }
    }

//stage 4 radix-4

    w50[0].real()=0;
    w50[0].imag()=-1;
    w60[0].real()=0;
    w60[0].imag()=1;
b=0;

    i=2;
    nb1=nb1/stg[0+i];
    nb2=nb1/stg[1+i];
    nb3=nb2/stg[2+i];
    nb4=nb4*stg[0+i];

    for(int b=0;b<stg[1+i];b=b+1)
    {
        for(int j=0;j<nb3;j=j+1)
        {
          w1[0].real()=cos(nb4*b*(0*nb3+j)*tmp5);
          w1[0].imag()=-sin(nb4*b*(0*nb3+j)*tmp5);
          w2[0].real()=cos(nb4*b*(1*nb3+j)*tmp5);
          w2[0].imag()=-sin(nb4*b*(1*nb3+j)*tmp5);
          w3[0].real()=cos(nb4*b*(2*nb3+j)*tmp5);
          w3[0].imag()=-sin(nb4*b*(2*nb3+j)*tmp5);
          w4[0].real()=cos(nb4*b*(3*nb3+j)*tmp5);
          w4[0].imag()=-sin(nb4*b*(3*nb3+j)*tmp5);

            for(int i=0;i<nb4;i=i+1)
            {
                tmp1=w1[0]*tab2[i*nb1+b*nb2+0*nb3+j];
                tmp2=w2[0]*tab2[i*nb1+b*nb2+1*nb3+j];
                tmp3=w3[0]*tab2[i*nb1+b*nb2+2*nb3+j];
                tmp4=w4[0]*tab2[i*nb1+b*nb2+3*nb3+j];
                //radix-4
                tab[i*nb1+b*nb2+0*nb3+j]   =tmp1+tmp2+tmp3+tmp4;
                tab[i*nb1+b*nb2+1*nb3+j]   =tmp1-tmp3+w50[0]*(tmp2-tmp4);
                tab[i*nb1+b*nb2+2*nb3+j]   =tmp1-tmp2+tmp3-tmp4;
                tab[i*nb1+b*nb2+3*nb3+j]   =tmp1-tmp3+w60[0]*(tmp2-tmp4);
            }
        }
    }

//stage 2 radix-3
b=0;
    i=3;
    nb1=nb1/stg[0+i];
    nb2=nb1/stg[1+i];
    nb3=nb2/stg[2+i];
    nb4=nb4*stg[0+i];

    for(int b=0;b<stg[1+i];b=b+1)
    {
        for(int j=0;j<nb3;j=j+1)
        {
          w1[0].real()=cos(nb4*b*(0*nb3+j)*tmp6);
          w1[0].imag()=-sin(nb4*b*(0*nb3+j)*tmp6);
          w2[0].real()=cos(nb4*b*(1*nb3+j)*tmp6);
          w2[0].imag()=-sin(nb4*b*(1*nb3+j)*tmp6);
          w3[0].real()=cos(nb4*b*(2*nb3+j)*tmp6);
          w3[0].imag()=-sin(nb4*b*(2*nb3+j)*tmp6);

            for(int i=0;i<nb4;i=i+1)
            {
              tmp1=w1[0]*tab[i*nb1+b*nb2+0*nb3+j];
              tmp2=w2[0]*tab[i*nb1+b*nb2+1*nb3+j];
              tmp3=w3[0]*tab[i*nb1+b*nb2+2*nb3+j];
              //radix-3
              tab2[i*nb1+b*nb2+0*nb3+j]   =w34[0]*tmp1+w34[0]*tmp2+w34[0]*tmp3;
              tab2[i*nb1+b*nb2+1*nb3+j]   =w34[0]*tmp1+w5[0]*tmp2+w6[0]*tmp3;
              tab2[i*nb1+b*nb2+2*nb3+j]   =w34[0]*tmp1+w6[0]*tmp2+w7[0]*tmp3;
            }
        }
    }

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