图像处理之离散傅里叶变换1-CFANZ编程社区

1 概述

离散傅里叶变换（Discrete Fourier Transform，缩写为DFT），是指傅里叶变换在时域和频域上都呈现离散的形式，将时域信号的采样变换为在离散时间傅里叶变换（DTFT）频域的采样。在形式上，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号做DFT，也应当对其经过周期延拓成为周期信号再进行变换。在实际应用中，通常采用快速傅里叶变换来高效计算DFT。

2 离散傅里叶变换的原理

简单来说，对一张图像使用傅里叶变换就是将它分解成正弦和余弦两部分，也就是将图像从空间域（spatial domain）转换到频域（frequency domain）。

这一转换的理论基础为：任一函数都可以表示成无数个正弦和余弦函数的和的形式。傅里叶变换就是一个用来将函数分解的工具。

二维图像的傅里叶变换可以用以下数学公式表达。

图像处理之离散傅里叶变换1_傅里叶变换

式中f是空间域（spatial domain）值，F是频域（frequency domain）值。转换之后的频域值是复数，因此，显示傅里叶变换之后的结果需要使用实数图像（real image）加虚数图像（complex image），或者幅度图像（magitude image）加相位图像（phase image）的形式。在实际的图像处理过程中，仅仅使用了幅度图像，因为幅度图像包含了原图像的几乎所有我们需要的几何信息。然而，如果想通过修改幅度图像或者相位图像的方法来间接修改原空间图像，需要使用逆傅里叶变换得到修改后的空间图像，这样就必须同时保留幅度图像和相位图像了。

在此示例中，我们将展示如何计算以及显示傅里叶变换后的幅度图像。由于数字图像的离散性，像素值的取值范围也是有限的。比如在一张灰度图像中，像素灰度值一般在0到255之间。因此，我们这里讨论的也仅仅是离散傅里叶变换（DFT）。如果需要得到图像中的几何结构信息，那你就要用到它了。

在频域里面，对于一幅图像，高频部分代表了图像的细节、纹理信息；低频部分代表了图像的轮廓信息。如果对一幅精细的图像使用低通滤波器，那么滤波后的结果就只剩下轮廓了。这与信号处理的基本思想是相通的。如果图像受到的噪声恰好位于某个特定的“频率”范围内，则可以通过滤波器来恢复原来的图像。傅里叶变换在图像处理中可以做到图像增强与图像去噪、图像分割之边缘检测、图像特征提取、图像压缩等。

3 dft()函数详解

dft函数的作用是对一维或二维浮点数数组进行正向或反向离散傅里叶变换。

C++:

void dft(InputArray src,OutputArray dst,int flags=0,int nonzeroRows=0)

第一个参数，InputArray类型的src。输入矩阵，可以为实数或者虚数。
第二个参数，OutputArray 类型的dst。函数调用后的运算结果存在这里，其尺寸和类型取决于标识符，也就是第三个参数 flags。
第三个参数，int类型的flags。转换的标识符，有默认值0，取值可以为下面中标识符的结合。

标识符取值列表标识符名称	意义
DFT_INVERSE	用一维或二维逆变换代替默认的正向变换
DFT_SCALE	缩放比例标识符，输出的结果都会以 1／N进行缩放，通常会结合DFT＿INVERSE一起使用
DFT_ROWS	对输入矩阵的每行进行正向或反向的变换，此标识符可以在处理多种矢量的时候用于减小资源开销，这些处理常常是三维或高维变换等复杂操作
DFT_COMPLEX_OUTPUT	行一维或二维实数数组正变换。这样的结果虽然是复数阵列，但拥有复数的共轭对称性（CCS），所以可以被写成一个拥有同样尺寸的实数阵列
DFT_REAL_OUTPUT	进行一维或二维复数数组反变换。这样的结果通常是一个大小相同的复矩阵。如果输入的矩阵有复数的共轭对称性（比如是一个带有DFT＿COMPLEX＿OUTPUT标识符的正变换结果），便会输出实矩阵

第四个参数，int 类型的 nonzeroRows，有默认值0。当此参数设为非零时（最好是取值为想要处理的那一行的值，比如C.rows），函数会假设只有输入矩阵的第一个非零行包含非零元素（没有设置DFT＿INVERSE标识符），或只有输出矩阵的第一个非零行包含非零元素（设置了DFT＿INVERSE标识符）。这样的话，函数就可对其他行进行更高效的处理，以节省时间开销。这项技术尤其是在采用DFT计算矩阵卷积时非常有效。

讲解完函数的参数含义，下面我们看一个用 dft 函数计算两个二维实矩阵卷积的示例核心片段。

void convolveDFT (InputArray A, InputArray B, OutputArray C) 
{
    ／／【1】初始化输出矩阵
    C.create (abs (A.rows - B.rows)+1, abs (A.cols - B.cols)+1,A.type()); 
    Size dftSize;
    ／／【2】计算DFT变换的尺寸
   dftSize.width =getOptimalDFTSize(A.cols+B.cols-1);
   dftSize.height = getOptimalDFTSize(A.rows + B.rows-1); 
   ／／【3】分配临时缓冲区并初始化置零
    Mat tempA(dftSize, A.type(), Scalar::al1(0)); 
    Mat tempB(dftSize, B.type(),Scalar::a11(0));
    ／／【4】分别复制A和B到tempA和tempB的左上角
    Mat roiA(tempA, Rect(0,0,A.cols,A.rows)); 
    A.copyTo(roiA);
    Mat roiB(tempB, Rect (0,0,B.cols,B.rows)); 
    B.copyTo(roiB);
    ／／【5】就地操作（in—place），进行快速傅里叶变换，并将nonzeroRows参数置为非零，以进行更加快速地处理
    dft(tempA, tempA, 0, A.rows);
    dft(tempB, tempB, 0, B.rows);
    ／／【6】将得到的频谱相乘，结果存放于tempA中
    mulSpectrums(tempA, tempB,tempA);
    ／／【7】将结果变换为频域，且尽管结果行（result rows）都为非零，我们只需要其中的C.rows的第一行，所以采用 nonzeroRows＝＝C.rows
    dft(tempA, tempA, DFT_INVERSE + DFT_SCALE,C.rows);
    ／／【8】将结果复制到C中
    tempA(Rect(0, 0, C.cols, C.rows)).copyTo(C);
    ／／所有的临时缓冲区将被自动释放，所以无须收尾操作
}

此示例中的注释已经十分详尽。其中出现了新的函数MulSpectrums，它的作用是计算两个傅里叶频谱的每个元素的乘法，前两个参数为输入的参加乘法运算的两个矩阵，第三个参数为得到的乘法结果矩阵。