鸟语天空

二维离散傅里叶变换(C#版)

作者：追风剑情 发布于：2021-3-23 12:44 分类：Algorithms

示例：图像处理

傅里叶变换类

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;
            FD2TD[i].im /= 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 centerH = 4;
        int centerW = 4;
        //在中间挖个洞(挖圆形效果更好,这里挖矩形方便写代码)
        for (int i = halfH - centerH; i < halfH + centerH; i++)
        {
            for (int j = halfW - centerW; j < halfW + centerW; 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 centerH = 32;
        int centerW = 32;
        for (int i=0; i < complex2D.Height; i++)
        {
            for (int j=0; j < complex2D.Width; j++)
            {
                if (i < halfH - centerH || i > halfH + centerH ||
                    j < halfW - centerW || j > halfW + centerW)
                {
                    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[height,width];
    }

    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[i,j];
        return row;
    }

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

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

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

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

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

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

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

    // 交换列
    public void SwapColumn(int i, int j)
    {
        for (int k = 0; k < m_height; k++)
        {
            Complex cpx0 = matrix[k,i];
            Complex cpx1 = matrix[k,j];
            matrix[k,i] = cpx1;
            matrix[k,j] = 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} ", Math.Floor(matrix[i,j].re).ToString().PadRight(5));
            sb.AppendLine();
        }
        System.IO.File.WriteAllText(string.Format("D://{0}.txt", fileName), sb.ToString());
    }
}

辅助类

using UnityEngine;
/// <summary>
/// 辅助类 For Unity
/// </summary>
public sealed class FourierHelper
{
    public static Texture2D ToTexture2D(Complex2D complex2D)
    {
        Texture2D tex = new Texture2D(complex2D.Width, complex2D.Height, TextureFormat.RGBA32, false);
        for (int i = 0; i < complex2D.Height; i++)
        {
            Complex[] cpxs = complex2D.GetRow(i);
            for (int j = 0; j < cpxs.Length; j++)
            {
                Complex cpx = cpxs[j];
                tex.SetPixel(j, i, new Color(cpx.re, cpx.re, cpx.re));
            }
        }
        tex.Apply();
        return tex;
    }

    // 转成频谱图
    public static Texture2D ToSpectrumTexture2D(Complex2D complex2D)
    {
        Texture2D tex = new Texture2D(complex2D.Width, complex2D.Height, TextureFormat.RGBA32, false);
        for (int i = 0; i < complex2D.Height; i++)
        {
            Complex[] cpxs = complex2D.GetRow(i);
            for (int j = 0; j < cpxs.Length; j++)
            {
                Complex cpx = cpxs[j];
                float p = (float)cpx.PixelAmplitude;
                tex.SetPixel(j, i, new Color(p, p, p));
            }
        }
        tex.Apply();
        return tex;
    }

    public static Complex2D ToComplex2D(Texture2D tex)
    {
        Complex2D complex2D = new Complex2D(tex.width, tex.height);
        for (int y = 0; y < tex.height; y++)
        {
            for (int x = 0; x < tex.width; x++)
            {
                Color c = tex.GetPixel(x, y);
                Complex cpx = new Complex();
                cpx.re = c.r * 0.3f + c.g * 0.59f + c.b * 0.11f; //灰度值
                cpx.im = 0;
                complex2D.SetComplex(y, x, cpx);
            }
        }
        return complex2D;
    }
}

测试——快速傅里叶变换与逆变换

using UnityEngine;
using UnityEngine.UI;
/// <summary>
/// 测试FFT->IFFT算法
/// </summary>
public class FFT_IFFT_Test : MonoBehaviour
{
    public RawImage rawImage;

    private void Start()
    {
        Texture2D tex = rawImage.texture as Texture2D;
        Complex2D complex2D = FourierHelper.ToComplex2D(tex);

        Fourier.FFT2(complex2D);
        Fourier.IFFT2(complex2D);

        Texture2D ifft_tex = FourierHelper.ToTexture2D(complex2D);
        rawImage.texture = ifft_tex;
    }
}

测试——频谱图

using UnityEngine;
using UnityEngine.UI;
/// <summary>
/// 显示频谱图
/// </summary>
public class SpectrumTexture : MonoBehaviour
{
    public RawImage rawImage;

    private void Start()
    {
        Texture2D tex = rawImage.texture as Texture2D;
        Complex2D complex2D = FourierHelper.ToComplex2D(tex);

        Fourier.FFT2(complex2D);
        Fourier.FFT2Shift(complex2D);

        Texture2D sp_tex = FourierHelper.ToSpectrumTexture2D(complex2D);
        rawImage.texture = sp_tex;
    }
}

测试——高通滤波

using UnityEngine;
using UnityEngine.UI;
/// <summary>
/// 高通滤波器
/// </summary>
public class HighPassFilter : MonoBehaviour
{
    public RawImage rawImage;

    void Start()
    {
        Texture2D tex = rawImage.texture as Texture2D;
        Complex2D complex2D = FourierHelper.ToComplex2D(tex);

        Fourier.FFT2(complex2D);
        Fourier.FFT2Shift(complex2D);
        Fourier.HighPassFilting(complex2D);
        Fourier.IFFT2(complex2D);

        Texture2D ifft_tex = FourierHelper.ToTexture2D(complex2D);
        rawImage.texture = ifft_tex;
    }
}

测试——低通滤波

using UnityEngine;
using UnityEngine.UI;
/// <summary>
/// 低通滤波器
/// </summary>
public class LowPassFilter : MonoBehaviour
{
    public RawImage rawImage;

    void Start()
    {
        Texture2D tex = rawImage.texture as Texture2D;
        Complex2D complex2D = FourierHelper.ToComplex2D(tex);

        Fourier.FFT2(complex2D);
        Fourier.FFT2Shift(complex2D);
        Fourier.LowPassFilting(complex2D);
        Fourier.IFFT2(complex2D);

        Texture2D ifft_tex = FourierHelper.ToTexture2D(complex2D);
        rawImage.texture = ifft_tex;
    }
}

运行效果(Unity中的效果)

标签: Algorithms