场景需求
之前分享了许多图像频域滤波器的功能，如理想滤波器、高斯滤波器、巴特沃斯滤波器等等，都是二维形式的写法；有粉丝提出，不知道一维数据该怎么进行类似的处理。其实搞懂原理，二维降一维也是非常方便的，为了方便大家使用，本文将用C++和OpenCV实现一维频域滤波器功能。
示例为巴特沃斯，如果想换成别的，只要把滤波核更换即可。
相关功能函数的C++实现代码
二维与一维的不同之处：
- 获取DFT变换最佳高度时，一维只需要考虑列或者行，二维需要都考虑。
- 频谱搬移操作，二维是四象限，一维是两部分，即上下或者左右交换即可。
- 频域中滤波器尺寸与距离中心的距离有关，二维的距离公式是一个圆表达式，一维的距离公式是差值绝对值。
// 图像边界处理（一维） cv::Mat image_make_border_onedim(cv::Mat &src) { int h = cv::getOptimalDFTSize(src.rows); // 获取DFT变换的最佳高度 cv::Mat padded; // 常量法扩充图像边界，常量 = 0 cv::copyMakeBorder(src, padded, 0, h - src.rows, 0, 0, cv::BORDER_CONSTANT, cv::Scalar::all(0)); padded.convertTo(padded, CV_32FC1); return padded; } // fft变换后进行频谱搬移（一维） void fftshift_onedim(cv::Mat &plane0, cv::Mat &plane1) { // 以下的操作是移动图像 (零频移到中心) int cy = plane0.rows / 2; cv::Mat part1_r(plane0, cv::Rect(0, 0, 1, cy)); // 元素坐标表示为(cx, cy) cv::Mat part2_r(plane0, cv::Rect(0, cy, 1, cy)); cv::Mat temp; part1_r.copyTo(temp); //上与下交换位置(实部) part2_r.copyTo(part1_r); temp.copyTo(part2_r); cv::Mat part1_i(plane1, cv::Rect(0, 0, 1, cy)); //元素坐标(cx,cy) cv::Mat part2_i(plane1, cv::Rect(0, cy, 1, cy)); part1_i.copyTo(temp); //上与下交换位置(虚部) part2_i.copyTo(part1_i); temp.copyTo(part2_i); } // pow操作 Mat powZ(cv::InputArray src, double power) { cv::Mat dst; cv::pow(src, power, dst); return dst; } // sqrt操作 Mat sqrtZ(cv::InputArray src) { cv::Mat dst; cv::sqrt(src, dst); return dst; } // 频率域滤波 cv::Mat frequency_filter(cv::Mat &scr, cv::Mat &blur) { cv::Mat mask = scr == scr; scr.setTo(0.0f, ~mask); //创建通道，存储dft后的实部与虚部（CV_32F，必须为单通道数） cv::Mat plane[] = { scr.clone(), cv::Mat::zeros(scr.size() , CV_32FC1) }; cv::Mat complexIm; cv::merge(plane, 2, complexIm); // 合并通道 （把两个矩阵合并为一个2通道的Mat类容器） cv::dft(complexIm, complexIm); // 进行傅立叶变换，结果保存在自身 // 分离通道（数组分离） cv::split(complexIm, plane); // 以下的操作是频域迁移 fftshift_onedim(plane[0], plane[1]); // *****************滤波器函数与DFT结果的乘积**************** cv::Mat blur_r, blur_i, BLUR; cv::multiply(plane[0], blur, blur_r); // 滤波（实部与滤波器模板对应元素相乘） cv::multiply(plane[1], blur, blur_i); // 滤波（虚部与滤波器模板对应元素相乘） cv::Mat plane1[] = { blur_r, blur_i }; // 再次搬移回来进行逆变换 fftshift_onedim(plane1[0], plane1[1]); cv::merge(plane1, 2, BLUR); // 实部与虚部合并 cv::idft(BLUR, BLUR); // idft结果也为复数 BLUR = BLUR / BLUR.rows / BLUR.cols; cv::split(BLUR, plane);//分离通道，主要获取通道 return plane[0]; } // 巴特沃斯低通滤波核函数（一维） cv::Mat butterworth_low_kernel_onedim(cv::Mat &scr, float sigma, int n) { cv::Mat butterworth_low_pass(scr.size(), CV_32FC1); //，CV_32FC1 float D0 = sigma;//半径D0越小，模糊越大；半径D0越大，模糊越小 for (int i = 0; i < scr.rows; i++) { float d = abs(float(i - scr.rows / 2)); butterworth_low_pass.at<float>(i, 0) = 1.0f / (1.0f + pow(d / D0, 2 * n)); } return butterworth_low_pass; } // 巴特沃斯低通滤波（一维） cv::Mat butterworth_low_pass_filter_onedim(cv::Mat &src, float d0, int n) { // H = 1 / (1+(D/D0)^2n) n表示巴特沃斯滤波器的次数 // 阶数n=1 无振铃和负值 阶数n=2 轻微振铃和负值 阶数n=5 明显振铃和负值 阶数n=20 与ILPF相似 cv::Mat padded = image_make_border_onedim(src); cv::Mat butterworth_kernel = butterworth_low_kernel_onedim(padded, d0, n); cv::Mat result = frequency_filter(padded, butterworth_kernel); return result; }
测试代码
#include <iostream> #include <fstream> #include <sstream> #include <vector> #include<opencv2/opencv.hpp> #include<ctime> using namespace cv; using namespace std; // 读取数据 cv::Mat ReadOneDimData(string &file) { cv::Mat result; std::ifstream in(file); std::string str; getline(in, str); while (str != "") { int a = str.find(','); string front = str.substr(0, a); string back = str.substr(a + 1); float wave = stof(front); float h = stof(back); cv::Mat temp = (cv::Mat_<float>(1, 2) << wave, h); result.push_back(temp); getline(in, str); } return result; } // 显示数据曲线 cv::Mat ShowDataCurve(cv::Mat input) { // 归一化 cv::Mat curve = input.clone(); cv::normalize(curve, curve, 0, 255, cv::NORM_MINMAX); curve.convertTo(curve, CV_8UC1); // 绘制简易曲线图 int size = input.rows; cv::Mat pic = cv::Mat::zeros(256, size, CV_8UC1); for (int i = 0; i < size; ++i) { uchar h = curve.at<uchar>(i, 0); pic.at<uchar>(255 - h, i) = 255; } return pic; } // 图像边界处理（一维） cv::Mat image_make_border_onedim(cv::Mat &src) { int h = cv::getOptimalDFTSize(src.rows); // 获取DFT变换的最佳高度 cv::Mat padded; // 常量法扩充图像边界，常量 = 0 cv::copyMakeBorder(src, padded, 0, h - src.rows, 0, 0, cv::BORDER_CONSTANT, cv::Scalar::all(0)); padded.convertTo(padded, CV_32FC1); return padded; } // fft变换后进行频谱搬移（一维） void fftshift_onedim(cv::Mat &plane0, cv::Mat &plane1) { // 以下的操作是移动图像 (零频移到中心) int cy = plane0.rows / 2; cv::Mat part1_r(plane0, cv::Rect(0, 0, 1, cy)); // 元素坐标表示为(cx, cy) cv::Mat part2_r(plane0, cv::Rect(0, cy, 1, cy)); cv::Mat temp; part1_r.copyTo(temp); //上与下交换位置(实部) part2_r.copyTo(part1_r); temp.copyTo(part2_r); cv::Mat part1_i(plane1, cv::Rect(0, 0, 1, cy)); //元素坐标(cx,cy) cv::Mat part2_i(plane1, cv::Rect(0, cy, 1, cy)); part1_i.copyTo(temp); //上与下交换位置(虚部) part2_i.copyTo(part1_i); temp.copyTo(part2_i); } // pow操作 Mat powZ(cv::InputArray src, double power) { cv::Mat dst; cv::pow(src, power, dst); return dst; } // sqrt操作 Mat sqrtZ(cv::InputArray src) { cv::Mat dst; cv::sqrt(src, dst); return dst; } // 频率域滤波 cv::Mat frequency_filter(cv::Mat &scr, cv::Mat &blur) { cv::Mat mask = scr == scr; scr.setTo(0.0f, ~mask); //创建通道，存储dft后的实部与虚部（CV_32F，必须为单通道数） cv::Mat plane[] = { scr.clone(), cv::Mat::zeros(scr.size() , CV_32FC1) }; cv::Mat complexIm; cv::merge(plane, 2, complexIm); // 合并通道 （把两个矩阵合并为一个2通道的Mat类容器） cv::dft(complexIm, complexIm); // 进行傅立叶变换，结果保存在自身 // 分离通道（数组分离） cv::split(complexIm, plane); // 以下的操作是频域迁移 fftshift_onedim(plane[0], plane[1]); // *****************滤波器函数与DFT结果的乘积**************** cv::Mat blur_r, blur_i, BLUR; cv::multiply(plane[0], blur, blur_r); // 滤波（实部与滤波器模板对应元素相乘） cv::multiply(plane[1], blur, blur_i); // 滤波（虚部与滤波器模板对应元素相乘） cv::Mat plane1[] = { blur_r, blur_i }; // 再次搬移回来进行逆变换 fftshift_onedim(plane1[0], plane1[1]); cv::merge(plane1, 2, BLUR); // 实部与虚部合并 cv::idft(BLUR, BLUR); // idft结果也为复数 BLUR = BLUR / BLUR.rows / BLUR.cols; cv::split(BLUR, plane);//分离通道，主要获取通道 return plane[0]; } // 巴特沃斯低通滤波核函数（一维） cv::Mat butterworth_low_kernel_onedim(cv::Mat &scr, float sigma, int n) { cv::Mat butterworth_low_pass(scr.size(), CV_32FC1); //，CV_32FC1 float D0 = sigma;//半径D0越小，模糊越大；半径D0越大，模糊越小 for (int i = 0; i < scr.rows; i++) { float d = abs(float(i - scr.rows / 2)); butterworth_low_pass.at<float>(i, 0) = 1.0f / (1.0f + pow(d / D0, 2 * n)); } return butterworth_low_pass; } // 巴特沃斯低通滤波（一维） cv::Mat butterworth_low_pass_filter_onedim(cv::Mat &src, float d0, int n) { // H = 1 / (1+(D/D0)^2n) n表示巴特沃斯滤波器的次数 // 阶数n=1 无振铃和负值 阶数n=2 轻微振铃和负值 阶数n=5 明显振铃和负值 阶数n=20 与ILPF相似 cv::Mat padded = image_make_border_onedim(src); cv::Mat butterworth_kernel = butterworth_low_kernel_onedim(padded, d0, n); cv::Mat result = frequency_filter(padded, butterworth_kernel); return result; } int main(void) { // 读取数据 string datafile = "data4.txt"; cv::Mat data = ReadOneDimData(datafile); // 获取原始曲线 cv::Mat test = data.col(1).clone(); cv::Mat show1 = ShowDataCurve(test); // 获取低通数据曲线 float D0 = 50.0f; Mat lowpass = butterworth_low_pass_filter_onedim(test, D0, 2); cv::Mat show2 = ShowDataCurve(lowpass); imshow("curve1", show1); imshow("curve2", show2); waitKey(0); system("pause"); return 0; }
测试效果
图1 原始数据 
图2 滤波处理后数据
为了方便大家看出差异，我写了一个简单的数据曲线显示函数，将就着看哈~上述示例数据为光谱仪的一组数据，原始数据是高频变化的信息，经过低通滤波器处理后，保留低频成分，过滤掉高频成分，因此曲线变得平滑了。
另外，如果我的代码有什么问题，欢迎大家提出异议批评指正，一同进步~
如果文章帮助到你了，可以点个赞让我知道，我会很快乐~加油！
