基于选权迭代法的Savitzky_Golay平滑算法

中心思想

在基于最小二乘原理对多项式进行拟合时，引入观测值的权值，有明显偏差的信号给予较小的权值，不含明显偏差的信号给予较大的权值，此时最小而成的优化函数有min(V^TV) 变为min(V^TPV)，其中P是根据某次平差的改正数V基于的权值，在下面的代码中P=1/|V|，迭代求解，知道两次结果的差小于给定的阈值

import numpy as np

class Robust(object):
    def __init__(self,X:np.array,Y:np.array,r:float) -> None:
        super().__init__()
        self.X=X
        self.Y=Y
        n,t=X.shape
        P=np.diag(np.ones(n))
        self.P=P
        self.n=n
        self.t=t
        # r迭代阈值
        self.r=r
        
    @staticmethod
    def ls(X:np.array,Y:np.array,P:np.array):
        N=X.T@P@X
        U=X.T@P@Y
        a=np.linalg.inv(N)@U
        return a
    def iterate(self):
        a0=None
        n=0
        while True:
            a1=Robust.ls(self.X,self.Y,self.P)
            if n and np.max(np.abs(a1-a0))<self.r:
                n+=1
                break
            else:
                a0=a1
                v=self.X@a1-self.Y
                p=np.ones(self.n)
                for i in range(self.n):
                    p[i]=1/(np.abs(v[i])+1e-13)
                self.P=np.diag(p)
                n+=1
        return a1
class SG():
    def __init__(self,x,y,r,k) -> None:
        # n 多项式模型的次数
        self.x=x
        self.y=y
        self.k=k
        self.r=r
    def sm(self):
        n=self.x.shape[0]
        if n<(2*self.r+1):
            self.smy=self.y
            return
        if (2*self.r+1)<(self.k+1):
            self.smy=self.y
            return
        smy=self.y[:]
        for i in range(self.r,n-self.r):
            start=i-self.r
            end=i+self.r
            xi=self.x[start:(end+1)]
            yi=self.y[start:(end+1)]
            Xi=list()
            for x in xi:
                xr=list()
                for j in range(self.k+1):
                    xr.append(pow(x,j))
                Xi.append(xr)
            Xi=np.array(Xi)
            yi=np.array(yi)
            robust=Robust(Xi,yi,0.001)
            a=robust.iterate()
            smy[i]=np.dot(a,Xi[self.r])
        self.smy=smy
        return
if __name__=="__main__":
    a=1
    b=2
    c=3
    x=list()
    X=list()
    y=list()
    for i in range(1,6):
        X.append([i**2,i,1])
        y.append(a*pow(i,2)+b*i+c+(np.random.rand()-0.5)/10)
    X=np.array(X)
    y=np.array(y)
    y[0]+=0.5
    P=np.diag(np.ones(5))
    a1=Robust.ls(X,y,P)
    print('最小二乘法:')
    print(a1)
    robust=Robust(X,y,0.1)
    a2=robust.iterate()
    print('选权迭代法:')
    print(a2)