piątek, 7 kwietnia 2017

fourier transform iFFT radix-4 for N=64 algorithm c++ source code


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

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


//fourier transform iFFT radix-4 for N=64 algorithm c++ source code
#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_64_official(int N,std::complex<double> tab[]);
void fun_inverse_fourier_transform_FFT_radix_4_N_64_official(int N,std::complex<double> tab[]);


void fun_fourier_transform_FFT_radix_4_N_64_method2(int N,std::complex<double> tab[]);
void fun_inverse_fourier_transform_FFT_radix_4_N_64_method2(int N,std::complex<double> tab[]);
void fun_fourier_transform_FFT_radix_4_N_64_method3(int N,std::complex<double> tab[]);
void fun_fourier_transform_FFT_radix_4_N_64_method4(int N,std::complex<double> tab[]);

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=64;
    std::complex<double> tab2[64]=
    {{1,65},{2,66},{3,67},{4,68},
    {5,69},{6,70},{7,71},{8,72},
    {9,73},{10,74},{11,75},{12,76},
    {13,77},{14,78},{15,79},{16,80},

    {17,81},{18,82},{19,83},{20,84},
    {21,85},{22,86},{23,87},{24,88},
    {25,89},{26,90},{27,91},{28,92},
    {29,93},{30,94},{31,95},{32,96},

    {33,97},{34,98},{35,99},{36,100},
    {37,101},{38,102},{39,103},{40,104},
    {41,105},{42,106},{43,107},{44,108},
    {45,109},{46,110},{47,111},{48,112},

    {49,113},{50,114},{51,115},{52,116},
    {53,117},{54,118},{55,119},{56,120},
    {57,121},{58,122},{59,123},{60,124},
    {61,125},{62,126},{63,127},{64,128}
    };


    double time2;
    double zmienna=0;
    /*
    std::fstream plik;
    plik.open("test.txt", std::ios::in | std::ios::out);
    if( plik.good() == false )
    {
        cout<<"nie otwarto pliku"<<endl;
        system("pause");
    }
   */
    cout<<"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;

    clock_t start = clock();

    //fun_fourier_transform_FFT_radix_4_N_64_method2(N,tab2);
    //fun_fourier_transform_FFT_radix_4_N_64_method3(N,tab2);
    //fun_fourier_transform_FFT_radix_4_N_64_method4(N,tab2);

    //////////////////////////////////////////////////////////
    fun_fourier_transform_FFT_radix_4_N_64_official(N,tab2);
    fun_inverse_bits_radix_4(N,tab2);
    ///////////////////////////////////////////////////////////

    time2=diffclock( start, clock() );

    cout<<"frequency Hz"<<endl;
    zmienna=0;
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
   cout<<round(tab2[j].real()*1000)/1000<<"  ";
    zmienna=zmienna+fabs(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<<"  ";
    zmienna=zmienna+fabs(round(tab2[j].imag()*1000)/1000);
    }
    cout<<endl;
    cout<<endl;

 // plik.close();
    system("pause");
    cout<<endl;
    //fun_inverse_fourier_transform_FFT_radix_4_N_64_method2(N,tab2);

    /////////////////////////////////////////////////////////////////
    fun_inverse_fourier_transform_FFT_radix_4_N_64_official(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;
    for(int j=0;j<N;j++)
    {
    cout.precision(4);
    cout<<round(tab2[j].imag()*1000)/1000<<"  ";
    }
    cout<<endl;
    cout<<endl;
    cout<<endl;

    cout<<endl;
    system("pause");
    return 0;
}


