my interpretation of FFT b*radix-4 for N=16 algorithm c++ source code

//my interpretation of FFT b*radix-4 for N=16 algorithm c++ source code




////https://www.beechwood.eu/fft-implementation-r2-dit-r4-dif-r8-dif/

//my interpretation of FFT 16*radix-4 for N=16 algorithm c++ source code it works witch
// FFT=+/-cos(w+fi)-isin(w+fi) and maybe witch other functions that are mirror inverse
//copyrights marcin matysek (r)ewertyn.PL
#include <iostream>
#include "conio.h"
#include <stdlib.h>
#include <math.h>
#include <complex>
#include <fstream>
using namespace std;

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

int b1,b2,b3;
double fi=-0.41;


int main()
{
    std::complex<double>tab2[16]={};
    std::complex<double>tab3[16]={};

    std::fstream plik;

    plik.open("test.txt", std::ios::in | std::ios::out);

/*
    if( plik.good() == false )
    {

        cout<<"create file : test.txt"<<endl;
        system("pause");

    }
*/

    int N=16;

    for(int i=0;i<N;i++)
    {
        tab2[i].real()=i*1.2+1;
        tab2[i].imag()=i*3.2+11;
        tab3[i].real()=i*1.2+1;
        tab3[i].imag()=i*3.2+11;
    }


    cout<<endl;
cout<<"signal g(x)"<<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 g(x)"<<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;


    for(int i1=0;i1<2;i1++)
    {
        if(i1==0){b1=0;}
        if(i1==1){b1=1;}
    for(int i2=0;i2<2;i2++)
    {
        if(i2==0){b2=0;}
        if(i2==1){b2=1;}
    for(int i3=0;i3<2;i3++)
    {
        if(i3==0){b3=0;}
        if(i3==1){b3=1;}


    for(int i=0;i<N;i++)
    {
        tab2[i].real()=i*1.2+1;
        tab2[i].imag()=i*3.2+11;
        tab3[i].real()=i*1.2+1;
        tab3[i].imag()=i*3.2+11;
    }


    fun_fourier_transform_DFT_method5_full_complex(N,tab2);
    fun_fourier_transform_FFT_radix_4_N_4(N,tab3);


    cout<<endl;
    cout<<endl;
    cout<<endl;
    cout<<endl;
    cout<<"frequency Hz DFT:"<<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;

    cout<<"frequency Hz radix4:"<<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<<endl;

        if((((fabs(tab3[0].real())-fabs(tab2[0].real())>=-0.03)&&(fabs(tab3[0].real())-fabs(tab2[0].real())<=0.03))&&
           ((fabs(tab3[1].real())-fabs(tab2[1].real())>=-0.03)&&(fabs(tab3[1].real())-fabs(tab2[1].real())<=0.03))&&
           ((fabs(tab3[2].real())-fabs(tab2[2].real())>=-0.03)&&(fabs(tab3[2].real())-fabs(tab2[2].real())<=0.03))&&
           ((fabs(tab3[3].real())-fabs(tab2[3].real())>=-0.03)&&(fabs(tab3[3].real())-fabs(tab2[3].real())<=0.03)))
          )
        {
            cout<<"DFT==FFT =>";
            cout<<endl;
        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<<" . ";
            }
        }
        cout<<endl<<endl;
        cout<<", b1="<<b1<<" b2="<<b2<<" b3="<<b3<<endl<<endl;system("pause");
        //    plik<<a1<<" "<<a2<<" "<<a3<<" "<<a4<<" "<<a5<<" "<<a11<<" , "<<a12<<" "<<a13<<" "<<a14<<" "<<a15<<" "<<a16<<" "<<endl;
       }
   }}}
    //cout<<"end:"<<endl;
    cout<<endl;
    cout<<endl;
    //plik.close();
    system("pause");
    return 0;
}
///////////////
void fun_fourier_transform_DFT_method5_full_complex(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double> tab2[16]={};    // 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_inverse_bits_radix_4(int N,std::complex<double> tab[])
{
    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_FFT_radix_4_N_4(int N,std::complex<double> tab[])
{
    const double pi=3.141592653589793238462;
    std::complex<double> tab2[64]={};    // tab2[]==N

    std::complex<double>  w0[1]={{1,0}};
    std::complex<double>  w1[1]={{1,0}};
    std::complex<double>  w2[1]={{1,0}};
    std::complex<double>  w3[1]={{1,0}};
    std::complex<double>  a0,a1,a2,a3,a4,a6,a9;

    double zm=2*pi/4;
    double wbm=2*pi/N;
    double fi1;

    int nr1=1;
    int nr2=1;
    int nr3=1;

        if(b1==0){nr1=-1;}
        else if(b1==1){nr1=1;}
        else{}
        if(b2==0){nr2=-1;}
        else if(b2==1){nr2=1;}
        else{}
        if(b3==0){nr3=-1;}
        else if(b3==1){nr3=1;}
        else{}
        //cout<<"b1="<<b1<<endl<<endl;

   // nr1=1;
   // nr2=-1;
   // nr3=1;

    fi1=fi;

          a0.real()=nr1*cos(0*zm+fi1);
          a0.imag()=-sin(0*zm+fi1);
          a1.real()=nr1*cos(1*zm+fi1);
          a1.imag()=-sin(1*zm+fi1);
          a2.real()=nr1*cos(2*zm+fi1);
          a2.imag()=-sin(2*zm+fi1);
          a3.real()=nr1*cos(3*zm+fi1);
          a3.imag()=-sin(3*zm+fi1);
          a4.real()=nr1*cos(4*zm+fi1);
          a4.imag()=-sin(4*zm+fi1);
          a6.real()=nr1*cos(6*zm+fi1);
          a6.imag()=-sin(6*zm+fi1);
          a9.real()=nr1*cos(9*zm+fi1);
          a9.imag()=-sin(9*zm+fi1);

        //cout.precision(3);
        //cout<<" ["<<round(a0.real()*1000)/1000<<" "<<round(a0.imag()*1000)/1000<<"i] ["<<round(a0.real()*1000)/1000<<" "<<round(a0.imag()*1000)/1000<<"i] ["<<round(a0.real()*1000)/1000<<" "<<round(a0.imag()*1000)/1000<<"i] ["<<round(a0.real()*1000)/1000<<" "<<round(a0.imag()*1000)/1000<<"i] "<<endl;
       //// cout<<" ["<<round(a0.real()*1000)/1000<<" "<<round(a0.imag()*1000)/1000<<"i] ["<<round(a1.real()*1000)/1000<<" "<<round(a1.imag()*1000)/1000<<"i] ["<<round(a2.real()*1000)/1000<<" "<<round(a2.imag()*1000)/1000<<"i] ["<<round(a3.real()*1000)/1000<<" "<<round(a3.imag()*1000)/1000<<"i] "<<endl;
       // cout<<" ["<<round(a0.real()*1000)/1000<<" "<<round(a0.imag()*1000)/1000<<"i] ["<<round(a2.real()*1000)/1000<<" "<<round(a2.imag()*1000)/1000<<"i] ["<<round(a4.real()*1000)/1000<<" "<<round(a4.imag()*1000)/1000<<"i] ["<<round(a6.real()*1000)/1000<<" "<<round(a6.imag()*1000)/1000<<"i] "<<endl;
        //cout<<" ["<<round(a0.real()*1000)/1000<<" "<<round(a0.imag()*1000)/1000<<"i] ["<<round(a3.real())<<" "<<round(a3.imag())<<"i] ["<<round(a6.real())<<" "<<round(a6.imag())<<"i] ["<<round(a9.real())<<" "<<round(a9.imag())<<"i] "<<endl;
       //system("pause");

    //radix-4

          w0[0].real()=nr3*cos(0+fi1);
          w0[0].imag()=-sin(0+fi1);

    for(int k=0;k<4;k++)
    {
        for(int n=0;n<4;n++)
        {

          w2[0].real()=nr2*cos(n*k*2*pi/(N/4)+fi1);
          w2[0].imag()=-sin(n*k*2*pi/(N/4)+fi1);

          tab2[4*k+0]=tab2[4*k+0]+(-a0*tab[n]-a0*tab[n+4]-a0*tab[n+8]-a0*tab[n+12])*w0[0]*w2[0];

          w1[0].real()=nr3*cos(n*wbm+fi1);
          w1[0].imag()=-sin(n*wbm+fi1);

          tab2[4*k+1]=tab2[4*k+1]+(-a0*tab[n]-a1*tab[n+4]-a2*tab[n+8]-a3*tab[n+12])*w1[0]*w2[0];

          w1[0].real()=nr3*cos(2*n*wbm+fi1);
          w1[0].imag()=-sin(2*n*wbm+fi1);
          tab2[4*k+2]=tab2[4*k+2]+(-a0*tab[n]-a2*tab[n+4]-a4*tab[n+8]-a6*tab[n+12])*w1[0]*w2[0];

          w1[0].real()=nr3*cos(3*n*wbm+fi1);
          w1[0].imag()=-sin(3*n*wbm+fi1);

          tab2[4*k+3]=tab2[4*k+3]+(-a0*tab[n]-a3*tab[n+4]-a6*tab[n+8]-a9*tab[n+12])*w1[0]*w2[0];
        }
    }
/////////////////////////////
    //inverse

    std::complex<double> tmp;
    tmp=tab2[4];
    tab2[4]=tab2[12];
    tab2[12]=tmp;

    tmp=tab2[4+1];
    tab2[4+1]=tab2[12+1];
    tab2[12+1]=tmp;

    tmp=tab2[4+2];
    tab2[4+2]=tab2[12+2];
    tab2[12+2]=tmp;

    tmp=tab2[4+3];
    tab2[4+3]=tab2[12+3];
    tab2[12+3]=tmp;


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

}