OpenCV 离散傅里叶变换

离散傅里叶变换（Discrete Fourier Transform，缩写为DFT），是傅里叶变换在时域和频域上都呈离散的形式，将信号的时域采样变换为其DFT的频域采样。

对于N点序列{X[n]}(0 <= n <= N),它的离散傅里叶变换为:

dft() 函数

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

函数原型

void dft(
    InoutArray src,  //输入矩阵
    OutputArray dst, ///输出矩阵
    int flags = 0,  //转换的标识符
    int onozeroRows,  
)
复制代码

第三个参数，转换的标识符分为:

• DFT_INVERSE 用一维或二维逆变换代替默认的正向变换。
• DFT_SCALE 缩放比例标识符，输出的结果都会以1 / N进行缩放。
• DFT_CMPLEX_OUTPUT、DFT_REAL_OUTPUT 进行一维或二维的数组反变换。

返回DFT最优尺寸大小：getOptimalDFTSize()函数

void int getOptimalDFTSize(int vecsize);
复制代码

扩充图像边界：copyMakeBorder()函数

void copyMakeBorder(
    InputArray src,  //输入图像
    OutputArray dst,  //输出图像
    int top,  //在图像上方扩充的像素值
    int bottom,  //在图像下方扩充的像素值
    int left,  //在图像左方扩充的像素值
    int right,  //在图像右方扩充的像素值
    int borderType, //边界类型·
    const Scalars,
)
复制代码

计算二维矢量的幅值：magnitude()函数

用于计算二维矢量的幅值

void magnitude(
    InputArray x,  //表示矢量的浮点型X坐标值,即实部
    InputArray y,  //表示矢量的浮点型Y坐标值，即虚部
    OutputArray magnitude,  //输出的幅值
)
复制代码

计算自然对数：log() 函数

log()函数的功能是计算每个数组元素绝对值的自然对数

void log(
    InputArray src,
    OutputArray dst,
);

原理即为：

if(src(I) != 0)
    log|src(I)|
else
    C
复制代码

矩阵归一化：normalize()函数

void normalize(
    InputArray src,
    OutputArray dst,
    double alpha = 1,  //归一化之后的最大值，有默认值1
    double beta = 0,  //归一化之后的最大值，有默认值0
    int norm_type = NORM_L2,  //归一化类型
    int dtype = -1,   //为负数时输出矩阵和src有同样的类型，否则，它和src有同样的通道数，深度为CV_MAT_DEPTH
    InputArray mask=noArray(),  //可选的操作掩模
)
复制代码

getOptimalDFTSize函数

该函数返回给定向量尺寸的傅里叶最优尺寸大小。为了提高离散傅里叶变换的运行速度，需要扩充图像。

使用dft函数计算两个二维实矩阵卷积的示例核心片段

• 其中MulSpectrums的作用是计算两个傅里叶频谱的每个元素的乘法

void convolveDft(InputArray A, InputArray B, OutputArray C)
{
    //初始化输出矩阵
    C.create(abs(A.rows - B.rows) + 1, abs(A.cols - B.cols) + 1, A.type);

    //计算DFT变换的尺寸
    dftSize.width = getOptimalDFTSize(A.cols + B.cols - 1);
    dftSize.height = getOptimalDFTSize(A.rows + B.rows - 1);

    //分配临时缓冲区并初始化置0
    Mat tempA(dftSize,A.type(),Scalar::all(0));
    Mat tempB(dftSize,B.type(),Scalar::all(0));

    //分别复制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);

    //就地操作，进行快速傅里叶变换，并将nonzeroRows 参数置为零，以进行更加快速的处理。
    dft(tempA,tempA,0,A.rows);
    dft(tempB,tempB,0,B.rows);

    //将得到的频谱相乘，结果存放于tempA当中
    mulSpectrums(tempA,tempB,tempA);  

    //将结果变换为频域，尽管结果行（result.rows）都为非零，我们只需其中的C.rows的第一行，所以采用nonzeroRows == C.rows
    dft(tempA,tempA,EFT_INVERSE + EFT_SCALE,C.rows);

    //将结果复制到C当中
    tempA(Rect(0,0,C.cols,C.rows)).copyTo(C);

    //所有的临时缓存区将被自动释放，所以无需收尾操作
}
复制代码

