场景需求
matlab的fftshift是进行傅里叶相关操作时常用的函数，傅里叶之后得到的结果频率范围是0到fs，为了便于进行频率域滤波，也便于观察频谱信息 ，通常将频率范围调整至-fs/2到fs/2，这样就将零频分量（直流分量）迁移到了图像中心，呈现的效果就是中心低频信息，四周外围是高频信息，这个实现我们就称为fftshift。C++中运用OpenCV也可以实现同matlab一致的效果，其实原理搞懂了，用数组、Eigen等其他数据结构或者运算库也都很方便实现。
功能函数代码
// fft变换后进行频谱搬移；分别将实部信息和虚部信息进行搬移 void fftshift(cv::Mat &plane0, cv::Mat &plane1) { // 以下的操作是移动图像 (零频移到中心) int cx = plane0.cols / 2; int cy = plane0.rows / 2; cv::Mat part1_r(plane0, cv::Rect(0, 0, cx, cy)); // 元素坐标表示为(cx, cy) cv::Mat part2_r(plane0, cv::Rect(cx, 0, cx, cy)); cv::Mat part3_r(plane0, cv::Rect(0, cy, cx, cy)); cv::Mat part4_r(plane0, cv::Rect(cx, cy, cx, cy)); cv::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); cv::Mat part1_i(plane1, cv::Rect(0, 0, cx, cy)); //元素坐标(cx,cy) cv::Mat part2_i(plane1, cv::Rect(cx, 0, cx, cy)); cv::Mat part3_i(plane1, cv::Rect(0, cy, cx, cy)); cv::Mat part4_i(plane1, cv::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); }
C++测试代码
#include<iostream> #include<opencv2/opencv.hpp> #include<ctime> using namespace std; using namespace cv; void fftshift(cv::Mat &plane0, cv::Mat &plane1); int main(void) { Mat test = imread("test.jpg", 0); test.convertTo(test, CV_32FC1); //创建通道，存储dft后的实部与虚部（CV_32F，必须为单通道数） cv::Mat plane[] = { test.clone(), cv::Mat::zeros(test.size() , CV_32FC1) }; cv::Mat complexIm; cv::merge(plane, 2, complexIm); // 合并通道 （把两个矩阵合并为一个2通道的Mat类容器） cv::dft(complexIm, complexIm); // 进行傅立叶变换，结果保存在自身 // 分离通道（数组分离） cv::split(complexIm, plane); // 以下的操作是频域迁移 fftshift(plane[0], plane[1]); // 创建一个简单的高斯低通滤波器 cv::Mat gaussianBlur(test.size(), CV_32FC1); //，CV_32FC1 float d0 = 5;//高斯函数参数，越小，频率高斯滤波器越窄，滤除高频成分越多，图像就越平滑 for (int i = 0; i < test.rows; i++) { for (int j = 0; j < test.cols; j++) { float d = pow(float(i - test.rows / 2), 2) + pow(float(j - test.cols / 2), 2);//分子,计算pow必须为float型 gaussianBlur.at<float>(i, j) = expf(-d / (2 * d0*d0));//expf为以e为底求幂（必须为float型） } } // *****************滤波器函数与DFT结果的乘积**************** cv::Mat blur_r, blur_i, BLUR; cv::multiply(plane[0], gaussianBlur, blur_r); // 滤波（实部与滤波器模板对应元素相乘） cv::multiply(plane[1], gaussianBlur, blur_i); // 滤波（虚部与滤波器模板对应元素相乘） cv::Mat plane1[] = { blur_r, blur_i }; // 再次搬移回来进行逆变换 fftshift(plane1[0], plane1[1]); cv::merge(plane1, 2, BLUR); // 实部与虚部合并 cv::idft(BLUR, BLUR); // idft结果也为复数 BLUR = BLUR / BLUR.rows / BLUR.cols; cv::split(BLUR, plane);//分离通道，主要获取通道 imshow("original", test/255); imshow("result", plane[0]/255); waitKey(0); system("pause"); return 0; } // fft变换后进行频谱搬移 void fftshift(cv::Mat &plane0, cv::Mat &plane1) { // 以下的操作是移动图像 (零频移到中心) int cx = plane0.cols / 2; int cy = plane0.rows / 2; cv::Mat part1_r(plane0, cv::Rect(0, 0, cx, cy)); // 元素坐标表示为(cx, cy) cv::Mat part2_r(plane0, cv::Rect(cx, 0, cx, cy)); cv::Mat part3_r(plane0, cv::Rect(0, cy, cx, cy)); cv::Mat part4_r(plane0, cv::Rect(cx, cy, cx, cy)); cv::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); cv::Mat part1_i(plane1, cv::Rect(0, 0, cx, cy)); //元素坐标(cx,cy) cv::Mat part2_i(plane1, cv::Rect(cx, 0, cx, cy)); cv::Mat part3_i(plane1, cv::Rect(0, cy, cx, cy)); cv::Mat part4_i(plane1, cv::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); }
测试效果
图1 处理前频谱图
图2 处理后频谱图
图3 滤波前后对比图
如上图所示，fftshift实现的效果就是将低频信息迁移到中心，高频信息在四周，此时用滤波器相乘，频谱图中心的信息保留，即保留了低频成分，四周的信息归0，即去除了高频成分，这就实现了我们常说的低通滤波~
如果文章帮助到你了，可以点个赞让我知道，我会很快乐~加油！
