diff options
Diffstat (limited to 'R/nonparam.R')
| -rw-r--r-- | R/nonparam.R | 88 |
1 files changed, 45 insertions, 43 deletions
diff --git a/R/nonparam.R b/R/nonparam.R index 3e9cf0b..540207a 100644 --- a/R/nonparam.R +++ b/R/nonparam.R @@ -152,20 +152,14 @@ step <- function(model){ #' Estimate frequency response #' -#' Estimates Frequency Response with fixed frequency resolution using -#' spectral analysis +#' Estimates frequency response and noise spectrum from data with +#' fixed resolution using spectral analysis #' -#' @param data an \code{idframe} object -#' @param 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) -#' -#' @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%. +#' @param x an \code{idframe} object +#' @param W lag size of the Hanning window (Default: \code{min +#' (length(x)/10,30)}) +#' @param freq frequency points at which the response is evaluated +#' (Default: \code{seq(1,128)/128*pi/Ts}) #' #' @return #' an \code{idfrd} object containing the estimated frequency response @@ -174,31 +168,55 @@ step <- function(model){ #' 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}} -#' #' @examples #' data(frf) #' frf <- spa(data) #' #' @import signal #' @export -spa <- function(data,npad=255){ - temp <- cbind(data$output,data$input) +spa <- function(x,winsize=NULL,freq=NULL){ + N <- dim(x$out)[1] + nout <- nOutputSeries(x); nin <- nInputSeries(x) - # Non-parametric Estimation of Spectral Densities - - # WOSA and Hanning window - gamma <- sapa::SDF(temp,method="wosa",sampling.interval = - deltat(data),npad=npad) - freq <- matrix(attributes(gamma)$frequency*2*pi) - resp <- Conj(gamma[,2])/Mod(gamma[,3]) + 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) - # power-spectrum - spec <- gamma[,2] - resp*gamma[,3] - - out <- idfrd(resp,freq,deltat(data),spec) + cov2spec <- function(omega,R,M){ + seq1 <- exp(-1i*(-M:M)*omega) + sum(R*signal::hanning(2*M+1)*seq1) + } + + G <- array(0,c(nout,nin,length(freq))) + for(i in 1:nout){ + 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 + } + } + out <- idfrd(G,matrix(freq),deltat(x)) return(out) } +mult_ccf <- function(X,Y=NULL,lag.max=30){ + N <- dim(X)[1]; nx <- dim(X)[2] + ny <- ifelse(is.null(Y),nx,dim(Y)[2]) + + ccvfij <- function(i,j,lag=30) ccf(X[,i],Y[,j],plot=F,lag=lag, + type="covariance") + Xindex <- matrix(sapply(1:nx,rep,nx),ncol=1)[,1] + temp <- mapply(ccvfij,i=Xindex,j=rep(1:ny,ny), + MoreArgs = list(lag=lag.max)) + + ccfextract <- function(i,l) l[,i]$acf + temp2 <- t(sapply(1:(nx*ny),ccfextract,l=temp)) + dim(temp2) <- c(nx,ny,2*lag.max+1) + return(temp2) +} + #' Estimate empirical transfer function #' #' Estimates the emperical transfer function from the data by taking the @@ -236,20 +254,4 @@ comdiv <- function(z1,z2){ phi1 <- Arg(z1); phi2 <- Arg(z2) complex(modulus=mag1/mag2,argument=signal::unwrap(phi1-phi2)) -} - -mult_ccf <- function(X,Y=NULL,lag.max=30){ - N <- dim(X)[1]; nx <- dim(X)[2] - ny <- ifelse(is.null(Y),nx,dim(Y)[2]) - - ccvfij <- function(i,j,lag=30) ccf(X[,i],Y[,j],plot=F,lag=lag, - type="covariance") - Xindex <- matrix(sapply(1:nx,rep,nx),ncol=1)[,1] - temp <- mapply(ccfij,i=Xindex,j=rep(1:ny,ny), - MoreArgs = list(lag=lag.max)) - - ccfextract <- function(i,l) l[,i]$acf - temp2 <- t(sapply(1:(nx*ny),ccfextract,l=temp)) - dim(temp2) <- c(nx,ny,2*lag.max+1) - return(temp2) }
\ No newline at end of file |