Быстрое преобразование Фурье

/**

* алгоритм взят с http://alglib.sources.ru/fft/

*

*/

/*

*************************************************************************

Быстрое преобразование Фурье

Алгоритм проводит быстрое преобразование Фурье вещественной

функции, заданной n отсчетами на действительной оси.

В зависимости от  переданных параметров, может выполняться

как прямое, так и обратное преобразование.

Входные параметры:

   p   -   Число значений функции. Степень двойки.

   a   -   array [0 .. tnn-1] of Real

           Значения функции.

   InverseFFT

       -   направление преобразования.

           True, если обратное, False, если прямое.

Выходные параметры:

   a   -   результат преобразования.

Выходные параметры:

   a   -   результат преобразования.

 четные индексы   - cos-инусные составляющие, re

 нечетные индексы - sin-инусные составляющие, im

Определение частоты:

 F = SampleRate * i / BuffSize / 2

 где:

         SampleRate - частота дискретизации

         i - номер отсчета

         BuffSize - размер буфера

Нормировка результатов:

 a[i] = 2 * a[i] / n

Определение амплитуды:

 ampl = 2 * sqrt(a[i] * a[i] + a[i+1] * a[i+1]) / n

*************************************************************************

*/

 public static void realFastFourierTransform(double[] a, int p, boolean inversefft)

 {

         int tnn = 1 << p;

         if( tnn < 2 ) return;

         double twpr = Math.cos( 2.0 * Math.PI / tnn ) - 1.0;

         double twpi = Math.sin( 2.0 * Math.PI / tnn );

         int nn = tnn >> 1;

         if( inversefft ) {

                 twpr = -twpr;

                 twpi = -twpi;

                 double twr = 1.0 + twpr;

                 double twi = twpi;

                 for( int i = 2; i < nn + 1; i += 2 ) {

                         int i2 = i + 1;

                         int i3 = tnn - i;

                         int i4 = i3 + 1;

                         double h1r = 0.5 * (a[i] + a[i3] );

                         double h1i = 0.5 * (a[i2] - a[i4]);

                         double h2r = -0.5 * (a[i2] + a[i4]);

                         double h2i = 0.5 * (a[i] - a[i3]);

                         a[i]  =  h1r + twr * h2r - twi * h2i;

                         a[i2] =  h1i + twr * h2i + twi * h2r;

                         a[i3] =  h1r - twr * h2r + twi * h2i;

                         a[i4] = -h1i + twr * h2i + twi * h2r;

                         h1i = twr;

                         twr = twr * twpr - twi * twpi + twr;

                         twi = twi * twpr + h1i * twpi + twi;

                 }

                 twr = a[0];

                 a[0] = 0.5 * (twr + a[1]);

                 a[1] = 0.5 * (twr - a[1]);

         }

         for( int i = 1, j = 1; i < tnn; i += 2 ) {

                 if( j > i ) {

                         double temp = a[j - 1];

                         a[j - 1] = a[i - 1];

                         a[i - 1] = temp;

                         temp = a[j];

                         a[j] = a[i];

                         a[i] = temp;

                 }

                 int m = nn;

                 while( m >= 2 && j > m ) {

                         j = j - m;

                         m = m >> 1;

                 }

                 j = j + m;

         }

         double isign = 2 * Math.PI;

         if( inversefft ) isign = -isign;

         int mmax = 2;

         while( tnn > mmax ) {

                 int istep = mmax << 1;

                 double theta = isign / mmax;

                 double wpr = Math.cos( theta ) - 1.0;

                 double wpi = Math.sin( theta );

                 double wr = 1.0;

                 double wi = 0.0;

                 for( int ii = 1; ii < mmax; ii += 2 ) {

                         for( int i = ii; i <= tnn; i += istep ) {

                                 int j = i + mmax;

                                 double tempr = wr * a[j - 1] - wi * a[j];

                                 double tempi = wr * a[j] + wi * a[j - 1];

                                 a[j - 1] = a[i - 1] - tempr;

                                 a[j] = a[i] - tempi;

                                 a[i - 1] = a[i - 1] + tempr;

                                 a[i] = a[i] + tempi;

                         }

                         double wtemp = wr;

                         wr = wr * wpr - wi * wpi + wr;

                         wi = wi * wpr + wtemp * wpi + wi;

                 }

                 mmax = istep;

         }

         if( inversefft ) {

                 for( int i = 0; i < (nn << 1); i++ ) a[i] /= nn;

         }

         else {

                 double twr = 1.0 + twpr;

                 double twi = twpi;

                 for( int i = 2; i <= nn; i += 2 ) {

                         int i2 = i + 1;

                         int i3 = tnn - i;

                         int i4 = i3 + 1;

                         double h1r =  0.5 * (a[i] + a[i3]);

                         double h1i =  0.5 * (a[i2] - a[i4]);

                         double h2r =  0.5 * (a[i2] + a[i4]);

                         double h2i = -0.5 * (a[i] - a[i3]);

                         a[i]  =  h1r + twr * h2r - twi * h2i;

                         a[i2] =  h1i + twr * h2i + twi * h2r;

                         a[i3] =  h1r - twr * h2r + twi * h2i;

                         a[i4] = -h1i + twr * h2i + twi * h2r;

                         double twtemp = twr;

                         twr = twr * twpr - twi * twpi + twr;

                         twi = twi * twpr + twtemp * twpi + twi;

                 }

                 twr = a[0];

                 a[0] = twr + a[1];

                 a[1] = twr - a[1];

         }

 }