Savitzky–Golay滤波 matlab/C++代码编写

       最近准备把matlab sgolayfilt函数转成opencv + C++代码，网上溜达了一圈，看到有配置GSL库，调用相关函数实现的C代码。GSL编译好了，觉得程序太复杂，决定看一看原理自己写吧。原理什么的参见下面的连接，话不多说，直接上代码。

相关参考连接

1、Savitzky-Golay 滤波器详解及C/matlab语言  https://blog.csdn.net/shenziheng1/article/details/53391422

2、C++编写S-G滤波  https://download.csdn.net/download/mr_wwww/10370036

3、Savitzky-Golay 卷积平滑算法 http://www.doc88.com/p-7788316275187.html

4、Windows下编译64为GSL  https://www.cnblogs.com/Ooman/p/11080694.html

matlab代码

function [test_data_sg] = Maxtrix_sg(test_data, M, N)
% test_data 待滤波数据
% M 滑动窗口大小 = 2 * M + 1
% N 拟合多项式最高阶数 = N - 1
% test_data_sg 滤波后的数据
test_data_sg = test_data;
[~, cols] = size(test_data);

if(cols <= 2 * M + 1)
    msgbox('待滤波数据长度太短', 'Warning:')
end

if(2 * M + 1 <= N)
    msgbox('滑动窗口比拟合多项式还小, 无法求解', 'Warning:')
end

move_box = zeros(2*M+1, N);
for i = -M : M
    for j = 0 : N - 1
        move_box(i+M+1, j+1) = i ^ j;
    end
end
I = ones(1, N);
I = diag(I);
move_box_move = move_box' * move_box;
mox_box_coef = (move_box * (I / move_box_move)) * move_box';

for i = M+1 : cols-M
    
    result = mox_box_coef * test_data(2, i-M : i+M)';
    test_data_sg(2, i) = result(M+1, 1);
    
end

matlab 自编写函数Maxtrix_sg(test_data, 7, 5) 与 sgolayfilt(test_data(2, :), 5, 15) 测试结果对比：

对比结果：自写Maxtrix_sg函数除了首尾端各M个数据没有处理之外，其余处理结果与matlab sgolayfilt处理结果一致。

opencv + C++代码

#include <iostream>
#include <opencv2/opencv.hpp>
#include <vector>
#include <cmath>

using namespace cv;
using namespace std;

// 矩阵X的广义逆
Mat Generalized_inverse(Mat &X)
{
	Mat Xt = X.t();

	Mat Xt_X = Xt * X;

	Mat inverse_Xt_X;
	invert(Xt_X, inverse_Xt_X, DECOMP_LU);		//方阵LU分解,求逆矩阵

	return inverse_Xt_X * Xt;
}

vector<vector<float>> opencv_sg(vector<vector<float>> &test_data, int M, int N)
{
	// SG平滑矩阵算子	
	Mat move_box(2 * M + 1, N, CV_32F);
	for (int i = -M; i <= M; ++i)
	{
		for (int j = 0; j < N; ++j)
		{
			move_box.at<float>(i + M, j) = pow(i, j);
		}
	}
	Mat gen_inv = Generalized_inverse(move_box);	// 矩阵move_box 的广义逆
	
	vector<vector<float>> test_data_clone = test_data;

	if (test_data[0].size() <= 2 * M + 1)
	{
		qDebug() << "the length of the input_data is too short!";
	}

	if (2 * M + 1 <= N)
	{
		qDebug() << "box size is too small!";
	}

	for (int i = M; i < test_data[0].size() - M; ++i)
	{
		Mat y_box_value(2 * M + 1, 1, CV_32F);
		for (int j = -M; j <= M; ++j)
		{
			y_box_value.at<float>(j + M, 0) = test_data[1][i + j];
		}

		Mat ls_sol = gen_inv * y_box_value;					// 最小二乘解
		Mat result = move_box * ls_sol;						// sg滤波值
		test_data_clone[1][i] = result.at<float>(M, 0);
	}

	return test_data_clone;
}

opencv + C++ 代码测试结果：

测试结果与matlab测试结果一致(注：C++下标从0开始，matlab下标从1开始。故这里从349开始截图)

完工~