//Implementation of Split Radix Algorithm for 6-Point inverse srFFT source code c++
#include <iostream>
#include "conio.h"
#include <stdlib.h>
#include <math.h>
#include <cmath>
#include <time.h>
#include <complex>
using namespace std;
//complex number method:
//void fun_inverse_bits(int N,std::complex<double> tab[]);
void fun_fourier_transform_srFFT_method_N_6(int N,std::complex<double> tab[]);
void fun_inverse_fourier_transform_srFFT_method_N_6(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
// period= 16 samples=16
/*
N=16;
std::complex<double> tab2[16]={{-0.923879533},{0.382683432},{1.03153013},{0.923879533},{0.923879533},{1.465075634},{1.796896994},{0.923879533},{-0.923879533},
{-2.230442498},{-1.796896994},{-0.158512669},{0.923879533},{0.382683432},{-1.03153013},{-1.689246397}};
*/
/*
N=12;
std::complex<double> tab2[12]={{0}, {2.366025404}, {1.732050808}, {0}, {0},
{0.633974596}, {0},{-0.633974596}, {0}, {0}, {-1.732050808}, {-2.366025404}
};
*/
N=6;
// std::complex<double> tab2[6]={{0}, {2.598076212},{0.866025404}, {0},
// {-0.866025404}, {-2.598076212}};
std::complex<double> tab2[6]={{0.707106781}, {1.473231763}, {-0.965925826},
{-0.707106781}, {0.258819045}, {-0.766124982}};
double time2;
cout<<"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;
clock_t start = clock();
//fun_inverse_bits(N,tab2);
fun_fourier_transform_srFFT_method_N_6(N,tab2);
time2=diffclock( start, clock() );
cout<<"frequency Hz"<<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;
system("pause");
//fun_inverse_bits(N,tab2);
fun_inverse_fourier_transform_srFFT_method_N_6(N,tab2);
cout<<"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;
system("pause");
return 0;
}
void fun_fourier_transform_srFFT_method_N_6(int N,std::complex<double> tab[])
{
const double pi=3.141592653589793238462;
std::complex<double> tab2[6]={}; // tab2[]==N
std::complex<double> w[1]={{1,1}};
/*
for(int j=0;j<N;j++)
{
cout<< tab[j]<<" ";
}
cout<<endl<<"........."<<endl;
*/
w[0]={{1,1}};
w[0].real()=cos(0*2*pi/(float)N);
w[0].imag()=-sin(0);
tab2[0]=tab[0]+tab[3]*w[0];
tab2[3]=tab[0]-tab[3]*w[0];
tab2[1]=tab[1]+tab[4];
tab2[4]=-tab[1]+tab[4];
tab2[2]=tab[2]+tab[5];
tab2[5]=tab[2]-tab[5];
///////////////////////
//3 point srFTT
tab[1]=tab2[1]+tab2[2];
tab[2]=tab2[1]-tab2[2];
tab[0]=tab2[0]+tab[1];
w[0].real()=cos(1*2*pi/3);
w[0].imag()=0;
tab2[1]=tab2[0]+tab[1]*w[0];
w[0].real()=0;
w[0].imag()=sin(1*2*pi/3);
tab[1]=tab2[1]+tab[2]*w[0];
w[0].real()=-0;
w[0].imag()=sin(1*2*pi/3);
tab[2]=tab2[1]-tab[2]*w[0];
////////////////////////
//3 point srFTT
tab[4]=tab2[4]+tab2[5];
tab[5]=tab2[4]-tab2[5];
tab[3]=tab2[3]+tab[4];
w[0].real()=cos(1*2*pi/3);
w[0].imag()=0;
tab2[4]=tab2[3]+tab[4]*w[0];
w[0].real()=0;
w[0].imag()=sin(1*2*pi/3);
tab[4]=tab2[4]+tab[5]*w[0];
w[0].real()=-0;
w[0].imag()=sin(1*2*pi/3);
tab[5]=tab2[4]-tab[5]*w[0];
/*
in official algorithm is:
tab2[0]=tab[0];
tab2[1]=tab[4];
tab2[2]=tab[1];
tab2[3]=tab[5];
tab2[4]=tab[2];
tab2[5]=tab[3];
*/
//should be like as in DFT n=6 that same result provides these combinations:
tab2[0]=tab[0];
tab2[1]=tab[4];
tab2[2]=tab[2];
tab2[3]=tab[3];
tab2[4]=tab[1];
tab2[5]=tab[5];
///////////////
tab[0]=tab2[0];
tab[1]=tab2[1];
tab[2]=tab2[2];
tab[3]=tab2[3];
tab[4]=tab2[4];
tab[5]=tab2[5];
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_srFFT_method_N_6(int N,std::complex<double> tab[])
{
const double pi=3.141592653589793238462;
std::complex<double> tab2[6]={}; // tab2[]==N
std::complex<double> w[1]={{1,1}};
/*
for(int j=0;j<N;j++)
{
cout<< tab[j]<<" ";
}
cout<<endl<<"........."<<endl;
*/
w[0]={{1,1}};
w[0].real()=cos(0*2*pi/(float)N);
w[0].imag()=sin(0);
tab2[0]=tab[0]+tab[3]*w[0];
tab2[3]=tab[0]-tab[3]*w[0];
tab2[1]=tab[1]+tab[4];
tab2[4]=-tab[1]+tab[4];
tab2[2]=tab[2]+tab[5];
tab2[5]=tab[2]-tab[5];
///////////////////////
//3 point srFTT
tab[1]=tab2[1]+tab2[2];
tab[2]=tab2[1]-tab2[2];
tab[0]=tab2[0]+tab[1];
w[0].real()=cos(1*2*pi/3);
w[0].imag()=0;
tab2[1]=tab2[0]+tab[1]*w[0];
w[0].real()=0;
w[0].imag()=-sin(1*2*pi/3);
tab[1]=tab2[1]+tab[2]*w[0];
w[0].real()=-0;
w[0].imag()=-sin(1*2*pi/3);
tab[2]=tab2[1]-tab[2]*w[0];
////////////////////////
//3 point srFTT
tab[4]=tab2[4]+tab2[5];
tab[5]=tab2[4]-tab2[5];
tab[3]=tab2[3]+tab[4];
w[0].real()=cos(1*2*pi/3);
w[0].imag()=0;
tab2[4]=tab2[3]+tab[4]*w[0];
w[0].real()=0;
w[0].imag()=-sin(1*2*pi/3);
tab[4]=tab2[4]+tab[5]*w[0];
w[0].real()=-0;
w[0].imag()=-sin(1*2*pi/3);
tab[5]=tab2[4]-tab[5]*w[0];
/*
in official algorithm is:
tab2[0]=tab[0];
tab2[1]=tab[4];
tab2[2]=tab[1];
tab2[3]=tab[5];
tab2[4]=tab[2];
tab2[5]=tab[3];
*/
//should be like as in DFT n=6 that same result provides these combinations:
tab2[0]=tab[0];
tab2[1]=tab[4];
tab2[2]=tab[2];
tab2[3]=tab[3];
tab2[4]=tab[1];
tab2[5]=tab[5];
///////////////
tab[0]=tab2[0];
tab[1]=tab2[1];
tab[2]=tab2[2];
tab[3]=tab2[3];
tab[4]=tab2[4];
tab[5]=tab2[5];
for(int j=0;j<N;j++)
{
tab[j].real() =tab[j].real()*0.5;
tab[j].imag() =tab[j].imag()*0.5;
}
}
Brak komentarzy:
Prześlij komentarz