本文共 2958 字，大约阅读时间需要 9 分钟。

最近由于作业原因，试着用OpenCV实现频率域滤波，但是OpenCV中并没有像MATLAB中fftshift这样的中心化操作，所以我写了一个频率域滤波的函数，以后用频率域滤波的时候在主函数中调用即可。当然，水平有限，编写的代码并不优美，有问题请大家留言批评指正。

在这里我不介绍傅里叶变换，频率域滤波和高斯低通滤波器的原理，想必大家已经有了大概了解，本文关注于OpenCV中的代码实现。

废话少说，先上一张例图：

这幅图像的噪声为周期性噪声，带用阻滤波器或者陷波滤波器去噪较好，这里用高斯低通去噪，效果一般，周期性没有完全去除。

下面是频率域滤波的流程：

1、计算原图像的DFT，得到F（U，V）；

2、将频谱F（U，V）的零频点移到中心位置；

3、计算滤波器函数H（U，V）与上步结果F（U，V）的乘积G（U，V）；

4、将上步的IDFT，得到g（X，Y）；

5、取g（X，Y）的模作为结果。

下面根据上述步骤编写代码如下：

#include 
   
    #include 
    
     #include 
     
      using namespace std;using namespace cv;Mat gaussianlbrf(Mat scr, float sigma) {    Mat gaussianBlur(scr.size(), CV_32FC1);    float d0 = 2 * sigma * sigma;    for (int i = 0; i < scr.rows; ++i) {        for (int j = 0; j < scr.cols; ++j) {            float d = powf((i - scr.rows / 2) * (i - scr.rows / 2) + (j - scr.cols / 2) * (j - scr.cols / 2));            gaussianBlur.at
      
       (i, j) = expf(-d / d0);        }    }    return gaussianBlur;}Mat freqfilt(Mat scr, Mat blur) {    Mat complexIm;    split(scr, complexIm);    dft(complexIm, complexIm);        // 将零频点移到中心    int cx = complexIm.cols / 2;    int cy = complexIm.rows / 2;        // 处理实部    Mat part1_r(complexIm, Rect(0, 0, cx, cy));    Mat part2_r(complexIm, Rect(cx, 0, cx, cy));    Mat part3_r(complexIm, Rect(0, cy, cx, cy));    Mat part4_r(complexIm, Rect(cx, cy, cx, cy));        Mat temp;    part1_r.copyTo(temp);    part4_r.copyTo(part1_r);    temp.copyTo(part4_r);    part2_r.copyTo(temp);    part3_r.copyTo(part2_r);    temp.copyTo(part3_r);        // 处理虚部    Mat part1_i(complexIm, Rect(0, 0, cx, cy));    Mat part2_i(complexIm, Rect(cx, 0, cx, cy));    Mat part3_i(complexIm, Rect(0, cy, cx, cy));    Mat part4_i(complexIm, Rect(cx, cy, cx, cy));        part1_i.copyTo(temp);    part4_i.copyTo(part1_i);    temp.copyTo(part4_i);    part2_i.copyTo(temp);    part3_i.copyTo(part2_i);    temp.copyTo(part3_i);        // 滤波    Mat blur_r, blur_i;    multiply(complexIm, blur, blur_r);    multiply(complexIm, blur, blur_i);        Mat BLUR;    merge({blur_r, blur_i}, BLUR);        magnitude(complexIm, complexIm, complexIm);    complexIm += Scalar::all(1);    log(complexIm, complexIm);    normalize(complexIm, complexIm, 1, 0, CV_MINMAX);        idft(BLUR, BLUR);    magnitude(BLUR, BUR, BUR);    normalize(BUR, BUR, 1, 0, CV_MINMAX);        return BUR;}
      
     
    
   

int main() {    Mat input = imread("imageHill.jpg", CV_LOAD_IMAGE_GRAYSCALE);    imshow("原图", input);        int w = getOptimalDFTSize(input.cols);    int h = getOptimalDFTSize(input.rows);        Mat padded;    copyMakeBorder(input, padded, 0, h - input.rows, 0, w - input.cols, BORDER_CONSTANT, Scalar::all(0));    padded.convertTo(padded, CV_32FC1);        Mat gaussianKernel = gaussianlbrf(padded, 45);    Mat out = freqfilt(padded, gaussianKernel);        imshow("结果图 sigma=40", out);    waitKey(0);    return 0;}

注：1，结果的黑边为填充效果。

2，务必注意矩阵数据类型，做DFT，必须为浮点型，计算滤波器函数H（U，V）与DFT结果的乘积G（U，V）时，数据类型务必一致，包括通道数，单通道类型不能与多通道类型相乘，关于数据类型，可以参看。

参考：