Opencv one dimensional frequency domain filter

author ： Zhai Tianbao Steven
Copyright notice ： The copyright belongs to the author , Commercial reprint please contact the author for authorization , Non-commercial reprint please indicate the source

Scenario requirements

        Previously, I shared the functions of many image frequency domain filters , Such as ideal filter 、 Gauss filter 、 Butterworth filter and so on , Are written in two-dimensional form ; Some fans proposed , I don't know how to do similar processing for one-dimensional data . In fact, understand the principle , It is also very convenient to reduce two dimensions to one dimension , For the convenience of everyone , This article will use C++ and OpenCV Realize the function of one-dimensional frequency domain filter .

        An example is Butterworth , If you want to change to something else , Just replace the filter core .

Related function C++ Implementation code

        The difference between two dimensions and one dimension ：

1. obtain DFT When changing the optimal height , One dimension only needs to consider columns or rows , Two dimensional needs to be considered .
2. Spectrum transfer operation , Two dimensions are four quadrants , One dimension is two parts , That is, up and down or left and right can be exchanged .
3. The size of the filter in the frequency domain is related to the distance from the center , The two-dimensional distance formula is a circular expression , The one-dimensional distance formula is the absolute value of the difference .

//  Image boundary processing （ A one-dimensional ）
cv::Mat image_make_border_onedim(cv::Mat &src)
{
	int h = cv::getOptimalDFTSize(src.rows); //  obtain DFT The optimal height of the transformation 

	cv::Mat padded;
	//  The constant method expands the image boundary , Constant  = 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 The spectrum is moved after transformation （ A one-dimensional ）
void fftshift_onedim(cv::Mat &plane0, cv::Mat &plane1)
{
	//  The following operation is to move the image   ( Zero frequency shift to the center )
	int cy = plane0.rows / 2;
	cv::Mat part1_r(plane0, cv::Rect(0, 0, 1, cy));  //  The element coordinates are expressed as (cx, cy)
	cv::Mat part2_r(plane0, cv::Rect(0, cy, 1, cy));

	cv::Mat temp;
	part1_r.copyTo(temp);  // Exchange positions up and down ( Real component )
	part2_r.copyTo(part1_r);
	temp.copyTo(part2_r);

	cv::Mat part1_i(plane1, cv::Rect(0, 0, 1, cy));  // Element coordinates (cx,cy)
	cv::Mat part2_i(plane1, cv::Rect(0, cy, 1, cy));

	part1_i.copyTo(temp);  // Exchange positions up and down ( Imaginary part )
	part2_i.copyTo(part1_i);
	temp.copyTo(part2_i);
}

// pow operation 
Mat powZ(cv::InputArray src, double power) {
	cv::Mat dst;
	cv::pow(src, power, dst);
	return dst;
}

// sqrt operation 
Mat sqrtZ(cv::InputArray src) {
	cv::Mat dst;
	cv::sqrt(src, dst);
	return dst;
}

//  Frequency domain filtering 
cv::Mat frequency_filter(cv::Mat &scr, cv::Mat &blur)
{
	cv::Mat mask = scr == scr;
	scr.setTo(0.0f, ~mask);

	// Create channels , Storage dft The real part and imaginary part after （CV_32F, Must be the number of single channels ）
	cv::Mat plane[] = { scr.clone(), cv::Mat::zeros(scr.size() , CV_32FC1) };

	cv::Mat complexIm;
	cv::merge(plane, 2, complexIm); //  Merge channel  （ Merge the two matrices into one 2 The tunnel Mat Class container ）
	cv::dft(complexIm, complexIm); //  Carry out Fourier transform , The results are saved in itself 

	//  The separation channel （ Array separation ）
	cv::split(complexIm, plane);

	//  The following operation is frequency domain migration 
	fftshift_onedim(plane[0], plane[1]);

	// ***************** Filter functions and DFT The product of the result ****************
	cv::Mat blur_r, blur_i, BLUR;
	cv::multiply(plane[0], blur, blur_r);  //  wave filtering （ The real part is multiplied by the corresponding element of the filter template ）
	cv::multiply(plane[1], blur, blur_i);  //  wave filtering （ The imaginary part is multiplied by the corresponding element of the filter template ）
	cv::Mat plane1[] = { blur_r, blur_i };

	//  Move back again for inverse transformation 
	fftshift_onedim(plane1[0], plane1[1]);
	cv::merge(plane1, 2, BLUR); //  The real part and the imaginary part merge 

	cv::idft(BLUR, BLUR);       // idft The result is also plural 
	BLUR = BLUR / BLUR.rows / BLUR.cols;

	cv::split(BLUR, plane);// The separation channel , The main access channel 

	return plane[0];
}

//  Butterworth low pass filter kernel function （ A one-dimensional ）
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;// radius D0 The smaller it is , The greater the blur ; radius D0 The bigger it is , The smaller the blur 
	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;
}

//  Butterworth low pass filter （ A one-dimensional ）
cv::Mat butterworth_low_pass_filter_onedim(cv::Mat &src, float d0, int n)
{
	// H = 1 / (1+(D/D0)^2n)   n Indicates the number of Butterworth filters 
	//  Order n=1  No ringing and negative values      Order n=2  Slight ringing and negative value    Order n=5  Obvious ringing and negative value     Order n=20  And ILPF be similar 
	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;
}

Test code

