"Fossies" - the Fresh Open Source Software Archive

Member "getdp-3.1.0-source/Numeric/kissfft.hh" (31 Jul 2018, 12203 Bytes) of package /linux/privat/getdp-3.1.0-source.tgz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "kissfft.hh" see the Fossies "Dox" file reference documentation.

    1 // KISS FFT 1.3.0
    2 //
    3 // Copyright (c) 2003-2010 Mark Borgerding
    4 //
    5 // All rights reserved.
    6 //
    7 // Redistribution and use in source and binary forms, with or without
    8 // modification, are permitted provided that the following conditions are met:
    9 //
   10 // * Redistributions of source code must retain the above copyright notice, this
   11 //   list of conditions and the following disclaimer.
   12 //
   13 // * Redistributions in binary form must reproduce the above copyright notice,
   14 //   this list of conditions and the following disclaimer in the documentation
   15 //   and/or other materials provided with the distribution.
   16 //
   17 // * Neither the author nor the names of any contributors may be used to endorse
   18 //   or promote products derived from this software without specific prior
   19 //   written permission.
   20 //
   21 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
   22 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
   23 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
   24 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
   25 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
   26 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
   27 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
   28 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
   29 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
   30 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
   31 // POSSIBILITY OF SUCH DAMAGE.
   32 
   33 #ifndef KISSFFT_CLASS_HH
   34 #include <complex>
   35 #include <vector>
   36 
   37 namespace kissfft_utils {
   38 
   39 template <typename T_scalar>
   40 struct traits
   41 {
   42     typedef T_scalar scalar_type;
   43     typedef std::complex<scalar_type> cpx_type;
   44     void fill_twiddles( std::complex<T_scalar> * dst ,int nfft,bool inverse)
   45     {
   46         T_scalar phinc =  (inverse?2:-2)* acos( (T_scalar) -1)  / nfft;
   47         for (int i=0;i<nfft;++i)
   48             dst[i] = exp( std::complex<T_scalar>(0,i*phinc) );
   49     }
   50 
   51     void prepare(
   52             std::vector< std::complex<T_scalar> > & dst,
   53             int nfft,bool inverse,
   54             std::vector<int> & stageRadix,
   55             std::vector<int> & stageRemainder )
   56     {
   57         _twiddles.resize(nfft);
   58         fill_twiddles( &_twiddles[0],nfft,inverse);
   59         dst = _twiddles;
   60 
   61         //factorize
   62         //start factoring out 4's, then 2's, then 3,5,7,9,...
   63         int n= nfft;
   64         int p=4;
   65         do {
   66             while (n % p) {
   67                 switch (p) {
   68                     case 4: p = 2; break;
   69                     case 2: p = 3; break;
   70                     default: p += 2; break;
   71                 }
   72                 if (p*p>n)
   73                     p=n;// no more factors
   74             }
   75             n /= p;
   76             stageRadix.push_back(p);
   77             stageRemainder.push_back(n);
   78         }while(n>1);
   79     }
   80     std::vector<cpx_type> _twiddles;
   81 
   82 
   83     const cpx_type twiddle(int i) { return _twiddles[i]; }
   84 };
   85 
   86 }
   87 
   88 template <typename T_Scalar,
   89          typename T_traits=kissfft_utils::traits<T_Scalar>
   90          >
   91 class kissfft
   92 {
   93     public:
   94         typedef T_traits traits_type;
   95         typedef typename traits_type::scalar_type scalar_type;
   96         typedef typename traits_type::cpx_type cpx_type;
   97 
   98         kissfft(int nfft,bool inverse,const traits_type & traits=traits_type() )
   99             :_nfft(nfft),_inverse(inverse),_traits(traits)
  100         {
  101             _traits.prepare(_twiddles, _nfft,_inverse ,_stageRadix, _stageRemainder);
  102         }
  103 
  104         void transform(const cpx_type * src , cpx_type * dst)
  105         {
  106             kf_work(0, dst, src, 1,1);
  107         }
  108 
  109     private:
  110         void kf_work( int stage,cpx_type * Fout, const cpx_type * f, size_t fstride,size_t in_stride)
  111         {
  112             int p = _stageRadix[stage];
  113             int m = _stageRemainder[stage];
  114             cpx_type * Fout_beg = Fout;
  115             cpx_type * Fout_end = Fout + p*m;
  116 
  117             if (m==1) {
  118                 do{
  119                     *Fout = *f;
  120                     f += fstride*in_stride;
  121                 }while(++Fout != Fout_end );
  122             }else{
  123                 do{
  124                     // recursive call:
  125                     // DFT of size m*p performed by doing
  126                     // p instances of smaller DFTs of size m,
  127                     // each one takes a decimated version of the input
  128                     kf_work(stage+1, Fout , f, fstride*p,in_stride);
  129                     f += fstride*in_stride;
  130                 }while( (Fout += m) != Fout_end );
  131             }
  132 
  133             Fout=Fout_beg;
  134 
  135             // recombine the p smaller DFTs
  136             switch (p) {
  137                 case 2: kf_bfly2(Fout,fstride,m); break;
  138                 case 3: kf_bfly3(Fout,fstride,m); break;
  139                 case 4: kf_bfly4(Fout,fstride,m); break;
  140                 case 5: kf_bfly5(Fout,fstride,m); break;
  141                 default: kf_bfly_generic(Fout,fstride,m,p); break;
  142             }
  143         }
  144 
  145         // these were #define macros in the original kiss_fft
  146         void C_ADD( cpx_type & c,const cpx_type & a,const cpx_type & b) { c=a+b;}
  147         void C_MUL( cpx_type & c,const cpx_type & a,const cpx_type & b) { c=a*b;}
  148         void C_SUB( cpx_type & c,const cpx_type & a,const cpx_type & b) { c=a-b;}
  149         void C_ADDTO( cpx_type & c,const cpx_type & a) { c+=a;}
  150         void C_FIXDIV( cpx_type & ,int ) {} // NO-OP for float types
  151         scalar_type S_MUL( const scalar_type & a,const scalar_type & b) { return a*b;}
  152         scalar_type HALF_OF( const scalar_type & a) { return a*.5;}
  153         void C_MULBYSCALAR(cpx_type & c,const scalar_type & a) {c*=a;}
  154 
  155         void kf_bfly2( cpx_type * Fout, const size_t fstride, int m)
  156         {
  157             for (int k=0;k<m;++k) {
  158                 cpx_type t = Fout[m+k] * _traits.twiddle(k*fstride);
  159                 Fout[m+k] = Fout[k] - t;
  160                 Fout[k] += t;
  161             }
  162         }
  163 
  164         void kf_bfly4( cpx_type * Fout, const size_t fstride, const size_t m)
  165         {
  166             cpx_type scratch[7];
  167             int negative_if_inverse = _inverse * -2 +1;
  168             for (size_t k=0;k<m;++k) {
  169                 scratch[0] = Fout[k+m] * _traits.twiddle(k*fstride);
  170                 scratch[1] = Fout[k+2*m] * _traits.twiddle(k*fstride*2);
  171                 scratch[2] = Fout[k+3*m] * _traits.twiddle(k*fstride*3);
  172                 scratch[5] = Fout[k] - scratch[1];
  173 
  174                 Fout[k] += scratch[1];
  175                 scratch[3] = scratch[0] + scratch[2];
  176                 scratch[4] = scratch[0] - scratch[2];
  177                 scratch[4] = cpx_type( scratch[4].imag()*negative_if_inverse , -scratch[4].real()* negative_if_inverse );
  178 
  179                 Fout[k+2*m]  = Fout[k] - scratch[3];
  180                 Fout[k] += scratch[3];
  181                 Fout[k+m] = scratch[5] + scratch[4];
  182                 Fout[k+3*m] = scratch[5] - scratch[4];
  183             }
  184         }
  185 
  186         void kf_bfly3( cpx_type * Fout, const size_t fstride, const size_t m)
  187         {
  188             size_t k=m;
  189             const size_t m2 = 2*m;
  190             cpx_type *tw1,*tw2;
  191             cpx_type scratch[5];
  192             cpx_type epi3;
  193             epi3 = _twiddles[fstride*m];
  194 
  195             tw1=tw2=&_twiddles[0];
  196 
  197             do{
  198                 C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3);
  199 
  200                 C_MUL(scratch[1],Fout[m] , *tw1);
  201                 C_MUL(scratch[2],Fout[m2] , *tw2);
  202 
  203                 C_ADD(scratch[3],scratch[1],scratch[2]);
  204                 C_SUB(scratch[0],scratch[1],scratch[2]);
  205                 tw1 += fstride;
  206                 tw2 += fstride*2;
  207 
  208                 Fout[m] = cpx_type( Fout->real() - HALF_OF(scratch[3].real() ) , Fout->imag() - HALF_OF(scratch[3].imag() ) );
  209 
  210                 C_MULBYSCALAR( scratch[0] , epi3.imag() );
  211 
  212                 C_ADDTO(*Fout,scratch[3]);
  213 
  214                 Fout[m2] = cpx_type(  Fout[m].real() + scratch[0].imag() , Fout[m].imag() - scratch[0].real() );
  215 
  216                 C_ADDTO( Fout[m] , cpx_type( -scratch[0].imag(),scratch[0].real() ) );
  217                 ++Fout;
  218             }while(--k);
  219         }
  220 
  221         void kf_bfly5( cpx_type * Fout, const size_t fstride, const size_t m)
  222         {
  223             cpx_type *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
  224             size_t u;
  225             cpx_type scratch[13];
  226             cpx_type * twiddles = &_twiddles[0];
  227             cpx_type *tw;
  228             cpx_type ya,yb;
  229             ya = twiddles[fstride*m];
  230             yb = twiddles[fstride*2*m];
  231 
  232             Fout0=Fout;
  233             Fout1=Fout0+m;
  234             Fout2=Fout0+2*m;
  235             Fout3=Fout0+3*m;
  236             Fout4=Fout0+4*m;
  237 
  238             tw=twiddles;
  239             for ( u=0; u<m; ++u ) {
  240                 C_FIXDIV( *Fout0,5); C_FIXDIV( *Fout1,5); C_FIXDIV( *Fout2,5); C_FIXDIV( *Fout3,5); C_FIXDIV( *Fout4,5);
  241                 scratch[0] = *Fout0;
  242 
  243                 C_MUL(scratch[1] ,*Fout1, tw[u*fstride]);
  244                 C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]);
  245                 C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]);
  246                 C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]);
  247 
  248                 C_ADD( scratch[7],scratch[1],scratch[4]);
  249                 C_SUB( scratch[10],scratch[1],scratch[4]);
  250                 C_ADD( scratch[8],scratch[2],scratch[3]);
  251                 C_SUB( scratch[9],scratch[2],scratch[3]);
  252 
  253                 C_ADDTO( *Fout0, scratch[7]);
  254                 C_ADDTO( *Fout0, scratch[8]);
  255 
  256                 scratch[5] = scratch[0] + cpx_type(
  257                         S_MUL(scratch[7].real(),ya.real() ) + S_MUL(scratch[8].real() ,yb.real() ),
  258                         S_MUL(scratch[7].imag(),ya.real()) + S_MUL(scratch[8].imag(),yb.real())
  259                         );
  260 
  261                 scratch[6] =  cpx_type(
  262                         S_MUL(scratch[10].imag(),ya.imag()) + S_MUL(scratch[9].imag(),yb.imag()),
  263                         -S_MUL(scratch[10].real(),ya.imag()) - S_MUL(scratch[9].real(),yb.imag())
  264                         );
  265 
  266                 C_SUB(*Fout1,scratch[5],scratch[6]);
  267                 C_ADD(*Fout4,scratch[5],scratch[6]);
  268 
  269                 scratch[11] = scratch[0] +
  270                     cpx_type(
  271                             S_MUL(scratch[7].real(),yb.real()) + S_MUL(scratch[8].real(),ya.real()),
  272                             S_MUL(scratch[7].imag(),yb.real()) + S_MUL(scratch[8].imag(),ya.real())
  273                             );
  274 
  275                 scratch[12] = cpx_type(
  276                         -S_MUL(scratch[10].imag(),yb.imag()) + S_MUL(scratch[9].imag(),ya.imag()),
  277                         S_MUL(scratch[10].real(),yb.imag()) - S_MUL(scratch[9].real(),ya.imag())
  278                         );
  279 
  280                 C_ADD(*Fout2,scratch[11],scratch[12]);
  281                 C_SUB(*Fout3,scratch[11],scratch[12]);
  282 
  283                 ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
  284             }
  285         }
  286 
  287         /* perform the butterfly for one stage of a mixed radix FFT */
  288         void kf_bfly_generic(
  289                 cpx_type * Fout,
  290                 const size_t fstride,
  291                 int m,
  292                 int p
  293                 )
  294         {
  295             int u,k,q1,q;
  296             cpx_type * twiddles = &_twiddles[0];
  297             cpx_type t;
  298             int Norig = _nfft;
  299             std::vector<cpx_type> scratchbuf(p);
  300 
  301             for ( u=0; u<m; ++u ) {
  302                 k=u;
  303                 for ( q1=0 ; q1<p ; ++q1 ) {
  304                     scratchbuf[q1] = Fout[ k  ];
  305                     C_FIXDIV(scratchbuf[q1],p);
  306                     k += m;
  307                 }
  308 
  309                 k=u;
  310                 for ( q1=0 ; q1<p ; ++q1 ) {
  311                     int twidx=0;
  312                     Fout[ k ] = scratchbuf[0];
  313                     for (q=1;q<p;++q ) {
  314                         twidx += fstride * k;
  315                         if (twidx>=Norig) twidx-=Norig;
  316                         C_MUL(t,scratchbuf[q] , twiddles[twidx] );
  317                         C_ADDTO( Fout[ k ] ,t);
  318                     }
  319                     k += m;
  320                 }
  321             }
  322         }
  323 
  324         int _nfft;
  325         bool _inverse;
  326         std::vector<cpx_type> _twiddles;
  327         std::vector<int> _stageRadix;
  328         std::vector<int> _stageRemainder;
  329         traits_type _traits;
  330 };
  331 #endif