C++实现二维离散傅里叶变换

在上一篇文章《C++实现一维离散傅里叶变换》中，我们介绍了一维信号傅立叶变换的公式和C++实现，并阐述了频域幅值的意义。

一维傅立叶变换只适用于一维信号，例如音频数据、心脑电图等。

在图像处理中，图像信号具有高度和宽度两个属性，属于二维空间信号。将图像信号从空间域转换到频域时，需使用二维离散傅立叶变换。因此，需要将傅立叶变换从一维推广至二维。

二维连续傅立叶变换公式如下：

经过离散化后，上述公式变为：

，其中u=0,1,2,…,M-1，v=0,1,2,…,N-1。              (式1-1)

相应地，其逆变换为：

对式1-1稍作变换，可得：

其中，刚好是一维傅立叶变换的离散形式。故二维傅立叶变换可理解为先对图像的每一行进行一维傅立叶变换，变换的结果构成一个新的矩阵，再对新矩阵的每一列进行一维傅立叶变换，这样经过两次变换的结果就是二维傅立叶变换的最终结果。

也就是说，二维傅里叶变换相当于先按行变换，再按列变换（也可以先按列变换，再按行变换）。最终生成的频域数据，就是由两个方向的频率强度信息叠加而成。

我们可以用伪代码来实现上述逻辑。假设一维傅里叶变换函数为dft1(vec)，二维信号高度为h，宽度为w，那么伪代码可写成如下形式：

for (int i=0; i<h; i++)

    dft1(row[i]);

for (int i=0; i<w; i++)

    dft1(column[i]);

在本文中，我们暂不采用上述伪代码思路。我们直接按式1-1来进行编程实现：

其中，e^{-i2π(ux/M+vy/N)}可转化为三角函数cos(-2π(ux/M+vy/N))+isin(-2π(ux/M+vy/N))，i为虚数单位。转化成三角函数形式，有利于计算机编程实现。只需要采用两个for循环，就可求出单个F(u,v)值，外面再嵌套两个for循环，即可求出所有F(u,v)值。

故得到有利于编程实现的离散傅里叶变换公式如下：

由于涉及复数运算，需要先定义一个复数类，并实现复数运算。

在上一篇文章《C++实现一维离散傅里叶变换》中，我们已经实现了一个复数类ComplexNumber。直接引用它即可。

（注：编译环境是VS2013，使用MFC对话框模板）

二维离散傅里叶变换的头文件如下：

Dft2.h

#pragma once
#include "ComplexNumber.h"

#define      MAX_MATRIX_SIZE                4194304             // 2048 * 2048
#define      PI                                           3.141592653

class CDft2
{
public:
	CDft2(void);
	~CDft2(void);

public:
	bool dft2(double IN const matrix[], int IN const width, int IN const height);       // 二维离散傅里叶变换
	bool idft2(LPVOID OUT *pMatrix, int OUT *width, int OUT *height);                 // 二维离散傅里叶逆变换

	void generate_spectrum();                   // 对变换结果求模，生成频谱/幅度谱
	void normalize_spectrum();                  // 对频谱进行归一化操作

	bool has_dft2_matrix();                         // 是否已存有变换结果
	bool is_shifted_to_center();                  // 是否已将频谱原点平移到图像中心

	void clear_dft2_matrix();                       // 清除已有的变换结果
	void print_matrix();                                // 打印变换结果
	void print_spectrum();                           // 打印频谱
	void shift_spectrum_to_center();          // 将频谱原点平移到图像中心

public:
	CComplexNumber *m_dft2_matrix;
	double                    *m_spectrum_data;

protected:
	bool      m_has_dft_matrix;
	bool      m_is_normalized;
	bool      m_is_spectrum_shifted;
	int         m_dft_matrix_height;
	int         m_dft_matrix_width;
};

由于二维傅里叶变换后得到的矩阵元素数值很大，并且包含实数和虚数。为便于观察分析，需要将变换后的结果进行求模，然后归一化到[0,25