void fun_fourier_transform_FFT_radix_4_N_64_method2(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double> tab2[64]={};    // tab2[]==N
    std::complex<double>  w1[1]={{1,1}};
    std::complex<double>  w2[1]={{1,1}};
    std::complex<double>  w3[1]={{1,1}};
/*
for(int j=0;j<N;j++)
{
   cout<< tab[j]<<"  ";
}
cout<<endl<<"........."<<endl;
*/
    //radix-4
    for(int k=0;k<16;k++)
    {
        for(int n=0;n<16;n++)
        {
          w2[0].real()=cos(n*k*2*pi/(N/4));
          w2[0].imag()=-sin(n*k*2*pi/(N/4));

          tab2[4*k+0]=tab2[4*k+0]+(tab[n]+tab[n+N/4]+tab[n+N/2]+tab[n+3*N/4])*w2[0];

          w1[0].real()=cos(n*2*pi/N);
          w1[0].imag()=-sin(n*2*pi/N);
          w3[0].real()=0;
          w3[0].imag()=-1;
          tab2[4*k+1]=tab2[4*k+1]+(tab[n]-tab[n+N/2]+w3[0]*(tab[n+N/4]-tab[n+3*N/4]))*w1[0]*w2[0];

          w1[0].real()=cos(2*n*2*pi/N);
          w1[0].imag()=-sin(2*n*2*pi/N);
          tab2[4*k+2]=tab2[4*k+2]+(tab[n]-tab[n+N/4]+tab[n+N/2]-tab[n+3*N/4])*w1[0]*w2[0];

          w1[0].real()=cos(3*n*2*pi/N);
          w1[0].imag()=-sin(3*n*2*pi/N);
          w3[0].real()=0;
          w3[0].imag()=1;
          tab2[4*k+3]=tab2[4*k+3]+(tab[n]-tab[n+N/2]+w3[0]*(tab[n+N/4]-tab[n+3*N/4]))*w1[0]*w2[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_inverse_fourier_transform_FFT_radix_4_N_64_method2(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double> tab2[64]={};    // tab2[]==N
    std::complex<double>  w1[1]={{1,1}};
    std::complex<double>  w2[1]={{1,1}};
    std::complex<double>  w3[1]={{1,1}};
/*
for(int j=0;j<N;j++)
{
   cout<< tab[j]<<"  ";
}
cout<<endl<<"........."<<endl;
*/
          //radix-4
    for(int k=0;k<16;k++)
    {
        for(int n=0;n<16;n++)
        {
          w2[0].real()=cos(n*k*2*pi/(N/4));
          w2[0].imag()=sin(n*k*2*pi/(N/4));

          tab2[4*k+0]=tab2[4*k+0]+(tab[n]+tab[n+N/4]+tab[n+N/2]+tab[n+3*N/4])*w2[0];

          w1[0].real()=cos(n*2*pi/N);
          w1[0].imag()=sin(n*2*pi/N);
          w3[0].real()=0;
          w3[0].imag()=1;
          tab2[4*k+1]=tab2[4*k+1]+(tab[n]-tab[n+N/2]+w3[0]*(tab[n+N/4]-tab[n+3*N/4]))*w1[0]*w2[0];

          w1[0].real()=cos(2*n*2*pi/N);
          w1[0].imag()=sin(2*n*2*pi/N);
          tab2[4*k+2]=tab2[4*k+2]+(tab[n]-tab[n+N/4]+tab[n+N/2]-tab[n+3*N/4])*w1[0]*w2[0];

          w1[0].real()=cos(3*n*2*pi/N);
          w1[0].imag()=sin(3*n*2*pi/N);
          w3[0].real()=0;
          w3[0].imag()=-1;
          tab2[4*k+3]=tab2[4*k+3]+(tab[n]-tab[n+N/2]+w3[0]*(tab[n+N/4]-tab[n+3*N/4]))*w1[0]*w2[0];
        }
    }

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


void fun_inverse_bits_radix_4(int N,std::complex<double> tab[])
{

}

void fun_fourier_transform_FFT_radix_4_N_64_method3(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double> tab2[64]={};    // tab2[]==N
    std::complex<double>  w1[1]={{1,1}};
    std::complex<double>  w2[1]={{1,1}};
    std::complex<double>  w3[1]={{1,1}};
/*
for(int j=0;j<N;j++)
{
   cout<< tab[j]<<"  ";
}
cout<<endl<<"........."<<endl;
*/
    //radix-4
    for(int k=0;k<16;k++)
    {
        for(int n=0;n<16;n++)
        {
          w2[0].real()=cos(n*k*2*pi/(N/4));
          w2[0].imag()=-sin(n*k*2*pi/(N/4));

          tab2[4*k+0]=tab2[4*k+0]+(tab[n]+tab[n+N/4]+tab[n+N/2]+tab[n+3*N/4])*w2[0];

          w1[0].real()=cos(n*2*pi/N);
          w1[0].imag()=-sin(n*2*pi/N);
          w3[0].real()=0;
          w3[0].imag()=-1;
          tab2[4*k+1]=tab2[4*k+1]+(w3[0]*(tab[n]-tab[n+N/2])+(tab[n+N/4]-tab[n+3*N/4]))*w1[0]*w2[0];

          w1[0].real()=cos(2*n*2*pi/N);
          w1[0].imag()=-sin(2*n*2*pi/N);
          tab2[4*k+2]=tab2[4*k+2]+(tab[n]-tab[n+N/4]+tab[n+N/2]-tab[n+3*N/4])*w1[0]*w2[0];

          w1[0].real()=cos(3*n*2*pi/N);
          w1[0].imag()=-sin(3*n*2*pi/N);
          w3[0].real()=0;
          w3[0].imag()=1;
          tab2[4*k+3]=tab2[4*k+3]+(w3[0]*(tab[n]-tab[n+N/2])+(tab[n+N/4]-tab[n+3*N/4]))*w1[0]*w2[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_64_method4(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double> tab2[64]={};    // tab2[]==N
    std::complex<double>  w1[1]={{1,1}};
    std::complex<double>  w2[1]={{1,1}};
    std::complex<double>  w3[1]={{1,1}};
/*
for(int j=0;j<N;j++)
{
   cout<< tab[j]<<"  ";
}
cout<<endl<<"........."<<endl;
*/

        for(int i=0;i<16;i++)
        {
          tab2[4*i]=(tab[0+i]+tab[N/4+i]+tab[N/2+i]+tab[3*N/4+i]);
        }

        for(int i=0;i<16;i++)
        {
          w1[0].real()=cos(i*2*pi/N);
          w1[0].imag()=-sin(i*2*pi/N);
          w2[0].real()=0;
          w2[0].imag()=-1;
          tab2[4*i+1]=(tab[0+i]-tab[i+N/2]+w2[0]*(tab[i+N/4]-tab[i+3*N/4]))*w1[0];
        }

        for(int i=0;i<16;i++)
        {
          w1[0].real()=cos(2*i*2*pi/N);
          w1[0].imag()=-sin(2*i*2*pi/N);
          tab2[4*i+2]=(tab[i]-tab[i+N/4]+tab[i+N/2]-tab[i+3*N/4])*w1[0];
        }

        for(int i=0;i<16;i++)
        {
          w1[0].real()=cos(3*i*2*pi/N);
          w1[0].imag()=-sin(3*i*2*pi/N);
          w3[0].real()=0;
          w3[0].imag()=1;
          tab2[4*i+3]=(tab[i]-tab[i+N/2]+w3[0]*(tab[i+N/4]-tab[i+3*N/4]))*w1[0];
        }

///////////////////////////////////
for(int k=0;k<N;k++)
{
    tab[k]={0,0};
}

        for(int i=0;i<16;i++)
        {
        for(int j=0;j<N;j=j+4)
        {
          w1[0].real()=cos(i*j/(4)*2*pi/(N/4));
          w1[0].imag()=-sin(i*j/(4)*2*pi/(N/4));
          tab[4*i]=tab[4*i]+tab2[j]*w1[0];
          tab[4*i+1]=tab[4*i+1]+tab2[j+1]*w1[0];
          tab[4*i+2]=tab[4*i+2]+tab2[j+2]*w1[0];
          tab[4*i+3]=tab[4*i+3]+tab2[j+3]*w1[0];
        }
        }

    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_64_official(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double> tab2[64]={};    // tab2[]==N
    std::complex<double>  w[1]={{1,1}};
    std::complex<double>  w1[1]={{1,1}};
    std::complex<double>  w2[1]={{1,1}};
    std::complex<double>  w3[1]={{1,1}};
    std::complex<double>  w4[1]={{1,1}};
    std::complex<double>  w5[1]={{1,1}};
    std::complex<double>  w6[1]={{1,1}};
/*
for(int j=0;j<N;j++)
{
   cout<< tab[j]<<"  ";
}
cout<<endl<<"........."<<endl;
*/
        w[0]={{0,1}};
        //radix-4
        for(int i=0;i<(N/4);i++)
        {
          tab2[i]     =(tab[i+0]+tab[i+N/4]+tab[i+N/2]+tab[i+3*N/4]);

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

          tab2[i+N/2] =(tab[i+0]-tab[i+N/4]+tab[i+N/2]-tab[i+3*N/4]);

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


for(int k=0;k<N;k++)
{
    tab[k]={0,0};
}

        for(int i=0;i<(N/16);i++)
        {
        for(int j=0;j<4;j++)
        {
          w1[0].real()=cos((0+i*j)*2*pi/(N/1));
          w1[0].imag()=sin((0+i*j)*2*pi/(N/1));
          w2[0].real()=cos((4*i+i*j)*2*pi/(N/1));
          w2[0].imag()=sin((4*i+i*j)*2*pi/(N/1));
          w3[0].real()=cos((8*i+i*j)*2*pi/(N/1));
          w3[0].imag()=sin((8*i+i*j)*2*pi/(N/1));
          w4[0].real()=cos((12*i+i*j)*2*pi/(N/1));
          w4[0].imag()=sin((12*i+i*j)*2*pi/(N/1));

          w5[0].real()=0;
          w5[0].imag()=1;
          tab[0+j+N/4*i]     =(w1[0]*tab2[0+j+N/4*i]+w2[0]*tab2[N/16+j+N/4*i]+       w3[0]*tab2[N/8 +j+N/4*i]+w4[0]*tab2[3*N/16+j+N/4*i]);
          tab[N/16+j+N/4*i]  =(w1[0]*tab2[0+j+N/4*i]-w3[0]*tab2[N/8 +j+N/4*i]+w5[0]*(w2[0]*tab2[N/16+j+N/4*i]-w4[0]*tab2[3*N/16+j+N/4*i]));
          tab[N/8+j+N/4*i]   =(w1[0]*tab2[0+j+N/4*i]-w2[0]*tab2[N/16+j+N/4*i]+       w3[0]*tab2[N/8 +j+N/4*i]-w4[0]*tab2[3*N/16+j+N/4*i]);
          w5[0].real()=0;
          w5[0].imag()=-1;
          tab[3*N/16+j+N/4*i]=(w1[0]*tab2[0+j+N/4*i]-w3[0]*tab2[N/8 +j+N/4*i]+w5[0]*(w2[0]*tab2[N/16+j+N/4*i]-w4[0]*tab2[3*N/16+j+N/4*i]));
        }
        }
///////////////////////////
for(int k=0;k<N;k++)
{
    tab2[k]={0,0};

}
        for(int i=0;i<(N/16);i++)
        {
        for(int j=0;j<4;j++)
        {
          w1[0].real()=cos(0*2*pi/(N/1));
          w1[0].imag()=sin(0*2*pi/(N/1));
          w2[0].real()=cos(1*4*j*2*pi/(N/1));
          w2[0].imag()=sin(1*4*j*2*pi/(N/1));
          w3[0].real()=cos(2*4*j*2*pi/(N/1));
          w3[0].imag()=sin(2*4*j*2*pi/(N/1));
          w4[0].real()=cos(3*4*j*2*pi/(N/1));
          w4[0].imag()=sin(3*4*j*2*pi/(N/1));

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

          tab2[N/4*i+0+4*j]=(w1[0]*tab[N/4*i+0+4*j]+w2[0]*tab[N/4*i+1+4*j]+       w3[0]*tab[N/4*i+2+4*j]+w4[0]*tab[N/4*i+3+4*j]);
          tab2[N/4*i+1+4*j]=(w1[0]*tab[N/4*i+0+4*j]-w3[0]*tab[N/4*i+2+4*j]+w5[0]*(w2[0]*tab[N/4*i+1+4*j]-w4[0]*tab[N/4*i+3+4*j]));
          tab2[N/4*i+2+4*j]=(w1[0]*tab[N/4*i+0+4*j]-w2[0]*tab[N/4*i+1+4*j]+       w3[0]*tab[N/4*i+2+4*j]-w4[0]*tab[N/4*i+3+4*j]);
          w5[0].real()=0;
          w5[0].imag()=-1;
          tab2[N/4*i+3+4*j]=(w1[0]*tab[N/4*i+0+4*j]-w3[0]*tab[N/4*i+2+4*j]+w5[0]*(w2[0]*tab[N/4*i+1+4*j]-w4[0]*tab[N/4*i+3+4*j]));
        }
        }

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







void fun_fourier_transform_FFT_radix_4_N_64_official(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double> tab2[64]={};    // tab2[]==N
    std::complex<double>  w[1]={{1,1}};
    std::complex<double>  w1[1]={{1,1}};
    std::complex<double>  w2[1]={{1,1}};
    std::complex<double>  w3[1]={{1,1}};
    std::complex<double>  w4[1]={{1,1}};
    std::complex<double>  w5[1]={{1,1}};
    std::complex<double>  w6[1]={{1,1}};
/*
for(int j=0;j<N;j++)
{
   cout<< tab[j]<<"  ";
}
cout<<endl<<"........."<<endl;
*/
        w[0]={{0,1}};
        //radix-4
        for(int i=0;i<(N/4);i++)
        {
          tab2[i]     =(tab[i+0]+tab[i+N/4]+tab[i+N/2]+tab[i+3*N/4]);

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

          tab2[i+N/2] =(tab[i+0]-tab[i+N/4]+tab[i+N/2]-tab[i+3*N/4]);

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


for(int k=0;k<N;k++)
{
    tab[k]={0,0};
}

        for(int i=0;i<(N/16);i++)
        {
        for(int j=0;j<4;j++)
        {
          w1[0].real()=cos((0+i*j)*2*pi/(N/1));
          w1[0].imag()=-sin((0+i*j)*2*pi/(N/1));
          w2[0].real()=cos((4*i+i*j)*2*pi/(N/1));
          w2[0].imag()=-sin((4*i+i*j)*2*pi/(N/1));
          w3[0].real()=cos((8*i+i*j)*2*pi/(N/1));
          w3[0].imag()=-sin((8*i+i*j)*2*pi/(N/1));
          w4[0].real()=cos((12*i+i*j)*2*pi/(N/1));
          w4[0].imag()=-sin((12*i+i*j)*2*pi/(N/1));

          w5[0].real()=0;
          w5[0].imag()=-1;
          tab[0+j+N/4*i]     =(w1[0]*tab2[0+j+N/4*i]+w2[0]*tab2[N/16+j+N/4*i]+       w3[0]*tab2[N/8 +j+N/4*i]+w4[0]*tab2[3*N/16+j+N/4*i]);
          tab[N/16+j+N/4*i]  =(w1[0]*tab2[0+j+N/4*i]-w3[0]*tab2[N/8 +j+N/4*i]+w5[0]*(w2[0]*tab2[N/16+j+N/4*i]-w4[0]*tab2[3*N/16+j+N/4*i]));
          tab[N/8+j+N/4*i]   =(w1[0]*tab2[0+j+N/4*i]-w2[0]*tab2[N/16+j+N/4*i]+       w3[0]*tab2[N/8 +j+N/4*i]-w4[0]*tab2[3*N/16+j+N/4*i]);
          w5[0].real()=0;
          w5[0].imag()=1;
          tab[3*N/16+j+N/4*i]=(w1[0]*tab2[0+j+N/4*i]-w3[0]*tab2[N/8 +j+N/4*i]+w5[0]*(w2[0]*tab2[N/16+j+N/4*i]-w4[0]*tab2[3*N/16+j+N/4*i]));
        }
        }
///////////////////////////
for(int k=0;k<N;k++)
{
    tab2[k]={0,0};

}

        for(int i=0;i<(N/16);i++)
        {
        for(int j=0;j<4;j++)
        {
          w1[0].real()=cos(0*2*pi/(N/1));
          w1[0].imag()=-sin(0*2*pi/(N/1));
          w2[0].real()=cos(1*4*j*2*pi/(N/1));
          w2[0].imag()=-sin(1*4*j*2*pi/(N/1));
          w3[0].real()=cos(2*4*j*2*pi/(N/1));
          w3[0].imag()=-sin(2*4*j*2*pi/(N/1));
          w4[0].real()=cos(3*4*j*2*pi/(N/1));
          w4[0].imag()=-sin(3*4*j*2*pi/(N/1));
          w5[0].real()=0;
          w5[0].imag()=-1;

          tab2[N/4*i+0+4*j]=(w1[0]*tab[N/4*i+0+4*j]+w2[0]*tab[N/4*i+1+4*j]+       w3[0]*tab[N/4*i+2+4*j]+w4[0]*tab[N/4*i+3+4*j]);
          tab2[N/4*i+1+4*j]=(w1[0]*tab[N/4*i+0+4*j]-w3[0]*tab[N/4*i+2+4*j]+w5[0]*(w2[0]*tab[N/4*i+1+4*j]-w4[0]*tab[N/4*i+3+4*j]));
          tab2[N/4*i+2+4*j]=(w1[0]*tab[N/4*i+0+4*j]-w2[0]*tab[N/4*i+1+4*j]+       w3[0]*tab[N/4*i+2+4*j]-w4[0]*tab[N/4*i+3+4*j]);
          w5[0].real()=0;
          w5[0].imag()=1;
          tab2[N/4*i+3+4*j]=(w1[0]*tab[N/4*i+0+4*j]-w3[0]*tab[N/4*i+2+4*j]+w5[0]*(w2[0]*tab[N/4*i+1+4*j]-w4[0]*tab[N/4*i+3+4*j]));
        }
        }

    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