FFT Chris Lomont

private const int A = 0;
        private const int B = 1;
        public void FFT(double[] data, bool forward)
        {
            var n = data.Length;
            if ((n & (n - 1)) != 0)
                throw new ArgumentException(
                    "data length " + n + " in FFT is not a power of 2");
            n /= 2;    // кол-во сэмплов

            Reverse(data, n); //перестановка элементов согласнно мод.БПФ

            // преобразование
            double sign = forward ? B : -B;
            var mmax = 1;
            while (n > mmax)
            {
                var istep = 2 * mmax;
                var theta = sign * Math.PI / mmax;
                double wr = 1, wi = 0;
                var wpr = Math.Cos(theta);
                var wpi = Math.Sin(theta);
                for (var m = 0; m < istep; m += 2)
                {
                    for (var k = m; k < 2 * n; k += 2 * istep)
                    {
                        var j = k + istep;
                        var tempr = wr * data[j] - wi * data[j + 1];
                        var tempi = wi * data[j] + wr * data[j + 1];
                        data[j] = data[k] - tempr;
                        data[j + 1] = data[k + 1] - tempi;
                        data[k] = data[k] + tempr;
                        data[k + 1] = data[k + 1] + tempi;
                    }
                    var t = wr;
                    wr = wr * wpr - wi * wpi;
                    wi = wi * wpr + t * wpi;
                }
                mmax = istep;
            }

            // маштабирование данных
            Scale(data, n, forward);
        }

        void Scale(double[] data, int n, bool forward)
        {
            // при прямом
            if ((forward) && (A != 1))
            {
                var scale = Math.Pow(n, (A - 1) / 2.0);
                for (var i = 0; i < data.Length; ++i)
                    data[i] *= scale;
            }
            // при обратном
            if ((!forward) && (A != -1))
            {
                var scale = Math.Pow(n, -(A + 1) / 2.0);
                for (var i = 0; i < data.Length; ++i)
                    data[i] *= scale;
            }
        }

        static void Reverse(double[] data, int n)
        {
            // битовая перестановка индексов. Взята из 5 упр. в секции
            // 7.2.1.1 в ИП Кнута
            int j = 0, k = 0; // Кнут R1: инициализация
            var top = n / 2;  // это 2^(n-1) Кнута
            // цикл перестановки битов
            while (true)
            {
                // Кнут R2: swap - swap j+1 и k+2^(n-1), по 2 записи каждая
                var t = data[j + 2];
                data[j + 2] = data[k + n];
                data[k + n] = t;
                t = data[j + 3];
                data[j + 3] = data[k + n + 1];
                data[k + n + 1] = t;
                if (j > k)
                { // поменять местами еще две
                    // j и k
                    t = data[j];
                    data[j] = data[k];
                    data[k] = t;
                    t = data[j + 1];
                    data[j + 1] = data[k + 1];
                    data[k + 1] = t;
                    // j + top + 1 иk+top + 1
                    t = data[j + n + 2];
                    data[j + n + 2] = data[k + n + 2];
                    data[k + n + 2] = t;
                    t = data[j + n + 3];
                    data[j + n + 3] = data[k + n + 3];
                    data[k + n + 3] = t;
                }
                // Кнут R3: увеличить k
                k += 4;
                if (k >= n)
                    break;
                // Кнут R4: увеличить j
                var h = top;
                while (j >= h)
                {
                    j -= h;
                    h /= 2;
                }
                j += h;
            }