Mam do sprawdzenia modyfikację algorytmu mieszany radix - szybka transformacja fouriera FFT i iFFT który jest uważany za jeden z 10 najważniejszych algorytmów XX wieku do analizy sygnałów,dźwięków.
Pierwsza modyfikajca to kazaine algorytmowi obliczać przesunięcie fi,
a druga modyfikacja to zamiana obiektu liczb zespolonych na tablice[a][b] co daje 2 krotne przyspieszenie,
oraz sprawdzić czy dobrze sortuje tablicę wyników (void fun_inverse_table_FFT(int M,std::complex<double> tab[])) czy nie ma tam jakiś nadmiarowych operacji,oraz mnożę w algorytmie przez inne modyfikatory dla FFT i iFFT na końcu obliczeń.
Algorytm jest duży więc daję link do strony i kilka najważniejszych funkcji tutaj do posta
tutaj linki do stron:
universal program for mixed radix fast fourier transform FFT radix-11 * radix-7 * radix-5 * radix-4 * radix-3 * radix-2 for N= points algorithm example:
http://mixedradixfastfouriertransformifft.blogspot.com/
inverse discrete fourier transform iDFT 1D algorithm c++ source code
https://inverse-fourier-transformdftalgorithm.blogspot.com/
void fun_fourier_transform_FFT_mixed_radix(int N,std::complex<double> tab[])
{
//open-source
const double pi=3.141592653589793238462;
std::complex<double> *tab2 = new std::complex<double>[N]; // tab2[]==N
double tmp2;
double tmp3;
double tmp6;
double tmp7;
double tmp8;
double tmp9;
double tmp10;
double tmp11;
double fi2=fi; //shift only for stage nr 1
int i=0;
int rx5=5,rx4=4,rx3=3,rx2=2,rx7=7,rx11=11;
int stg[100]={};
int nb_stages=0;
int nb1,nb2,nb3,nb4,nb5_stg_previous,stg_first;
tmp6=2*pi/(N/1);
tmp2=2*pi/(2/1);
tmp3=2*pi/(3/1);
tmp10=2*pi/(4/1);
tmp8=2*pi/(5/1);
tmp7=2*pi/(7/1);
tmp11=2*pi/(11/1);
std::complex<double> z_rx2[2]={{1,0}};
std::complex<double> z_rx3[3]={{1,0}};
std::complex<double> z_rx4[2]={{1,0}};
std::complex<double> z_rx5[5]={{1,0}};
std::complex<double> z_rx7[7]={{1,0}};
std::complex<double> z_rx11[11]={{1,0}};
//radix 2 fundament
z_rx2[0].real(cos(0*tmp2));
z_rx2[0].imag(-sin(0*tmp2));
z_rx2[1].real(cos(1*tmp2));
z_rx2[1].imag(-sin(1*tmp2));
//radix 3 fundament
z_rx3[0].real(cos(0*tmp3));
z_rx3[0].imag(-sin(0*tmp3));
z_rx3[1].real(cos(1*tmp3));
z_rx3[1].imag(-sin(1*tmp3));
z_rx3[2].real(cos(2*tmp3));
z_rx3[2].imag(-sin(2*tmp3));
//radix 4 fundament
z_rx4[0].real(cos(0*tmp10));
z_rx4[0].imag(-sin(0*tmp10));
z_rx4[1].real(cos(1*tmp10));
z_rx4[1].imag(-sin(1*tmp10));
//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));
//radix 7 fundament
z_rx7[0].real(cos(0*tmp7));
z_rx7[0].imag(-sin(0*tmp7));
z_rx7[1].real(cos(1*tmp7));
z_rx7[1].imag(-sin(1*tmp7));
z_rx7[2].real(cos(2*tmp7));
z_rx7[2].imag(-sin(2*tmp7));
z_rx7[3].real(cos(3*tmp7));
z_rx7[3].imag(-sin(3*tmp7));
z_rx7[4].real(cos(4*tmp7));
z_rx7[4].imag(-sin(4*tmp7));
z_rx7[5].real(cos(5*tmp7));
z_rx7[5].imag(-sin(5*tmp7));
z_rx7[6].real(cos(6*tmp7));
z_rx7[6].imag(-sin(6*tmp7));
//radix 11 fundament
z_rx11[0].real(cos(0*tmp11));
z_rx11[0].imag(-sin(0*tmp11));
z_rx11[1].real(cos(1*tmp11));
z_rx11[1].imag(-sin(1*tmp11));
z_rx11[2].real(cos(2*tmp11));
z_rx11[2].imag(-sin(2*tmp11));
z_rx11[3].real(cos(3*tmp11));
z_rx11[3].imag(-sin(3*tmp11));
z_rx11[4].real(cos(4*tmp11));
z_rx11[4].imag(-sin(4*tmp11));
z_rx11[5].real(cos(5*tmp11));
z_rx11[5].imag(-sin(5*tmp11));
z_rx11[6].real(cos(6*tmp11));
z_rx11[6].imag(-sin(6*tmp11));
z_rx11[7].real(cos(7*tmp11));
z_rx11[7].imag(-sin(7*tmp11));
z_rx11[8].real(cos(8*tmp11));
z_rx11[8].imag(-sin(8*tmp11));
z_rx11[9].real(cos(9*tmp11));
z_rx11[9].imag(-sin(9*tmp11));
z_rx11[10].real(cos(10*tmp11));
z_rx11[10].imag(-sin(10*tmp11));
/*
for(int j=0;j<102;j++)
{
for(int i=0;i<102;i++)
{
if(((fabs(round(z_rx11[j].imag()*1000)/1000-round(z_rx11[i].imag()*1000)/1000)<0.001)
&&(fabs(round(z_rx11[j].real()*1000)/1000-round(z_rx11[i].real()*1000)/1000)<0.001))
||((fabs(round(z_rx11[j].imag()*1000)/1000+round(z_rx11[i].imag()*1000)/1000)<0.001)
&&(fabs(round(z_rx11[j].real()*1000)/1000+round(z_rx11[i].real()*1000)/1000)<0.001)))
{
cout<<j<<" "<<round(z_rx11[j].real()*1000)/1000<<" "<<round(z_rx11[j].imag()*1000)/1000<<" ";
cout<<i<<" "<<round(z_rx11[i].real()*1000)/1000<<" "<<round(z_rx11[i].imag()*1000)/1000<<endl;
}
//cout<<endl;
}
//system("pause");
}
*/
nb_stages=radix_base(N,stg);
if(nb_stages>=1){cout<<"N= "<<N<<endl;}
for(int m=1;m<=nb_stages;m++)
{
cout <<"stage:"<<m<<" = radix-"<<stg[m] << endl;
}
cout << endl;
if(nb_stages>=1)
{
stg_first=N/stg[1];
if(stg[1]==2)
{
fun_fourier_transform_FFT_radix_2_stg_first(tab,stg_first,fi2,z_rx2);
}
else if(stg[1]==3)
{
fun_fourier_transform_FFT_radix_3_stg_first(tab,stg_first,fi2,z_rx3);
}
else if(stg[1]==4)
{
fun_fourier_transform_FFT_radix_4_stg_first(tab,stg_first,fi2,z_rx4);
}
else if(stg[1]==5)
{
fun_fourier_transform_FFT_radix_5_stg_first(tab,stg_first,fi2,z_rx5);
}
else if(stg[1]==7)
{
fun_fourier_transform_FFT_radix_7_stg_first(tab,stg_first,fi2,z_rx7);
}
else if(stg[1]==11)
{
fun_fourier_transform_FFT_radix_11_stg_first(tab,stg_first,fi2,z_rx11);
}
else{}
nb1=N;
nb4=1;
for(int i=0;i<nb_stages-1;i++)
{
nb1=nb1/stg[0+i];
nb2=nb1/stg[1+i];
nb3=nb2/stg[2+i];
nb4=nb4*stg[0+i];
nb5_stg_previous=stg[1+i];
if(stg[i+2]==2)
{
fun_fourier_transform_FFT_radix_2_stg_rest(tab,nb1,nb2,nb3,nb4,nb5_stg_previous,tmp6,z_rx2);
}
else if(stg[i+2]==3)
{
fun_fourier_transform_FFT_radix_3_stg_rest(tab,nb1,nb2,nb3,nb4,nb5_stg_previous,tmp6,z_rx3);
}
else if(stg[i+2]==4)
{
fun_fourier_transform_FFT_radix_4_stg_rest(tab,nb1,nb2,nb3,nb4,nb5_stg_previous,tmp6,z_rx4);
}
else if(stg[i+2]==5)
{
fun_fourier_transform_FFT_radix_5_stg_rest(tab,nb1,nb2,nb3,nb4,nb5_stg_previous,tmp6,z_rx5);
}
else if(stg[i+2]==7)
{
fun_fourier_transform_FFT_radix_7_stg_rest(tab,nb1,nb2,nb3,nb4,nb5_stg_previous,tmp6,z_rx7);
}
else if(stg[i+2]==11)
{
fun_fourier_transform_FFT_radix_11_stg_rest(tab,nb1,nb2,nb3,nb4,nb5_stg_previous,tmp6,z_rx11);
}
else{}
}
}
//new:
for(int j=0;j<N;j++)
{
tab[j].real(tab[j].real()*2/N);
tab[j].imag(tab[j].imag()*2/N);
}
delete [] tab2;
}
///////////////////////////////////////////////////
void fun_fourier_transform_FFT_radix_2_stg_rest(std::complex<double> tab[],int nb1,int nb2,int nb3,int nb4,int nb5_stg_previous,double tmp6,std::complex<double> z[])
{
std::complex<double> tmp1,tmp2,tmp40,tmp50;
std::complex<double> w[2]={{1,0}};
double tmp100=0.0;
double tmp200=0.0;
double tmp300=0.0;
int nb_tmp1=0.0;
int nb_tmp2=0.0;
int nb_tmp3=0.0;
for(int b=0;b<nb5_stg_previous;b=b+1)
{
tmp300=nb4*b*tmp6;
tmp100=nb3*tmp300;
for(int j=0;j<nb3;j=j+1)
{
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(nb4*b*(1*nb3+j)*tmp5));
w[1].imag(-sin(1*tmp100+tmp200));
nb_tmp1=b*nb2+j;
for(int i=0;i<nb4;i=i+1)
{
nb_tmp2=i*nb1;
nb_tmp3=nb_tmp1+nb_tmp2;
tmp1=w[0]*tab[nb_tmp3+0*nb3];
tmp2=w[1]*tab[nb_tmp3+1*nb3];
tmp40=z[0]*(tmp1+tmp2);
tab[nb_tmp3+0*nb3]=tmp40;
tab[nb_tmp3+1*nb3]=z[0]*tmp1+z[1]*tmp2;
}
}
}
}
///////////////////////////////////////////////////////////////////////
void fun_inverse_table_FFT(int M,std::complex<double> tab[])
{
int rx5=5,rx4=4,rx3=3,rx2=2,rx7=7,rx11=11;
int stg[100]={};
int *tab8 = new int[M];
int *tab9 = new int[M];
std::complex<double> *tab11 = new std::complex<double>[M];
int nb_stages=0;
int nb1=0;
int nb2=0;
int nb3=0;
nb_stages=6;
nb_stages=radix_base(M,stg);
for(int i=0;i<M;i++)
{
//tab9[i]=tab2[i];
//tab8[i]=tab2[i];
tab9[i]=i;
tab8[i]=i;
}
nb3=1;
for(int op=nb_stages;op>=2;op--)
{
nb1=stg[op];
nb3=nb3*stg[op];
if(op==nb_stages)
{
nb2=stg[0];
}
else
{
nb2=nb2*stg[op+1];
}
for(int i=0,n=0,p=0;i<M;i=i+M/nb3,n++)
{
if(n>=nb1)
{
n=0,p=p+M/nb2;
}
for(int j=0,k=0;j<M/nb3;j++,k=k+nb1)
{
if(op%2==0)
{
tab8[i+j]=tab9[k+n+p];
}
else
{
tab9[i+j]=tab8[k+n+p];
}
}
}
}
for(int i=0;i<M;i++)
{
tab11[i]=tab[tab8[i]];
}
for(int i=0;i<M;i++)
{
tab[i]=tab11[i];
}
delete [] tab8;
delete [] tab9;
delete [] tab11;
}
///////////////////////////////////////
int radix_base(int N,int stg[])
{
int k=0;
double M=(double)N;
double epsilon1;
stg[0]=1;
//*flg2=0;
//cout<<"M= "<<M<<endl;
epsilon1=0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001;
for(int j=0;j<200;j++)
{
if(fmod(M,11)<=epsilon1)
{
k++;
M=M/11.0;
stg[k]=11;
}
else if(fmod(M,7)<=epsilon1)
{
k++;
M=M/7.0;
stg[k]=7;
}
else if(fmod(M,5)<=epsilon1)
{
k++;
M=M/5.0;
stg[k]=5;
}
else if(fmod(M,4)<=epsilon1)
{
k++;
M=M/4.0;
stg[k]=4;
}
else if(fmod(M,3)<=epsilon1)
{
k++;
M=M/3.0;
stg[k]=3;
}
else if(fmod(M,2)<=epsilon1)
{
k++;
M=M/2.0;
stg[k]=2;
}
else if(M>=1.0-epsilon1&&M<=1.0+epsilon1)
{
//*flag=*flag+1;
//cout<<"*flag= "<<*flag<<" N= "<<N<<endl;
break;
}
else
{
cout<<endl<< "Unsupported signal: N= "<<N<< endl;
if(k>0)
{
for(int m=1;m<=k;m++)
{
cout <<"stage:"<<m<<" = radix-"<<stg[m] << endl;
}
}
cout <<"stage:"<<k+1<<" = radix-??" << endl;
k=0;
break;
}
//*flg2=*flg2+1;
}
return k;
}
/////////////////////////////////////////////////////////////
void fun_fourier_transform_FFT_radix_2_stg_first(std::complex<double> tab[],int stg_first,double fi2,std::complex<double> z[])
{
std::complex<double> tmp1,tmp2,tmp40,tmp50;
std::complex<double> w[2]={{1,0}};
w[0].real(cos(0+fi2));
w[0].imag(-sin(0+fi2));
w[1].real(cos(0+fi2));
w[1].imag(-sin(0+fi2));
for(int i=0;i<stg_first;i=i+1)
{
tmp1=w[0]*tab[i+0*stg_first];
tmp2=w[1]*tab[i+1*stg_first];
tmp40=z[0]*(tmp1+tmp2);
tab[i+0*stg_first]=tmp40;
tab[i+1*stg_first]=z[0]*tmp1+z[1]*tmp2;
}
}
///////////////////////////////////////////////////////////////