niedziela, 18 czerwca 2017

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


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

void fun_fourier_transform_FFT_mixed_radix_4(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>  w50[1]={{1,0}};
    std::complex<double>  w60[1]={{1,0}};
    std::complex<double> tmp1,tmp2,tmp3,tmp4;
    double tmp5;
    double tmp6;
    double tmp7;

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

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


          w4[0].real()=cos(0*tmp7);
          w4[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);

//stage 1  radix-4
          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);

        for(int i=0;i<(N/4);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];
         //radix-4
          tab2[i]       =tmp1+tmp2+tmp3+tmp4;
          tab2[i+N/4]   =tmp1-tmp3+w50[0]*(tmp2-tmp4);
          tab2[i+N/2]   =tmp1-tmp2+tmp3-tmp4;
          tab2[i+3*N/4] =tmp1-tmp3+w60[0]*(tmp2-tmp4);
        }
//stage 2 radix-3
    int i=0;
    int b=0;

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

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

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

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

     i=6;
     b=1;

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

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

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

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

    i=12;
    b=2;

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

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

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

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

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

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

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

          tmp1=w1[0]*tab2[i+1];
          tmp2=w2[0]*tab2[i+3];
          tmp3=w3[0]*tab2[i+5];
         //radix-3
          tab[i+1]   =w4[0]*tmp1+w4[0]*tmp2+w4[0]*tmp3;
          tab[i+3]   =w4[0]*tmp1+w5[0]*tmp2+w6[0]*tmp3;
          tab[i+5]   =w4[0]*tmp1+w6[0]*tmp2+w7[0]*tmp3;
//stage 3  radix-2
 //////////////////////////////////////
        i=0;
                w1[0].real()= cos(0*tmp6);
                w1[0].imag()=-sin(0*tmp6);
                w2[0].real()= cos(0*tmp6);
                w2[0].imag()=-sin(0*tmp6);
                tmp1=w1[0]*tab[0+i];
                tmp2=w2[0]*tab[1+i];
                tab2[0+i]=tmp1+tmp2;
                tab2[1+i]=tmp1-tmp2;

                w1[0].real()= cos(0*tmp6);
                w1[0].imag()=-sin(0*tmp6);
                w2[0].real()= cos(4*tmp6);
                w2[0].imag()=-sin(4*tmp6);
                tmp1=w1[0]*tab[2+i];
                tmp2=w2[0]*tab[3+i];
                tab2[2+i]=tmp1+tmp2;
                tab2[3+i]=tmp1-tmp2;

                w1[0].real()= cos(0);
                w1[0].imag()=-sin(0);
                w2[0].real()= cos(8*tmp6);
                w2[0].imag()=-sin(8*tmp6);
                tmp1=w1[0]*tab[4+i];
                tmp2=w2[0]*tab[5+i];
                tab2[4+i]=tmp1+tmp2;
                tab2[5+i]=tmp1-tmp2;

        i=6;
                w1[0].real()= cos(0*tmp6);
                w1[0].imag()=-sin(0*tmp6);
                w2[0].real()= cos(0*tmp6);
                w2[0].imag()=-sin(0*tmp6);
                tmp1=w1[0]*tab[0+i];
                tmp2=w2[0]*tab[1+i];
                tab2[0+i]=tmp1+tmp2;
                tab2[1+i]=tmp1-tmp2;

                w1[0].real()= cos(0*tmp6);
                w1[0].imag()=-sin(0*tmp6);
                w2[0].real()= cos(4*tmp6);
                w2[0].imag()=-sin(4*tmp6);
                tmp1=w1[0]*tab[2+i];
                tmp2=w2[0]*tab[3+i];
                tab2[2+i]=tmp1+tmp2;
                tab2[3+i]=tmp1-tmp2;

                w1[0].real()= cos(0);
                w1[0].imag()=-sin(0);
                w2[0].real()= cos(8*tmp6);
                w2[0].imag()=-sin(8*tmp6);
                tmp1=w1[0]*tab[4+i];
                tmp2=w2[0]*tab[5+i];
                tab2[4+i]=tmp1+tmp2;
                tab2[5+i]=tmp1-tmp2;

        i=12;
                w1[0].real()= cos(0*tmp6);
                w1[0].imag()=-sin(0*tmp6);
                w2[0].real()= cos(0*tmp6);
                w2[0].imag()=-sin(0*tmp6);
                tmp1=w1[0]*tab[0+i];
                tmp2=w2[0]*tab[1+i];
                tab2[0+i]=tmp1+tmp2;
                tab2[1+i]=tmp1-tmp2;

                w1[0].real()= cos(0*tmp6);
                w1[0].imag()=-sin(0*tmp6);
                w2[0].real()= cos(4*tmp6);
                w2[0].imag()=-sin(4*tmp6);
                tmp1=w1[0]*tab[2+i];
                tmp2=w2[0]*tab[3+i];
                tab2[2+i]=tmp1+tmp2;
                tab2[3+i]=tmp1-tmp2;

                w1[0].real()= cos(0);
                w1[0].imag()=-sin(0);
                w2[0].real()= cos(8*tmp6);
                w2[0].imag()=-sin(8*tmp6);
                tmp1=w1[0]*tab[4+i];
                tmp2=w2[0]*tab[5+i];
                tab2[4+i]=tmp1+tmp2;
                tab2[5+i]=tmp1-tmp2;

        i=18;
                w1[0].real()= cos(0*tmp6);
                w1[0].imag()=-sin(0*tmp6);
                w2[0].real()= cos(0*tmp6);
                w2[0].imag()=-sin(0*tmp6);
                tmp1=w1[0]*tab[0+i];
                tmp2=w2[0]*tab[1+i];
                tab2[0+i]=tmp1+tmp2;
                tab2[1+i]=tmp1-tmp2;

                w1[0].real()= cos(0*tmp6);
                w1[0].imag()=-sin(0*tmp6);
                w2[0].real()= cos(4*tmp6);
                w2[0].imag()=-sin(4*tmp6);
                tmp1=w1[0]*tab[2+i];
                tmp2=w2[0]*tab[3+i];
                tab2[2+i]=tmp1+tmp2;
                tab2[3+i]=tmp1-tmp2;

                w1[0].real()= cos(0);
                w1[0].imag()=-sin(0);
                w2[0].real()= cos(8*tmp6);
                w2[0].imag()=-sin(8*tmp6);
                tmp1=w1[0]*tab[4+i];
                tmp2=w2[0]*tab[5+i];
                tab2[4+i]=tmp1+tmp2;
                tab2[5+i]=tmp1-tmp2;

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


Brak komentarzy:

Prześlij komentarz