R/internal_psdsPlot_s3.R

Defines functions image.psds

Documented in image.psds

#' Plot method for psds class
#' @description This function produces a spectogram from a psds object.
#' @export
#' @param x an object of class psds
#' @param func statistics to be applied to \code{x}, e.g., \code{median} or \code{mean}
#' @param logZ logical value (default is TRUE) to determine if the PSD should be on log scale
#' @param zoomFreq interval which specifies the frequency axis to be plotted, for instance, c(0.2, 0.7) leaves out the 20\% and 30\% of the low and high frequencies, respectively.  The default value is c(0,1), i.e., all the frequencies
#' @param fs sampling frequency
#' @param ... other graphical parameters from the image.plot function
#' @return plot of the estimate of the (log) spectrum
#' @seealso \link{gibbs_pspline}
#' @importFrom fields image.plot
#'
#' @examples
#' \dontrun{
#'
#' set.seed(1)
#'
#' # Generate AR(1) data with rho = 0.9
#' n = 128
#' data = arima.sim(n, model = list(ar = 0.9));
#' data = data - mean(data);
#'
#' # Spectrum estimate via p-splines
#' spec = spec_pspline(data, l=30, p=20, Ntotal1=3000, burnin1=1000, thin1=10,
#'                     Ntotal=3000, burnin=1000, thin=10, k=30);
#'
#' image(spec) # Plot log PSD (see documentation of iplot.psds)
#' }
image.psds = function(d, time = NULL, func = 'median', logZ = TRUE, zoomFreq = c(0,1), fs = 16384, ...) {  # Plot method for "psds" class

  if((zoomFreq[1]<0) || (zoomFreq[2]>1) || (zoomFreq[1] > zoomFreq[2])){
    stop("zoomFreq must be a vector c(a,b) with values 0 <= a < b <= 1")
  }

  N   = length(d$psds); # number of fpsd.sample elements
  po  = d$info$p/100;
  pc  = 1-po;
  aux = NULL;

  for(i in 1:N){
    aux = cbind(aux, apply(d$psds[[i]], 1, get(func))); # MEDIAN
  }

  aux = sqrt(aux);

  if(logZ == TRUE){

    aux = log10(aux);

  }

  i <- 1:length(d$psds);
  x <- round(d$info$l / fs * (i-1) * pc * 1000, 2); # number of intervals

  if(!is.null(time)){

    if(length(time)==2){

      x <- seq(from= time[1], to=time[2], length = length(x));

    }else{

      x = seq(from= time[1], to=time[length(time)], length = length(x));

    }
  }

  #y <- seq(0,1,length=dim(aux)[1]) # number of frequencies
  y <- seq(1,fs/2,length=dim(aux)[1]);

  if(all(zoomFreq == c(0,1))){

    z <- t(aux);
    z <- z[, -c(1, dim(z)[2])]; #deleting first and last row in image
    y <- y[-c(1, dim(z)[2])]; # adjustment

  }else{

    index    = round(dim(aux)[1] * zoomFreq);
    index    = seq(index[1] + 1, index[2]);

    y <- y[index]
    z <- t(aux[index, ])

  }

  fields::image.plot(x, y, z, xlab = "Time [ms]", ylab = "Frequency [Hz]",
                     main = paste(func), ...); # yaxt='n'

}
pmat747/psplinePsd documentation built on July 7, 2023, 9:06 p.m.