Traduire la fonction Savitzky-Golay de Matlab à R

Question

J'essaie de traduire cette fonction de Savitzky-Golay de Matlab à R. Cependant cela ne fonctionne pas dans R. Comment faire fonctionner la fonction R?

x exemple peut être téléchargé ici: https://drive.google.com/open?id=0B5AOSYBy_josMUgtRi1wLW4tZEE

entrées de fonction sont:

• x (mxn) de données pour traiter
• largeur (1 x 1) nombre de points
• commande (1 x 1) ordre polynomial
• dérivé (1 x 1) ordre dérivé

sortie de la fonction est la suivante:

• données traitées de XSG (mxn)

Atacched sont des exemples de Matlab et 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') 

parcelles Adresser XSG:

Answer 1

Probablement pas une réponse complète, mais regardez autour de:

--- --- Matlab Y = round (X) arrondit chaque élément de X à la l'entier le plus proche. Dans le cas d'un lien, où un élément a une partie fractionnaire d'exactement 0.5, la fonction arrondit de zéro à l'entier de plus grande amplitude.

--- R --- Notez que pour arrondir un 5, la norme CEI 60559 devrait être utilisée, 'aller au chiffre pair'.

ronde (0,5) se traduira par 0 à R et 1 sur Matlab

Answer 2

je l'ai résolu à R par:

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 
}