C#中实现FFT的两种方法

最近工作中有个需求，在C#环境中实现FFT算法，在网上找了些资料，最后实现了下面的两种方式，实际应用任选其一就好。

第一种方法：

不依赖C#中的Complex，需要实现计算过程的每一步详细步骤。

输入序列长度为2的N次幂，使用前需先定义序列长度：

FFT filter = new FFT(256);

filter.fft(x,y) 其中x为实部y为虚部，计算后x为FFT后的实部，y为FFT后的虚部。

 class FFT
    {
        int n, m;
        // Lookup tables. Only need to recompute when size of FFT changes.
        double[] cos;
        double[] sin;

        public FFT(int n)
        {
            this.n = n;
            this.m = (int)(Math.Log(n) / Math.Log(2));

            // Make sure n is a power of 2
            if (n != (1 << m))
                Console.Out.WriteLine("FFT length must be power of 2");

            // precompute tables
            cos = new double[n / 2];
            sin = new double[n / 2];

            for (int i = 0; i < n / 2; i++)
            {
                cos[i] = Math.Cos(-2 * Math.PI * i / n);
                sin[i] = Math.Sin(-2 * Math.PI * i / n);
            }

        }
///////x为实部y为虚部///////////
        private void fft(double[] x, double[] y)
        {
            int i, j, k, n1, n2, a;
            double c, s, t1, t2;

            // Bit-reverse
            j = 0;
            n2 = n / 2;
            for (i = 1; i < n - 1; i++)
            {
                n1 = n2;
                while (j >= n1)
                {
                    j = j - n1;
                    n1 = n1 / 2;
                }
                j = j + n1;

                if (i < j)
                {
                    t1 = x[i];
                    x[i] = x[j];
                    x[j] = t1;
                    t1 = y[i];
                    y[i] = y[j];
                    y[j] = t1;
                }
            }

            // FFT
            n1 = 0;
            n2 = 1;

            for (i = 0; i < m; i++)
            {
                n1 = n2;
                n2 = n2 + n2;
                a = 0;

                for (j = 0; j < n1; j++)
                {
                    c = cos[a];
                    s = sin[a];
                    a += 1 << (m - i - 1);

                    for (k = j; k < n; k = k + n2)
                    {
                        t1 = c * x[k + n1] - s * y[k + n1];
                        t2 = s * x[k + n1] + c * y[k + n1];
                        x[k + n1] = x[k] - t1;
                        y[k + n1] = y[k] - t2;
                        x[k] = x[k] + t1;
                        y[k] = y[k] + t2;
                    }
                }
            }
        }


    }

第二种方式：

傅里叶变换的计算中需要用到复数，无意中发现C#中有复数的定义，考虑使用复数直接计算的方式实现，这样可以利用C#中的复数计算。C#中使用复数需要引用using System.Numerics;

FFT filter = new FFT(256);

filter.DoFFT(x) ; x为输入参数(Complex[]类型)，计算后x为FFT后的结果。

突然发现C#好好用啊。。。。。。

class FFT
    {
       
        private Complex[] W;//定义旋转因子

        /// <summary>
        /// 构造函数用于初始化旋转因子。
        /// </summary>
        public FFT(int N)
        {                
            W = initW(N);
        }


        /// <summary>
        /// 快速傅立叶正变换函数
        /// W：输入旋转因子序列数组名
        /// x: 输入信号序列数组名
        /// </summary>
        public void DoFFT(Complex[] x)
        {
            int i = 0, j = 0, k = 0, l = 0;
            Complex up, down, product;
            ReverseOder(x);         ///对输入序列进行倒序
            for (i = 0; i < Math.Log(N, 2); i++)            ///外循环，变换级数循环
            {
                l = 1 << i;
                for (j = 0; j < N; j += 2 * l)             ///中间循环，旋转因子循环
                {
                    for (k = 0; k < l; k++)                     ///内循环，序列循环
                    {
                        product = x[j + k + l] * W[N * k / 2 / l];
                        up = x[j + k] + product;
                        down = x[j + k] - product;
                        x[j + k] = up;
                        x[j + k + l] = down;
                    }
                }
            }
        }

       

        /// <summary>
        /// 倒序函数
        /// W：输入旋转因子序列长度
        /// </summary>
        private Complex[] initW(int N)
        {
            Complex[] W = new Complex[N];
            for (int i = 0; i < N; i++)
            {
                double real = Math.Cos(2 * Math.PI * i / N);        //旋转因子实部展开
                double imag = -1 * Math.Sin(2 * Math.PI * i / N);   //旋转因子虚部展开
                W[i] = new Complex(real, imag);
            }
            return W;
        }


        /// <summary>
        /// ReverseOder()倒序变换
        /// </summary>
        /// <param name="x"></param>
        private void ReverseOder(Complex[] x)
        {
            Complex temp;
            int i = 0, j = 0, k = 0;
            double t;
            for (i = 0; i < N; i++)
            {
                k = i; j = 0;
                t = Math.Log(N, 2);
                while ((t--) > 0)
                {
                    j = j << 1;
                    j |= (k & 1);
                    k = k >> 1;
                }
                if (j > i)
                {
                    temp = x[i];
                    x[i] = x[j];
                    x[j] = temp;
                }
            }
        }


    }