//inverse fourier transform DFT and 1 stage FFT
//author marcin matysek (r)ewertyn.PL
#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_fourier_transform_DFT_method5_full_complex(int N,std::complex<double> tab[]);
void fun_inverse_fourier_transform_DFT_method5_full_complex(int N,std::complex<double> tab[]);
void DFT_FFT_1_stage(int N,std::complex<double> tab[]);
void inverse_DFT_FFT_1_stage(int N,std::complex<double> tab[]);
int a2,a3,a4,a5,a6,a7;
int fi=0;
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=10;
std::complex<double> tab2[4096]={};
std::complex<double> tab3[4096]={};
std::complex<double> tab4[4096]={};
for(int i=0;i<N;i++)
{
tab2[i].real()=i*1.2+1;
tab2[i].imag()=-N-i-1;
tab3[i].real()=i*1.2+1;
tab3[i].imag()=-N-i-1;
tab4[i].real()=i*1.2+1;
tab4[i].imag()=-N-i-1;
}
double time2;
double zmienna=0;
cout<<"signal 1="<<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 2="<<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;
system("pause");
clock_t start = clock();
for(int i=0;i<N;i++)
{
tab2[i].real()=i*1.2+1;
tab2[i].imag()=-N-i-1;
tab3[i].real()=i*1.2+1;
tab3[i].imag()=-N-i-1;
tab4[i].real()=i*1.2+1;
tab4[i].imag()=-N-i-1;
}
for(int i2=0;i2<9;i2++)
{
if(i2==0)
{
a2=0;
}
else if(i2==1)
{
a2=1;
}
else if(i2==2)
{
a2=2;
}
else if(i2==3)
{
a2=3;
}
else if(i2==4)
{
a2=4;
}
else if(i2==5)
{
a2=5;
}
else if(i2==6)
{
a2=6;
}
else if(i2==7)
{
a2=7;
}
else if(i2==8)
{
a2=8;
}
for(int i3=0;i3<2;i3++)
{
if(i3==0)
{
a3=0;
}
else{a3=1;}
for(int i4=0;i4<2;i4++)
{
if(i4==0)
{
a4=0;
}
else{a4=1;}
for(int i5=0;i5<2;i5++)
{
if(i5==0)
{
a5=0;
}
else{a5=1;}
for(int i6=0;i6<2;i6++)
{
if(i6==0)
{
a6=0;
}
else{a6=1;}
for(int i7=0;i7<2;i7++)
{
if(i7==0)
{
a7=0;
}
else{a7=1;}
for(int i=0;i<N;i++)
{
tab2[i].real()=i*1.2+1;
tab2[i].imag()=-N-i-1;
tab3[i].real()=i*1.2+1;
tab3[i].imag()=-N-i-1;
tab4[i].real()=i*1.2+1;
tab4[i].imag()=-N-i-1;
}
fun_fourier_transform_DFT_method5_full_complex(N,tab3);
DFT_FFT_1_stage(N,tab2);
time2=diffclock( start, clock() );
/*
cout<<"frequency Hz radix-4"<<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 DFT"<<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;
*/
// system("pause");
// cout<<"if radix-4 == DFT tab2[j].real(): "<<endl;system("pause");
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<<"j= "<<j<<" "<<round(tab2[j].real()*1000)/1000<<" ";
//system("pause");
}
else {
cout<<-99<<" .. ";
}
}
//system("pause");
// cout<<endl;
/////////////////////////////////////////////////////////////////
fun_inverse_fourier_transform_DFT_method5_full_complex(N,tab3);
inverse_DFT_FFT_1_stage(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;
cout<<endl;
for(int j=0;j<N;j++)
{
cout.precision(4);
cout<<round(tab2[j].imag()*1000)/1000<<" ";
}
cout<<endl;
cout<<endl;
*/
//fun_inverse_fourier_transform_DFT_method5_full_complex(N,tab3);
for(int j=0;j<N;j++)
{
if(((tab4[j].imag()-tab3[j].imag()>=-0.03)&&(tab4[j].imag()-tab3[j].imag()<=0.03))&&
((tab4[j].real()-tab3[j].real()>=-0.03)&&(tab4[j].real()-tab3[j].real()<=0.03))
)
{
cout.precision(4);
// cout<<"j= "<<j<<" "<<a2<<" "<<a3<<" "<<a4<<" "<<a5<<" "<<a6<<" "<<a7<<" "<<tab4[j].real()<<" = "<<round(tab3[j].real()*1000)/1000<<" , "<<tab4[j].imag()<<" = "<<round(tab3[j].imag()*1000)/1000<<" ";
//system("pause");
}
else {
//cout<<-99<<" .. ";
}
}
//inverse_DFT_FFT_1_stage(N,tab2);
if(((tab4[1].imag()-tab2[1].imag()>=-0.03)&&(tab4[1].imag()-tab2[1].imag()<=0.03))&&
((tab4[3].imag()-tab2[3].imag()>=-0.03)&&(tab4[3].imag()-tab2[3].imag()<=0.03))
)
{
for(int j=0;j<N;j++)
{
if(((tab4[j].imag()-tab2[j].imag()>=-0.03)&&(tab4[j].imag()-tab2[j].imag()<=0.03))&&
((tab4[j].real()-tab2[j].real()>=-0.03)&&(tab4[j].real()-tab2[j].real()<=0.03))
)
{
cout<<"OK ";
cout.precision(4);
//cout<<"j= "<<j<<" "<<a2<<" "<<a3<<" "<<a4<<" "<<a5<<" "<<a6<<" "<<a7<<" "<<tab4[j].real()<<" = "<<round(tab2[j].real()*1000)/1000<<" , "<<tab4[j].imag()<<" = "<<round(tab2[j].imag()*1000)/1000<<" ";
//system("pause");
}
else {
//cout<<-99<<" .. ";
}
if(((tab4[j].imag()-tab2[1].imag()>=-0.03)&&(tab4[j].imag()-tab2[1].imag()<=0.03)))
{
cout<<"j= "<<j<<" "<<a2<<" "<<a3<<" "<<a4<<" "<<a5<<" "<<a6<<" "<<a7<<" | "<<tab4[1].real()<<" = "<<round(tab2[1].real()*1000)/1000<<" , "<<tab4[1].imag()<<" = "<<round(tab2[1].imag()*1000)/1000<<" "<<endl;
cout<<tab4[3].real()<<" = "<<round(tab2[3].real()*1000)/1000<<" , "<<tab4[3].imag()<<" = "<<round(tab2[3].imag()*1000)/1000<<" "<<endl;
system("pause");
}
}
}
}}}}}}
cout<<endl;
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[4096]={}; // tab2[]==N
std::complex<double> w[1]={{1,0}};
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];
}
//cout<<tab2[1]<<endl;
}
for(int j=0;j<N;j++)
{
tab[j].real() =tab2[j].real()*2/N;
tab[j].imag() =tab2[j].imag()*2/N;
}
/*
for(int j=0;j<N;j++)
{
tab[j].real() =tab2[j].real();
tab[j].imag() =tab2[j].imag();
}
*/
}
///////////////////////////////////////////////////
void fun_inverse_fourier_transform_DFT_method5_full_complex(int N,std::complex<double> tab[])
{
const double pi=3.141592653589793238462;
std::complex<double> tab2[4096]={}; // tab2[]==N
std::complex<double> w[1]={{1,0}};
std::complex<double> w2[1]={{1,0}};
std::complex<double> w3[1]={{1,0}};
int nr1=0,nr2=0,nr3=0,nr4=0;
std::complex<double> nr5=0;
double zmienna1=2*pi/(float)N;
double fi2=fi;
//////////////////////////
if(a2==0)
{
w3[0].real()=-1;
w3[0].imag()=-1;
}
else if(a2==1)
{
w3[0].real()=-1;
w3[0].imag()=0;
}
else if(a2==2)
{
w3[0].real()=-1;
w3[0].imag()=1;
}
else if(a2==3)
{
w3[0].real()=0;
w3[0].imag()=-1;
}
else if(a2==4)
{
w3[0].real()=0;
w3[0].imag()=0;
}
else if(a2==5)
{
w3[0].real()=0;
w3[0].imag()=1;
}
else if(a2==6)
{
w3[0].real()=1;
w3[0].imag()=-1;
}
else if(a2==7)
{
w3[0].real()=1;
w3[0].imag()=0;
}
else if(a2==8)
{
w3[0].real()=1;
w3[0].imag()=1;
}
////////////////
if(a3==0){nr1=-1;}
else{nr1=1;}
if(a4==0){nr2=-1;}
else{nr2=1;}
if(a5==0){nr3=-1;}
else{nr3=1;}
if(a6==0){nr4=-1;}
else{nr4=1;}
if(a7==0){nr5.real()=-1,nr5.imag()=0;}
else{nr5.real()=1,nr5.imag()=0;}
for (int i=0;i<N;i++)
{
for(int j=0;j<N;j++)
{
//complex number method:
w[0].real()=nr1*cos(nr3*i*j*zmienna1+fi2);
w[0].imag()=nr2*sin(nr4*i*j*zmienna1+fi2);
tab2[i]=tab2[i]+nr5*tab[j]*w[0]*w3[0];
}
//cout<<tab2[1]<<endl;
}
for(int j=0;j<N;j++)
{
tab[j].real() =tab2[j].real()*0.5;
tab[j].imag() =tab2[j].imag()*0.5;
}
/*
for(int j=0;j<N;j++)
{
tab[j].real() =tab2[j].real();
tab[j].imag() =tab2[j].imag();
}
*/
}
void inverse_DFT_FFT_1_stage(int N,std::complex<double> tab[])
{
const double pi=3.141592653589793238462;
std::complex<double> tab2[4096]={}; // tab2[]==N
std::complex<double> w[8]={{1,0}};
std::complex<double> w2[1]={{1,0}};
std::complex<double> w3[1]={{1,0}};
std::complex<double> w4[1]={{1,0}};
int nr1=0,nr2=0,nr3=0,nr4=0,nr6=0,nr7=0,nr8=0,nr9=0;
std::complex<double> nr5=0;
std::complex<double> nr10=0;
double zmienna1=2*pi/(float)N;
double fi2=fi;
//////////////////////////
if(a2==0)
{
w3[0].real()=-1;
w3[0].imag()=-1;
}
else if(a2==1)
{
w3[0].real()=-1;
w3[0].imag()=0;
}
else if(a2==2)
{
w3[0].real()=-1;
w3[0].imag()=1;
}
else if(a2==3)
{
w3[0].real()=0;
w3[0].imag()=-1;
}
else if(a2==4)
{
w3[0].real()=0;
w3[0].imag()=0;
}
else if(a2==5)
{
w3[0].real()=0;
w3[0].imag()=1;
}
else if(a2==6)
{
w3[0].real()=1;
w3[0].imag()=-1;
}
else if(a2==7)
{
w3[0].real()=1;
w3[0].imag()=0;
}
else if(a2==8)
{
w3[0].real()=1;
w3[0].imag()=1;
}
////////////////
if(a3==0){nr1=-1;}
else{nr1=1;}
if(a4==0){nr2=-1;}
else{nr2=1;}
if(a5==0){nr3=-1;}
else{nr3=1;}
if(a6==0){nr4=-1;}
else{nr4=1;}
if(a7==0){nr5.real()=-1,nr5.imag()=0;}
else{nr5.real()=1,nr5.imag()=0;}
/////////////////////////
for(int k=0;k<N;k=k+2)
{
for(int n=0;n<(N/2);n++)
{
//complex number method:
w[0].real()=nr1*cos(nr3*k*n*zmienna1);
w[0].imag()=nr2*sin(nr4*k*n*zmienna1);
w[1].real()=nr1*cos(nr3*(k+1)*n*zmienna1);
w[1].imag()=nr2*sin(nr4*(k+1)*n*zmienna1);
//parzyste:
tab2[k]=tab2[k] +nr5*w[0]*(tab[n]+tab[N/2+n])*w3[0];
//nieparzyste
tab2[k+1]=tab2[k+1]+nr5*w[1]*(tab[n]-tab[N/2+n])*w3[0];
}
}
for(int j=0;j<N;j++)
{
tab[j].real() =tab2[j].real()*0.5;
tab[j].imag() =tab2[j].imag()*0.5;
}
/*
for(int i=0;i<N;i++)
{
}
*/
}
void DFT_FFT_1_stage(int N,std::complex<double> tab[])
{
const double pi=3.141592653589793238462;
std::complex<double> tab2[4096]={}; // tab2[]==N
std::complex<double> w[8]={{1,0},{1,0},{1,0},{1,0}};
double zmienna1=2*pi/(float)N;
double fi2=fi;
for(int k=0;k<N;k=k+2)
{
for(int n=0;n<(N/2);n++)
{
//complex number method:
w[0].real()= cos(k*n*zmienna1);
w[0].imag()=-sin(k*n*zmienna1);
w[1].real()= cos((k+1)*n*zmienna1);
w[1].imag()=-sin((k+1)*n*zmienna1);
//parzyste:
tab2[k]=tab2[k]+w[0]*(tab[n]+tab[N/2+n]);
//nieparzyste
tab2[k+1]=tab2[k+1]+w[1]*(tab[n]-tab[N/2+n]);
}
}
for(int j=0;j<N;j++)
{
tab[j].real() =tab2[j].real()*2/N;
tab[j].imag() =tab2[j].imag()*2/N;
}
/*
for(int i=0;i<N;i++)
{
tab[i]=tab2[i];
}
*/
}
Brak komentarzy:
Prześlij komentarz