信号分析——S-G卷积平滑(Java/Matlab)

Savitzky-Golay滤波，即S-G滤波平滑，利用多项式进行数据平滑，基于最小二乘法，能够保留分析信号中的有用信息，消除随机噪声，在信号图谱中最直接的结果就是将图谱的“毛刺”去掉，整个图谱更加平滑。

平滑窗口宽度：2m+1，即为该窗口原始数据点个数n（n=2m+1），一般要求窗口宽度为奇数，从测量原始数据中选取一子集作为窗口，而非整个数据集。

假设窗口内的原始数据点可以用一个k-1次多项式进行拟合，即

其中：i=（-m，-m+1，···，0,1，m-1，m）。因此，该窗口内的n个原始数据点均可以得到1个上述的多项式，n个这样的多项式构2成了k元线性方程组，需要求解其中的k个拟合参数aj（j=0,1,2，···，k-1），一般选取的滤波窗口宽度n应大于或至少等于k。当n=k时，可用线性代数方法求解拟合参数；若n>k，则采用最小二乘法求解。

将n个上述多项式并列，可得到下列矩阵运算：

简化为下列超定方程组：

A的最小二乘解为：

求解该方程组的最小二乘解方法参考：http://blog.csdn.net/acdreamers/article/details/44662633

则Y的滤波值计算公式如下：

其中：，B为滤波系数矩阵，由且仅由X矩阵决定，B矩阵为(2m+1)×(2m+1)阶矩阵，根据系数矩阵即可获得S-G平滑拟合方程。

下面分别通过Matlab、Java实现该方法：

Matlab：

function [ filter_matrix ] = S_Gtest(  n,k )
%S-G 平滑一般方法-n点k-1次平滑 
%Input: -n，滤波带宽，即原始数据点个数，一般要求为奇数
%       -k，拟合多项式最高次幂+1
%Output:-S-G滤波系数矩阵
m=(n-1)/2;
X=[];
for i=0:(n-1)
    for j=0:(k-1)
        X(i+1,j+1)=power(i-m,j);
    end
end
filter_matrix=X*inv(X'*X)*X';
end

Java：

public static RealMatrix getXMatrix(int n,int k) {
   int step = (n - 1) / 2;
   RealMatrix realMatrix = MatrixUtils.createRealMatrix(n, k );
   for (int i = 0; i < realMatrix.getRowDimension(); i++) {
        for (int j = 0; j < k ; j++) {
            realMatrix.setEntry(i, j, Math.pow(i - step, j));
        }
   }
   return realMatrix;
}
public static RealMatrix filterCoeff(int n,int k) {
   RealMatrix realMatrix = getXMatrix(n,k);
   RealMatrix B = realMatrix.multiply(inverseMatrix(realMatrix.transpose().
   multiply(realMatrix))).multiply(realMatrix.transpose());
   return B;
}

——参考《化学计量学方法及MATLAB实现》史永刚等编著