我写的二维傅里叶逆变换未还原数据，求帮忙调试下!

Question

先说下问题：

我用8x8的数据测试FFT2和IFFT2两个算法，发现数据未还原，求帮忙检查下代码哪写错了？

(注：如果只对每一行或者每一列做FFT和IFFT，可以还原数据。)

打印的日志下如

封装的傅里叶变换类

using System;

/// <summary>
/// 傅里叶变换
/// </summary>
public sealed class Fourier
{
    //快速傅里叶变换
    public static void FFT(Complex[] TD2FD)
    {
        FFT_Core(TD2FD, WT_LUT(TD2FD.Length, 1));
    }

    //快速傅里叶变换 (二维)
    public static void FFT2(Complex2D TD2FD)
    {
        //对每一行做FFT
        for (int i = 0; i < TD2FD.Height; i++)
        {
            Complex[] row = TD2FD.GetRow(i);
            FFT(row);
            TD2FD.SetRow(i, row);
        }

        //对每一列做FFT
        for (int i = 0; i < TD2FD.Width; i++)
        {
            Complex[] column = TD2FD.GetColumn(i);
            FFT(column);
            TD2FD.SetColumn(i, column);
        }
    }

    //快速傅里叶逆变换
    public static void IFFT(Complex[] FD2TD)
    {
        //做FFT变换
        Complex[] WT = WT_LUT(FD2TD.Length, -1);
        FFT_Core(FD2TD, WT);
        //实部除以N
        for (int i = 0; i < FD2TD.Length; i++)
            FD2TD[i].re /= FD2TD.Length;
    }

    //快速傅里叶逆变换 (二维)
    public static void IFFT2(Complex2D FD2TD)
    {
        //对每一行做IFFT
        for (int i = 0; i < FD2TD.Height; i++)
        {
            Complex[] row = FD2TD.GetRow(i);
            IFFT(row);
            FD2TD.SetRow(i, row);
        }

        //对每一列做IFFT
        for (int i = 0; i < FD2TD.Width; i++)
        {
            Complex[] column = FD2TD.GetColumn(i);
            IFFT(column);
            FD2TD.SetColumn(i, column);
        }
    }

    // 将直流分量移到频谱图的中心
    public static void FFT2Shift(Complex2D complex2D)
    {
        int halfH = complex2D.Height / 2;
        int halfW = complex2D.Width / 2;
        //将图像四个象限区域按对角线交换
        for (int i=0; i<halfH; i++)
        {
            for (int j=0; j<complex2D.Width; j++)
            {
                if (j < halfW)
                    complex2D.SwapComplex(i, j, i + halfH, j + halfW);
                else
                    complex2D.SwapComplex(i, j, i + halfH, j - halfW);
            }
        }
    }

    // 高通滤波
    public static void HighPassFilting(Complex2D complex2D)
    {
        int halfH = complex2D.Height / 2;
        int halfW = complex2D.Width / 2;
        int H4 = complex2D.Height / 8;
        int W4 = complex2D.Width / 8;
        for (int i = halfH - H4; i < halfH + H4; i++)
        {
            for (int j = halfW - W4; j < halfW + W4; j++)
            {
                Complex cpx = complex2D.GetComplex(i, j);
                cpx.re = 0;
                cpx.im = 0;
                complex2D.SetComplex(i, j, cpx);
            }
        }
    }

    // 低通滤波
    public static void LowPassFilting(Complex2D complex2D)
    {
        int halfH = complex2D.Height / 2;
        int halfW = complex2D.Width / 2;
        int H4 = complex2D.Height / 8;
        int W4 = complex2D.Width / 8;
        for (int i=0; i < complex2D.Height; i++)
        {
            for (int j=0; j < complex2D.Width; j++)
            {
                if (i < halfH - H4 || i > halfH + H4 ||
                    j < halfW - W4 || j > halfW + W4)
                {
                    Complex cpx = complex2D.GetComplex(i, j);
                    cpx.re = 0;
                    cpx.im = 0;
                    complex2D.SetComplex(i, j, cpx);
                }
            }    
        }
    }

    // 返回旋转因子查询表(twiddle factor lookup table)
    private static Complex[] WT_LUT(int N, int flag = 1)
    {
        Complex[] WT = new Complex[N];
        for (int i = 0; i < N; i++)
        {
            Complex cpx_wt = new Complex();
            float angle = (float)(-i * Math.PI * 2 / N);
            cpx_wt.re = (float)Math.Cos(angle);
            //IFFT flag=-1, FFT flag=1
            cpx_wt.im = flag * (float)Math.Sin(angle);
            WT[i] = cpx_wt;
        }
        return WT;
    }

