fast fourier transform FFT radix-3 for N=9 points example

//fast fourier transform FFT radix-3 for N=9 points example


void fun_fourier_transform_FFT_radix_3_N_9(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>  w7[1]={{1,0}};
    std::complex<double> tmp1,tmp2,tmp3,tmp4;
    double tmp5;
    double tmp6;
    tmp5=2*pi/(N/1);
    tmp6=2*pi/(3/1);


          w4[0].real()=cos(0*tmp6);
          w4[0].imag()=-sin(0*tmp6);
          w5[0].real()=cos(1*tmp6);
          w5[0].imag()=-sin(1*tmp6);
          w6[0].real()=cos(2*tmp6);
          w6[0].imag()=-sin(2*tmp6);
          w7[0].real()=cos(4*tmp6);
          w7[0].imag()=-sin(4*tmp6);



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

          tmp1=w1[0]*tab[i+0];
          tmp2=w2[0]*tab[i+3];
          tmp3=w3[0]*tab[i+6];
         //radix-4
          tab2[i+0]  =w4[0]*tmp1+w4[0]*tmp2+w4[0]*tmp3;
          tab2[i+3]  =w4[0]*tmp1+w5[0]*tmp2+w6[0]*tmp3;
          tab2[i+6]  =w4[0]*tmp1+w6[0]*tmp2+w7[0]*tmp3;
        }
//stage 2 radix-3

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

    int i=0;
          tmp1=w1[0]*tab2[i+0];
          tmp2=w2[0]*tab2[i+1];
          tmp3=w3[0]*tab2[i+2];
         //radix-4
          tab[i+0]   =w4[0]*tmp1+w4[0]*tmp2+w4[0]*tmp3;
          tab[i+1]   =w4[0]*tmp1+w5[0]*tmp2+w6[0]*tmp3;
          tab[i+2]   =w4[0]*tmp1+w6[0]*tmp2+w7[0]*tmp3;

          w1[0].real()=cos(0*tmp5);
          w1[0].imag()=-sin(0*tmp5);
          w2[0].real()=cos(1*tmp5);
          w2[0].imag()=-sin(1*tmp5);
          w3[0].real()=cos(2*tmp5);
          w3[0].imag()=-sin(2*tmp5);

     i=3;
          tmp1=w1[0]*tab2[i+0];
          tmp2=w2[0]*tab2[i+1];
          tmp3=w3[0]*tab2[i+2];
         //radix-4
          tab[i+0]   =w4[0]*tmp1+w4[0]*tmp2+w4[0]*tmp3;
          tab[i+1]   =w4[0]*tmp1+w5[0]*tmp2+w6[0]*tmp3;
          tab[i+2]   =w4[0]*tmp1+w6[0]*tmp2+w7[0]*tmp3;

          w1[0].real()=cos(0*tmp5);
          w1[0].imag()=-sin(0*tmp5);
          w2[0].real()=cos(2*tmp5);
          w2[0].imag()=-sin(2*tmp5);
          w3[0].real()=cos(4*tmp5);
          w3[0].imag()=-sin(4*tmp5);

     i=6;
          tmp1=w1[0]*tab2[i+0];
          tmp2=w2[0]*tab2[i+1];
          tmp3=w3[0]*tab2[i+2];
         //radix-4
          tab[i+0]   =w4[0]*tmp1+w4[0]*tmp2+w4[0]*tmp3;
          tab[i+1]   =w4[0]*tmp1+w5[0]*tmp2+w6[0]*tmp3;
          tab[i+2]   =w4[0]*tmp1+w6[0]*tmp2+w7[0]*tmp3;

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

}