diff options
author | Suraj Yerramilli | 2016-02-27 14:39:29 +0530 |
---|---|---|
committer | Suraj Yerramilli | 2016-02-27 14:40:05 +0530 |
commit | ae603683420f2782407b798a66de700717aad1af (patch) | |
tree | a409b8479ece7bbdbab6f70fdd8ea80ba16931f5 /R/sim.R | |
parent | 36c9ca307cc178643f1aa6233e89a838bb3526b2 (diff) | |
download | SysID-R-code-ae603683420f2782407b798a66de700717aad1af.tar.gz SysID-R-code-ae603683420f2782407b798a66de700717aad1af.tar.bz2 SysID-R-code-ae603683420f2782407b798a66de700717aad1af.zip |
Implemented the addNoise attribute in the sim routine
Diffstat (limited to 'R/sim.R')
-rw-r--r-- | R/sim.R | 89 |
1 files changed, 43 insertions, 46 deletions
@@ -23,64 +23,61 @@ #' y <- sim(model,u,sigma=0.1) #' #' @export -sim <- function(model,input,innov=NULL,sigma=0,seed=NULL) UseMethod("sim") +sim <- function(model,input,addNoise = F,innov=NULL,seed=NULL) UseMethod("sim") #' @export -sim.default <- function(model,input,sigma=0,seed=NULL){ +sim.default <- function(model,input,addNoise = F,innov=NULL,seed=NULL){ print("The sim method is not developed for the current class of the object") } #' @import signal polynom #' @export -sim.idpoly <- function(model,input,innov=NULL,sigma=0,seed=NULL){ - n <- length(input)[1] - if(!is.null(innov)){ - ek <- innov - } else{ - if(!is.null(seed)) set.seed(seed) - ek <- rnorm(n,sd=sigma) - } +sim.idpoly <- function(model,input,addNoise = F,innov=NULL,seed=NULL){ + B <- c(rep(0,model$ioDelay),model$B) + Gden <- as.numeric(polynomial(model$A)*polynomial(model$F1)) + G <- signal::Arma(b=B,a=Gden) + ufk <- signal::filter(G,input) + yk <- as.numeric(ufk) - if(model$type=="arx"){ - sim_arx(model,input,ek) - } else{ - + if(addNoise){ + n <- ifelse(is.numeric(input),length(input),dim(input)[1]) + if(!is.null(innov)){ + ek <- innov + } else{ + if(!is.null(seed)) set.seed(seed) + ek <- rnorm(n,sd=model$noiseVar) + } den1 <- as.numeric(polynomial(model$A)*polynomial(model$D)) filt1 <- Arma(b=model$C,a=den1) vk <- signal::filter(filt1,ek) - B <- c(rep(0,model$ioDelay),model$B) - den2 <- as.numeric(polynomial(model$A)*polynomial(model$F1)) - filt2 <- signal::Arma(b=B,a=den2) - ufk <- signal::filter(filt2,input) - - yk <- as.numeric(ufk) + as.numeric(vk) - - return(as.numeric(yk)) + yk <- yk + as.numeric(vk) } + + return(yk) } -sim_arx <- function(model,input,ek){ - na <- length(model$A) - 1; nk <- model$ioDelay; - nb <- length(model$B) - 1; nb1 <- nb+nk - n <- max(na,nb1) - coef <- matrix(c(model$A[-1],model$B),nrow=na+(nb+1)) - - if(class(input)=="idframe"){ - uk <- input$input[,1,drop=T] - } else if(class(input) %in% c("matrix","data.frame")){ - uk <- input[,1,drop=T] - } else if(is.numeric(input)){ - uk <- input - } - - y <- rep(0,length(input)+n) - u <- c(rep(0,n),uk) - - for(i in n+1:length(uk)){ - if(nk==0) v <- u[i-0:nb] else v <- u[i-nk:nb1] - reg <- matrix(c(-(y[i-1:na]),v),ncol=na+(nb+1)) - y[i] <- reg%*%coef + ek[i] - } - return(y[n+1:length(uk)]) -}
\ No newline at end of file +# sim_arx <- function(model,input,ek){ +# na <- length(model$A) - 1; nk <- model$ioDelay; +# nb <- length(model$B) - 1; nb1 <- nb+nk +# n <- max(na,nb1) +# coef <- matrix(c(model$A[-1],model$B),nrow=na+(nb+1)) +# +# if(class(input)=="idframe"){ +# uk <- input$input[,1,drop=T] +# } else if(class(input) %in% c("matrix","data.frame")){ +# uk <- input[,1,drop=T] +# } else if(is.numeric(input)){ +# uk <- input +# } +# +# y <- rep(0,length(input)+n) +# u <- c(rep(0,n),uk) +# +# for(i in n+1:length(uk)){ +# if(nk==0) v <- u[i-0:nb] else v <- u[i-nk:nb1] +# reg <- matrix(c(-(y[i-1:na]),v),ncol=na+(nb+1)) +# y[i] <- reg%*%coef + ek[i] +# } +# return(y[n+1:length(uk)]) +# }
\ No newline at end of file |