自适应过滤法代码

music813

#自适应过滤法
weight<-function(x,p,mse)   
{
  biaozhun<-function(x,p)  ##标准化函数
  {
    n<-length(x)
    a=b=rep(NA,n)
    y=matrix(data=NA,nrow=n,ncol=p+1)
    for(t in (p+1):n)
    {
      a[t]=sum(x[(t-1)t-p)]^2)
      b[t]=(a[t])^(1/2)
      for(j in 0:p)
      {
        y[t,j+1]=as.numeric(x[t-j])/b[t]
      }
    }
    out<-data.frame(x,a,b,y)
    return(out)
  }
  y=biaozhun(x,2)[44+p)]  ##提取标准化后值
  mse1<-function(x,y)     ##计算mse函数
  {      
    e=x-y
    mse=mean(e^2,na.rm=TRUE)
    return(mse)
  }
  fai=rep(1/p,p)        ##权数的初始值
  k=min(1/p,ceiling(1/sum(sort(-y[,1])[1:p]^2)))   
  #sort默认是升序排列
  MSE=10     
  #初始均方误差，可以任意给值，适当大些
  n=0  
  mm=list()  #用来存放迭代的结果
  proc=e=matrix(NA,nrow=nrow(y))  
  #用来存放每一步迭代中的预测值和误差   
  while(MSE >= mse)
  {
    for(t in (p+1):nrow(y))
    {
      proc[t]=sum(fai[1:p]*y[t,2p+1)])
      e[t]=y[t,1]-proc[t]
      for(j in 1:p)
      {
        fai[j]=fai[j]+2*k*e[t]*y[t,j+1]
      }
    }
    n=n+1  #用来表示迭代次数
    MSE=mse1(y[,1],proc)  #更新均方误差
    mm[[n]]=data.frame(y,proc,e,rep(MSE,nrow(y)))
    names(mm)=paste(rep("第",n),1:n,rep("次迭代",n))
    #拼接字符串
  }
  pred=sum(fai[1:p]*x[length(x)length(x)-p+1)])  
  #计算预测值
  m=list(mm,fai,k,n,pred)  #所有结果放在列表m里
  names(m)=c('迭代摘要','fai参数值','k值','迭代次数','下期预测值')
  return(m)
}