1921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol#ifndef KISSFFT_CLASS_HH 2921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol#include <complex> 3921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol#include <vector> 4921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 5921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignolnamespace kissfft_utils { 6921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 7921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignoltemplate <typename T_scalar> 8921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignolstruct traits 9921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol{ 10921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol typedef T_scalar scalar_type; 11921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol typedef std::complex<scalar_type> cpx_type; 12921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol void fill_twiddles( std::complex<T_scalar> * dst ,int nfft,bool inverse) 13921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol { 14921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol T_scalar phinc = (inverse?2:-2)* acos( (T_scalar) -1) / nfft; 15921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol for (int i=0;i<nfft;++i) 16921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol dst[i] = exp( std::complex<T_scalar>(0,i*phinc) ); 17921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 18921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 19921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol void prepare( 20921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol std::vector< std::complex<T_scalar> > & dst, 21921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol int nfft,bool inverse, 22921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol std::vector<int> & stageRadix, 23921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol std::vector<int> & stageRemainder ) 24921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol { 25921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol _twiddles.resize(nfft); 26921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol fill_twiddles( &_twiddles[0],nfft,inverse); 27921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol dst = _twiddles; 28921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 29921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol //factorize 30921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol //start factoring out 4's, then 2's, then 3,5,7,9,... 31921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol int n= nfft; 32921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol int p=4; 33921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol do { 34921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol while (n % p) { 35921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol switch (p) { 36921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol case 4: p = 2; break; 37921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol case 2: p = 3; break; 38921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol default: p += 2; break; 39921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 40921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol if (p*p>n) 41921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol p=n;// no more factors 42921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 43921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol n /= p; 44921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol stageRadix.push_back(p); 45921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol stageRemainder.push_back(n); 46921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol }while(n>1); 47921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 48921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol std::vector<cpx_type> _twiddles; 49921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 50921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 51921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol const cpx_type twiddle(int i) { return _twiddles[i]; } 52921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol}; 53921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 54921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol} 55921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 56921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignoltemplate <typename T_Scalar, 57921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol typename T_traits=kissfft_utils::traits<T_Scalar> 58921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol > 59921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignolclass kissfft 60921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol{ 61921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol public: 62921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol typedef T_traits traits_type; 63921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol typedef typename traits_type::scalar_type scalar_type; 64921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol typedef typename traits_type::cpx_type cpx_type; 65921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 66921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol kissfft(int nfft,bool inverse,const traits_type & traits=traits_type() ) 67921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol :_nfft(nfft),_inverse(inverse),_traits(traits) 68921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol { 69921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol _traits.prepare(_twiddles, _nfft,_inverse ,_stageRadix, _stageRemainder); 70921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 71921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 72921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol void transform(const cpx_type * src , cpx_type * dst) 73921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol { 74921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol kf_work(0, dst, src, 1,1); 75921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 76921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 77921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol private: 78921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol void kf_work( int stage,cpx_type * Fout, const cpx_type * f, size_t fstride,size_t in_stride) 79921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol { 80921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol int p = _stageRadix[stage]; 81921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol int m = _stageRemainder[stage]; 82921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol cpx_type * Fout_beg = Fout; 83921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol cpx_type * Fout_end = Fout + p*m; 84921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 85921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol if (m==1) { 86921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol do{ 87921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol *Fout = *f; 88921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol f += fstride*in_stride; 89921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol }while(++Fout != Fout_end ); 90921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol }else{ 91921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol do{ 92921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol // recursive call: 93921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol // DFT of size m*p performed by doing 94921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol // p instances of smaller DFTs of size m, 95921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol // each one takes a decimated version of the input 96921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol kf_work(stage+1, Fout , f, fstride*p,in_stride); 97921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol f += fstride*in_stride; 98921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol }while( (Fout += m) != Fout_end ); 99921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 100921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 101921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol Fout=Fout_beg; 102921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 103921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol // recombine the p smaller DFTs 104921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol switch (p) { 105921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol case 2: kf_bfly2(Fout,fstride,m); break; 106921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol case 3: kf_bfly3(Fout,fstride,m); break; 107921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol case 4: kf_bfly4(Fout,fstride,m); break; 108921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol case 5: kf_bfly5(Fout,fstride,m); break; 109921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol default: kf_bfly_generic(Fout,fstride,m,p); break; 110921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 111921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 112921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 113921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol // these were #define macros in the original kiss_fft 114921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol void C_ADD( cpx_type & c,const cpx_type & a,const cpx_type & b) { c=a+b;} 115921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol void C_MUL( cpx_type & c,const cpx_type & a,const cpx_type & b) { c=a*b;} 116921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol void C_SUB( cpx_type & c,const cpx_type & a,const cpx_type & b) { c=a-b;} 117921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol void C_ADDTO( cpx_type & c,const cpx_type & a) { c+=a;} 118921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol void C_FIXDIV( cpx_type & ,int ) {} // NO-OP for float types 119921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol scalar_type S_MUL( const scalar_type & a,const scalar_type & b) { return a*b;} 120921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol scalar_type HALF_OF( const scalar_type & a) { return a*.5;} 121921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol void C_MULBYSCALAR(cpx_type & c,const scalar_type & a) {c*=a;} 122921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 123921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol void kf_bfly2( cpx_type * Fout, const size_t fstride, int m) 124921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol { 125921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol for (int k=0;k<m;++k) { 126921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol cpx_type t = Fout[m+k] * _traits.twiddle(k*fstride); 127921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol Fout[m+k] = Fout[k] - t; 128921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol Fout[k] += t; 129921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 130921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 131921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 132921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol void kf_bfly4( cpx_type * Fout, const size_t fstride, const size_t m) 133921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol { 134921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol cpx_type scratch[7]; 135921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol int negative_if_inverse = _inverse * -2 +1; 136921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol for (size_t k=0;k<m;++k) { 137921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol scratch[0] = Fout[k+m] * _traits.twiddle(k*fstride); 138921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol scratch[1] = Fout[k+2*m] * _traits.twiddle(k*fstride*2); 139921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol scratch[2] = Fout[k+3*m] * _traits.twiddle(k*fstride*3); 140921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol scratch[5] = Fout[k] - scratch[1]; 141921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 142921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol Fout[k] += scratch[1]; 143921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol scratch[3] = scratch[0] + scratch[2]; 144921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol scratch[4] = scratch[0] - scratch[2]; 145921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol scratch[4] = cpx_type( scratch[4].imag()*negative_if_inverse , -scratch[4].real()* negative_if_inverse ); 146921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 147921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol Fout[k+2*m] = Fout[k] - scratch[3]; 148921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol Fout[k] += scratch[3]; 149921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol Fout[k+m] = scratch[5] + scratch[4]; 150921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol Fout[k+3*m] = scratch[5] - scratch[4]; 151921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 152921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 153921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 154921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol void kf_bfly3( cpx_type * Fout, const size_t fstride, const size_t m) 155921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol { 156921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol size_t k=m; 157921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol const size_t m2 = 2*m; 158921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol cpx_type *tw1,*tw2; 159921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol cpx_type scratch[5]; 160921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol cpx_type epi3; 161921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol epi3 = _twiddles[fstride*m]; 162921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 163921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol tw1=tw2=&_twiddles[0]; 164921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 165921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol do{ 166921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3); 167921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 168921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_MUL(scratch[1],Fout[m] , *tw1); 169921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_MUL(scratch[2],Fout[m2] , *tw2); 170921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 171921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_ADD(scratch[3],scratch[1],scratch[2]); 172921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_SUB(scratch[0],scratch[1],scratch[2]); 173921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol tw1 += fstride; 174921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol tw2 += fstride*2; 175921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 176921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol Fout[m] = cpx_type( Fout->real() - HALF_OF(scratch[3].real() ) , Fout->imag() - HALF_OF(scratch[3].imag() ) ); 177921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 178921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_MULBYSCALAR( scratch[0] , epi3.imag() ); 179921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 180921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_ADDTO(*Fout,scratch[3]); 181921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 182921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol Fout[m2] = cpx_type( Fout[m].real() + scratch[0].imag() , Fout[m].imag() - scratch[0].real() ); 183921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 184921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_ADDTO( Fout[m] , cpx_type( -scratch[0].imag(),scratch[0].real() ) ); 185921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol ++Fout; 186921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol }while(--k); 187921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 188921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 189921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol void kf_bfly5( cpx_type * Fout, const size_t fstride, const size_t m) 190921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol { 191921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol cpx_type *Fout0,*Fout1,*Fout2,*Fout3,*Fout4; 192921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol size_t u; 193921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol cpx_type scratch[13]; 194921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol cpx_type * twiddles = &_twiddles[0]; 195921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol cpx_type *tw; 196921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol cpx_type ya,yb; 197921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol ya = twiddles[fstride*m]; 198921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol yb = twiddles[fstride*2*m]; 199921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 200921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol Fout0=Fout; 201921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol Fout1=Fout0+m; 202921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol Fout2=Fout0+2*m; 203921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol Fout3=Fout0+3*m; 204921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol Fout4=Fout0+4*m; 205921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 206921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol tw=twiddles; 207921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol for ( u=0; u<m; ++u ) { 208921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_FIXDIV( *Fout0,5); C_FIXDIV( *Fout1,5); C_FIXDIV( *Fout2,5); C_FIXDIV( *Fout3,5); C_FIXDIV( *Fout4,5); 209921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol scratch[0] = *Fout0; 210921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 211921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_MUL(scratch[1] ,*Fout1, tw[u*fstride]); 212921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]); 213921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]); 214921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]); 215921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 216921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_ADD( scratch[7],scratch[1],scratch[4]); 217921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_SUB( scratch[10],scratch[1],scratch[4]); 218921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_ADD( scratch[8],scratch[2],scratch[3]); 219921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_SUB( scratch[9],scratch[2],scratch[3]); 220921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 221921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_ADDTO( *Fout0, scratch[7]); 222921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_ADDTO( *Fout0, scratch[8]); 223921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 224921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol scratch[5] = scratch[0] + cpx_type( 225921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol S_MUL(scratch[7].real(),ya.real() ) + S_MUL(scratch[8].real() ,yb.real() ), 226921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol S_MUL(scratch[7].imag(),ya.real()) + S_MUL(scratch[8].imag(),yb.real()) 227921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol ); 228921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 229921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol scratch[6] = cpx_type( 230921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol S_MUL(scratch[10].imag(),ya.imag()) + S_MUL(scratch[9].imag(),yb.imag()), 231921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol -S_MUL(scratch[10].real(),ya.imag()) - S_MUL(scratch[9].real(),yb.imag()) 232921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol ); 233921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 234921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_SUB(*Fout1,scratch[5],scratch[6]); 235921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_ADD(*Fout4,scratch[5],scratch[6]); 236921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 237921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol scratch[11] = scratch[0] + 238921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol cpx_type( 239921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol S_MUL(scratch[7].real(),yb.real()) + S_MUL(scratch[8].real(),ya.real()), 240921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol S_MUL(scratch[7].imag(),yb.real()) + S_MUL(scratch[8].imag(),ya.real()) 241921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol ); 242921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 243921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol scratch[12] = cpx_type( 244921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol -S_MUL(scratch[10].imag(),yb.imag()) + S_MUL(scratch[9].imag(),ya.imag()), 245921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol S_MUL(scratch[10].real(),yb.imag()) - S_MUL(scratch[9].real(),ya.imag()) 246921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol ); 247921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 248921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_ADD(*Fout2,scratch[11],scratch[12]); 249921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_SUB(*Fout3,scratch[11],scratch[12]); 250921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 251921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4; 252921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 253921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 254921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 255921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol /* perform the butterfly for one stage of a mixed radix FFT */ 256921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol void kf_bfly_generic( 257921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol cpx_type * Fout, 258921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol const size_t fstride, 259921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol int m, 260921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol int p 261921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol ) 262921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol { 263921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol int u,k,q1,q; 264921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol cpx_type * twiddles = &_twiddles[0]; 265921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol cpx_type t; 266921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol int Norig = _nfft; 267921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol cpx_type scratchbuf[p]; 268921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 269921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol for ( u=0; u<m; ++u ) { 270921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol k=u; 271921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol for ( q1=0 ; q1<p ; ++q1 ) { 272921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol scratchbuf[q1] = Fout[ k ]; 273921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_FIXDIV(scratchbuf[q1],p); 274921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol k += m; 275921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 276921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 277921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol k=u; 278921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol for ( q1=0 ; q1<p ; ++q1 ) { 279921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol int twidx=0; 280921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol Fout[ k ] = scratchbuf[0]; 281921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol for (q=1;q<p;++q ) { 282921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol twidx += fstride * k; 283921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol if (twidx>=Norig) twidx-=Norig; 284921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_MUL(t,scratchbuf[q] , twiddles[twidx] ); 285921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol C_ADDTO( Fout[ k ] ,t); 286921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 287921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol k += m; 288921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 289921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 290921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol } 291921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol 292921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol int _nfft; 293921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol bool _inverse; 294921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol std::vector<cpx_type> _twiddles; 295921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol std::vector<int> _stageRadix; 296921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol std::vector<int> _stageRemainder; 297921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol traits_type _traits; 298921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol}; 299921cb97acb79e5096d5441f35027e401468d6cf5Andrew Rossignol#endif 300