下面是频域滤波示例程序：
在本程序中，共有五个自定义函数，分别是：
1. myMagnitude（），在该函数中封装了Opencv中的magnitude函数，实现对于复数图像的幅值计算。
2. dftshift()，该函数实现对图像四个象限的对角互换，相当于MatLab中 fftshift（），将频谱的原点（0，0）移到图像中心。示例1中采用了该函数实现了频谱图中心化。
3. srcCentralized()用于傅里叶变换前的预处理，以便得到傅里叶频谱的原点（0，0）位于图像的中心。
该函数与dftshift（）目的一致，实现方法不同，一个是变换前预处理，一个是变换后处理。
示例2中采用了该函数实现了频谱图中心化。
4. displayDftSpectrum()，该函数用于显示复数图像。
5. makeFilter（），用于制作频域滤波器，该函数利用了ptr（）
指针遍历图像的方法，实现了圆等滤波函数。
下面是两个示例，分别采用了dftshift()、srcCentralized()实现频谱图的中心化。示例2与冈萨雷斯的《数字图像处理（第3版）》的4.6.7小结中的流程一致。
在主程序内容，参见注释。
示例1：在该程序中采用了第2个函数dftshift()，实现了频谱图中心化。
#include <iostream>
#include<opencv2/opencv.hpp>
using namespace cv;
using namespace std; void myMagnitude(Mat & complexImg,Mat & mI)
{
Mat planes[];
split(complexImg,planes);
magnitude(planes[],planes[],mI);
}
void dftshift(Mat& ds)
{
int cx=ds.cols/;//图像的中心点x 坐标
int cy=ds.rows/;//图像的中心点y 坐标
Mat q0=ds(Rect(,,cx,cy));//左上
Mat q1=ds(Rect(cx,,cx,cy));//右上
Mat q2=ds(Rect(,cy,cx,cy));//左下
Mat q3=ds(Rect(cx,cy,cx,cy));//右下
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
}
void displayDftSpectrum(Mat&dftDst,String winName,bool inverseSpectrum)
{
Mat magI;
myMagnitude(dftDst,magI);
if(!inverseSpectrum)//如果是正向傅里叶变换谱
{
magI+=Scalar::all();
log(magI,magI);
}
normalize(magI,magI,,,NORM_MINMAX);
imshow(winName,magI);
}
void makeFilter(Mat&filter,int R,bool isLowPassFilter)
{
int cx=filter.cols/,cy=filter.rows/;
int R2=R*R; for(int i=;i<filter.rows;i++)
{
Vec2d* pf=filter.ptr<Vec2d>(i);
for(int j=;j<filter.cols;j++)
{
int r2=(j-cx)*(j-cx)+(i-cy)*(i-cy);
if(r2<R2)
{
pf[j][]=;
pf[j][]=pf[j][];
}
}
}
if(!isLowPassFilter)
{
filter =Scalar::all()-filter;
}
Mat displayFilter;
extractChannel(filter,displayFilter,);
imshow("filter image",displayFilter);
} int main()
{
///读入灰度图像
Mat src=imread("D:\\Qt\\MyImage\\baboon.jpg",);
int ny=src.rows,nx=src.cols;
cv::copyMakeBorder(src,src,,ny,,nx,BORDER_CONSTANT);
imshow("original image",src);
/***************傅里叶变换*************/
src.convertTo(src,CV_64FC1);//转成浮点数据类型
Mat dftDst(src.size(),CV_64FC2);//预设dft的输出结果矩阵
dft(src,dftDst,DFT_COMPLEX_OUTPUT,);//离散傅立叶变换
dftshift(dftDst);//将傅里叶变换的结果，四象限对角互换
displayDftSpectrum(dftDst,"dftSpectrum display",false);//显示傅里叶频谱图 /***************频域滤波*************/
Mat filter(src.size(),CV_64FC2,Scalar::all());
makeFilter(filter,,false);
Mat fftTemp=dftDst.mul(filter);//复数频谱图与滤波器相乘，实现频域滤波
Mat srcFiltered;
dft(fftTemp,srcFiltered,DFT_INVERSE);//逆傅里叶变换
//显示经过滤波后的图像
displayDftSpectrum(srcFiltered,"filtered image",true);//显示逆傅里叶频谱图
waitKey();
return ;
}
示例2：在该示例中采用了第3个函数srcCentralized()，实现频谱图中心化。程序流程即
在本教材的6.2.3频域滤波步骤小结中，或者在冈萨雷斯的《数字图像处理》4.6.7节中，频域滤波流程总结如下：
- 给定一幅大小为M×N的输入图像f(x,y)，从式（6.1-25）和式（6.1-26）得到填充参数P和Q。典型地，我们选择P=2M和Q=2N；
- 对f(x,y)添加必要数量的0，形成大小为P×Q填充后的图像
- 用（-1）乘以f(x,y)，进行频谱中心化的预处理；
- 计算中心化预处理过的f(x,y)的傅里叶变换，得到F(u,v);
- 生成一个实的、对称的滤波函数H(u,v)，其大小为P×Q，频谱零点位于(P/2,Q/2)处。用阵列相乘形成乘积G(u,v)=F(u,v)H(u,v)；
- 经过频域滤波后的图像：
g(x,y)={real[fft[G(u,v)]]}（-1）
其中，为忽略由于计算不准确导致的寄生复分量，选择了实部，下标P指出我们处理的是填充后的阵列； - 通过从 g(x,y)的左上象限提取M×N区域，得到最终处理结果g(x,y)。
/*求取复数矩阵的幅值*/
void myMagnitude(Mat & complexImg,Mat & mI)
{
Mat planes[];
split(complexImg,planes);
magnitude(planes[],planes[],mI);
}
/*傅里叶变换后的频谱图后处理，将傅里叶普的原点(0,0)平移到图像的中心*/
void dftshift(Mat& ds)
{
int cx=ds.cols/;//图像的中心点x 坐标
int cy=ds.rows/;//图像的中心点y 坐标
Mat q0=ds(Rect(,,cx,cy));//左上
Mat q1=ds(Rect(cx,,cx,cy));//右上
Mat q2=ds(Rect(,cy,cx,cy));//左下
Mat q3=ds(Rect(cx,cy,cx,cy));//右下
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
}
/*傅里叶变换前的预处理，以便频谱图的原点（0，0）移动到图像的中心*/
void srcCentralized(Mat& src)
{
for(int i=;i<src.rows;i++)
{
for(int j=;j<src.cols;j++)
{
float* mv=CV_MAT_ELEM2(src,float,i,j);
if((i+j)%!=)
mv[]=-mv[];//如果i+j为奇数，该像素值取负值
}
}
}
/*在窗口中显示复数图像，如果是正向傅里叶矩阵，需要取log才能显示更多频谱信息
*如果是逆傅里叶变换，通过normalize归一化后，显示频谱图*/
void displayDftSpectrum(Mat&dftDst,String winName,bool inverseSpectrum)
{
Mat magI;
myMagnitude(dftDst,magI);
if(!inverseSpectrum)//如果是正向傅里叶变换谱
{
magI+=Scalar::all();
log(magI,magI);
}
normalize(magI,magI,,,NORM_MINMAX);
imshow(winName,magI);
}
void makeFilter(Mat&filter,int R,bool isLowPassFilter)
{
int cx=filter.cols/,cy=filter.rows/;
int R2=R*R; for(int i=;i<filter.rows;i++)
{
Vec2d* pf=filter.ptr<Vec2d>(i);
for(int j=;j<filter.cols;j++)
{
int r2=(j-cx)*(j-cx)+(i-cy)*(i-cy);
if(r2<R2)
{
pf[j][]=;
pf[j][]=pf[j][];
}
}
}
if(!isLowPassFilter)
{
filter =Scalar::all()-filter;
}
Mat displayFilter;
extractChannel(filter,displayFilter,);
imshow("filter image",displayFilter);
} int main()
{
///读入灰度图像
Mat src=imread("D:\\Qt\\MyImage\\interferometer2.jpg",);
int ny=src.rows,nx=src.cols;
imshow("original image",src);
src.convertTo(src,CV_32FC1);
srcCentralized(src);
cv::copyMakeBorder(src,src,,ny,,nx,BORDER_CONSTANT); /***************傅里叶变换*************/
src.convertTo(src,CV_64FC1);//转成浮点数据类型
Mat dftDst(src.size(),CV_64FC2);//预设dft的输出结果矩阵
dft(src,dftDst,DFT_COMPLEX_OUTPUT,);//离散傅立叶变换
// dftshift(dftDst);//将傅里叶变换的结果，四象限对角互换
displayDftSpectrum(dftDst,"dftSpectrum display",false);//显示傅里叶频谱图 /***************频域滤波*************/
Mat filter(src.size(),CV_64FC2,Scalar::all());
makeFilter(filter,,true);
Mat fftTemp=dftDst.mul(filter);//复数频谱图与滤波器相乘，实现频域滤波
srcCentralized(fftTemp); Mat srcFiltered;
dft(fftTemp,srcFiltered,DFT_INVERSE);//逆傅里叶变换
// 显示经过滤波后的图像
displayDftSpectrum(srcFiltered,"filtered image",true);//显示逆傅里叶频谱图
waitKey();
return ;
}