Simplifying sim functions

1 files changed, 34 insertions, 54 deletions

diff --git a/R/sim.R b/R/sim.R
index d849c15..8426722 100644
--- a/R/sim.R
+++ b/R/sim.R
@@ -6,12 +6,13 @@ sim.default <- function(model,input,sigma=0,seed=NULL){
   print("The sim method is not developed for the current class of the object")
 }
 
-#' Simulate from an ARX Model
+#' Simulate from a Polynomial Model
 #' 
-#' Simulate the response of an ARX system, given the input
+#' Simulate the response of a system system governed by a polynomial model
+#' , given the input
 #' 
-#' @param model an object of class \code{arx} containing the coefficients
-#' @param input a vector/matrix/idframe containing the input
+#' @param model an object of class \code{idpoly} containing the coefficients
+#' @param input a vector/matrix containing the input
 #' @param sigma standard deviation of the innovations (Default= \code{0})
 #' @param seed integer indicating the seed value of the random number generator
 #' 
@@ -23,15 +24,40 @@ sim.default <- function(model,input,sigma=0,seed=NULL){
 #' include support for MIMO systems
 #' 
 #' @seealso 
-#' \code{\link{arx}} for defining ARX models
+#' \code{\link{idpoly}} for defining polynomial models
 #' 
 #' @examples
+#' # ARX Model
 #' u <- rnorm(200,sd=1)
-#' model <- arx(A=c(1,-1.5,0.7),B=c(0.8,-0.25),ioDelay=1)
+#' model <- idpoly(A=c(1,-1.5,0.7),B=c(0.8,-0.25),ioDelay=1)
 #' y <- sim(model,u,sigma=0.1)
-#'  
+#' 
 #' @export
-sim.arx <- function(model,input,sigma=0,seed=NULL){
+sim.idpoly <- function(model,input,sigma=0,seed=NULL){
+  if(model$type=="arx"){
+    sim_arx(model,input,sigma,seed)
+  } else{
+    require(signal);require(polynom)
+    
+    n <- length(input)[1]
+    if(!is.null(seed)) set.seed(seed)
+    ek <- rnorm(n,sd=sigma)
+    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 <- Arma(b=B,a=den2)
+    ufk <- signal::filter(filt2,input)
+    
+    yk <- as.numeric(ufk) + as.numeric(vk)
+    
+    return(as.numeric(yk)) 
+  }
+}
+
+sim_arx <- function(model,input,sigma=0,seed=NULL){
   na <- length(model$A) - 1; nk <- model$ioDelay; 
   nb <- length(model$B) - 1; nb1 <- nb+nk
   n <- max(na,nb1)
@@ -60,50 +86,4 @@ sim.arx <- function(model,input,sigma=0,seed=NULL){
     y[i] <- reg%*%coef + ek[i-n]
   }
   return(y[n+1:length(uk)])
-}
-
-#' Simulate from a Polynomial Model
-#' 
-#' Simulate the response of a system system governed by a polynomial model
-#' , given the input
-#' 
-#' @param model an object of class \code{idpoly} containing the coefficients
-#' @param input a vector/matrix containing the input
-#' @param sigma standard deviation of the innovations (Default= \code{0})
-#' @param seed integer indicating the seed value of the random number generator
-#' 
-#' @return
-#' a vector containing the output
-#' 
-#' @details
-#' The routine is currently built only for SISO systems. Future Versions will
-#' include support for MIMO systems
-#' 
-#' @seealso 
-#' \code{\link{idpoly}} for defining polynomial models
-#' 
-#' @examples
-#' u <- rnorm(200,sd=1)
-#' model <- idpoly(A=c(1,-1.5,0.7),B=c(0.8,-0.25),C=1,D=1,F1=1,ioDelay=1)
-#' y <- sim(model,u,sigma=0.1)
-#'  
-#' @export
-sim.idpoly <- function(model,input,sigma=0,seed=NULL){
-  require(signal);require(polynom)
-  
-  n <- length(input)[1]
-  if(!is.null(seed)) set.seed(seed)
-  ek <- rnorm(n,sd=sigma)
-  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 <- Arma(b=B,a=den2)
-  ufk <- signal::filter(filt2,input)
-  
-  yk <- as.numeric(ufk) + as.numeric(vk)
-  
-  return(as.numeric(yk))
 }
\ No newline at end of file