例示程序：离散傅里叶变换

【1】载入原始图像

//【1】ui灰度模式读取原始图像并显示
Mat srcImage = imread("D:\\Desktop\\lena.jpg",0);
if (!srcImage.data)
{
    cout << "lena图片读取错误" << endl;
    return false;
}
imshow("原始图像", srcImage);
复制代码

【2】将图像扩大到合适的尺寸

//【2】将输入的图像延扩到最佳的尺寸，边界用0补充
int m = getOptimalDFTSize(I.rows);
int n = getOptimalDFTSize(I.cols);
//将添加的像素初始化为0
Mat padded;
copyMakeBorder(I,padded,0,m - I,rows,0,n - I.cols,BORDER_CONSTANT,Scalar::all(0));
复制代码

【3】为傅里叶变换的结果分配储存空间

傅里叶变换的结果是复数，每个原图像值，结果会有两个图像值。

为傅里叶变换的结果（实部和虚部）分配储存空间
//将planes数组组合合并成一个多通道的数组complexI
Mat planes[] = (Mat_<float>(padded),Mat::zeros(padded.size(),CV_32F));
Mat complexI;
merge(planes,2,complexI);
复制代码

【4】进行离散傅里叶变换

def(complexI,complexI);
复制代码

【5】将复数转换为幅值

离散傅里叶变换的结果是复数，对应的幅值可表示为：

//将复数转换为幅值，即=> log(1 + sqrt(Re(DFT(I)) ^ 2 + Im(DFT(I)) ^ 2))
split(complexI,planes);  //将多通道complexI分离成几个单通道数组
planes[0] = Re(DFT(I),planes[1] = Im(DFT(I))
magnitude(planes[0],planes[1],planes[0]);
Mat magnitudeImage = planes[0];
复制代码

【6】进行对数尺度缩放

傅里叶变换的幅度之范围大到不合适在屏幕上显示。高值在屏幕上显示为白点，而低值为黑点，高低值的变化无法有效判断。为了在屏幕上凸显出高低变化的连续性，可以用对数尺度来代替线性尺度，公式如下：

//进行对数尺度缩放
magnitudeImage += Scalar::all(1);
log(magnitudeImage,magnitudeImage);  //求自然对数
复制代码

【7】剪切和重分布幅度图像象限

因为在第二步中延扩了图像，那现在是时候将新添加的像素剔除了。为了方便显示，可以重新分布图像象限位置。

//若有奇数行或技术列，进行谱写裁剪
magnitudeImage = magnitudeImage(Rect(0,0,magnitudeImage.cols & -2,magnitudeImage.rows & -2));
//重新排列傅里叶图像中的象限，使得原点位于图像中心
int cx = magnitudeImage.cols / 2;
int cy = magnitudeImage.rows / 2;
Mat q0(magnitudeImage,Rect(0,0,cx,cy)); //ROI区域的左上
Mat q1(magnitudeImage,Rect(cx,0,cx,cy));  //ROI区域的右上
Mat q2(magnitudeImage,Rect(0,cy,cx,cy));  //ROI区域的左下
Mat q3(magnitudeImage,Rect(cx,cy,cx,cy));  //ROI区域的右下

//交换象限(左上与右下)
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);

//交换象限(右上与左下)
q1.copyTo(tmp);
q4.copyTo(q1);
tmp.copyTo(q4);
复制代码

【8】归一化

现在有了重分布后的幅度图，但是幅度值仍然超过可显示范围[0,1],我们可以使用 normalize()函数归一化到可显示范围。

normalize(magnitudeImage,magnitudeImage,0,1,NORM_MINMAX);
复制代码

【9】显示效果

imshow("频谱幅值"，magnitudeImage);
复制代码