blob: a586cb115413f7b85bfea69da5762e61a338272f [file] [log] [blame]
 #ifndef KISSFFT_CLASS_HH #include #include namespace kissfft_utils { template struct traits { typedef T_scalar scalar_type; typedef std::complex cpx_type; void fill_twiddles( std::complex * dst ,int nfft,bool inverse) { T_scalar phinc = (inverse?2:-2)* acos( (T_scalar) -1) / nfft; for (int i=0;i(0,i*phinc) ); } void prepare( std::vector< std::complex > & dst, int nfft,bool inverse, std::vector & stageRadix, std::vector & stageRemainder ) { _twiddles.resize(nfft); fill_twiddles( &_twiddles[0],nfft,inverse); dst = _twiddles; //factorize //start factoring out 4's, then 2's, then 3,5,7,9,... int n= nfft; int p=4; do { while (n % p) { switch (p) { case 4: p = 2; break; case 2: p = 3; break; default: p += 2; break; } if (p*p>n) p=n;// no more factors } n /= p; stageRadix.push_back(p); stageRemainder.push_back(n); }while(n>1); } std::vector _twiddles; const cpx_type twiddle(int i) { return _twiddles[i]; } }; } template > class kissfft { public: typedef T_traits traits_type; typedef typename traits_type::scalar_type scalar_type; typedef typename traits_type::cpx_type cpx_type; kissfft(int nfft,bool inverse,const traits_type & traits=traits_type() ) :_nfft(nfft),_inverse(inverse),_traits(traits) { _traits.prepare(_twiddles, _nfft,_inverse ,_stageRadix, _stageRemainder); } void transform(const cpx_type * src , cpx_type * dst) { kf_work(0, dst, src, 1,1); } private: void kf_work( int stage,cpx_type * Fout, const cpx_type * f, size_t fstride,size_t in_stride) { int p = _stageRadix[stage]; int m = _stageRemainder[stage]; cpx_type * Fout_beg = Fout; cpx_type * Fout_end = Fout + p*m; if (m==1) { do{ *Fout = *f; f += fstride*in_stride; }while(++Fout != Fout_end ); }else{ do{ // recursive call: // DFT of size m*p performed by doing // p instances of smaller DFTs of size m, // each one takes a decimated version of the input kf_work(stage+1, Fout , f, fstride*p,in_stride); f += fstride*in_stride; }while( (Fout += m) != Fout_end ); } Fout=Fout_beg; // recombine the p smaller DFTs switch (p) { case 2: kf_bfly2(Fout,fstride,m); break; case 3: kf_bfly3(Fout,fstride,m); break; case 4: kf_bfly4(Fout,fstride,m); break; case 5: kf_bfly5(Fout,fstride,m); break; default: kf_bfly_generic(Fout,fstride,m,p); break; } } // these were #define macros in the original kiss_fft void C_ADD( cpx_type & c,const cpx_type & a,const cpx_type & b) { c=a+b;} void C_MUL( cpx_type & c,const cpx_type & a,const cpx_type & b) { c=a*b;} void C_SUB( cpx_type & c,const cpx_type & a,const cpx_type & b) { c=a-b;} void C_ADDTO( cpx_type & c,const cpx_type & a) { c+=a;} void C_FIXDIV( cpx_type & ,int ) {} // NO-OP for float types scalar_type S_MUL( const scalar_type & a,const scalar_type & b) { return a*b;} scalar_type HALF_OF( const scalar_type & a) { return a*.5;} void C_MULBYSCALAR(cpx_type & c,const scalar_type & a) {c*=a;} void kf_bfly2( cpx_type * Fout, const size_t fstride, int m) { for (int k=0;kreal() - HALF_OF(scratch[3].real() ) , Fout->imag() - HALF_OF(scratch[3].imag() ) ); C_MULBYSCALAR( scratch[0] , epi3.imag() ); C_ADDTO(*Fout,scratch[3]); Fout[m2] = cpx_type( Fout[m].real() + scratch[0].imag() , Fout[m].imag() - scratch[0].real() ); C_ADDTO( Fout[m] , cpx_type( -scratch[0].imag(),scratch[0].real() ) ); ++Fout; }while(--k); } void kf_bfly5( cpx_type * Fout, const size_t fstride, const size_t m) { cpx_type *Fout0,*Fout1,*Fout2,*Fout3,*Fout4; size_t u; cpx_type scratch[13]; cpx_type * twiddles = &_twiddles[0]; cpx_type *tw; cpx_type ya,yb; ya = twiddles[fstride*m]; yb = twiddles[fstride*2*m]; Fout0=Fout; Fout1=Fout0+m; Fout2=Fout0+2*m; Fout3=Fout0+3*m; Fout4=Fout0+4*m; tw=twiddles; for ( u=0; u=Norig) twidx-=Norig; C_MUL(t,scratchbuf[q] , twiddles[twidx] ); C_ADDTO( Fout[ k ] ,t); } k += m; } } } int _nfft; bool _inverse; std::vector _twiddles; std::vector _stageRadix; std::vector _stageRemainder; traits_type _traits; }; #endif