FFT與二維FFT的C#實現

　　傅里葉變換在信息處理應用中具有很實用的價值，而快速傅里葉變換，即FFT，是實用的計算算法。

　　本文介紹FFT和2維FFT的C#實現方法。

1. FFT編程依據

　　FFT是按照如圖結構（也稱蝶形結構）進行運算（圖片來源於網絡）。

　　圖中，箭頭代表數據流向，箭頭與箭頭的合並點代表相加，箭頭下面的常數代表相乘，WN(注：N為下標) = exp(-j * 2 * PI / N)，j為虛數單位，DIT的輸入（或DIF的輸出）序列重排點y[m]與初始輸入點x[n]的關系為：m與n互為倒位序（例：001 -- 100，011 -- 110）。

　　該圖給出的是8點FFT例程，其余2^N點FFT可按照這個圖進行類推，類推后可以發現，2^N點FFT只需用三個循環套嵌即可實現。

　　對於idft的算法，根據公式X[k] = 對n求和{x[n] * (WN) ^ kn}以及x[n] = 對k求和{X[k] * (WN) ^ (-kn)} /N可知，只需將FFT算法中的WN換為(1 / WN)並將最后結果除以N即可得到idft。

2. FFT算法實現

　　由於作者並不擅長面向對象編程，僅僅是在功能上實現了本文需求，代碼在閱讀、效率上有很多缺陷。

　　下面是FFT算法的實現，完整的實現在文章后面。

        // 快速傅里葉變換
        // 作者：死貓
        // ComplexList src：數據輸入
        // ref int lenght：欲變換數據的長度，函數調用后此值更改為實際變換長度
        // int flag：區分fft或dtft，為1為fft，為-1為idft
        private static ComplexList fft_core(ComplexList src, ref int lenght, int flag)
        {
            //按時間抽取FFT方法(DIT)

            //補零后長度
            int relog2N = ReLog2N(lenght);

            int bitlenghth = relog2N;
            int N = 0x01 << bitlenghth;

            //重新復制數據，同時進行
            //    逆序排放，並補零
            int index;
            ComplexList Register = new ComplexList(N);
            for (int i = 0; i < N; i++)
            {
                index = ReArrange(i, bitlenghth);
                Register[i] = (index < src.Lenght) ? src[index] : new Complex(0);
            }

            //生成WN表，以免運行時進行重復計算
            ComplexList WN = new ComplexList(N / 2);
            for (int i = 0; i < N / 2; i++)
            {
                WN[i] = new Complex(Math.Cos(2 * Math.PI / N * i), -flag * Math.Sin(2 * Math.PI / N * i));
            }

            //蝶形運算
            int Index0, Index1;
            Complex temp;
            for (int steplenght = 2; steplenght <= N; steplenght *= 2)
            {
                for (int step = 0; step < N / steplenght; step++)
                {
                    for (int i = 0; i < steplenght / 2; i++)
                    {
                        Index0 = steplenght * step + i;
                        Index1 = steplenght * step + i + steplenght / 2;

                        temp = Register[Index1] * WN[N / steplenght * i];
                        Register[Index1] = Register[Index0] - temp;
                        Register[Index0] = Register[Index0] + temp;
                    }
                }
            }

            //若為idft
            if (-1 == flag)
            {
                for (int i = 0; i < Register.Lenght; i++)
                {
                    Register[i] /= N;
                }
            }

            //賦值
            lenght = N;

            /*
            //清理內存
            WN = null;
            temp = null;
            GC.Collect();
            */

            //返回
            return Register;
        }

 

3. 二維FFT變換

　　2維傅里葉變換是在1維傅里葉變換基礎上實現的，實現方法為：

　　1. 對2維輸入數據的每一行進行FFT變換，變換結果仍然按行存入一個二維數組中；

　　2. 對按行FFT變換后的結果，對每一列進行FFT變換，變換結果仍然按列存入一個二維數組中，該數組即為2維FFT變換結果。

