czwartek, 13 kwietnia 2017
fast inverse fourier transform iFFT radix-4 for N=1024 algorithm c++ source code implementation
//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
//author marcin matysek (r)ewertyn.PL
//fast inverse fourier transform iFFT radix-4 for N=1024 algorithm c++ source code implementation
//szybka odwrotna transformacja fouriera iFFT radix-4 dla N=1024 algorytm c++ kod źródłowy implementacja
#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_1024(int N,std::complex<double> tab[]);
void fun_inverse_fourier_transform_FFT_radix_4_N_1024(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=1024;
std::complex<double> tab2[4096]={};
std::complex<double> tab3[4096]={};
for(int i=0;i<N;i++)
{
tab2[i].real()=i+1;
tab2[i].imag()=N+i+1;
tab3[i].real()=i+1;
tab3[i].imag()=N+i+1;
}
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;
system("pause");
clock_t start = clock();
fun_fourier_transform_DFT_method5_full_complex(N,tab3);
//////////////////////////////////////////////////////////
fun_fourier_transform_FFT_radix_4_N_1024(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_1024(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
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;
}
}
///////////////////////////////////////////////////
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_fourier_transform_FFT_radix_4_N_1024(int N,std::complex<double> tab[])
{
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> w5[1]={{1,0}};
std::complex<double> w6[1]={{1,0}};
std::complex<double> tmp1,tmp2,tmp3,tmp4;
double tmp5;
w5[0].real()=0;
w5[0].imag()=-1;
w6[0].real()=0;
w6[0].imag()=1;
tmp5=2*pi/(N/1);
//stage 1
for(int i=0;i<(N/4);i++)//256
{
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];
//radix-4
tab2[i] =tmp1+tmp2+tmp3+tmp4;
tab2[i+N/4] =tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab2[i+N/2] =tmp1-tmp2+tmp3-tmp4;
tab2[i+3*N/4] =tmp1-tmp3+w6[0]*(tmp2-tmp4);
}
/////////////////////////////////////////////////////////////////////////////////
for(int m=0;m<1;m++)//stage 2
{
for(int i=0;i<4;i++)
{
for(int j=0;j<4*4*4;j++)
{
w1[0].real()= cos((1*0*i*(j+4*4*4*0)+1*1*i*j)*tmp5);
w1[0].imag()=-sin((1*0*i*(j+4*4*4*0)+1*1*i*j)*tmp5);
w2[0].real()= cos((1*1*i*(j+4*4*4*1)+1*0*i*j)*tmp5);
w2[0].imag()=-sin((1*1*i*(j+4*4*4*1)+1*0*i*j)*tmp5);
w3[0].real()= cos((1*2*i*(j+4*4*4*1)-1*1*i*j)*tmp5);
w3[0].imag()=-sin((1*2*i*(j+4*4*4*1)-1*1*i*j)*tmp5);
w4[0].real()= cos((1*3*i*(j+4*4*4*1)-1*2*i*j)*tmp5);
w4[0].imag()=-sin((1*3*i*(j+4*4*4*1)-1*2*i*j)*tmp5);
tmp1=w1[0]*tab2[0*N/(4*4)+N/(4)*i+N/(1)*m+j];
tmp2=w2[0]*tab2[1*N/(4*4)+N/(4)*i+N/(1)*m+j];
tmp3=w3[0]*tab2[2*N/(4*4)+N/(4)*i+N/(1)*m+j];
tmp4=w4[0]*tab2[3*N/(4*4)+N/(4)*i+N/(1)*m+j];
//radix-4
tab[0*N/(4*4)+N/(4)*i+N/(1)*m+j]=tmp1+tmp2+ tmp3+tmp4;
tab[1*N/(4*4)+N/(4)*i+N/(1)*m+j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab[2*N/(4*4)+N/(4)*i+N/(1)*m+j]=tmp1-tmp2+ tmp3-tmp4;
tab[3*N/(4*4)+N/(4)*i+N/(1)*m+j]=tmp1-tmp3+w6[0]*(tmp2-tmp4);
}
}
}
////////////////////////////////////////////////////////////////////////////////
for(int m=0;m<4;m++)//stage 3
{
for(int i=0;i<4;i++)
{
for(int j=0;j<4*4;j++)
{
w1[0].real()= cos((4*0*i*(j+4*4*0)+4*1*i*j)*tmp5);
w1[0].imag()=-sin((4*0*i*(j+4*4*0)+4*1*i*j)*tmp5);
w2[0].real()= cos((4*1*i*(j+4*4*1)+4*0*i*j)*tmp5);
w2[0].imag()=-sin((4*1*i*(j+4*4*1)+4*0*i*j)*tmp5);
w3[0].real()= cos((4*2*i*(j+4*4*1)-4*1*i*j)*tmp5);
w3[0].imag()=-sin((4*2*i*(j+4*4*1)-4*1*i*j)*tmp5);
w4[0].real()= cos((4*3*i*(j+4*4*1)-4*2*i*j)*tmp5);
w4[0].imag()=-sin((4*3*i*(j+4*4*1)-4*2*i*j)*tmp5);
tmp1=w1[0]*tab[0*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j];
tmp2=w2[0]*tab[1*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j];
tmp3=w3[0]*tab[2*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j];
tmp4=w4[0]*tab[3*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j];
//radix-4
tab2[0*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=tmp1+tmp2+ tmp3+tmp4;
tab2[1*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab2[2*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=tmp1-tmp2+ tmp3-tmp4;
tab2[3*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=tmp1-tmp3+w6[0]*(tmp2-tmp4);
}
}
}
////////////////////////////////////////////////////////////////////////////////
for(int m=0;m<4*4;m++)//stage 4
{
for(int i=0;i<4;i++)
{
for(int j=0;j<4;j++)
{
w1[0].real()= cos((4*4*0*i*(j+4*0)+4*4*1*i*j)*tmp5);
w1[0].imag()=-sin((4*4*0*i*(j+4*0)+4*4*1*i*j)*tmp5);
w2[0].real()= cos((4*4*1*i*(j+4*1)+4*4*0*i*j)*tmp5);
w2[0].imag()=-sin((4*4*1*i*(j+4*1)+4*4*0*i*j)*tmp5);
w3[0].real()= cos((4*4*2*i*(j+4*1)-4*4*1*i*j)*tmp5);
w3[0].imag()=-sin((4*4*2*i*(j+4*1)-4*4*1*i*j)*tmp5);
w4[0].real()= cos((4*4*3*i*(j+4*1)-4*4*2*i*j)*tmp5);
w4[0].imag()=-sin((4*4*3*i*(j+4*1)-4*4*2*i*j)*tmp5);
tmp1=w1[0]*tab2[0*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j];
tmp2=w2[0]*tab2[1*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j];
tmp3=w3[0]*tab2[2*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j];
tmp4=w4[0]*tab2[3*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j];
//radix-4
tab[0*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=tmp1+tmp2+ tmp3+tmp4;
tab[1*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab[2*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=tmp1-tmp2+ tmp3-tmp4;
tab[3*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=tmp1-tmp3+w6[0]*(tmp2-tmp4);
}
}
}
///////////////////////////////////////////////////////////////////////////////
for(int m=0;m<4*4*4;m++)//stage 5
{
for(int i=0;i<4;i++)
{
for(int j=0;j<1;j++)
{
w1[0].real()= cos((4*4*4*0*i*(j+1*0)+4*4*4*1*i*j)*tmp5);
w1[0].imag()=-sin((4*4*4*0*i*(j+1*0)+4*4*4*1*i*j)*tmp5);
w2[0].real()= cos((4*4*4*1*i*(j+1*1)+4*4*4*0*i*j)*tmp5);
w2[0].imag()=-sin((4*4*4*1*i*(j+1*1)+4*4*4*0*i*j)*tmp5);
w3[0].real()= cos((4*4*4*2*i*(j+1*1)-4*4*4*1*i*j)*tmp5);
w3[0].imag()=-sin((4*4*4*2*i*(j+1*1)-4*4*4*1*i*j)*tmp5);
w4[0].real()= cos((4*4*4*3*i*(j+1*1)-4*4*4*2*i*j)*tmp5);
w4[0].imag()=-sin((4*4*4*3*i*(j+1*1)-4*4*4*2*i*j)*tmp5);
tmp1=w1[0]*tab[0*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j];
tmp2=w2[0]*tab[1*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j];
tmp3=w3[0]*tab[2*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j];
tmp4=w4[0]*tab[3*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j];
//radix-4
tab2[0*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j]=tmp1+tmp2+ tmp3+tmp4;
tab2[1*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab2[2*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j]=tmp1-tmp2+ tmp3-tmp4;
tab2[3*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j]=tmp1-tmp3+w6[0]*(tmp2-tmp4);
}
}
}
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_1024(int N,std::complex<double> tab[])
{
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> w5[1]={{1,0}};
std::complex<double> w6[1]={{1,0}};
std::complex<double> tmp1,tmp2,tmp3,tmp4;
double tmp5;
w5[0].real()=0;
w5[0].imag()=1;
w6[0].real()=0;
w6[0].imag()=-1;
tmp5=2*pi/(N/1);
//stage 1
for(int i=0;i<(N/4);i++)//256
{
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];
//radix-4
tab2[i] =tmp1+tmp2+tmp3+tmp4;
tab2[i+N/4] =tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab2[i+N/2] =tmp1-tmp2+tmp3-tmp4;
tab2[i+3*N/4] =tmp1-tmp3+w6[0]*(tmp2-tmp4);
}
/////////////////////////////////////////////////////////////////////////////////
for(int m=0;m<1;m++)//stage 2
{
for(int i=0;i<4;i++)
{
for(int j=0;j<4*4*4;j++)
{
w1[0].real()= cos((1*0*i*(j+4*4*4*0)+1*1*i*j)*tmp5);
w1[0].imag()=sin((1*0*i*(j+4*4*4*0)+1*1*i*j)*tmp5);
w2[0].real()= cos((1*1*i*(j+4*4*4*1)+1*0*i*j)*tmp5);
w2[0].imag()=sin((1*1*i*(j+4*4*4*1)+1*0*i*j)*tmp5);
w3[0].real()= cos((1*2*i*(j+4*4*4*1)-1*1*i*j)*tmp5);
w3[0].imag()=sin((1*2*i*(j+4*4*4*1)-1*1*i*j)*tmp5);
w4[0].real()= cos((1*3*i*(j+4*4*4*1)-1*2*i*j)*tmp5);
w4[0].imag()=sin((1*3*i*(j+4*4*4*1)-1*2*i*j)*tmp5);
tmp1=w1[0]*tab2[0*N/(4*4)+N/(4)*i+N/(1)*m+j];
tmp2=w2[0]*tab2[1*N/(4*4)+N/(4)*i+N/(1)*m+j];
tmp3=w3[0]*tab2[2*N/(4*4)+N/(4)*i+N/(1)*m+j];
tmp4=w4[0]*tab2[3*N/(4*4)+N/(4)*i+N/(1)*m+j];
//radix-4
tab[0*N/(4*4)+N/(4)*i+N/(1)*m+j]=tmp1+tmp2+ tmp3+tmp4;
tab[1*N/(4*4)+N/(4)*i+N/(1)*m+j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab[2*N/(4*4)+N/(4)*i+N/(1)*m+j]=tmp1-tmp2+ tmp3-tmp4;
tab[3*N/(4*4)+N/(4)*i+N/(1)*m+j]=tmp1-tmp3+w6[0]*(tmp2-tmp4);
}
}
}
////////////////////////////////////////////////////////////////////////////////
for(int m=0;m<4;m++)//stage 3
{
for(int i=0;i<4;i++)
{
for(int j=0;j<4*4;j++)
{
w1[0].real()= cos((4*0*i*(j+4*4*0)+4*1*i*j)*tmp5);
w1[0].imag()=sin((4*0*i*(j+4*4*0)+4*1*i*j)*tmp5);
w2[0].real()= cos((4*1*i*(j+4*4*1)+4*0*i*j)*tmp5);
w2[0].imag()=sin((4*1*i*(j+4*4*1)+4*0*i*j)*tmp5);
w3[0].real()= cos((4*2*i*(j+4*4*1)-4*1*i*j)*tmp5);
w3[0].imag()=sin((4*2*i*(j+4*4*1)-4*1*i*j)*tmp5);
w4[0].real()= cos((4*3*i*(j+4*4*1)-4*2*i*j)*tmp5);
w4[0].imag()=sin((4*3*i*(j+4*4*1)-4*2*i*j)*tmp5);
tmp1=w1[0]*tab[0*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j];
tmp2=w2[0]*tab[1*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j];
tmp3=w3[0]*tab[2*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j];
tmp4=w4[0]*tab[3*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j];
//radix-4
tab2[0*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=tmp1+tmp2+ tmp3+tmp4;
tab2[1*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab2[2*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=tmp1-tmp2+ tmp3-tmp4;
tab2[3*N/(4*4*4)+N/(4*4)*i+N/(4)*m+j]=tmp1-tmp3+w6[0]*(tmp2-tmp4);
}
}
}
////////////////////////////////////////////////////////////////////////////////
for(int m=0;m<4*4;m++)//stage 4
{
for(int i=0;i<4;i++)
{
for(int j=0;j<4;j++)
{
w1[0].real()= cos((4*4*0*i*(j+4*0)+4*4*1*i*j)*tmp5);
w1[0].imag()=sin((4*4*0*i*(j+4*0)+4*4*1*i*j)*tmp5);
w2[0].real()= cos((4*4*1*i*(j+4*1)+4*4*0*i*j)*tmp5);
w2[0].imag()=sin((4*4*1*i*(j+4*1)+4*4*0*i*j)*tmp5);
w3[0].real()=cos((4*4*2*i*(j+4*1)-4*4*1*i*j)*tmp5);
w3[0].imag()=sin((4*4*2*i*(j+4*1)-4*4*1*i*j)*tmp5);
w4[0].real()= cos((4*4*3*i*(j+4*1)-4*4*2*i*j)*tmp5);
w4[0].imag()=sin((4*4*3*i*(j+4*1)-4*4*2*i*j)*tmp5);
tmp1=w1[0]*tab2[0*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j];
tmp2=w2[0]*tab2[1*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j];
tmp3=w3[0]*tab2[2*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j];
tmp4=w4[0]*tab2[3*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j];
//radix-4
tab[0*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=tmp1+tmp2+ tmp3+tmp4;
tab[1*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab[2*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=tmp1-tmp2+ tmp3-tmp4;
tab[3*N/(4*4*4*4)+N/(4*4*4)*i+N/(4*4)*m+j]=tmp1-tmp3+w6[0]*(tmp2-tmp4);
}
}
}
///////////////////////////////////////////////////////////////////////////////
for(int m=0;m<4*4*4;m++)//stage 5
{
for(int i=0;i<4;i++)
{
for(int j=0;j<1;j++)
{
w1[0].real()= cos((4*4*4*0*i*(j+1*0)+4*4*4*1*i*j)*tmp5);
w1[0].imag()=sin((4*4*4*0*i*(j+1*0)+4*4*4*1*i*j)*tmp5);
w2[0].real()= cos((4*4*4*1*i*(j+1*1)+4*4*4*0*i*j)*tmp5);
w2[0].imag()=sin((4*4*4*1*i*(j+1*1)+4*4*4*0*i*j)*tmp5);
w3[0].real()= cos((4*4*4*2*i*(j+1*1)-4*4*4*1*i*j)*tmp5);
w3[0].imag()=sin((4*4*4*2*i*(j+1*1)-4*4*4*1*i*j)*tmp5);
w4[0].real()= cos((4*4*4*3*i*(j+1*1)-4*4*4*2*i*j)*tmp5);
w4[0].imag()=sin((4*4*4*3*i*(j+1*1)-4*4*4*2*i*j)*tmp5);
tmp1=w1[0]*tab[0*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j];
tmp2=w2[0]*tab[1*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j];
tmp3=w3[0]*tab[2*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j];
tmp4=w4[0]*tab[3*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j];
//radix-4
tab2[0*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j]=tmp1+tmp2+ tmp3+tmp4;
tab2[1*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j]=tmp1-tmp3+w5[0]*(tmp2-tmp4);
tab2[2*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j]=tmp1-tmp2+ tmp3-tmp4;
tab2[3*N/(4*4*4*4*4)+N/(4*4*4*4)*i+N/(4*4*4)*m+j]=tmp1-tmp3+w6[0]*(tmp2-tmp4);
}
}
}
for(int j=0;j<N;j++)
{
tab[j].real() =tab2[j].real()*0.5;
tab[j].imag() =tab2[j].imag()*0.5;
}
}
//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)
Subskrybuj:
Komentarze do posta (Atom)
Brak komentarzy:
Prześlij komentarz