#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include<opencv2/opencv.hpp>
#include<ctime>

using namespace cv;
using namespace std;

//  Reading data 
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;
}

//  Display data curve 
cv::Mat ShowDataCurve(cv::Mat input)
{
	//  normalization 
	cv::Mat curve = input.clone();
	cv::normalize(curve, curve, 0, 255, cv::NORM_MINMAX);
	curve.convertTo(curve, CV_8UC1);

	//  Draw a simple curve 
	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;
}

//  Image boundary processing （ A one-dimensional ）
cv::Mat image_make_border_onedim(cv::Mat &src)
{
	int h = cv::getOptimalDFTSize(src.rows); //  obtain DFT The optimal height of the transformation 

	cv::Mat padded;
	//  The constant method expands the image boundary , Constant  = 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 The spectrum is moved after transformation （ A one-dimensional ）
void fftshift_onedim(cv::Mat &plane0, cv::Mat &plane1)
{
	//  The following operation is to move the image   ( Zero frequency shift to the center )
	int cy = plane0.rows / 2;
	cv::Mat part1_r(plane0, cv::Rect(0, 0, 1, cy));  //  The element coordinates are expressed as (cx, cy)
	cv::Mat part2_r(plane0, cv::Rect(0, cy, 1, cy));

	cv::Mat temp;
	part1_r.copyTo(temp);  // Exchange positions up and down ( Real component )
	part2_r.copyTo(part1_r);
	temp.copyTo(part2_r);

	cv::Mat part1_i(plane1, cv::Rect(0, 0, 1, cy));  // Element coordinates (cx,cy)
	cv::Mat part2_i(plane1, cv::Rect(0, cy, 1, cy));

	part1_i.copyTo(temp);  // Exchange positions up and down ( Imaginary part )
	part2_i.copyTo(part1_i);
	temp.copyTo(part2_i);
}

// pow operation 
Mat powZ(cv::InputArray src, double power) {
	cv::Mat dst;
	cv::pow(src, power, dst);
	return dst;
}

// sqrt operation 
Mat sqrtZ(cv::InputArray src) {
	cv::Mat dst;
	cv::sqrt(src, dst);
	return dst;
}

//  Frequency domain filtering 
cv::Mat frequency_filter(cv::Mat &scr, cv::Mat &blur)
{
	cv::Mat mask = scr == scr;
	scr.setTo(0.0f, ~mask);

	// Create channels , Storage dft The real part and imaginary part after （CV_32F, Must be the number of single channels ）
	cv::Mat plane[] = { scr.clone(), cv::Mat::zeros(scr.size() , CV_32FC1) };

	cv::Mat complexIm;
	cv::merge(plane, 2, complexIm); //  Merge channel  （ Merge the two matrices into one 2 The tunnel Mat Class container ）
	cv::dft(complexIm, complexIm); //  Carry out Fourier transform , The results are saved in itself 

	//  The separation channel （ Array separation ）
	cv::split(complexIm, plane);

	//  The following operation is frequency domain migration 
	fftshift_onedim(plane[0], plane[1]);

	// ***************** Filter functions and DFT The product of the result ****************
	cv::Mat blur_r, blur_i, BLUR;
	cv::multiply(plane[0], blur, blur_r);  //  wave filtering （ The real part is multiplied by the corresponding element of the filter template ）
	cv::multiply(plane[1], blur, blur_i);  //  wave filtering （ The imaginary part is multiplied by the corresponding element of the filter template ）
	cv::Mat plane1[] = { blur_r, blur_i };

	//  Move back again for inverse transformation 
	fftshift_onedim(plane1[0], plane1[1]);
	cv::merge(plane1, 2, BLUR); //  The real part and the imaginary part merge 

	cv::idft(BLUR, BLUR);       // idft The result is also plural 
	BLUR = BLUR / BLUR.rows / BLUR.cols;

	cv::split(BLUR, plane);// The separation channel , The main access channel 

	return plane[0];
}

//  Butterworth low pass filter kernel function （ A one-dimensional ）
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;// radius D0 The smaller it is , The greater the blur ; radius D0 The bigger it is , The smaller the blur 
	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;
}

//  Butterworth low pass filter （ A one-dimensional ）
cv::Mat butterworth_low_pass_filter_onedim(cv::Mat &src, float d0, int n)
{
	// H = 1 / (1+(D/D0)^2n)   n Indicates the number of Butterworth filters 
	//  Order n=1  No ringing and negative values      Order n=2  Slight ringing and negative value    Order n=5  Obvious ringing and negative value     Order n=20  And ILPF be similar 
	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)
{
	//  Reading data 
	string datafile = "data4.txt";
	cv::Mat data = ReadOneDimData(datafile);

	//  Get original curve 
	cv::Mat test = data.col(1).clone();
	cv::Mat show1 = ShowDataCurve(test);

	//  Get low-pass data curve 
	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;
}

The test results

chart 1 Raw data

chart 2 Filtered data

        In order to facilitate everyone to see the difference , I wrote a simple data curve display function , Make do with it ~ The above example data is a set of data of the spectrometer , Raw data is information with high frequency changes , After low-pass filter processing , Retain low-frequency components , Filter out high-frequency components , So the curve becomes smooth .

        in addition , If there is something wrong with my code , Welcome to raise objections, criticisms and corrections , Progress together ~

        If the article helps you , You can like it and let me know , I'll be happy ~ come on. ！