//Implementation of Split Radix Algorithm for 4-Point inverse srFFT source code c++
//results are equal as in discrete fourier DFT
#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(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
N=4;
std::complex<double> tab2[4]={{1,5}, {2,11} ,{3,18}, {4,6}};
//std::complex<double> tab2[4]={{0.707106781}, {0.292893219} ,{-0.707106781}, {-0.292893219}};
//std::complex<double> tab2[4]={{0,1}, {0,2} ,{0,3}, {0,4}};
//std::complex<double> tab2[4]={{0,0}, {0,0} ,{0,0}, {0,3}};
//std::complex<double> tab2[4]={{7,0}, {0,0} ,{0,0}, {0,0}};
//std::complex<double> tab2[4]={{1}, {2} ,{3}, {4}};
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;
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;
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;
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;
// plik.close();
system("pause");
cout<<endl;
//fun_inverse_bits(N,tab2);
fun_inverse_fourier_transform_srFFT_method_N_6(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;
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}};
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;
*/
//srFFT 4-point
w[0]={{1,1}};
//w[0].real()=cos(0*2*pi/4);
//w[0].imag()=-sin(0*2*pi/4);
tab2[0]=(tab[0]+tab[2]+tab[1]+tab[3]);
//tab2[0]=(tab[0]+tab[2]+tab[1]+tab[3])*w[0];
w[0].real()=cos(1*2*pi/4);
w[0].imag()=-sin(1*2*pi/4);
w2[0].real()=0;
w2[0].imag()=1;
//w3[0].real()=-1;
//w3[0].imag()=0;
tab2[1]=(w2[0]*(tab[0]-tab[2])+(tab[1]-tab[3]))*w[0];
//tab2[1]=(w2[0]*(tab[0]+w3[0]*tab[2])+(tab[1]+w3[0]*tab[3]))*w[0];
//w[0].real()=cos(2*2*pi/4);
//w[0].imag()=-sin(2*2*pi/4);
tab2[2]=(tab[0]+tab[2]-tab[1]-tab[3]);
//tab2[2]=(-tab[0]-tab[2]+(tab[1]+tab[3]))*w[0];
w[0].real()=cos(3*2*pi/4);
w[0].imag()=-sin(3*2*pi/4);
w2[0].real()=0;
w2[0].imag()=-1;
//w3[0].real()=-1;
//w3[0].imag()=0;
tab2[3]=(w2[0]*(tab[0]-tab[2])+(tab[1]-tab[3]))*w[0];
//tab2[3]=(w2[0]*(tab[0]+w3[0]*tab[2])+(tab[1]+w3[0]*tab[3]))*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_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}};
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;
*/
//srFFT 4-point
w[0]={{1,1}};
//w[0].real()=cos(0*2*pi/4);
//w[0].imag()=sin(0*2*pi/4);
tab2[0]=(tab[0]+tab[2]+tab[1]+tab[3]);
//tab2[0]=(tab[0]+tab[2]+tab[1]+tab[3])*w[0];
w[0].real()=cos(1*2*pi/4);
w[0].imag()=sin(1*2*pi/4);
w2[0].real()=0;
w2[0].imag()=-1;
//w3[0].real()=-1;
//w3[0].imag()=0;
tab2[1]=(w2[0]*(tab[0]-tab[2])+(tab[1]-tab[3]))*w[0];
//tab2[1]=(w2[0]*(tab[0]+w3[0]*tab[2])+(tab[1]+w3[0]*tab[3]))*w[0];
//w[0].real()=cos(2*2*pi/4);
//w[0].imag()=sin(2*2*pi/4);
tab2[2]=(tab[0]+tab[2]-tab[1]-tab[3]);
//tab2[2]=(-tab[0]-tab[2]+(tab[1]+tab[3]))*w[0];
w[0].real()=cos(3*2*pi/4);
w[0].imag()=sin(3*2*pi/4);
w2[0].real()=0;
w2[0].imag()=1;
//w3[0].real()=-1;
//w3[0].imag()=0;
tab2[3]=(w2[0]*(tab[0]-tab[2])+(tab[1]-tab[3]))*w[0];
//tab2[3]=(w2[0]*(tab[0]+w3[0]*tab[2])+(tab[1]+w3[0]*tab[3]))*w[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[])
{
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;
}
}
it is something like radix-4 algorlithm because orginal algorlithm from that pictures dont work for me for all signals
algorithm for DFT: http://inverse-fourier-transformdftalgorithm.blogspot.com/
Brak komentarzy:
Prześlij komentarz