fourier.cs

用c#编写的任意项傅立叶变换和2的n次项快速傅立叶变换的算法。

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;


namespace Mymath
{
    public static class Fourier
    {
        public enum TransformDirection
        {
            Forward = 1, Backward = -1
        };

        public static Complex[,] exp(int Len,int direction)
        {
            int i,j,n,t,L;
            Complex[,] Result=new Complex[Len,Len];
            Complex[] expn = new Complex[2*Len];
            double ang;
            L = 2 * Len;
            for (i = 0; i < L; i++)
            {
                ang = direction * Math.PI * i / Len;
                expn[i] =new Complex(Math.Cos(ang), Math.Sin(ang));
            }
            for (j = 0; j < Len; j++)
            {
                Result[0, j] = expn[0];
                Result[j, 0] = expn[0];
            }
            for (i = 1; i < Len; i++)
            {   
                n=2*i*i;
                t=2*i;
                while (n >= L)
                {
                    n -= L;
                }
                Result[i, i] = expn[n];
                for (j = i+1; j < Len; j++)
                {
                    n += t;
                    while (n >= L)
                    {
                        n -= L;
                    }
                    Result[i, j] = expn[n];
                    Result[j, i] = expn[n];
                }
            }
            return(Result);
        }

        public static Complex[] DFt(Complex[] date, int direction)
        {
            int n=date.GetLength(0);
            Complex[] Result=new Complex[n];
            Complex[,] matrix = exp(n, direction);
            for (int i = 0; i < n; i++)
            {
                Result[i] = Complex.Zero;
                for (int j = 0; j < n; j++)
                {
                    Result[i] += matrix[i, j] * date[j];
                }
            }
            return (Result);
        }

        public static Complex[,] DFt2(Complex[,] date, int direction)
        {
            int m = date.GetLength(0);
            int n = date.GetLength(1);
            Complex[,] expm = exp(m, direction);
            Complex[,] expn = exp(n, direction);
            Complex[,] MiddleResult =new Complex[m, n];
            Complex[,] Result = new Complex[m, n];
            for (int k = 0; k < m; k++)
            {
                for (int i = 0; i < n; i++)
                {
                    MiddleResult[k, i] = Complex.Zero;
                    for (int j = 0; j < n; j++)
                    {
                        MiddleResult[k, i] += expn[i, j] * date[k, j];
                    }
                }
            }
            for (int k = 0; k < n; k++)
            {
                for (int i = 0; i < m; i++)
                {
                    Result[i, k] = Complex.Zero;
                    for (int j = 0; j < m; j++)
                    {
                        Result[i, k] =Result[i,k]+MiddleResult[j, k] * expm[i, j];
                    }
                }
            }
            if (direction == 1)
            {
                double d =(double) m * n;
                for (int k = 0; k < m; k++)
                    for (int t = 0; t < n; t++)
                    {
                        Result[k, t].Re = Result[k, t].Re / d;
                        Result[k, t].Im = Result[k, t].Im / d;
                    }
            }
            return (Result);
        }
        public static int BitReverse(int src, int size)
        {
            int des = 0;
            for (int i = size - 1; i >= 0; i--)
            {
                des=((src & 0x1)<<i)|des;
                src = src >> 1;
            }
            return des;
        }
        public static Complex[] Reorder(Complex[] x,int n,int m)
        {
            Complex[] new_x =(Complex[]) x.Clone();
            for (int i = 0; i < n; i++)
                new_x[i] = x[BitReverse(i, m)];
            return new_x;
        }
        public static Complex[] Calcw(int n, int flag)
        {
            double ang;
            Complex[] w = new Complex[n];
            for (int i = 0; i < n / 2; i++)
            {
                ang = flag * 2 * Math.PI * i / n;
                w[i] = new Complex(Math.Cos(ang), Math.Sin(ang));
            }
            return w;
        }
        /// <summary>
        /// One dimensional Fast Fourier Transform
        /// </summary>
        /// <param name="data">Data to transform</param>
        /// <param name="direction">Transformation direction</param>
        /// <returns>The result array</returns>
        public static Complex[] FFT(this Complex[] data, TransformDirection direction)
        {
            int n = data.Length;
            int m = MathHelper.Log2(n);

            Complex[] result = (Complex[])data.Clone();

            // reorder data first
            ReorderData(result);

            // compute FFT
            int tn = 1, tm;

            for (int k = 1; k <= m; k++)
            {
                Complex[] rotation =GetComplexRotation(k, direction);

                tm = tn;
                tn <<= 1;

                for (int i = 0; i < tm; i++)
                {
                    Complex t = rotation[i];

                    for (int even = i; even < n; even += tn)
                    {
                        int odd = even + tm;
                        Complex ce = result[even];
                        Complex co = result[odd];

                        double tr = co.Re * t.Re - co.Im * t.Im;
                        double ti = co.Re * t.Im + co.Im * t.Re;

                        result[even].Re += tr;
                        result[even].Im += ti;

                        result[odd].Re = ce.Re - tr;
                        result[odd].Im = ce.Im - ti;
                    }
                }
            }

            if (direction == TransformDirection.Forward)
            {
                for (int i = 0; i < n; i++)
                {
                    result[i].Re /= (double)n;
                    result[i].Im /= (double)n;
                }
            }

            return result;
        }

