void fun_fourier_transform_FFT_N_6(int N,std::complex<double> tab[])
{
//N=6
//author marcin matysek (r)ewertyn.PL
const double pi=3.141592653589793238462;
std::complex<double> tab2[6]={}; // tab2[]==N
std::complex<double> w1[1]={{1,0}};
std::complex<double> w2[1]={{1,0}};
std::complex<double> w3[1]={{1,0}};
std::complex<double> w4[1]={{1,0}};
std::complex<double> w5[1]={{1,0}};
std::complex<double> w6[1]={{1,0}};
std::complex<double> tmp=0;
double tmp5;
int k=0;
int n=0;
tmp5=2*pi/(N/1);
//stage 1
k=0,n=0;
w2[0].real()= cos(k*n*tmp5);
w2[0].imag()=-sin(k*n*tmp5);
w1[0].real()= cos(k*n*tmp5);
w1[0].imag()=-sin(k*n*tmp5);
tab2[0]=w1[0]*tab[0]+w2[0]*tab[3];
tab2[3]=w1[0]*tab[0]-w2[0]*tab[3];
w2[0].real()= cos(k*n*tmp5);
w2[0].imag()=-sin(k*n*tmp5);
w1[0].real()= cos(k*n*tmp5);
w1[0].imag()=-sin(k*n*tmp5);
tab2[1]=w1[0]*tab[1]+w2[0]*tab[4];
tab2[4]=w1[0]*tab[1]-w2[0]*tab[4];
w2[0].real()= cos(k*n*tmp5);
w2[0].imag()=-sin(k*n*tmp5);
w1[0].real()= cos(k*n*tmp5);
w1[0].imag()=-sin(k*n*tmp5);
tab2[2]=w1[0]*tab[2]+w2[0]*tab[5];
tab2[5]=w1[0]*tab[2]-w2[0]*tab[5];
////////////////////////////////////////////////
//stage 2 nb1
w1[0].real()= cos(0*tmp5);//0
w1[0].imag()=-sin(0*tmp5);//0
w2[0].real()= cos(0*tmp5);//0
w2[0].imag()=-sin(0*tmp5);//0
w5[0].real()= cos(0*tmp5);//0
w5[0].imag()=-sin(0*tmp5);//0
tab[0]=w1[0]*tab2[0]+w2[0]*tab2[1]+w5[0]*tab2[2];
w3[0].real()= cos(0*tmp5);//0
w3[0].imag()=-sin(0*tmp5);//0
w4[0].real()= cos(1*tmp5);//1
w4[0].imag()=-sin(1*tmp5);//1
w6[0].real()= cos(2*tmp5);//2
w6[0].imag()=-sin(2*tmp5);//2
tab[1]=w3[0]*tab2[0]-w4[0]*tab2[1]+w6[0]*tab2[2];
//or
//w3[0].real()= cos(0*tmp5);//0
//w3 [0].imag()=-sin(0*tmp5);//0
//w4[0].real()= cos(4*tmp5);//
//w4[0].imag()=-sin(4*tmp5);//
//w6[0].real()= cos(2*tmp5);//
//w6[0].imag()=-sin(2*tmp5);//
//tab[1]=w3[0]*tab2[0]+w4[0]*tab2[1]+w6[0]*tab2[2];
w1[0].real()= cos(0*tmp5);//0
w1[0].imag()=-sin(0*tmp5);//0
w2[0].real()= cos(2*tmp5);//2
w2[0].imag()=-sin(2*tmp5);//2
w5[0].real()= cos(4*tmp5);//4
w5[0].imag()=-sin(4*tmp5);//4
tab[2]=w1[0]*tab2[0]+w2[0]*tab2[1]+w5[0]*tab2[2];
//or
//w1[0].real()= cos(0*tmp5);//0
//w1[0].imag()=-sin(0*tmp5);//0
//w2[0].real()= cos(2*tmp5);//
//w2[0].imag()=-sin(2*tmp5);//
//w5[0].real()= cos(1*tmp5);//
//w5[0].imag()=-sin(1*tmp5);//
//tab[2]=w1[0]*tab2[0]+w2[0]*tab2[1]-w5[0]*tab2[2];
//stage 2 nb2
int n1=3;
w1[0].real()= cos(0*tmp5);//0
w1[0].imag()=-sin(0*tmp5);//0
w2[0].real()= cos(0*tmp5);//0
w2[0].imag()=-sin(0*tmp5);//0
w3[0].real()= cos(0*tmp5);//0
w3[0].imag()=-sin(0*tmp5);//0
tab[0+n1]=w1[0]*tab2[0+n1]+w2[0]*tab2[1+n1]-w3[0]*tab2[2+n1];
w3[0].real()= cos(0*tmp5);//0
w3[0].imag()=-sin(0*tmp5);//0
w4[0].real()= cos(1*tmp5);//1
w4[0].imag()=-sin(1*tmp5);//1
w6[0].real()= cos(2*tmp5);//2
w6[0].imag()=-sin(2*tmp5);//2
tab[1+n1]=w3[0]*tab2[0+n1]+w4[0]*tab2[1+n1]+w6[0]*tab2[2+n1];
w1[0].real()= cos(0*tmp5);//0
w1[0].imag()=-sin(0*tmp5);//0
w2[0].real()= cos(2*tmp5);//2
w2[0].imag()=-sin(2*tmp5);//2
w5[0].real()= cos(4*tmp5);//4
w5[0].imag()=-sin(4*tmp5);//4
tab[2+n1]=w1[0]*tab2[0+n1]-w2[0]*tab2[1+n1]+w5[0]*tab2[2+n1];
//inverse
tmp=tab[1];
tab[1]=tab[4];
tab[4]=tmp;
for(int j=0;j<N;j++)
{
tab[j].real() =tab[j].real()*2/N;
tab[j].imag() =tab[j].imag()*2/N;
}
}
void void fun_fourier_transform_FFT_N_6v2(int N,std::complex<double> tab[])
{
//N=6
//author marcin matysek (r)ewertyn.PL
const double pi=3.141592653589793238462;
std::complex<double> tab2[6]={}; // tab2[]==N
std::complex<double> w1[1]={{1,0}};
std::complex<double> w2[1]={{1,0}};
std::complex<double> w3[1]={{1,0}};
std::complex<double> w4[1]={{1,0}};
std::complex<double> w5[1]={{1,0}};
std::complex<double> w6[1]={{1,0}};
std::complex<double> tmp=0;
double tmp5;
int k=0;
int n=0;
tmp5=2*pi/(N/1);
//stage 1
k=0,n=0;
w2[0].real()= cos(k*n*tmp5);
w2[0].imag()=-sin(k*n*tmp5);
w1[0].real()= cos(k*n*tmp5);
w1[0].imag()=-sin(k*n*tmp5);
tab2[0]=w1[0]*tab[0]+w2[0]*tab[3];
tab2[3]=w1[0]*tab[0]-w2[0]*tab[3];
w2[0].real()= cos(k*n*tmp5);
w2[0].imag()=-sin(k*n*tmp5);
w1[0].real()= cos(k*n*tmp5);
w1[0].imag()=-sin(k*n*tmp5);
tab2[1]=w1[0]*tab[1]+w2[0]*tab[4];
tab2[4]=w1[0]*tab[1]-w2[0]*tab[4];
w2[0].real()= cos(k*n*tmp5);
w2[0].imag()=-sin(k*n*tmp5);
w1[0].real()= cos(k*n*tmp5);
w1[0].imag()=-sin(k*n*tmp5);
tab2[2]=w1[0]*tab[2]+w2[0]*tab[5];
tab2[5]=w1[0]*tab[2]-w2[0]*tab[5];
//stage 2 nb1
w1[0].real()= cos(0*tmp5);//0
w1[0].imag()=-sin(0*tmp5);//0
w2[0].real()= cos(0*tmp5);//0
w2[0].imag()=-sin(0*tmp5);//0
w5[0].real()= cos(0*tmp5);//0
w5[0].imag()=-sin(0*tmp5);//0
tab[0]=w1[0]*tab2[0]+w2[0]*tab2[1]+w5[0]*tab2[2];
w3[0].real()= cos(0*tmp5);//0
w3[0].imag()=-sin(0*tmp5);//0
w4[0].real()= cos(1*tmp5);//1
w4[0].imag()=-sin(1*tmp5);//1
w6[0].real()= cos(2*tmp5);//2
w6[0].imag()=-sin(2*tmp5);//2
tab[1]=w3[0]*tab2[0]-w4[0]*tab2[1]+w6[0]*tab2[2];
w1[0].real()= cos(0*tmp5);//0
w1[0].imag()=-sin(0*tmp5);//0
w2[0].real()= cos(2*tmp5);//2
w2[0].imag()=-sin(2*tmp5);//2
w5[0].real()= cos(4*tmp5);//4
w5[0].imag()=-sin(4*tmp5);//4
tab[2]=w1[0]*tab2[0]+w2[0]*tab2[1]+w5[0]*tab2[2];
int n1=3;
w1[0].real()= cos(0*tmp5);//0
w1[0].imag()=-sin(0*tmp5);//0
w2[0].real()= cos(0*tmp5);//0
w2[0].imag()=-sin(0*tmp5);//0
w3[0].real()= cos(0*tmp5);//0
w3[0].imag()=-sin(0*tmp5);//0
tab[0+n1]=w1[0]*tab2[0+n1]+w2[0]*tab2[1+n1]-w3[0]*tab2[2+n1];
w3[0].real()= cos(0*tmp5);//0
w3[0].imag()=-sin(0*tmp5);//0
w4[0].real()= cos(4*tmp5);//
w4[0].imag()=-sin(4*tmp5);//
w6[0].real()= cos(2*tmp5);//
w6[0].imag()=-sin(2*tmp5);//
tab[1+n1]=w3[0]*tab2[0+n1]-w4[0]*tab2[1+n1]+w6[0]*tab2[2+n1];
//or
// w3[0].real()= cos(0*tmp5);//0
// w3[0].imag()=-sin(0*tmp5);//0
// w4[0].real()= cos(2*tmp5);//
// w4[0].imag()=-sin(2*tmp5);//
// w6[0].real()= cos(4*tmp5);//
// w6[0].imag()=-sin(4*tmp5);//
// tab[1+n1]=w3[0]*tab2[0+n1]+w4[0]*tab2[1+n1]-w6[0]*tab2[2+n1];
w1[0].real()= cos(0*tmp5);//0
w1[0].imag()=-sin(0*tmp5);//0
w2[0].real()= cos(4*tmp5);//
w2[0].imag()=-sin(4*tmp5);//
w5[0].real()= cos(2*tmp5);//
w5[0].imag()=-sin(2*tmp5);//
tab[2+n1]=w1[0]*tab2[0+n1]+w2[0]*tab2[1+n1]-w5[0]*tab2[2+n1];
//inverse
tmp=tab[1];
tab[1]=tab[4];
tab[4]=tmp;
for(int j=0;j<N;j++)
{
tab[j].real() =tab[j].real()*2/N;
tab[j].imag() =tab[j].imag()*2/N;
}
}
niedziela, 30 kwietnia 2017
sobota, 15 kwietnia 2017
współdzielenie pamięci przez 2 zmienne zamiast wskaźników, 2 zmienne w 1 pamięci
//współdzielenie pamięci przez 2 zmienne zamiast wskaźników, 2 zmienne w 1 pamięci
//author marcin matysek (r)ewertyn.PL
//source:
//http://stackoverflow.com/questions/4629317/what-does-int-mean
// http://www.cplusplus.com/forum/windows/17153/
#include <iostream>
#include <conio.h>
#include <stdlib.h>
using namespace std;
void funkcja(int & nb3);
int main()
{
int nb1=0;
int &nb2=nb1;
cout<<"wartosc zmiennej1= "<<nb1<<" jej adres to: "<<&nb1<<endl;
cout<<"wartosc zmiennej2= "<<nb2<<" jej adres to: "<<&nb2<<endl;
nb2=10;
cout<<endl;
cout<<"po modyfikacji zmiennej 2 zmienne maja wartosci:"<<endl;
cout<<"wartosc zmiennej1= "<<nb1<<" jej adres to: "<<&nb1<<endl;
cout<<"wartosc zmiennej2= "<<nb2<<" jej adres to: "<<&nb2<<endl;
funkcja(nb1);
cout<<endl;
cout<<"po wyjsciu z funkcji:"<<endl;
cout<<"wartosc zmiennej1= "<<nb1<<" jej adres to: "<<&nb1<<endl;
cout<<"wartosc zmiennej2= "<<nb2<<" jej adres to: "<<&nb2<<endl;
system("pause");
return 0;
}
void funkcja(int & nb3)
{
cout<<endl;
cout<<"w funcki zmienna 3 przyjmuje wartosc zmiennej 1:"<<endl;
cout<<"wartosc zmiennej3= "<<nb3<<" jej adres to: "<<&nb3<<endl;
cout<<"po modyfikacji zmiennej 3 ,zmienna 3 w funkcji ma wartosc:"<<endl;
nb3=20;
cout<<endl;
cout<<"wartosc zmiennej3= "<<nb3<<" jej adres to: "<<&nb3<<endl;
}
poniedziałek, 10 kwietnia 2017
moja interpretacja algorytmu FFT 16*radix-4 dla N=16 c++ source code
////https://www.beechwood.eu/fft-implementation-r2-dit-r4-dif-r8-dif/
//copyrights marcin matysek (r)ewertyn.PL
//moja interpretacja algorytmu FFT 16*radix-4 dla N=16 c++ source code działa z
// FFT=+/-cos(w+fi)-isin(w+fi) and maybe witch other functions that are mirror inverse
#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[])
{
//code by Sidney Burrus
//http://dsp.stackexchange.com/questions/3481/radix-4-fft-implementation
// Radix-4 bit-reverse
}
///////////////////
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;
}
}
moja interpretacja algorytmu FFT radix-4 c++ source code
//https://www.beechwood.eu/fft-implementation-r2-dit-r4-dif-r8-dif/
//moja interpretacja algorytmu FFT radix-4 c++ source code działa z
// 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[])
{
//code by Sidney Burrus
//http://dsp.stackexchange.com/questions/3481/radix-4-fft-implementation
}
///////////////////
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;
}
}
sobota, 8 kwietnia 2017
szybka transformacja fouriera iFFT radix-4 dla N=256 algorytm c++ kod źródłowy implementacja v2
//source:
//https://www.google.ch/patents/US6957241
//http://www.google.ch/patents/US20020083107
//https://www.beechwood.eu/fft-implementation-r2-dit-r4-dif-r8-dif/
//http://www.cmlab.csie.ntu.edu.tw/cml/dsp/training/coding/transform/fft.html
//author marcin matysek (r)ewertyn.PL
//szybka transformacja fouriera iFFT radix-4 dla N=256 algorytm c++ kod źródłowy implementacja v2
//this is new in that method:
//when you want to have equal results that are in false modificator in normal FFT then change this:
/*
fun_fourier_transform_FFT_radix_4_N_256_official
{
for(int j=0;j<N;j++)
{
tab[j].real() =tab[j].real()*2/N;
tab[j].imag() =tab[j].imag()*2/N;
}
}
//and:
fun_inverse_fourier_transform_FFT_radix_4_N_256_official
{
for(int j=0;j<N;j++)
{
tab[j].real() =tab[j].real()*0.5;
tab[j].imag() =tab[j].imag()*0.5;
}
}
//for official modificator that is only in inverse FFT:
fun_fourier_transform_FFT_radix_4_N_256_official
{
}
fun_inverse_fourier_transform_FFT_radix_4_N_256_official
{
for(int i=0;i<N;i++)
{
tablica1[0][i]=tablica1[0][i]*1/(float)N;
tablica1[1][i]=tablica1[1][i]*1/(float)N;
}
}
*/
//haven't try it with other function that cos(x)+jsin(x)=sin(x+pi/2)+jsin(x)
//fourier transform iFFT radix-4 for N=256 algorithm c++ source code implementation v2
#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_radix_4(int N,std::complex<double> tab[]);
void fun_fourier_transform_FFT_radix_4_N_256_official_v2(int N,std::complex<double> tab[]);
void fun_inverse_fourier_transform_FFT_radix_4_N_256_official_v2(int N,std::complex<double> tab[]);
void fun_fourier_transform_DFT_method5_full_complex(int N,std::complex<double> tab[]);
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=256;
std::complex<double> tab2[256]={};
std::complex<double> tab3[256]={};
for(int i=0;i<N;i++)
{
tab2[i].real()=i;
tab2[i].imag()=N+i;
tab3[i].real()=i;
tab3[i].imag()=N+i;
}
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;
clock_t start = clock();
fun_fourier_transform_DFT_method5_full_complex(N,tab3);
//////////////////////////////////////////////////////////
fun_fourier_transform_FFT_radix_4_N_256_official_v2(N,tab2);
fun_inverse_bits_radix_4(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;
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<<round(tab2[j].real()*1000)/1000<<" ";//system("pause");
}
else {
cout<<-1<<" . ";
}
}
system("pause");
cout<<endl;
/////////////////////////////////////////////////////////////////
fun_inverse_fourier_transform_FFT_radix_4_N_256_official_v2(N,tab2);
fun_inverse_bits_radix_4(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;
cout<<endl;
system("pause");
return 0;
}
void fun_inverse_bits_radix_4(int N,std::complex<double> tab[])
{
//code by Sidney Burrus
//http://dsp.stackexchange.com/questions/3481/radix-4-fft-implementation
}
///////////////
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,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_fourier_transform_FFT_radix_4_N_256_official_v2(int N,std::complex<double> tab[])
{
const double pi=3.141592653589793238462;
std::complex<double> tab2[256]={}; // tab2[]==N
std::complex<double> w1[1]={{1,1}};
std::complex<double> w2[1]={{1,1}};
std::complex<double> w3[1]={{1,1}};
std::complex<double> w4[1]={{1,1}};
std::complex<double> w5[1]={{1,1}};
std::complex<double> tmp1,tmp2,tmp3,tmp4;
double tmp5;
int flag=0;
//radix-4
//stage 1
tmp5=2*pi/(N/1);
for(int i=0;i<(N/4);i++)//64
{
w1[0].real()=cos(0*tmp5);
w1[0].imag()=sin(0*tmp5);
w2[0].real()=cos(0*tmp5);
w2[0].imag()=sin(0*tmp5);
w3[0].real()=cos(0*tmp5);
w3[0].imag()=sin(0*tmp5);
w4[0].real()=cos(0*tmp5);
w4[0].imag()=sin(0*tmp5);
tmp1=w1[0]*tab[i+0];
tmp2=w2[0]*tab[i+N/4];
tmp3=w3[0]*tab[i+N/2];
tmp4=w4[0]*tab[i+3*N/4];
tab2[i] =tmp1+tmp2+tmp3+tmp4;
w5[0].real()=0;
w5[0].imag()=1;
tab2[i+N/4] =tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab2[i+N/2] =tmp1-tmp2+tmp3-tmp4;
w5[0].real()=0;
w5[0].imag()=-1;
tab2[i+3*N/4] =tmp1-tmp3+w5[0]*(tmp2-tmp4);
}
///////////////////////////////////
//stage 2
for(int k=0;k<N;k++)
{
tab[k]={0,0};
}
tmp5=2*pi/(N/1);
for(int i=0;i<(N/64);i++)//4
{
for(int j=0;j<(N/16);j++)//16
{
w1[0].real()=cos((0*i+i*j)*tmp5);
w1[0].imag()=sin((0*i+i*j)*tmp5);
w2[0].real()=cos((1*16*i+i*j)*tmp5);
w2[0].imag()=sin((1*16*i+i*j)*tmp5);
w3[0].real()=cos((2*16*i+i*j)*tmp5);
w3[0].imag()=sin((2*16*i+i*j)*tmp5);
w4[0].real()=cos((3*16*i+i*j)*tmp5);
w4[0].imag()=sin((3*16*i+i*j)*tmp5);
tmp1=w1[0]*tab2[0+j+N/(4*1)*i];
tmp2=w2[0]*tab2[N/(16*1)+j+N/(4*1)*i];
tmp3=w3[0]*tab2[N/(8*1) +j+N/(4*1)*i];
tmp4=w4[0]*tab2[3*N/(16*1)+j+N/(4*1)*i];
tab[0+j+N/(4*1)*i] =tmp1+tmp2+ tmp3+tmp4;
w5[0].real()=0;
w5[0].imag()=1;
tab[N/(16*1)+j+N/(4*1)*i] =tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab[N/(8 *1)+j+N/(4*1)*i] =tmp1-tmp2+ tmp3-tmp4;
w5[0].real()=0;
w5[0].imag()=-1;
tab[3*N/(16*1)+j+N/(4*1)*i]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
}
}
//////////////////////////////
//stage 3
int tab4[256]={0};
for(int k=0;k<N;k++)
{
tab2[k]={0,0};
}
tmp5=2*pi/(N/1);
for(int i=0;i<(N/16);i++)//16
{
for(int j=0;j<4;j++)
{
w1[0].real()=cos((i-4*flag)*(16*0+4*j)*tmp5);
w1[0].imag()=sin((i-4*flag)*(16*0+4*j)*tmp5);
w2[0].real()=cos((i-4*flag)*(1*16+4*j)*tmp5);
w2[0].imag()=sin((i-4*flag)*(1*16+4*j)*tmp5);
w3[0].real()=cos((i-4*flag)*(2*16+4*j)*tmp5);
w3[0].imag()=sin((i-4*flag)*(2*16+4*j)*tmp5);
w4[0].real()=cos((i-4*flag)*(3*16+4*j)*tmp5);
w4[0].imag()=sin((i-4*flag)*(3*16+4*j)*tmp5);
tmp1=w1[0]*tab[0+j+N/(4*4)*i];
tmp2=w2[0]*tab[N/(16*4)+j+N/(4*4)*i];
tmp3=w3[0]*tab[N/(8*4) +j+N/(4*4)*i];
tmp4=w4[0]*tab[3*N/(16*4)+j+N/(4*4)*i];
tab2[0+j+N/(4*4)*i] =tmp1+tmp2+ tmp3+tmp4;
w5[0].real()=0;
w5[0].imag()=1;
tab2[N/(16*4)+j+N/(4*4)*i] =tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab2[N/(8*4)+j+N/(4*4)*i] =tmp1-tmp2+ tmp3-tmp4;
w5[0].real()=0;
w5[0].imag()=-1;
tab2[3*N/(16*4)+j+N/(4*4)*i]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
//tab4[0+j+N/(4*4)*i] =(i-4*flag)*(16*0+4*j);
//tab4[N/(16*4)+j+N/(4*4)*i] =(i-4*flag)*(1*16+4*j);
//tab4[N/(8*4)+j+N/(4*4)*i]=(i-4*flag)*(2*16+4*j);
//tab4[3*N/(16*4)+j+N/(4*4)*i]=(i-4*flag)*(3*16+4*j);
}
if((i+1)%4==0){flag++;}
}
/*
for(int i2=0;i2<256;i2++)
{
cout<<tab4[i2]<<endl;
if((i2+1)%16==0){system("pause");}
}
*/
///////////////////////////
//stage 4
for(int k=0;k<N;k++)
{
tab[k]={0,0};
}
tmp5=2*pi/(N/1);
for(int i=0;i<(N/16);i++)//16
{
for(int j=0;j<4;j++)
{
w1[0].real()=cos(0*16*j*tmp5);
w1[0].imag()=sin(0*16*j*tmp5);
w2[0].real()=cos(1*16*j*tmp5);
w2[0].imag()=sin(1*16*j*tmp5);
w3[0].real()=cos(2*16*j*tmp5);
w3[0].imag()=sin(2*16*j*tmp5);
w4[0].real()=cos(3*16*j*tmp5);
w4[0].imag()=sin(3*16*j*tmp5);
tmp1=w1[0]*tab2[N/(4*4)*i+0+4*j];
tmp2=w2[0]*tab2[N/(4*4)*i+1+4*j];
tmp3=w3[0]*tab2[N/(4*4)*i+2+4*j];
tmp4=w4[0]*tab2[N/(4*4)*i+3+4*j];
tab[N/(4*4)*i+0+4*j]=tmp1+tmp2+ tmp3+tmp4;
w5[0].real()=0;
w5[0].imag()=1;
tab[N/(4*4)*i+1+4*j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab[N/(4*4)*i+2+4*j]=tmp1-tmp2+ tmp3-tmp4;
w5[0].real()=0;
w5[0].imag()=-1;
tab[N/(4*4)*i+3+4*j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
//tab4[N/(4*4)*i+0+4*j] =0*16*j;
//tab4[N/(4*4)*i+1+4*j] =1*16*j;
//tab4[N/(4*4)*i+2+4*j]=2*16*j;
//tab4[N/(4*4)*i+3+4*j]=3*16*j;
}
}
/*
for(int i2=0;i2<256;i2++)
{
cout<<tab4[i2]<<endl;
if((i2+1)%16==0){system("pause");}
}
*/
for(int j=0;j<N;j++)
{
tab[j].real() =tab[j].real()*0.5;
tab[j].imag() =tab[j].imag()*0.5;
}
}
//////////////////////////////
//////////////////////////////
void fun_fourier_transform_FFT_radix_4_N_256_official_v2(int N,std::complex<double> tab[])
{
const double pi=3.141592653589793238462;
std::complex<double> tab2[256]={}; // tab2[]==N
std::complex<double> w1[1]={{1,1}};
std::complex<double> w2[1]={{1,1}};
std::complex<double> w3[1]={{1,1}};
std::complex<double> w4[1]={{1,1}};
std::complex<double> w5[1]={{1,1}};
std::complex<double> tmp1,tmp2,tmp3,tmp4;
double tmp5;
int flag=0;
int flag2=0;
int flag3=0;
int flag4=0;
int flg=0;
int tab4[256]={0};
int tab5[256]={0};
//radix-4
//stage 1
tmp5=2*pi/(N/1);
for(int i=0;i<(N/4);i++)//64
{
w1[0].real()=cos(0*tmp5);
w1[0].imag()=-sin(0*tmp5);
w2[0].real()=cos(0*tmp5);
w2[0].imag()=-sin(0*tmp5);
w3[0].real()=cos(0*tmp5);
w3[0].imag()=-sin(0*tmp5);
w4[0].real()=cos(0*tmp5);
w4[0].imag()=-sin(0*tmp5);
tmp1=w1[0]*tab[i+0];
tmp2=w2[0]*tab[i+N/4];
tmp3=w3[0]*tab[i+N/2];
tmp4=w4[0]*tab[i+3*N/4];
tab2[i] =tmp1+tmp2+tmp3+tmp4;
w5[0].real()=0;
w5[0].imag()=-1;
tab2[i+N/4] =tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab2[i+N/2] =tmp1-tmp2+tmp3-tmp4;
w5[0].real()=0;
w5[0].imag()=1;
tab2[i+3*N/4] =tmp1-tmp3+w5[0]*(tmp2-tmp4);
}
///////////////////////////////////
flag2=0;
flag3=0;
flag=0;
flg=0;
//stage 2
for(int k=0;k<N;k++)
{
tab[k]={0,0};
}
tmp5=2*pi/(N/1);
for(int i=0;i<(N/64);i++)//4
{
for(int j=0;j<(N/16);j++)//16
{
w1[0].real()=cos((0*i+i*j)*tmp5);
w1[0].imag()=-sin((0*i+i*j)*tmp5);
w2[0].real()=cos((1*16*i+i*j)*tmp5);
w2[0].imag()=-sin((1*16*i+i*j)*tmp5);
w3[0].real()=cos((2*16*i+i*j)*tmp5);
w3[0].imag()=-sin((2*16*i+i*j)*tmp5);
w4[0].real()=cos((3*16*i+i*j)*tmp5);
w4[0].imag()=-sin((3*16*i+i*j)*tmp5);
tmp1=w1[0]*tab2[0+j+N/(4*1)*i];
tmp2=w2[0]*tab2[N/(16*1)+j+N/(4*1)*i];
tmp3=w3[0]*tab2[N/(8*1) +j+N/(4*1)*i];
tmp4=w4[0]*tab2[3*N/(16*1)+j+N/(4*1)*i];
tab[0+j+N/(4*1)*i] =tmp1+tmp2+ tmp3+tmp4;
w5[0].real()=0;
w5[0].imag()=-1;
tab[N/(16*1)+j+N/(4*1)*i] =tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab[N/(8 *1)+j+N/(4*1)*i] =tmp1-tmp2+ tmp3-tmp4;
w5[0].real()=0;
w5[0].imag()=1;
tab[3*N/(16*1)+j+N/(4*1)*i]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
/*
tab4[0+j+N/(4*1)*i]=(0*i+i*j);
tab4[N/(16*1)+j+N/(4*1)*i]=(1*16*i+i*j);
tab4[N/(8 *1)+j+N/(4*1)*i]=(2*16*i+i*j);
tab4[3*N/(16*1)+j+N/(4*1)*i]=(3*16*i+i*j);
*/
if((j+1)%16==0){flag++;}
if((j+1)%16==0){flag2++;}
}
if((i+1)%4==0){flag3++;}
}
/*
for(int m=0;m<1;m++)//stage 2
{
for(int i=0;i<4;i++)
{
for(int j=0;j<16;j++)
{
tab5[0*N/(4*4)+N/(4)*i+j]=1*0*i*(j+16*0)+1*1*i*j;
tab5[1*N/(4*4)+N/(4)*i+j]=1*1*i*(j+16*1)+1*0*i*j;
tab5[2*N/(4*4)+N/(4)*i+j]=1*2*i*(j+16*1)-1*1*i*j;
tab5[3*N/(4*4)+N/(4)*i+j]=1*3*i*(j+16*1)-1*2*i*j;
}
}
}
for(int i2=0;i2<256;i2++)
{
cout<<"i2= "<<i2<<" tb4= "<<tab4[i2]<<" =tb5 "<<tab5[i2];
if(tab4[i2]==tab5[i2])cout<<" O..K ";
cout<<endl;
if((i2+1)%16==0){
system("pause");
}
}
*/
//////////////////////////////
flag=0;
flag2=0;
//stage 3
for(int k=0;k<N;k++)
{
tab2[k]={0,0};
}
tmp5=2*pi/(N/1);
for(int i=0;i<(N/16);i++)//16
{
for(int j=0;j<4;j++)
{
w1[0].real()=cos((i-4*flag)*(16*0+4*j)*tmp5);
w1[0].imag()=-sin((i-4*flag)*(16*0+4*j)*tmp5);
w2[0].real()=cos((i-4*flag)*(1*16+4*j)*tmp5);
w2[0].imag()=-sin((i-4*flag)*(1*16+4*j)*tmp5);
w3[0].real()=cos((i-4*flag)*(2*16+4*j)*tmp5);
w3[0].imag()=-sin((i-4*flag)*(2*16+4*j)*tmp5);
w4[0].real()=cos((i-4*flag)*(3*16+4*j)*tmp5);
w4[0].imag()=-sin((i-4*flag)*(3*16+4*j)*tmp5);
tmp1=w1[0]*tab[0+j+N/(4*4)*i];
tmp2=w2[0]*tab[N/(16*4)+j+N/(4*4)*i];
tmp3=w3[0]*tab[N/(8*4) +j+N/(4*4)*i];
tmp4=w4[0]*tab[3*N/(16*4)+j+N/(4*4)*i];
tab2[0+j+N/(4*4)*i] =tmp1+tmp2+ tmp3+tmp4;
w5[0].real()=0;
w5[0].imag()=-1;
tab2[N/(16*4)+j+N/(4*4)*i] =tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab2[N/(8*4)+j+N/(4*4)*i] =tmp1-tmp2+ tmp3-tmp4;
w5[0].real()=0;
w5[0].imag()=1;
tab2[3*N/(16*4)+j+N/(4*4)*i]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
/*
tab4[0+j+N/(4*4)*i] =(i-4*flag)*(16*0+4*j);
tab4[N/(16*4)+j+N/(4*4)*i] =(i-4*flag)*(1*16+4*j);
tab4[N/(8*4)+j+N/(4*4)*i]=(i-4*flag)*(2*16+4*j);
tab4[3*N/(16*4)+j+N/(4*4)*i]=(i-4*flag)*(3*16+4*j);
tab5[0*N/(16*4)+N/(4*4)*i+j]=4*0*(N/(16*4)*(i-4*flag))+4*(i-4*flag)*j;
tab5[1*N/(16*4)+N/(4*4)*i+j]=4*1*(N/(16*4)*(i-4*flag))+4*(i-4*flag)*j;
tab5[2*N/(16*4)+N/(4*4)*i+j]=4*2*(N/(16*4)*(i-4*flag))+4*(i-4*flag)*j;
tab5[3*N/(16*4)+N/(4*4)*i+j]=4*3*(N/(16*4)*(i-4*flag))+4*(i-4*flag)*j;
tab5[0*N/(4*4*4)+N/(4*4)*i+j]=(4*0*(N/(4*4*4)*i)+4*i*j)-0*N/(4)*flag-1*N/(4*4)*j*flag;
tab5[1*N/(4*4*4)+N/(4*4)*i+j]=(4*1*(N/(4*4*4)*i)+4*i*j)-1*N/(4)*flag-1*N/(4*4)*j*flag;
tab5[2*N/(4*4*4)+N/(4*4)*i+j]=(4*2*(N/(4*4*4)*i)+4*i*j)-2*N/(4)*flag-1*N/(4*4)*j*flag;
tab5[3*N/(4*4*4)+N/(4*4)*i+j]=(4*3*(N/(4*4*4)*i)+4*i*j)-3*N/(4)*flag-1*N/(4*4)*j*flag;
tab5[0*N/(4*4*4)+N/(4*4)*i+j]=(4*0*(N/(4*4*4)*(i-4*flag))+4*N/(4*4)*j*(i-4*flag))-1*4*15*(i-4*flag)*j;
tab5[1*N/(4*4*4)+N/(4*4)*i+j]=(4*1*(N/(4*4*4)*(i-4*flag))+4*N/(4*4)*j*(i-4*flag))-1*4*15*(i-4*flag)*j;
tab5[2*N/(4*4*4)+N/(4*4)*i+j]=(4*2*(N/(4*4*4)*(i-4*flag))+4*N/(4*4)*j*(i-4*flag))-1*4*15*(i-4*flag)*j;
tab5[3*N/(4*4*4)+N/(4*4)*i+j]=(4*3*(N/(4*4*4)*(i-4*flag))+4*N/(4*4)*j*(i-4*flag))-1*4*15*(i-4*flag)*j;
tab5[0*N/(4*4*4)+N/(4*4)*i+j]=4*0*(i-4*flag)*(j+4*0)+4*(i-4*flag)*j;
tab5[1*N/(4*4*4)+N/(4*4)*i+j]=4*1*(i-4*flag)*(j+4*1);
tab5[2*N/(4*4*4)+N/(4*4)*i+j]=4*2*(i-4*flag)*(j+4*1)-4*(i-4*flag)*j;
tab5[3*N/(4*4*4)+N/(4*4)*i+j]=4*3*(i-4*flag)*(j+4*1)-4*2*(i-4*flag)*j;
*/
}
if((i+1)%4==0){flag++;}
}
/*
for(int m=0;m<4;m++)//stage 3
{
for(int i=0;i<4;i++)
{
for(int j=0;j<4;j++)
{
tab5[0*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=4*0*i*(j+4*0)+4*1*i*j;
tab5[1*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=4*1*i*(j+4*1)+4*0*i*j;
tab5[2*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=4*2*i*(j+4*1)-4*1*i*j;
tab5[3*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=4*3*i*(j+4*1)-4*2*i*j;
}
}
}
for(int i2=0;i2<256;i2++)
{
cout<<i2<<" "<<tab4[i2]<<"="<<tab5[i2];
if(tab4[i2]==tab5[i2])cout<<" OK ";
cout<<endl;
if((i2+1)%16==0){
system("pause");
}
}
*/
///////////////////////////
flag=0;
//stage 4
for(int k=0;k<N;k++)
{
tab[k]={0,0};
}
tmp5=2*pi/(N/1);
for(int i=0;i<(N/16);i++)//16
{
for(int j=0;j<4;j++)
{
w1[0].real()=cos(0*16*j*tmp5);
w1[0].imag()=-sin(0*16*j*tmp5);
w2[0].real()=cos(1*16*j*tmp5);
w2[0].imag()=-sin(1*16*j*tmp5);
w3[0].real()=cos(2*16*j*tmp5);
w3[0].imag()=-sin(2*16*j*tmp5);
w4[0].real()=cos(3*16*j*tmp5);
w4[0].imag()=-sin(3*16*j*tmp5);
tmp1=w1[0]*tab2[N/(4*4)*i+0+4*j];
tmp2=w2[0]*tab2[N/(4*4)*i+1+4*j];
tmp3=w3[0]*tab2[N/(4*4)*i+2+4*j];
tmp4=w4[0]*tab2[N/(4*4)*i+3+4*j];
tab[N/(4*4)*i+0+4*j]=tmp1+tmp2+ tmp3+tmp4;
w5[0].real()=0;
w5[0].imag()=-1;
tab[N/(4*4)*i+1+4*j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab[N/(4*4)*i+2+4*j]=tmp1-tmp2+ tmp3-tmp4;
w5[0].real()=0;
w5[0].imag()=1;
tab[N/(4*4)*i+3+4*j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab4[N/(4*4)*i+0+4*j] =0*16*j;
tab4[N/(4*4)*i+1+4*j] =1*16*j;
tab4[N/(4*4)*i+2+4*j]=2*16*j;
tab4[N/(4*4)*i+3+4*j]=3*16*j;
//tab5[N/(4*4)*i+0+4*j]=16*0*(i-4*flag+j)-0*16*flag2;
//tab5[N/(4*4)*i+1+4*j]=16*1*(i-4*flag+j)-1*16*flag2;
//tab5[N/(4*4)*i+2+4*j]=16*2*(i-4*flag+j)-2*16*flag2;
//tab5[N/(4*4)*i+3+4*j]=16*3*(i-4*flag+j)-3*16*flag2;
//tab99[N/(4*4)*i+0+4*j]=4*0*(N/(4*4)*(i-4*flag)+0+4*(j-4*flag2));
//tab99[N/(4*4)*i+1+4*j]=4*1*(N/(4*4)*(i-4*flag)+0+4*(j-4*flag2));
//tab99[N/(4*4)*i+2+4*j]=4*2*(N/(4*4)*(i-4*flag)+0+4*(j-4*flag2));
//tab99[N/(4*4)*i+3+4*j]=4*3*(N/(4*4)*(i-4*flag)+0+4*(j-4*flag2));
}
if((i+1)%4==0){flag++;}
flag2++;
if((i+1)%4==0){flag2=0;}
}
/////////////////////////
//stage 4 bis
/*
flag=0;
flag3=0;
for(int i=0;i<(N/16);i++)//64
{
for(int j=0;j<1;j++)
{
tab5[0*N/(4*4*4*4)+N/(4*4*4)*i+j]=16*0*(N/(4*4*4*4)*(i-4*flag))+16*(i-4*flag)*j;
tab5[1*N/(4*4*4*4)+N/(4*4*4)*i+j]=16*1*(N/(4*4*4*4)*(i-4*flag))+16*(i-4*flag)*j;
tab5[2*N/(4*4*4*4)+N/(4*4*4)*i+j]=16*2*(N/(4*4*4*4)*(i-4*flag))+16*(i-4*flag)*j;
tab5[3*N/(4*4*4*4)+N/(4*4*4)*i+j]=16*3*(N/(4*4*4*4)*(i-4*flag))+16*(i-4*flag)*j;
tab5[0*N/(4*4*4*4)+N/(4*4*4)*i+j]=(16*0*(N/(4*4*4*4)*(i-4*flag))+16*N/(4*4*4)*j*(i-4*flag))-1*16*15*(i-4*flag)*j;
tab5[1*N/(4*4*4*4)+N/(4*4*4)*i+j]=(16*1*(N/(4*4*4*4)*(i-4*flag))+16*N/(4*4*4)*j*(i-4*flag))-1*16*15*(i-4*flag)*j;
tab5[2*N/(4*4*4*4)+N/(4*4*4)*i+j]=(16*2*(N/(4*4*4*4)*(i-4*flag))+16*N/(4*4*4)*j*(i-4*flag))-1*16*15*(i-4*flag)*j;
tab5[3*N/(4*4*4*4)+N/(4*4*4)*i+j]=(16*3*(N/(4*4*4*4)*(i-4*flag))+16*N/(4*4*4)*j*(i-4*flag))-1*16*15*(i-4*flag)*j;
//tab5[0*N/(4*4*4*4)+N/(4*4*4)*i+j]=16*0*(i-4*flag3)*(j+16*0)+16*1*(i-4*flag3)*j;
//tab5[1*N/(4*4*4*4)+N/(4*4*4)*i+j]=16*1*(i-4*flag3)*(j+16*1);
//tab5[2*N/(4*4*4*4)+N/(4*4*4)*i+j]=16*2*(i-4*flag3)*(j+16*1)-16*1*(i-4*flag3)*j;
//tab5[3*N/(4*4*4*4)+N/(4*4*4)*i+j]=16*3*(i-4*flag3)*(j+16*1)-16*2*(i-4*flag3)*j;
//cout<<0*N/(4*4*4*4)+N/(4*4*4)*i+j<<endl;
//system("pause");
}
if((i+1)%4==0){flag++;}
}
*/
/*
for(int m=0;m<16;m++)//stage 4
{
for(int i=0;i<4;i++)
{
for(int j=0;j<1;j++)
{
tab5[0*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=16*0*i*(j+1*0)+16*1*i*j;
tab5[1*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=16*1*i*(j+1*1)+16*0*i*j;
tab5[2*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=16*2*i*(j+1*1)-16*1*i*j;
tab5[3*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=16*3*i*(j+1*1)-16*2*i*j;
}
}
}
for(int i2=0;i2<256;i2++)
{
cout<<tab4[i2]<<"="<<tab5[i2];
if(tab4[i2]==tab5[i2])cout<<" O1K ";
cout<<endl;
if((i2+1)%16==0){
system("pause");
}
}
*/
for(int j=0;j<N;j++)
{
tab[j].real() =tab[j].real()*2/N;
tab[j].imag() =tab[j].imag()*2/N;
}
}
Subskrybuj:
Komentarze (Atom)