poniedziałek, 10 kwietnia 2017
my interpretation of FFT radix-4 N=4 algorithm c++ source code
//https://www.beechwood.eu/fft-implementation-r2-dit-r4-dif-r8-dif/
//my interpretation of FFT radix-4 algorithm c++ source code it works witch
// FFT=+/-cos(w+fi)-isin(w+fi) and maybe witch other functions that are mirror inverse
//author 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[]);
int b1,b2,b3;
int main()
{
std::complex<double>tab2[4]={};
std::complex<double>tab3[4]={};
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=4;
for(int i=0;i<4;i++)
{
tab2[i].real()=i*1.2+1;
tab2[i].imag()=i*3.2+1;
tab3[i].real()=i*1.2+1;
tab3[i].imag()=i*3.2+1;
}
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<4;i++)
{
tab2[i].real()=i*1.2+1;
tab2[i].imag()=i*3.2+1;
tab3[i].real()=i*1.2+1;
tab3[i].imag()=i*3.2+1;
}
fun_fourier_transform_DFT_method5_full_complex(N,tab2);
fun_fourier_transform_FFT_radix_4_N_4(N,tab3);
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<<"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;
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 =>";
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<<", b1="<<b1<<" b2="<<b2<<" b3="<<b3<<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;
for (int i=0;i<N;i++)
{
for(int j=0;j<N;j++)
{
//complex number method:
w[0].real()=-cos(i*j*zmienna1+0.11);
w[0].imag()=(-sin(i*j*zmienna1+0.11));
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[])
{
}
///////////////////
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 fi;
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;
fi=0.11;
a0.real()=nr1*cos(0*zm+fi);
a0.imag()=-sin(0*zm+fi);
a1.real()=nr1*cos(1*zm+fi);
a1.imag()=-sin(1*zm+fi);
a2.real()=nr1*cos(2*zm+fi);
a2.imag()=-sin(2*zm+fi);
a3.real()=nr1*cos(3*zm+fi);
a3.imag()=-sin(3*zm+fi);
a4.real()=nr1*cos(4*zm+fi);
a4.imag()=-sin(4*zm+fi);
a6.real()=nr1*cos(6*zm+fi);
a6.imag()=-sin(6*zm+fi);
a9.real()=nr1*cos(9*zm+fi);
a9.imag()=-sin(9*zm+fi);
//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+fi);
w0[0].imag()=-sin(0+fi);
w2[0].real()=nr2*cos(0*0*2*pi/(N/4)+fi);
w2[0].imag()=-sin(0*0*2*pi/(N/4)+fi);
tab2[0]=(-a0*tab[0]-a0*tab[1]-a0*tab[2]-a0*tab[3])*w0[0]*w2[0];
w1[0].real()=nr3*cos(0*wbm+fi);
w1[0].imag()=-sin(0*wbm+fi);
tab2[1]=(-a0*tab[0]-a1*tab[1]-a2*tab[2]-a3*tab[3])*w1[0]*w2[0];
w1[0].real()=nr3*cos(2*0*wbm+fi);
w1[0].imag()=-sin(2*0*wbm+fi);
tab2[2]=(-a0*tab[0]-a2*tab[1]-a4*tab[2]-a6*tab[3])*w1[0]*w2[0];
w1[0].real()=nr3*cos(3*0*wbm+fi);
w1[0].imag()=-sin(3*0*wbm+fi);
tab2[3]=(-a0*tab[0]-a3*tab[1]-a6*tab[2]-a9*tab[3])*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;
}
}
Subskrybuj:
Komentarze do posta (Atom)
Brak komentarzy:
Prześlij komentarz