code to calculate the noise spectrum

2 files changed, 24 insertions, 23 deletions

diff --git a/R/nonparam.R b/R/nonparam.R
index 540207a..d1e56a9 100644
--- a/R/nonparam.R
+++ b/R/nonparam.R
@@ -162,7 +162,8 @@ step <- function(model){
 #' (Default: \code{seq(1,128)/128*pi/Ts})
 #' 
 #' @return
-#' an \code{idfrd} object containing the estimated frequency response
+#' an \code{idfrd} object containing the estimated frequency response 
+#' and the noise spectrum
 #' 
 #' @references
 #' Arun K. Tangirala (2015), \emph{Principles of System Identification: 
@@ -181,8 +182,10 @@ spa <- function(x,winsize=NULL,freq=NULL){
   if(is.null(winsize)) winsize <- min(N/10,30)
   if(is.null(freq)) freq <- (1:128)/128*pi/deltat(x)
   M <- winsize
+  
   Ryu <- mult_ccf(x$out,x$input,lag.max = M)
   Ruu <- mult_ccf(x$input,x$input,lag.max=M)
+  Ryy <- mult_ccf(x$out,x$out,lag.max = M)
   
   cov2spec <- function(omega,R,M){
     seq1 <- exp(-1i*(-M:M)*omega)
@@ -190,14 +193,19 @@ spa <- function(x,winsize=NULL,freq=NULL){
   }
   
   G <- array(0,c(nout,nin,length(freq)))
+  spec <- array(0,c(nout,nout,length(freq)))
   for(i in 1:nout){
+    phi_y <- sapply(freq,cov2spec,Ryy[i,i,],M)
+    temp <- phi_y
     for(j in 1:nin){
-      num <- sapply(freq,cov2spec,Ryu[i,j,],M)
-      den <- sapply(freq,cov2spec,Ruu[,j,],M)
-      G[i,j,] <- num/den
+      phi_yu <- sapply(freq,cov2spec,Ryu[i,j,],M)
+      phi_u <- sapply(freq,cov2spec,Ruu[j,j,],M)
+      G[i,j,] <- phi_yu/phi_u
+      temp <- phi_y - phi_yu*Conj(phi_yu)/phi_u
     }
+    spec[i,i,] <- temp
   }
-  out <- idfrd(G,matrix(freq),deltat(x))
+  out <- idfrd(G,matrix(freq),deltat(x),spec)
   return(out)
 }
 
diff --git a/man/spa.Rd b/man/spa.Rd
index 3979b7a..53dffd8 100644
--- a/man/spa.Rd
+++ b/man/spa.Rd
@@ -4,28 +4,24 @@
 \alias{spa}
 \title{Estimate frequency response}
 \usage{
-spa(data, npad = 255)
+spa(x, winsize = NULL, freq = NULL)
 }
 \arguments{
-\item{data}{an \code{idframe} object}
+\item{x}{an \code{idframe} object}
 
-\item{npad}{an integer representing the total length of each time series 
-to analyze after padding with zeros. This argument allows the user to 
-control the spectral resolution of the SDF estimates: the normalized 
-frequency interval is deltaf=1/npad. (Default: 255)}
+\item{freq}{frequency points at which the response is evaluated
+(Default: \code{seq(1,128)/128*pi/Ts})}
+
+\item{W}{lag size of the Hanning window (Default: \code{min
+(length(x)/10,30)})}
 }
 \value{
-an \code{idfrd} object containing the estimated frequency response
+an \code{idfrd} object containing the estimated frequency response 
+and the noise spectrum
 }
 \description{
-Estimates Frequency Response with fixed frequency resolution using 
-spectral analysis
-}
-\details{
-The function calls the \code{SDF} function in the \code{sapa} package to
-compute the cross-spectral densities. The method used is \strong{Welch's 
-Overlapped Segment Averaging} with a normalized \strong{Hanning} window.
-The overlap used is 50%.
+Estimates frequency response and noise spectrum from data with 
+fixed resolution using spectral analysis
 }
 \examples{
 data(frf)
@@ -36,7 +32,4 @@ frf <- spa(data)
 Arun K. Tangirala (2015), \emph{Principles of System Identification: 
 Theory and Practice}, CRC Press, Boca Raton. Sections 16.5 and 20.4
 }
-\seealso{
-\code{\link[sapa]{SDF}}
-}