    private static void FFT_Core(Complex[] TD2FD, Complex[] WT)
    {
        int power = (int)Math.Log(TD2FD.Length, 2);
        int butterfly;
        int p, s;
        Complex x1, x2, wt;

        //倒位排序
        BitReverse(TD2FD);

        //蝶形运算
        for (int k = 0; k < power; k++) //级数
        {
            for (int j = 0; j < 1 << (power - k - 1); j++) //组数
            {
                //每组有几个元素
                butterfly = 1 << k + 1;
                //计算参与蝶形运算的两个元素的索引
                p = j * butterfly;
                s = p + butterfly / 2;
                for (int i = 0; i < butterfly / 2; i++) //蝶形运算次数
                {
                    x1 = TD2FD[i + p];
                    x2 = TD2FD[i + s];
                    wt = WT[i * TD2FD.Length / butterfly];
                    TD2FD[i + p] = x1 + x2 * wt;
                    TD2FD[i + s] = x1 - x2 * wt;
                }
            }
        }
    }

    // 倒位排序——雷德算法
    private static void BitReverse(Complex[] array)
    {
        //倒位排序原理
        //0   1   2   3   4   5   6   7   十进制
        //000 001 010 011 100 101 110 111 十进制对应的二进制
        //000 100 010 110 001 101 011 111 码位反转
        //0   4   2   6   1   5   3   7   码位反转后对应的十进制

        int i, j, k;
        int N = array.Length;
        Complex temp;
        j = 0;

        for (i = 0; i < N - 1; i++)
        {
            if (i < j)
            {
                temp = array[i];
                array[i] = array[j];
                array[j] = temp;
            }
            // 求j的下一个倒序位
            // N/2的二进制最高位为1，其他位都为0
            // 判断最高位是否为1，可与N/2进行比较
            // 判断次高位是否为1，可与N/4进行比较
            k = N >> 1;
            //如果k<=j,表示j的最高位为1
            while (k <= j)
            {
                //当k<=j时，需要将最高位变为0
                j = j - k;
                //判断次高位是否为1,依次类推，逐个比较，直到j某个位为0
                k >>= 1;
            }
            j = j + k;//将0变为1
        }
    }

    // 打印
    public static void Print(Complex[] TD2FD)
    {
        for (int i = 0; i < TD2FD.Length; i++)
        {
            Console.WriteLine(TD2FD[i].ToString());
        }
        Console.WriteLine();
    }
}

//定义复数
public class Complex
{
    public float re;//实数部
    public float im;//虚数部

    // 幅值
    public double Amplitude
    {
        get
        {
            //测试发现取值范围为
            //min=0.0009918213, max=412.4615
            return Math.Sqrt(re * re + im * im);
        }
    }

    // 频谱图像素值
    public double PixelAmplitude
    {
        get
        {
            //幅值范围很大,需要做以下处理:
            //1. 将幅值范围调到 [1, ?]
            //2. 利用Log函数压缩范围
            //3. 将范围映射到颜色值[0,1]
            double p = Math.Log(Amplitude * 10000) / 16f;
            return p;
        }
    }

    // 相位
    public double Phase
    {
        get
        {
            return Math.Atan2(im, re);
        }
    }

    public override string ToString()
    {
        return string.Format("re={0}, im={1}", re, im);
    }

    public static Complex operator +(Complex lhs, Complex rhs)
    {
        Complex result = new Complex();
        result.re = lhs.re + rhs.re;
        result.im = lhs.im + rhs.im;
        return result;
    }

    public static Complex operator -(Complex lhs, Complex rhs)
    {
        Complex result = new Complex();
        result.re = lhs.re - rhs.re;
        result.im = lhs.im - rhs.im;
        return result;
    }

    public static Complex operator *(Complex lhs, Complex rhs)
    {
        Complex result = new Complex();
        result.re = lhs.re * rhs.re - lhs.im * rhs.im;
        result.im = lhs.re * rhs.im + lhs.im * rhs.re;
        return result;
    }

    public static Complex operator *(float lhs, Complex rhs)
    {
        Complex result = new Complex();
        result.re = lhs * rhs.re;
        result.im = lhs * rhs.im;
        return result;
    }

    public static Complex operator *(Complex lhs, float rhs)
    {
        Complex result = new Complex();
        result.re = lhs.re * rhs;
        result.im = lhs.im * rhs;
        return result;
    }
}

public class Complex2D
{
    private Complex[,] matrix;
    private int m_width;
    private int m_height;

