common/alcomplex.c

   1
   2 #include "config.h"
   3
   4 #include "alcomplex.h"
   5 #include "math_defs.h"
   6
   7
   8 extern inline ALcomplex complex_add(ALcomplex a, ALcomplex b);
   9 extern inline ALcomplex complex_sub(ALcomplex a, ALcomplex b);
  10 extern inline ALcomplex complex_mult(ALcomplex a, ALcomplex b);
  11
  12 void complex_fft(ALcomplex *FFTBuffer, ALsizei FFTSize, ALdouble Sign)
  13 {
  14     ALsizei i, j, k, mask, step, step2;
  15     ALcomplex temp, u, w;
  16     ALdouble arg;
  17
  18     /* Bit-reversal permutation applied to a sequence of FFTSize items */
  19     for(i = 1;i < FFTSize-1;i++)
  20     {
  21         for(mask = 0x1, j = 0;mask < FFTSize;mask <<= 1)
  22         {
  23             if((i&mask) != 0)
  24                 j++;
  25             j <<= 1;
  26         }
  27         j >>= 1;
  28
  29         if(i < j)
  30         {
  31             temp         = FFTBuffer[i];
  32             FFTBuffer[i] = FFTBuffer[j];
  33             FFTBuffer[j] = temp;
  34         }
  35     }
  36
  37     /* Iterative form of DanielsonLanczos lemma */
  38     for(i = 1, step = 2;i < FFTSize;i<<=1, step<<=1)
  39     {
  40         step2 = step >> 1;
  41         arg   = M_PI / step2;
  42
  43         w.Real = cos(arg);
  44         w.Imag = sin(arg) * Sign;
  45
  46         u.Real = 1.0;
  47         u.Imag = 0.0;
  48
  49         for(j = 0;j < step2;j++)
  50         {
  51             for(k = j;k < FFTSize;k+=step)
  52             {
  53                 temp               = complex_mult(FFTBuffer[k+step2], u);
  54                 FFTBuffer[k+step2] = complex_sub(FFTBuffer[k], temp);
  55                 FFTBuffer[k]       = complex_add(FFTBuffer[k], temp);
  56             }
  57
  58             u = complex_mult(u, w);
  59         }
  60     }
  61 }
  62
  63 void complex_hilbert(ALcomplex *Buffer, ALsizei size)
  64 {
  65     const ALdouble inverse_size = 1.0/(ALdouble)size;
  66     ALsizei todo, i;
  67
  68     for(i = 0;i < size;i++)
  69         Buffer[i].Imag = 0.0;
  70
  71     complex_fft(Buffer, size, 1.0);
  72
  73     todo = size >> 1;
  74     Buffer[0].Real *= inverse_size;
  75     Buffer[0].Imag *= inverse_size;
  76     for(i = 1;i < todo;i++)
  77     {
  78         Buffer[i].Real *= 2.0*inverse_size;
  79         Buffer[i].Imag *= 2.0*inverse_size;
  80     }
  81     Buffer[i].Real *= inverse_size;
  82     Buffer[i].Imag *= inverse_size;
  83     i++;
  84
  85     for(;i < size;i++)
  86     {
  87         Buffer[i].Real = 0.0;
  88         Buffer[i].Imag = 0.0;
  89     }
  90
  91     complex_fft(Buffer, size, -1.0);
  92 }