        /// <summary>
        /// Two dimensional Fast Fourier Transform
        /// </summary>
        /// <param name="data">Data to transform</param>
        /// <param name="direction">Transformation direction</param>
        /// <returns>The result matrix</returns>
        public static Complex[,] FFt2(this Complex[,] data, TransformDirection direction)
        {
            Complex[,] result = (Complex[,])data.Clone();

            int k = result.GetLength(0);
            int n = result.GetLength(1);

            // check data size
            if (
                (!MathHelper.IsPowerOf2(k)) || (!MathHelper.IsPowerOf2(n)) ||
                (k < minLength) || (k > maxLength) ||
                (n < minLength) || (n > maxLength)
                )
            {
                throw new ArgumentException();
            }

            // process rows
            Complex[] row = new Complex[n];

            for (int i = 0; i < k; i++)
            {
                // copy row
                for (int j = 0; j < n; j++)
                    row[j] = result[i, j];
                // transform it
                row = FFT(row, direction);
                // copy back
                for (int j = 0; j < n; j++)
                    result[i, j] = row[j];
            }

            // process columns
            Complex[] col = new Complex[k];

            for (int j = 0; j < n; j++)
            {
                // copy column
                for (int i = 0; i < k; i++)
                    col[i] = result[i, j];
                // transform it
                col = FFT(col, direction);
                // copy back
                for (int i = 0; i < k; i++)
                    result[i, j] = col[i];
            }

            return result;
        }

        #region Private Region

        private const int minLength = 2;
        private const int maxLength = 16384;
        private const int minBits = 1;
        private const int maxBits = 14;
        private static int[][] reversedBits = new int[maxBits][];
        private static Complex[,][] complexRotation = new Complex[maxBits, 2][];

        // Get array, indicating which data members should be swapped before FFT
        private static int[] GetReversedBits(int numberOfBits)
        {
            if ((numberOfBits < minBits) || (numberOfBits > maxBits))
                throw new ArgumentOutOfRangeException();

            // check if the array is already calculated
            if (reversedBits[numberOfBits - 1] == null)
            {
                int n = MathHelper.Pow2(numberOfBits);
                int[] rBits = new int[n];

                // calculate the array
                for (int i = 0; i < n; i++)
                {
                    int oldBits = i;
                    int newBits = 0;

                    for (int j = 0; j < numberOfBits; j++)
                    {
                        newBits = (newBits << 1) | (oldBits & 1);
                        oldBits = (oldBits >> 1);
                    }
                    rBits[i] = newBits;
                }
                reversedBits[numberOfBits - 1] = rBits;
            }
            return reversedBits[numberOfBits - 1];
        }

        // Get rotation of complex number
        private static Complex[] GetComplexRotation(int numberOfBits, TransformDirection direction)
        {
            int directionIndex = (direction == TransformDirection.Forward) ? 0 : 1;

            // check if the array is already calculated
            if (complexRotation[numberOfBits - 1, directionIndex] == null)
            {
                int n = 1 << (numberOfBits - 1);
                double uR = 1.0;
                double uI = 0.0;
                double angle = System.Math.PI / n * (int)direction;
                double wR = System.Math.Cos(angle);
                double wI = System.Math.Sin(angle);
                double t;
                Complex[] rotation = new Complex[n];

                for (int i = 0; i < n; i++)
                {
                    rotation[i] = new Complex(uR, uI);
                    t = uR * wI + uI * wR;
                    uR = uR * wR - uI * wI;
                    uI = t;
                }

                complexRotation[numberOfBits - 1, directionIndex] = rotation;
            }
            return complexRotation[numberOfBits - 1, directionIndex];
        }

        // Reorder data for FFT using
        private static void ReorderData(Complex[] data)
        {
            int len = data.Length;

            // check data length
            if ((len < minLength) || (len > maxLength))
                throw new ArgumentException();

            int[] rBits = GetReversedBits(MathHelper.Log2(len));
            for (int i = 0; i < len; i++)
            {
                int s = rBits[i];

                if (s > i)
                {
                    Complex t = data[i];
                    data[i] = data[s];
                    data[s] = t;
                }
            }
        }

        #endregion
    }
}