wtorek, 20 czerwca 2017
mixed radix fast fourier transform FFT radix-5 * radix-3 * radix-2 * radix-4 * radix-5 for N=600 points example c++ + inverse table
//mixed radix fast fourier transform FFT radix-5 * radix-3 * radix-2 * radix-4 * radix-5 for N=600 points example c++ + inverse table
#include <complex>
using namespace std;
//complex number method:
void fun_fourier_transform_FFT_radix_3_stg_rest(std::complex<double> tab[],int j,int b,int nb1,int nb2,int nb3,int nb4,double tmp100,double
tmp300,std::complex<double> z[3]);
void fun_fourier_transform_FFT_radix_4_stg_rest(std::complex<double> tab[],int j,int b,int nb1,int nb2,int nb3,int nb4,double tmp100,double
tmp300,std::complex<double> z[2]);
void fun_fourier_transform_FFT_radix_5_stg_rest(std::complex<double> tab[],int j,int b,int nb1,int nb2,int nb3,int nb4,double tmp100,double
tmp300,std::complex<double> z[5]);
void fun_inverse_table_FFT_beta(int M,std::complex<double> tab[]);
void fun_fourier_transform_FFT_mixed_radix_N_600(int N,std::complex<double> tab[]);
void fun_fourier_transform_DFT(int N,std::complex<double> tab[]);
/////////////////////////////////////////////////////////////////////////////////////////////
void fun_fourier_transform_FFT_mixed_radix_N_600(int N,std::complex<double> tab[])
{
//author marcin matysek (r)ewertyn.PL
//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
//http://dsp.stackexchange.com/questions/3481/radix-4-fft-implementation
//https://community.arm.com/graphics/b/blog/posts/speeding-up-fast-fourier-transform-mixed-radix-on-mobile-arm-mali-gpu-by-means-of-
opencl---part-1
//book: "Cyfrowe przetwarzanie sygnalow" - Tomasz Zielinski it has quick version of radix-2 because it calculates less sin() and cos()
//Rabiner L.R., Gold B. Theory and application of digital signal processing p 378
const double pi=3.141592653589793238462;
std::complex<double> tab2[4096]={}; // 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> w34[1]={{1,0}};
std::complex<double> w40[1]={{1,0}};
std::complex<double> w5[1]={{1,0}};
std::complex<double> w6[1]={{1,0}};
std::complex<double> w50[1]={{1,0}};
std::complex<double> w60[1]={{1,0}};
std::complex<double> w7[1]={{1,0}};
std::complex<double> w11[1]={{1,0}};
std::complex<double> w12[1]={{1,0}};
std::complex<double> w13[1]={{1,0}};
std::complex<double> w14[1]={{1,0}};
std::complex<double> w15[1]={{1,0}};
std::complex<double> w16[1]={{1,0}};
std::complex<double> w17[1]={{1,0}};
std::complex<double> w18[1]={{1,0}};
std::complex<double> w19[1]={{1,0}};
std::complex<double> w20[1]={{1,0}};
std::complex<double> tmp1,tmp2,tmp3,tmp4,tmp11;
std::complex<double> tmp29,tmp30,tmp31,tmp32,tmp33,tmp34,tmp35,tmp36;
double tmp5;
double tmp6;
double tmp7;
double tmp8;
double tmp118;
double tmp100;
double tmp200;
double tmp300;
double fi2=fi;
int rx5=5,rx4=4,rx3=3,rx2=2;
int stg[100]={};
int nb1,nb2,nb3,nb4;
stg[0]=1;
stg[1]=rx5;
stg[2]=rx3;
stg[3]=rx2;
stg[4]=rx4;
stg[5]=rx5;
tmp5=2*pi/(N/1);
tmp6=2*pi/(N/1);
tmp7=2*pi/(3/1);
tmp8=2*pi/(5/1);
tmp118=2*pi/(4/1);
//radix 3 fundament
w34[0].real()=cos(0*tmp7);
w34[0].imag()=-sin(0*tmp7);
w5[0].real()=cos(1*tmp7);
w5[0].imag()=-sin(1*tmp7);
w6[0].real()=cos(2*tmp7);
w6[0].imag()=-sin(2*tmp7);
w7[0].real()=cos(4*tmp7);
w7[0].imag()=-sin(4*tmp7);
std::complex<double> z_rx3[3]={{1,0}};
std::complex<double> z_rx4[2]={{1,0}};
std::complex<double> z[2]={{1,0}};
std::complex<double> z_rx5[10]={{1,0}};
std::complex<double> w[11]={{1,0}};
//radix 3 fundament
z_rx3[0].real()=cos(0*tmp7);
z_rx3[0].imag()=-sin(0*tmp7);
z_rx3[1].real()=cos(1*tmp7);
z_rx3[1].imag()=-sin(1*tmp7);
z_rx3[2].real()=cos(2*tmp7);
z_rx3[2].imag()=-sin(2*tmp7);
//radix 4 fundament
z_rx4[0].real()=cos(0*tmp118);
z_rx4[0].imag()=-sin(0*tmp118);
z_rx4[1].real()=cos(1*tmp118);
z_rx4[1].imag()=-sin(1*tmp118);
//radix 5 fundament
z_rx5[0].real()=cos(0*tmp8);
z_rx5[0].imag()=-sin(0*tmp8);
z_rx5[1].real()=cos(1*tmp8);
z_rx5[1].imag()=-sin(1*tmp8);
z_rx5[2].real()=cos(2*tmp8);
z_rx5[2].imag()=-sin(2*tmp8);
z_rx5[3].real()=cos(3*tmp8);
z_rx5[3].imag()=-sin(3*tmp8);
z_rx5[4].real()=cos(4*tmp8);
z_rx5[4].imag()=-sin(4*tmp8);
/*
z_rx5[5].real()=cos(6*tmp8);
z_rx5[5].imag()=-sin(6*tmp8);
z_rx5[6].real()=cos(8*tmp8);
z_rx5[6].imag()=-sin(8*tmp8);
z_rx5[7].real()=cos(9*tmp8);
z_rx5[7].imag()=-sin(9*tmp8);
z_rx5[8].real()=cos(12*tmp8);
z_rx5[8].imag()=-sin(12*tmp8);
z_rx5[9].real()=cos(16*tmp8);
z_rx5[9].imag()=-sin(16*tmp8);
*/
for(int j=0;j<2;j++)
{
for(int i=0;i<2;i++)
{
if(((fabs(round(z_rx4[j].imag()*1000)/1000-round(z_rx4[i].imag()*1000)/1000)<0.001)
&&(fabs(round(z_rx4[j].real()*1000)/1000-round(z_rx4[i].real()*1000)/1000)<0.001))
||((fabs(round(z_rx4[j].imag()*1000)/1000+round(z_rx4[i].imag()*1000)/1000)<0.001)
&&(fabs(round(z_rx4[j].real()*1000)/1000+round(z_rx4[i].real()*1000)/1000)<0.001)))
{
cout<<j<<" "<<round(z_rx4[j].real()*1000)/1000<<" "<<round(z_rx4[j].imag()*1000)/1000<<" ";
cout<<i<<" "<<round(z_rx4[i].real()*1000)/1000<<" "<<round(z_rx4[i].imag()*1000)/1000<<endl;
}
cout<<endl;
}
}
system("pause");
w1[0].real()=cos(0+fi2);
w1[0].imag()=-sin(0+fi2);
w2[0].real()=cos(0+fi2);
w2[0].imag()=-sin(0+fi2);
w3[0].real()=cos(0+fi2);
w3[0].imag()=-sin(0+fi2);
w4[0].real()=cos(0+fi2);
w4[0].imag()=-sin(0+fi2);
w40[0].real()=cos(0+fi2);
w40[0].imag()=-sin(0+fi2);
for(int i=0;i<N/(stg[1]);i++)
{
tmp1=w1[0]*tab[i+0];
tmp2=w2[0]*tab[i+1*N/(stg[1])];
tmp3=w3[0]*tab[i+2*N/(stg[1])];
tmp4=w4[0]*tab[i+3*N/(stg[1])];
tmp11=w40[0]*tab[i+4*N/(stg[1])];
tmp33=z_rx5[0]*tmp1;
tmp34=z_rx5[0]*(tmp2+tmp3+tmp4+tmp11);
//radix-5
tab2[i] =tmp33+tmp34;
tab2[i+1*N/(stg[1])] =tmp33+z_rx5[1]*tmp2+z_rx5[2]*tmp3+z_rx5[3]*tmp4+z_rx5[4]*tmp11;
tab2[i+2*N/(stg[1])] =tmp33+z_rx5[2]*tmp2+z_rx5[4]*tmp3+z_rx5[1]*tmp4+z_rx5[3]*tmp11;
tab2[i+3*N/(stg[1])] =tmp33+z_rx5[3]*tmp2+z_rx5[1]*tmp3+z_rx5[4]*tmp4+z_rx5[2]*tmp11;
tab2[i+4*N/(stg[1])] =tmp33+z_rx5[4]*tmp2+z_rx5[3]*tmp3+z_rx5[2]*tmp4+z_rx5[1]*tmp11;
}
//stage 2 radix-3
int b=0;
int i=0;
nb1=N;
nb4=1;
nb1=nb1/stg[0+i];
nb2=nb1/stg[1+i];
nb3=nb2/stg[2+i];
nb4=nb4*stg[0+i];
for(int b=0;b<stg[1+i];b=b+1)
{
tmp300=nb4*b*tmp6;
tmp100=nb3*tmp300;
for(int j=0;j<nb3;j=j+1)
{
fun_fourier_transform_FFT_radix_3_stg_rest(tab2,j,b,nb1,nb2,nb3,nb4,tmp100,tmp300,z_rx3);
}
}
//stage 3 radix-2
b=0;
i=1;
nb1=nb1/stg[0+i];
nb2=nb1/stg[1+i];
nb3=nb2/stg[2+i];
nb4=nb4*stg[0+i];
for(int b=0;b<stg[1+i];b=b+1)
{
tmp100=nb4*b*nb3*tmp6;
for(int j=0;j<nb3;j=j+1)
{
tmp200=nb4*b*j*tmp6;
w1[0].real()= cos(0*tmp100+tmp200);
w1[0].imag()=-sin(0*tmp100+tmp200);
w2[0].real()= cos(1*tmp100+tmp200);
//w2[0].imag()=-sin(nb4*b*(1*nb3+j)*tmp5);
w2[0].imag()=-sin(1*tmp100+tmp200);
for(int i=0;i<nb4;i=i+1)
{
tmp1=w1[0]*tab2[i*nb1+b*nb2+0*nb3+j];
tmp2=w2[0]*tab2[i*nb1+b*nb2+1*nb3+j];
tab2[i*nb1+b*nb2+0*nb3+j]=tmp1+tmp2;
tab2[i*nb1+b*nb2+1*nb3+j]=tmp1-tmp2;
}
}
}
//stage 4 radix-4
w50[0].real()=0;
w50[0].imag()=-1;
w60[0].real()=0;
w60[0].imag()=1;
b=0;
i=2;
nb1=nb1/stg[0+i];
nb2=nb1/stg[1+i];
nb3=nb2/stg[2+i];
nb4=nb4*stg[0+i];
std::complex<double> tmp101,tmp102,tmp103,tmp104;
std::complex<double> tmp111,tmp112,tmp113,tmp114;
// std::complex<double> w[4]={{1,0}};
for(int b=0;b<stg[1+i];b=b+1)
{
tmp300=nb4*b*tmp6;
tmp100=nb4*b*nb3*tmp6;
for(int j=0;j<nb3;j=j+1)
{
fun_fourier_transform_FFT_radix_4_stg_rest(tab2,j,b,nb1,nb2,nb3,nb4,tmp100,tmp300,z_rx4);
}
}
//stage 5 radix-3
b=0;
i=3;
nb1=nb1/stg[0+i];
nb2=nb1/stg[1+i];
nb3=nb2/stg[2+i];
nb4=nb4*stg[0+i];
for(int b=0;b<stg[1+i];b=b+1)
{
tmp300=nb4*b*tmp6;
tmp100=nb4*b*nb3*tmp6;
for(int j=0;j<nb3;j=j+1)
{
fun_fourier_transform_FFT_radix_5_stg_rest(tab2,j,b,nb1,nb2,nb3,nb4,tmp100,tmp300,z_rx5);
}
}
//new:
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_fourier_transform_FFT_radix_3_stg_rest(std::complex<double> tab[],int j,int b,int nb1,int nb2,int nb3,int nb4,double
tmp100,double tmp300,std::complex<double> z[3])
{
std::complex<double> tmp1,tmp2,tmp3,tmp4,tmp5;
std::complex<double> w[3]={{1,0}};
double tmp200=0.0;
std::complex<double> tab_tmp[3]={{}};
tmp200=j*tmp300;
w[0].real()=cos(0*tmp100+tmp200);
w[0].imag()=-sin(0*tmp100+tmp200);
w[1].real()=cos(1*tmp100+tmp200);
w[1].imag()=-sin(1*tmp100+tmp200);
//w[2].real()=cos(nb4*b*(2*nb3+j)*tmp6);
w[2].real()=cos(2*tmp100+tmp200);
w[2].imag()=-sin(2*tmp100+tmp200);
for(int i=0;i<nb4;i=i+1)
{
tmp1=w[0]*tab[i*nb1+b*nb2+0*nb3+j];
tmp2=w[1]*tab[i*nb1+b*nb2+1*nb3+j];
tmp3=w[2]*tab[i*nb1+b*nb2+2*nb3+j];
tmp4=z[0]*tmp1;
tmp5=z[0]*(tmp2+tmp3);
//radix-3
tab_tmp[0] = tmp4+tmp5;
tab_tmp[1] = tmp4+z[1]*tmp2+z[2]*tmp3;
tab_tmp[2] = tmp4+z[2]*tmp2+z[1]*tmp3;
tab[i*nb1+b*nb2+0*nb3+j] = tab_tmp[0];
tab[i*nb1+b*nb2+1*nb3+j] = tab_tmp[1];
tab[i*nb1+b*nb2+2*nb3+j] = tab_tmp[2];
}
}
///////////////////////////////////////////////////////////////
void fun_fourier_transform_FFT_radix_5_stg_rest(std::complex<double> tab[],int j,int b,int nb1,int nb2,int nb3,int nb4,double
tmp100,double tmp300,std::complex<double> z[5])
{
std::complex<double> tmp1,tmp2,tmp3,tmp4,tmp5,tmp10,tmp20;
std::complex<double> w[5]={{1,0}};
double tmp200=0.0;
std::complex<double> tab_tmp[5]={{}};
tmp200=j*tmp300;
w[0].real()=cos(0*tmp100+tmp200);
w[0].imag()=-sin(0*tmp100+tmp200);
w[1].real()=cos(1*tmp100+tmp200);
w[1].imag()=-sin(1*tmp100+tmp200);
//w[2].real()=cos(nb4*b*(2*nb3+j)*tmp6);
w[2].real()=cos(2*tmp100+tmp200);
w[2].imag()=-sin(2*tmp100+tmp200);
w[3].real()=cos(3*tmp100+tmp200);
w[3].imag()=-sin(3*tmp100+tmp200);
w[4].real()=cos(4*tmp100+tmp200);
w[4].imag()=-sin(4*tmp100+tmp200);
for(int i=0;i<nb4;i=i+1)
{
tmp1=w[0]*tab[i*nb1+b*nb2+0*nb3+j];
tmp2=w[1]*tab[i*nb1+b*nb2+1*nb3+j];
tmp3=w[2]*tab[i*nb1+b*nb2+2*nb3+j];
tmp4=w[3]*tab[i*nb1+b*nb2+3*nb3+j];
tmp5=w[4]*tab[i*nb1+b*nb2+4*nb3+j];
tmp10=z[0]*tmp1;
tmp20=z[0]*(tmp1+tmp2+tmp3+tmp4+tmp5);
//radix-5
tab[i*nb1+b*nb2+0*nb3+j] =tmp20;
tab[i*nb1+b*nb2+1*nb3+j] =tmp10+z[1]*tmp2+z[2]*tmp3+z[3]*tmp4+z[4]*tmp5;
tab[i*nb1+b*nb2+2*nb3+j] =tmp10+z[2]*tmp2+z[4]*tmp3+z[1]*tmp4+z[3]*tmp5;
tab[i*nb1+b*nb2+3*nb3+j] =tmp10+z[3]*tmp2+z[1]*tmp3+z[4]*tmp4+z[2]*tmp5;
tab[i*nb1+b*nb2+4*nb3+j] =tmp10+z[4]*tmp2+z[3]*tmp3+z[2]*tmp4+z[1]*tmp5;
}
}
///////////////////////////////////////////////////////////////////////
void fun_fourier_transform_FFT_radix_4_stg_rest(std::complex<double> tab[],int j,int b,int nb1,int nb2,int nb3,int nb4,double
tmp100,double tmp300,std::complex<double> z[2])
{
std::complex<double> tmp1,tmp2,tmp3,tmp4;
std::complex<double> tmp101,tmp102,tmp103,tmp104;
std::complex<double> tmp111,tmp112,tmp113,tmp114;
std::complex<double> w[4]={{1,0}};
double tmp200=0.0;
tmp200=j*tmp300;
w[0].real()=cos(0*tmp100+tmp200);
w[0].imag()=-sin(0*tmp100+tmp200);
w[1].real()=cos(1*tmp100+tmp200);
w[1].imag()=-sin(1*tmp100+tmp200);
w[2].real()=cos(2*tmp100+tmp200);
w[2].imag()=-sin(2*tmp100+tmp200);
//w[3].real()=cos(nb4*b*(3*nb3+j)*tmp5);
w[3].real()=cos(3*tmp100+tmp200);
w[3].imag()=-sin(3*tmp100+tmp200);
for(int i=0;i<nb4;i=i+1)
{
tmp1=w[0]*tab[i*nb1+b*nb2+0*nb3+j];
tmp2=w[1]*tab[i*nb1+b*nb2+1*nb3+j];
tmp3=w[2]*tab[i*nb1+b*nb2+2*nb3+j];
tmp4=w[3]*tab[i*nb1+b*nb2+3*nb3+j];
tmp101=tmp2-tmp4;
tmp102=tmp2+tmp4;
tmp103=tmp1-tmp3;
tmp104=tmp1+tmp3;
tmp111=z[0]*(tmp104+tmp102);
tmp112=z[0]*tmp103;
tmp113=z[0]*(tmp104-tmp102);
tmp114=z[0]*tmp103;
//radix-4
tab[i*nb1+b*nb2+0*nb3+j] =tmp111;
tab[i*nb1+b*nb2+1*nb3+j] =tmp112+z[1]*tmp101;
tab[i*nb1+b*nb2+2*nb3+j] =tmp113;
tab[i*nb1+b*nb2+3*nb3+j] =tmp114-z[1]*tmp101;
}
}
/////////////////////////////////////////////////////////////////
void fun_inverse_table_FFT_beta(int M,std::complex<double> tab[])
{
int rx5=5,rx4=4,rx3=3,rx2=2;
int stg[100]={};
int tab9[1024]={};
int tab8[1024]={};
std::complex<double> tab11[1024]={};
stg[0]=1;
stg[1]=rx5;
stg[2]=rx3;
stg[3]=rx2;
stg[4]=rx4;
stg[5]=rx5;
for(int i=0;i<M;i++)
{
//tab9[i]=tab2[i];
tab9[i]=i;
}
//dla radix5*radix3*radix2*radix4*radix5 N=600 points "bit/digit" inverse
//stage 5
for(int i=0,n=0,p=0;i<M;i=i+M/(stg[5]),n++)
{
if(n>=stg[5]){n=0,p=p+M/(stg[0]);}
for(int j=0,k=0;j<M/(stg[5]);j++,k=k+stg[5])
{
tab8[i+j]=tab9[k+n+p];
}
}
//stage 4
for(int i=0,n=0,p=0;i<M;i=i+M/(stg[5]*stg[4]),n++)
{
if(n>=stg[4]){n=0,p=p+M/(stg[0]*stg[5]);}
for(int j=0,k=0;j<M/(stg[5]*stg[4]);j++,k=k+stg[4])
{
tab9[i+j]=tab8[k+n+p];
}
}
//stage 3
for(int i=0,n=0,p=0;i<M;i=i+M/(stg[5]*stg[4]*stg[3]),n++)
{
if(n>=stg[3]){n=0,p=p+M/(stg[0]*stg[5]*stg[4]);}
for(int j=0,k=0;j<M/(stg[5]*stg[4]*stg[3]);j++,k=k+stg[3])
{
tab8[i+j]=tab9[k+n+p];
}
}
//stage 2
for(int i=0,n=0,p=0;i<M;i=i+M/(stg[5]*stg[4]*stg[3]*stg[2]),n++)
{
if(n>=stg[2]){n=0,p=p+M/(stg[0]*stg[5]*stg[4]*stg[3]);}
for(int j=0,k=0;j<M/(stg[5]*stg[4]*stg[3]*stg[2]);j++,k=k+stg[2])
{
tab9[i+j]=tab8[k+n+p];
}
}
//system("pause");
for(int i=0;i<M;i++)
{
tab11[i]=tab[tab9[i]];
}
for(int i=0;i<M;i++)
{
tab[i]=tab11[i];
}
}
////////////////////////////////////////
void fun_fourier_transform_DFT(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;
}
}
//////////////////
Subskrybuj:
Komentarze do posta (Atom)
Brak komentarzy:
Prześlij komentarz