将 Savitzky-Golay 函数从 Matlab 转换为 R

Question

我正在尝试将这个 Savitzky-Golay 函数从 Matlab 转换为 R。但是它在 R 中不起作用。如何使 R 函数起作用？

x 示例可以在这里下载：https://drive.google.com/open?id=0B5AOSYBy_josMUgtRi1wLW4tZEE

函数输入是：

• 要处理的 x (m x n) 个数据
• 宽度 (1 x 1) 点数
• 阶 (1 x 1) 多项式阶
• deriv (1 x 1) 导数顺序

函数输出为：

• xsg (m x n) 处理数据

附上来自 Matlab 和 R 的示例：

---MATLAB---

function [xsg]= savgol(x,width,order,deriv)
[m,n]=size(x);
w=max( 3, 1+2*round((width-1)/2) );
o=min([max(0,round(order)),5,w-1]);
d=min(max(0,round(deriv)),o);
p=(w-1)/2;
xc=((-p:p)'*ones(1,1+o)).^(ones(size(1:w))'*(0:o));
we=xc\eye(w);
b=prod(ones(d,1)*[1:o+1-d]+[0:d-1]'*ones(1,o+1-d,1),1);
di=spdiags(ones(n,1)*we(d+1,:)*b(1),p:-1:-p,n,n);
w1=diag(b)*we(d+1:o+1,:);
di(1:w,1:p+1)=[xc(1:p+1,1:1+o-d)*w1]'; 
di(n-w+1:n,n-p:n)=[xc(p+1:w,1:1+o-d)*w1]';
xsg=x*di;

plot(xsg')

---R---

savgol=function(x,width,order,deriv)
{
  m=nrow(x)
  n=ncol(x)
  w=max(3,1+2*round((width-1)/2) )
  o=min(c(max(0,round(order)),5,w-1))
  d=min(max(0,round(deriv)),o)
  p=(w-1)/2
  xc=((-p:p)%*%matrix(1,1,1+o))^(t(matrix(1,1,w))%*%(0:o))
  we=qr.solve(xc,diag(w))
  b=apply((matrix(1,d,1)%*%matrix(1:(o+1-d),1,(o+1-d))+t(matrix(0:(d-1),1,(d)))%*%matrix(1,1,o+1-d)),2,prod)
  gg=matrix(1,n,1)%*%we[(d+1),]*b[1]
  library(Matrix)
  di=sparseMatrix(i=1:n,j=1:n,x=gg)
  w1=diag(b,nrow=length(b))%*%we[(d+1):(o+1),]
  di[1:w,1:(p+1)]=t(xc[1:(p+1),1:(1+o-d)]%*%w1)
  di[(n-w+1):n,(n-p):n]=t(xc[(p+1):w,1:(1+o-d)]%*%w1)
  xsg=x%*%di
    }

matplot(t(xsg),type='l')

获取到 xsg 的图：

Answer 1

可能不是一个完整的答案，但看看圆形：

---matlab--- Y = round(X) 将 X 的每个元素四舍五入为最接近的整数。在平局的情况下，元素的小数部分正好是 0.5，round 函数会从零四舍五入到更大的整数。

---R--- 请注意，对于 5 的四舍五入，预计将使用 IEC 60559 标准，“转到偶数位”。

round(0.5) 将导致 R 为 0，matlab 为 1

Answer 2

我在 R 中通过以下方式解决了它：

savgol=function(x,width,order,deriv)
{
  ##insert spdiags function
  spdiags <-function (arg1,arg2,arg3,arg4){
    B <- arg1 
    if (is.matrix(arg2))
      d <- matrix(arg2,dim(arg2)[1]*dim(arg2)[2],1)
    else
      d <- arg2
    p <- length(d) 
    A <- sparseMatrix(i = 1:arg3, j = 1:arg3, x = 0, dims= c(arg3,arg4))
    m <- dim(A)[1] 
    n <- dim(A)[2] 
    len<-matrix(0,p+1,1)
    for (k in 1:p)
      len[k+1] <- len[k]+length(max(1,1-d[k]):min(m,n-d[k])) 
    a <- matrix(0, len[p+1],3) 
    for (k in 1:p)
    {
      i <- t(max(1,1-d[k]):min(m,n-d[k])) 
      a[(len[k]+1):len[k+1],] <- c(i, i+d[k], B[(i+(m>=n)*d[k]),k]) 
    }
    res1 <- sparseMatrix(i = a[,1], j = a[,2], x = a[,3], dims = c(m,n)) 
    return (res1)
  }
  ##end spdiags function

  m=nrow(x)
  n=ncol(x)
  w=max(3,1+2*round((width-1)/2) )
  o=min(c(max(0,round(order)),5,w-1))
  d=min(max(0,round(deriv)),o)
  p=(w-1)/2
  xc=((-p:p)%*%matrix(1,1,1+o))^(t(matrix(1,1,w))%*%(0:o))
  we=qr.solve(xc,diag(w))
    options(warn=-1)
  b=apply((matrix(1,d,1)%*%matrix(1:(o+1-d),1,(o+1-d))+t(matrix(0:(d-1),1,(d)))%*%matrix(1,1,o+1-d)),2,prod)
  gg=matrix(1,n,1)%*%we[(d+1),]*b[1]
    library(Matrix)
  di=spdiags(gg,p:(-p),n,n)
    options(warn=0)
  w1=diag(b,nrow=length(b))%*%we[(d+1):(o+1),]
  di[1:w,1:(p+1)]=t(xc[1:(p+1),1:(1+o-d)]%*%w1)
  di[(n-w+1):n,(n-p):n]=t(xc[(p+1):w,1:(1+o-d)]%*%w1)
  result=x%*%di
}