FFT與二維FFT的C#實現


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

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

1. FFT編程依據

  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


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM