src/calf/fft.h

   1 /* Calf DSP Library
   2  * FFT class
   3  *
   4  * Copyright (C) 2007 Krzysztof Foltman
   5  *
   6  * This program is free software; you can redistribute it and/or
   7  * modify it under the terms of the GNU Lesser General Public
   8  * License as published by the Free Software Foundation; either
   9  * version 2 of the License, or (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General
  17  * Public License along with this program; if not, write to the
  18  * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  19  * Boston, MA 02111-1307, USA.
  20  */
  21
  22 #ifndef __CALF_FFT_H
  23 #define __CALF_FFT_H
  24
  25 #include <complex>
  26
  27 namespace dsp {
  28
  29 /// FFT routine copied from my old OneSignal library, modified to
  30 /// match Calf's style. It's not fast at all, just a straightforward
  31 /// implementation.
  32 template<class T, int O>
  33 class fft
  34 {
  35 public:
  36     typedef typename std::complex<T> complex;
  37 private:
  38     int scramble[1<<O];
  39     complex sines[1<<O];
  40 public:
  41     fft()
  42     {
  43         int N=1<<O;
  44         assert(N >= 4);
  45         for (int i=0; i<N; i++)
  46         {
  47             int v=0;
  48             for (int j=0; j<O; j++)
  49                 if (i&(1<<j))
  50                     v+=(N>>(j+1));
  51             scramble[i]=v;
  52         }
  53         int N90 = N >> 2;
  54         T divN = 2 * M_PI / N;
  55         // use symmetry
  56         for (int i=0; i<N90; i++)
  57         {
  58             T angle = divN * i;
  59             T c = cos(angle), s = sin(angle);
  60             sines[i + 3 * N90] = -(sines[i + N90] = complex(-s, c));
  61             sines[i + 2 * N90] = -(sines[i] = complex(c, s));
  62         }
  63     }
  64     void calculate(complex *input, complex *output, bool inverse) const
  65     {
  66         int N=1<<O;
  67         int N1=N-1;
  68         int i;
  69         // Scramble the input data
  70         if (inverse)
  71         {
  72             float mf=1.0/N;
  73             for (i=0; i<N; i++)
  74             {
  75                 complex &c=input[scramble[i]];
  76                 output[i]=mf*complex(c.imag(),c.real());
  77             }
  78         }
  79         else
  80             for (i=0; i<N; i++)
  81                 output[i]=input[scramble[i]];
  82
  83         // O butterfiles
  84         for (i=0; i<O; i++)
  85         {
  86             int PO=1<<i, PNO=1<<(O-i-1);
  87             int j,k;
  88             for (j=0; j<PNO; j++)
  89             {
  90                 int base=j<<(i+1);
  91                 for (k=0; k<PO; k++)
  92                 {
  93                     int B1=base+k;
  94                     int B2=base+k+(1<<i);
  95                     complex r1=output[B1];
  96                     complex r2=output[B2];
  97                     output[B1]=r1+r2*sines[(B1<<(O-i-1))&N1];
  98                     output[B2]=r1+r2*sines[(B2<<(O-i-1))&N1];
  99                 }
 100             }
 101         }
 102         if (inverse)
 103         {
 104             for (i=0; i<N; i++)
 105             {
 106                 const complex &c=output[i];
 107                 output[i]=complex(c.imag(),c.real());
 108             }
 109         }
 110     }
 111     template<class InType>
 112     void calculateN(InType *input, complex *output, bool inverse, int order) const
 113     {
 114         assert(order <= O);
 115         int N=1<<order;
 116         int rsh=O - order;
 117         int N1=(N-1) << rsh;
 118         int i;
 119         // Scramble the input data
 120         if (inverse)
 121         {
 122             float mf=1.0/N;
 123             for (i=0; i<N; i++)
 124             {
 125                 const complex &c=input[scramble[i] >> rsh];
 126                 output[i]=mf*complex(c.imag(),c.real());
 127             }
 128         }
 129         else
 130             for (i=0; i<N; i++)
 131                 output[i]=input[scramble[i] >> rsh];
 132
 133         // O butterfiles
 134         for (i=0; i<order; i++)
 135         {
 136             int PO=1<<i, PNO=1<<(order-i-1);
 137             int j,k;
 138             for (j=0; j<PNO; j++)
 139             {
 140                 int base=j<<(i+1);
 141                 for (k=0; k<PO; k++)
 142                 {
 143                     int B1=base+k;
 144                     int B2=base+k+(1<<i);
 145                     complex r1=output[B1];
 146                     complex r2=output[B2];
 147                     output[B1]=r1+r2*sines[(B1<<(O-i-1))&N1];
 148                     output[B2]=r1+r2*sines[(B2<<(O-i-1))&N1];
 149                 }
 150             }
 151         }
 152         if (inverse)
 153         {
 154             for (i=0; i<N; i++)
 155             {
 156                 const complex &c=output[i];
 157                 output[i]=complex(c.imag(),c.real());
 158             }
 159         }
 160     }
 161     void execute_r2r(int order, float *input, float *output, complex *tmp, bool inverse = false) const
 162     {
 163         calculateN<float>(input, tmp, inverse, order);
 164         size_t s = 1 << order;
 165         size_t s2 = 1 << (order - 1);
 166         output[0] = tmp[0].real();
 167         output[s2] = tmp[0].imag();
 168         for (size_t i = 1; i < s2; ++i)
 169         {
 170             output[i] = tmp[i].real();
 171             output[s - 1 - i] = tmp[i].imag();
 172         }
 173     }
 174 };
 175
 176 };
 177
 178 #endif