图像二维FFT

类型声明

#include <vector>

#define PI 3.1415926

typedef struct {
    double real;
    double imag;
}complex;

using std::vector;
using uchar = unsigned char;
using Complex1D = vector<complex>;
using Complex2D = vector<vector<complex>>;
using Double1D = vector<double>;
using Double2D = vector<vector<double>>;

复数运算及相关工具函数

// a * b
complex Mul(const complex& a, const complex& b)
{
    complex res;
    res.real = a.real * b.real - a.imag * b.imag;
    res.imag = a.imag * b.real + a.real * b.imag;
    return res;
}

// a + b
complex Add(const complex& a, const complex& b)
{
    return complex{ a.real + b.real, a.imag + b.imag };
}

// a - b
complex Sub(const complex& a, const complex& b)
{
    return complex{ a.real - b.real, a.imag - b.imag };
}

// e^(ix)
complex Euler(const double& x)
{
    return complex{ cos(x), sin(x) };
}

// |a|
double Magnitude(const complex& a)
{
    return sqrt(a.real * a.real + a.imag * a.imag);
}

// 二维复数矩阵转置
Complex2D Transform(const Complex2D& a)
{
    Complex2D res(a);
    for (int i{ 0 }; i < a.size(); i++) {
        for (int j{ 0 }; j < a[i].size(); j++) {
            res[i][j] = a[j][i];
        }
    }
    return res;
}

// 二维频谱图平移(一、三象限对换，二、四象限对换)
void Shift(Complex2D& src)
{
    Complex2D tmp(src);
    int height = src.size();
    int width = src[0].size();
    int width_2 = width / 2;
    int height_2 = height / 2;
    for (int i{ 0 }; i < height; i++) {
        for (int j{ 0 }; j < width; j++) {
            int i1, j1;
            if (i < height_2)	
                i1 = i + height_2;
            else 
                i1 = i - height_2;
            if (j < width_2)	
                j1 = j + width_2;
            else 
                j1 = j - width_2;
            tmp[i1][j1] = src[i][j];
        }
    }
    for (int i{ 0 }; i < width; i++) {
        for (int j{ 0 }; j < height; j++) {
            src[i][j] = tmp[i][j];
        }
    }
}

// 返回最近（不小于）num的2的幂
int GetClosest2Power(int num)
{
    int e = 1;
    while (e < num) {
        e <<= 1;
    }
    return e;
}

一维FFT

// 一维的处理过程
// flag=1 IFFT
// flag=-1 FFT
Complex1D FFTProgress1D(const Complex1D& src, int len, int flag)
{
    // 递归出口
    if (len == 1) {
        return Complex1D{ src };
    }
    else if (len == 2) {
        return Complex1D{
            Add(src[0], src[1]), Sub(src[0], src[1])
        };
    }

    int len_2 = len / 2;
    Complex1D a(len_2);
    Complex1D b(len_2);
    for (int i{ 0 }; i < len_2; i++) {
        a[i] = src[2 * i];      // 偶数
        b[i] = src[2 * i + 1];  // 奇数
    }
    // 分治
    Complex1D&& A = FFTProgress1D(a, len_2, flag);
    Complex1D&& B = FFTProgress1D(b, len_2, flag);
    Complex1D dst(len);
    // DFT为-2PI*k/N, IDFT为2PI*k/N
    double exp = flag * 2 * PI / len;
    for (int i{ 0 }; i < len_2; i++) {
        // e^(2PI*k/N) * P(k)
        complex&& tmp = Mul(Euler(exp * i), B[i]);
        // F(k) = G(k) + e^(2PI*k/N) * P(k)
        dst[i] = Add(A[i], tmp);
        // F(k+N/2) = G(k)+e^(2PI*(k+N/2)/N) * P(k+N/2)
        dst[i + len_2] = Sub(A[i], tmp);
    }
    return dst;
}

二维FFT

// 二维的处理过程
// flag = -1 DFT
// flag = 1 IDFT
Complex2D FFTProgress2D(const Complex2D& src, int len, int flag)
{
    Complex2D src1(src);

    // 每行进行一次一维FFT
    for (int i{ 0 }; i < len; i++) {
        src1[i] = FFTProgress1D(src[i], len, flag);
    }

    // 转置
    Complex2D&& dst = Transform(src1);

    // 再对每一行进行一维DFT（实际就是转置之前的每一列）
    for (int i{ 0 }; i < len; i++) {
        dst[i] = FFTProgress1D(dst[i], len, flag);
    }

    // 再转置回来
    dst = Transform(dst);
    return dst;
}

Complex2D FFT2D(const Double2D& src)
{
    int row = src.size();
    int column = src[0].size();
    
    // 扩展到2的幂次
    int len = GetClosest2Power(row > column ? row : column);
    Complex2D src1(len);

    // 并用0填充
    for (int i{ 0 }; i < len; i++) {
        src1[i].resize(len);
        for (int j{ 0 }; j < len; j++) {
            src1[i][j].imag = 0;
            if (i < row && j < column) {
                src1[i][j].real = src[i][j];
            }
            else {
                src1[i][j].real = 0;
            }
        }
    }

    // 进入计算过程
    return FFTProgress2D(src1, len, -1);
}

Double2D IFFT2D(const Complex2D& src, int len, int realRow, int realColumn)
{
    int size = len * len;
    Complex2D&& src1 = FFTProgress2D(src, len, 1);
    Double2D dst(realRow);
    int min{ 255 }, max{ 0 };
    for (int i{ 0 }; i < realRow; i++) {
        dst[i].resize(realColumn);
        for (int j{ 0 }; j < realColumn; j++) {
            dst[i][j] = Magnitude(src1[i][j]);
            min = dst[i][j] < min ? dst[i][j] : min;
            max = dst[i][j] > max ? dst[i][j] : max;
        }
    }
    // 缩放处理
    double scale = 255.0 / (max - min);
    for (auto& items : dst) {
        for (auto& item : items) {
            item = (item - min) * scale;
        }
    }
    return dst;
}

// 将FFT2D的结果转化成频谱图用于显示
void Complex2Image(const Complex2D& src, uchar* image)
{
    Complex2D src1(src);
    int len = src.size();
    Shift(src1);
    int* tmp = new int[len * len];
    for (int i{ 0 }; i < len; i++) {
        for (int j{ 0 }; j < len; j++) {
            int index = i * len + j;
            tmp[index] = (int)Magnitude(src1[i][j]);
            if (tmp[index] < 0)	tmp[index] = 0;
        }
    }
    // scale
    {
        for (int i{ len * len - 1 }; i >= 0; i--) {
            tmp[i] /= 500;
            image[i] = tmp[i] < 255 ? tmp[i] : 255;
        }
    }
    delete[] tmp;
}

计算BMP灰度图的频谱

对256色（一个像素占一个字节）的BMP灰度图像进行处理

uchar* image = m_dib->m_BMPdata;
int height = m_dib->m_height;
int width = m_dib->m_width;

int len = GetClosest2Power(height > width ? height : width);
uchar* newImage = new uchar[len * len]{ 0 };

Double2D src(height);
for (int i{ 0 }; i < height; i++) {
    src[i].resize(width);
    for (int j{ 0 }; j < width; j++) {
        src[i][j] = (double)image[i * width + j];
    }
}
Complex2D&& dst = FFT2D(src);
Complex2Image(dst, newImage);

m_dib->m_BMPdata = newImage;
m_dib->m_height = len;
m_dib->m_width = len;

平移之前的频谱图：

（将Complex2Image中Shift那行注释掉）

平移之后的频谱图：