Перевести функцию Savitzky-Golay от Matlab до R

Question

Я пытаюсь перевести эту функцию Savitzky-Golay из Matlab в R. Однако это не работает в R. Как заставить функцию R работать?

х пример можно скачать здесь: https://drive.google.com/open?id=0B5AOSYBy_josMUgtRi1wLW4tZEE

Функциональные входы: данные

• х (MXN) для обработки
• ширину (1 х 1) число точек
• порядка (1 x 1) полиномиальный порядок
• производный (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 = круглый (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 
}