    // width:图像宽度 height:图像高度
    public Complex2D(int width, int height)
    {
        m_width = width;
        m_height = height;
        matrix = new Complex[width, height];
    }

    public int Width { get { return m_width; } }
    public int Height { get { return m_height; } }

    public Complex[] GetRow(int i)
    {
        Complex[] row = new Complex[m_width];
        for (int j = 0; j < m_width; j++)
            row[j] = matrix[j,i];
        return row;
    }

    public void SetRow(int i, Complex[] array)
    {
        for (int j = 0; j < m_width; j++)
            matrix[j, i] = array[j];
    }

    public Complex[] GetColumn(int i)
    {
        Complex[] column = new Complex[m_height];
        for (int j = 0; j < m_height; j++)
            column[j] = matrix[i,j];
        return column;
    }

    public void SetColumn(int i, Complex[] array)
    {
        for (int j = 0; j < m_width; j++)
            matrix[i, j] = array[j];
    }

    //i: 第几行  j: 第几列
    public Complex GetComplex(int i, int j)
    {
        return matrix[j,i];
    }

    //i: 第几行  j: 第几列
    public void SetComplex(int i, int j, Complex src)
    {
        matrix[j, i] = src;
    }

    // 交换两个元素
    // i: 第几行  j: 第几列
    public void SwapComplex(int i0, int j0, int i1, int j1)
    {
        Complex tmp = matrix[j0,i0];
        matrix[j0, i0] = matrix[j1, i1];
        matrix[j1, i1] = tmp;
    }

    // 交换行
    public void SwapRow(int i, int j)
    {
        for (int k=0; k<m_width; k++)
        {
            Complex cpx0 = matrix[k,i];
            Complex cpx1 = matrix[k,j];
            matrix[k,i] = cpx1;
            matrix[k,j] = cpx0;
        }
    }

    // 交换列
    public void SwapColumn(int i, int j)
    {
        for (int k = 0; k < m_height; k++)
        {
            Complex cpx0 = matrix[i,k];
            Complex cpx1 = matrix[j,k];
            matrix[i,k] = cpx1;
            matrix[j,k] = cpx0;
        }
    }

    public void Print(string fileName)
    {
        System.Text.StringBuilder sb = new System.Text.StringBuilder();
        for (int i = 0; i < m_height; i++)
        {
            for (int j = 0; j < m_width; j++)
                sb.AppendFormat("{0:G} ", matrix[j,i].re.ToString().PadRight(5));
            sb.AppendLine();
        }
        System.IO.File.WriteAllText(string.Format("D://{0}.txt", fileName), sb.ToString());
    }
}

测试代码

//测试
        Complex2D complex2D = new Complex2D(8, 8);
        for (int i = 0; i < complex2D.Height; i++)
        {
            for (int j = 0; j < complex2D.Width; j++)
            {
                Complex cpx = new Complex();
                cpx.re = i * complex2D.Width + j;
                cpx.im = 0;
                complex2D.SetComplex(i, j, cpx);
            }
        }

        complex2D.Print("8x8");
        Fourier.FFT2(complex2D);
        complex2D.Print("fft2");
        Fourier.IFFT2(complex2D);
        complex2D.Print("ifft2");

Answer 1

ok,搞定了；

引发错误的地方在于：傅里叶逆变换，实数部分和虚数部分都需要除N，你只给实数部分除了N，一维数据没有问题是因为，你的一维数据全是实数，逆变换后虚数部分为0，所以不影响；而二维傅里叶变换需要两次逆变换（行和列的逆变换），因此就有了较大影响；

public static void IFFT(Complex[] FD2TD)
    {
        //做FFT变换
        Complex[] WT = WT_LUT(FD2TD.Length, -1);
        FFT_Core(FD2TD, WT);
        //实部除以N
        for (int i = 0; i < FD2TD.Length; i++)
        {
            FD2TD[i].re /= FD2TD.Length;
            FD2TD[i].im /= FD2TD.Length;
        }
    }

这样改一下就好了，我不太会输出txt，我用控制台输出试了一下，这是逆变换以后的数据：

ok,误差在可忽略的范围内；

 

还有一个规范性的问题，一般矩阵（row,columns）的行是row(Height),列是columns(Width),你刚开始初始化时确实是这么定义的，但后面你就全部反着来了，因为这次输入的是8*8的数据，行列数相等，没有报错，如果下个数据行列数不相等，就有错了，还是规范起来写的好。（这个问题不影响你的这次的结果，但还是改掉的好）