　　代碼實現如下：

        // 2維快速傅里葉變換
        // 作者：死貓
        // int width, int height：欲變換數據的長度和寬度，函數調用后此值更改為實際變換長度
        // int flag：區分fft或dtft，為1為fft，為-1為idft
        private static Complex2D fft_2D_core(Complex2D src, ref int width, ref int height, int flag)
        {
            //補零后長度
            int width_Log2N = ReLog2N(width);
            int height_Log2N = ReLog2N(height);
            int Relog2N = Math.Max(width_Log2N, height_Log2N);
            int ReWidth = 0x01 << Relog2N;
            int ReHeight = 0x01 << Relog2N;

            //重新復制數據，清零
            Complex2D ReList2D = new Complex2D(ReWidth, ReHeight);
            ReList2D.Clear();
            int width_temp = Math.Min(src.Width, width);
            int height_temp = Math.Min(src.Height, height);
            for (int i = 0; i < height_temp; i++)
            {
                for (int j = 0; j < width_temp; j++)
                {
                    ReList2D[i, j] = new Complex(src[i, j].Re, src[i, j].Im);
                }
            }

            ComplexList Xn;
            ComplexList Xk;
            int Lenght_temp;

            //第1遍fft
            for (int i = 0; i < ReHeight; i++)
            {
                Xn = ReList2D.GetRow(i);
                Lenght_temp = Xn.Lenght;
                Xk = fft_core(Xn, ref Lenght_temp, flag);
                ReList2D.SetRow(i, Xk);
            }

            //第2遍fft
            for (int i = 0; i < ReWidth; i++)
            {
                Xn = ReList2D.GetColumn(i);
                Lenght_temp = Xn.Lenght;
                Xk = fft_core(Xn, ref Lenght_temp, flag);
                ReList2D.SetColumn(i, Xk);
            }

            //賦值
            width = ReWidth;
            height = ReHeight;

            //清理內存
            Xn = null;
            Xk = null;
            GC.Collect();

            //返回
            return ReList2D;
        }

3. 二維FFT測試結果

　　我對一張圖片進行低通、高通濾波，其傳輸函數分別為：

　　低通：Hlp(z) = (1 - a) / 2 * (1 + 1 / z) / (1 - a * 1 / z)，a = 0.8；

　　高通：Hhp(z) = (1 + a) / 2 * (1 - 1 / z) / (1 - a * 1 / z)，a = 0.7。

　　低通、高通濾波測試結果分別如下：

4. 其余相關代碼實現

　　由於作者並不擅長面向對象編程，僅僅是在功能上實現了本文需求，代碼在閱讀、效率上有很多缺陷。

　　代碼1：

        // 快速傅里葉變換
        // ComplexList src：數據輸入，若長度非2的n次冪，
        //     則函數對末尾補零后的數據進行運算
        public static ComplexList fft(ComplexList src)
        {
            int lenght = src.Lenght;
            return fft_core(src, ref lenght, 1);
        }
        // 快速傅里葉變換
        // ComplexList src：數據輸入
        // ref int lenght：欲變換的長度，函數調用后此值更改為實際變換長度
        public static ComplexList fft(ComplexList src, ref int lenght)
        {
            return fft_core(src, ref lenght, 1);
        }

        // 快速傅里葉變換的逆變換
        // ComplexList src：數據輸入
        public static ComplexList idft(ComplexList src)
        {
            int lenght = src.Lenght;
            return idft(src, ref lenght);
        }

        // 快速傅里葉變換的逆變換
        // ComplexList src：數據輸入
        // ref int lenght：欲變換數據的長度，函數調用后此值更改為實際變換長度
        public static ComplexList idft(ComplexList src, ref int lenght)
        {
            return fft_core(src, ref lenght, -1);
        }

        public static Complex2D fft_2D(Complex2D src)
        {
            int width = src.Width;
            int height = src.Height;
            return fft_2D_core(src, ref width, ref height, 1);
        }

        public static Complex2D idft_2D(Complex2D src)
        {
            int width = src.Width;
            int height = src.Height;
            return fft_2D_core(src, ref width, ref height, -1);
        }

