將Savitzky-Golay函數從Matlab轉換爲R

Question

我正在試圖將這個Savitzky-Golay函數從Matlab轉換爲R.但是它在R.中不起作用。如何使R函數工作？

X例如可以在這裏下載：https://drive.google.com/open?id=0B5AOSYBy_josMUgtRi1wLW4tZEE

功能輸入是：

• X（MXN）的數據來處理
• 寬度（1×1）個
• 順序的號碼（ 1 x 1）多項式order
• deriv（1 x 1）衍生訂單

功能的輸出是：

• XSG（MXN）的處理過的數據

Atacched是從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，則舍入函數從零到大的整數。

--- R --- 請注意，爲了舍入一個5，IEC 60559標準預計會被使用，'轉到偶數位'。

輪（0.5）將導致0上R和1上的matlab

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 
}