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     typedef typename std::complex<T> complex;
  36     int scramble[1<<O];
  37     complex sines[1<<O];
  38 public:
  39     fft()
  40     {
  41         int N=1<<O;
  42         assert(N >= 4);
  43         for (int i=0; i<N; i++)
  44         {
  45             int v=0;
  46             for (int j=0; j<O; j++)
  47                 if (i&(1<<j))
  48                     v+=(N>>(j+1));
  49             scramble[i]=v;
  50         }
  51         int N90 = N >> 2;
  52         T divN = 2 * M_PI / N;
  53         // use symmetry
  54         for (int i=0; i<N90; i++)
  55         {
  56             T angle = divN * i;
  57             T c = cos(angle), s = sin(angle);
  58             sines[i + 3 * N90] = -(sines[i + N90] = complex(-s, c));
  59             sines[i + 2 * N90] = -(sines[i] = complex(c, s));
  60         }
  61     }
  62     void calculate(complex *input, complex *output, bool inverse)
  63     {
  64         int N=1<<O;
  65         int N1=N-1;
  66         int i;
  67         // Scramble the input data
  68         if (inverse)
  69         {
  70             float mf=1.0/N;
  71             for (i=0; i<N; i++)
  72             {
  73                 complex &c=input[scramble[i]];
  74                 output[i]=mf*complex(c.imag(),c.real());
  75             }
  76         }
  77         else
  78             for (i=0; i<N; i++)
  79                 output[i]=input[scramble[i]];
  80
  81         // O butterfiles
  82         for (i=0; i<O; i++)
  83         {
  84             int PO=1<<i, PNO=1<<(O-i-1);
  85             int j,k;
  86             for (j=0; j<PNO; j++)
  87             {
  88                 int base=j<<(i+1);
  89                 for (k=0; k<PO; k++)
  90                 {
  91                     int B1=base+k;
  92                     int B2=base+k+(1<<i);
  93                     complex r1=output[B1];
  94                     complex r2=output[B2];
  95                     output[B1]=r1+r2*sines[(B1<<(O-i-1))&N1];
  96                     output[B2]=r1+r2*sines[(B2<<(O-i-1))&N1];
  97                 }
  98             }
  99         }
 100         if (inverse)
 101         {
 102             for (i=0; i<N; i++)
 103             {
 104                 const complex &c=output[i];
 105                 output[i]=complex(c.imag(),c.real());
 106             }
 107         }
 108     }
 109 };
 110
 111 };
 112
 113 #endif