　　代碼2：

        // 獲取按位逆序，bitlenght為數據長度
        // fft函數內使用
        private static int ReArrange(int dat, int bitlenght)
        {
            int ret = 0;
            for (int i = 0; i < bitlenght; i++)
            {
                if (0 != (dat & (0x01 << i))) ret |= ((0x01 << (bitlenght - 1)) >> i);
            }
            return ret;
        }

        // 獲取擴展長度后的冪次
        // 由於fft要求長度為2^n，所以用此函數來獲取所需長度
        public static int ReLog2N(int count)
        {
            int log2N = 0;
            uint mask = 0x80000000;
            for (int i = 0; i < 32; i++)
            {
                if (0 != ((mask >> i) & count))
                {
                    if ((mask >> i) == count) log2N = 31 - i;
                    else log2N = 31 - i + 1;
                    break;
                }
            }
            return log2N;
        }

　　代碼3：

    public class Complex2D
    {
        double[] _Complex2D_Re;
        double[] _Complex2D_Im;
        public int Width { get; private set; }
        public int Height { get; private set; }
        public Complex this[int Row, int Column]
        {
            get 
            {
                return new Complex(_Complex2D_Re[Row * Width + Column], _Complex2D_Im[Row * Width + Column]); 
            }
            set 
            {
                _Complex2D_Re[Row * Width + Column] = ((Complex)value).Re;
                _Complex2D_Im[Row * Width + Column] = ((Complex)value).Im; 
            }
        }

        public Complex2D(int width, int height)
        {
            Width = width;
            Height = height;
            int lenght = Width * Height;
            _Complex2D_Re = new double[lenght];
            _Complex2D_Im = new double[lenght];
        }

        public void Clear()
        {
            Array.Clear(_Complex2D_Re, 0, _Complex2D_Re.Length);
            Array.Clear(_Complex2D_Im, 0, _Complex2D_Im.Length);
        }

        public ComplexList GetColumn(int index)
        {
            ComplexList ret = new ComplexList(Height);
            for(int i = 0; i < Height; i++)
            {
                ret[i] = this[i, index];
            }
            return ret;
        }
        public int SetColumn(int index, ComplexList src)
        {
            for (int i = 0; i < Height; i++)
            {
                this[i, index] = (i < src.Lenght) ? src[i] : new Complex(0);
            }
            return 0;
        }
        public ComplexList GetRow(int index)
        {
            ComplexList ret = new ComplexList(Width);
            for(int i = 0; i < Width; i++)
            {
                ret[i] = this[index, i];
            }
            return ret;
        }
        public int SetRow(int index, ComplexList src)
        {
            for (int i = 0; i < Width; i++)
            {
                this[index, i] = (i < src.Lenght) ? src[i] : new Complex(0);
            }
            return 0;
        }
        public Complex2D GetAmplitude()
        {
            Complex2D ret = new Complex2D(Width, Height);
            for (int i = 0; i < Height; i++)
            {
                for (int j = 0; j < Width; j++)
                {
                    ret[i, j] = new Complex(this[i, j].Modulus());
                }
            }
            return ret;
        }
    }
    public class ComplexList
    {
        double[] _ComplexList_Re;
        double[] _ComplexList_Im;
        public int Lenght { get; private set; }
        public Complex this[int Index]
        {
            get
            {
                return new Complex(_ComplexList_Re[Index], _ComplexList_Im[Index]);
            }
            set
            {
                _ComplexList_Re[Index] = ((Complex)value).Re;
                _ComplexList_Im[Index] = ((Complex)value).Im;
            }
        }

        public ComplexList(int lenght)
        {
            Lenght = lenght;
            _ComplexList_Re = new double[Lenght];
            _ComplexList_Im = new double[Lenght];
        }
        public ComplexList(double[] re)
        {
            Lenght = re.Length;
            _ComplexList_Re = re;
            _ComplexList_Im = new double[Lenght];
        }
        public ComplexList(double[] re, double[] im)
        {
            Lenght = Math.Max(re.Length, im.Length);
            if (re.Length == im.Length)
            {
                _ComplexList_Re = re;
                _ComplexList_Im = im;
            }
            else
            {
                _ComplexList_Re = new double[Lenght];
                _ComplexList_Im = new double[Lenght];
                for (int i = 0; i < re.Length; i++) _ComplexList_Re[i] = re[i];
                for (int i = 0; i < im.Length; i++) _ComplexList_Im[i] = im[i];
            }
        }

        public void Clear()
        {
            Array.Clear(_ComplexList_Re, 0, _ComplexList_Re.Length);
            Array.Clear(_ComplexList_Im, 0, _ComplexList_Im.Length);
        }

        public double[] GetRePtr()
        {
            return _ComplexList_Re;
        }
        public double[] GetImPtr()
        {
            return _ComplexList_Im;
        }
        public ComplexList Clone()
        {
            return new ComplexList((double[])(_ComplexList_Re.Clone()), (double[])(_ComplexList_Im.Clone()));
        }
        public ComplexList GetAmplitude()
        {
            double[] amp = new double[Lenght];
            for (int i = 0; i < Lenght; i++)
            {
                amp[i] = this[i].Modulus();
            }
            return new ComplexList(amp);
        }
    }
    public class Complex
    {
        public double Re;
        public double Im;
        public Complex()
        {
            Re = 0;
            Im = 0;
        }
        public Complex(double re)
        {
            Re = re;
            Im = 0;
        }
        public Complex(double re, double im)
        {
            Re = re;
            Im = im;
        }

        public double Modulus()
        {
            return Math.Sqrt(Re * Re + Im * Im);
        }

        public override string ToString()
        {
            string retStr;
            if (Math.Abs(Im) < 0.0001) retStr = Re.ToString("f4");
            else if (Math.Abs(Re) < 0.0001)
            {
                if (Im > 0) retStr = "j" + Im.ToString("f4");
                else retStr = "- j" + (0 - Im).ToString("f4");
            }
            else
            {
                if (Im > 0) retStr = Re.ToString("f4") + "+ j" + Im.ToString("f4");
                else retStr = Re.ToString("f4") + "- j" + (0 - Im).ToString("f4");
            }
            retStr += " ";
            return retStr;
        }

        //操作符重載
        public static Complex operator +(Complex c1, Complex c2)
        {
            return new Complex(c1.Re + c2.Re, c1.Im + c2.Im);
        }
        public static Complex operator +(double d, Complex c)
        {
            return new Complex(d + c.Re, c.Im);
        }
        public static Complex operator -(Complex c1, Complex c2)
        {
            return new Complex(c1.Re - c2.Re, c1.Im - c2.Im);
        }
        public static Complex operator -(double d, Complex c)
        {
            return new Complex(d - c.Re, -c.Im);
        }
        public static Complex operator *(Complex c1, Complex c2)
        {
            return new Complex(c1.Re * c2.Re - c1.Im * c2.Im, c1.Re * c2.Im + c2.Re * c1.Im);
        }
        public static Complex operator *(Complex c, double d)
        {
            return new Complex(c.Re * d, c.Im * d);
        }
        public static Complex operator *(double d, Complex c)
        {
            return new Complex(c.Re * d, c.Im * d);
        }
        public static Complex operator /(Complex c, double d)
        {
            return new Complex(c.Re / d, c.Im / d);
        }
        public static Complex operator /(double d, Complex c)
        {
            double temp = d / (c.Re * c.Re + c.Im * c.Im);
            return new Complex(c.Re * temp, -c.Im * temp);
        }
        public static Complex operator /(Complex c1, Complex c2)
        {
            double temp = 1 / (c2.Re * c2.Re + c2.Im * c2.Im);
            return new Complex((c1.Re * c2.Re + c1.Im * c2.Im) * temp, (-c1.Re * c2.Im + c2.Re * c1.Im) * temp);
        }
    }

 ==================

20200208補充源工程，下載地址：

https://files.cnblogs.com/files/Si-Mao/FFT.rar