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