本文主要是介绍基于选权迭代法的Savitzky_Golay平滑算法，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
中心思想
在基于最小二乘原理对多项式进行拟合时，引入观测值的权值，有明显偏差的信号给予较小的权值，不含明显偏差的信号给予较大的权值，此时最小而成的优化函数有min(VTV) 变为min(VTPV)，其中P是根据某次平差的改正数V基于的权值，在下面的代码中P=1/|V|，迭代求解，知道两次结果的差小于给定的阈值
import numpy as npclass Robust(object):def __init__(self,X:np.array,Y:np.array,r:float) -> None:super().__init__()self.X=Xself.Y=Yn,t=X.shapeP=np.diag(np.ones(n))self.P=Pself.n=nself.t=t# r迭代阈值self.r=r@staticmethoddef ls(X:np.array,Y:np.array,P:np.array):N=X.T@P@XU=X.T@P@Ya=np.linalg.inv(N)@Ureturn adef iterate(self):a0=Nonen=0while True:a1=Robust.ls(self.X,self.Y,self.P)if n and np.max(np.abs(a1-a0))<self.r:n+=1breakelse:a0=a1v=self.X@a1-self.Yp=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+=1return a1
class SG():def __init__(self,x,y,r,k) -> None:# n 多项式模型的次数self.x=xself.y=yself.k=kself.r=rdef sm(self):n=self.x.shape[0]if n<(2*self.r+1):self.smy=self.yreturnif (2*self.r+1)<(self.k+1):self.smy=self.yreturnsmy=self.y[:]for i in range(self.r,n-self.r):start=i-self.rend=i+self.rxi=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=smyreturn
if __name__=="__main__":a=1b=2c=3x=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.5P=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)
这篇关于基于选权迭代法的Savitzky_Golay平滑算法的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
