# R/windowed_scalogram.R In wavScalogram: Wavelet Scalogram Tools for Time Series Analysis

#### Documented in windowed_scalogram

#' @title Windowed scalograms of a signal
#'
#' @description This function computes the normalized windowed scalograms of a signal for
#' the scales given. It is computed using time windows with radius \code{windowrad}
#' centered at a vector of central times with increment of time \code{delta_t}. It is
#' important to note that the notion of scalogram here is analogous to the spectrum of the
#' Fourier transform. It gives the contribution of each scale to the total energy of the
#' signal. For each scale \eqn{s} and central time \eqn{tc}, it is defined as the square
#' root of the integral of the squared modulus of the wavelet transform w.r.t the time
#' variable \eqn{t}, i.e.
#'
#'
#' "Normalized" means that the windowed scalograms are divided by the square root of the
#' length of the respective time windows in order to be comparable between them.
#'
#'
#' @usage windowed_scalogram(signal,
#'                           dt = 1,
#'                           scales = NULL,
#'                           powerscales = TRUE,
#'                           delta_t = NULL,
#'                           wname = c("MORLET", "DOG", "PAUL", "HAAR", "HAAR2"),
#'                           wparam = NULL,
#'                           border_effects = c("BE", "INNER", "PER", "SYM"),
#'                           energy_density = TRUE,
#'                           makefigure = TRUE,
#'                           time_values = NULL,
#'                           figureperiod = TRUE,
#'                           xlab = "Time",
#'                           ylab = NULL,
#'                           main = "Windowed Scalogram",
#'                           zlim = NULL)
#'
#' @param signal A vector containing the signal whose windowed scalogram is wanted.
#' @param dt Numeric. The time step of the signal.
#' @param scales A vector containing the wavelet scales at wich the windowed scalograms
#' are computed. This can be either a vector with all the scales or, following Torrence
#' and Compo 1998, a vector of 3 elements with the minimum scale, the maximum scale and
#' the number of suboctaves per octave. In the first case, \code{powerscales} must be
#' \code{FALSE} if the given scales are not power 2 scales. In the second case,
#' \code{powerscales} must be \code{TRUE} in order to construct power 2 scales using a
#' base 2 logarithmic scale). If \code{scales} is NULL, they are automatically constructed.
#' @param powerscales Logical. It must be TRUE (default) only in these cases:
#' \itemize{
#' \item If \code{scales} are power 2 scales, i.e. they use a base 2 logarithmic scale.
#' \item If we want to construct power 2 scales automatically. In this case, \code{scales}
#' must be \code{NULL}.
#' \item If we want to construct power 2 scales from \code{scales}. In this case,
#' \code{length(scales)} must be 3.
#' }
#' Otherwise, it must be \code{FALSE}.
#' @param windowrad Integer. Time radius for the windows, measured in \code{dt}. By
#' default, it is set to \eqn{ceiling(length(signal) / 20)}.
#' @param delta_t Integer. Increment of time for the construction of windows central
#' times, measured in \code{dt}. By default, it is set to
#' \eqn{ceiling(length(signal) / 256)}.
#' @param wname A string, equal to "MORLET", "DOG", "PAUL", "HAAR" or "HAAR2". The
#' difference between "HAAR" and "HAAR2" is that "HAAR2" is more accurate but slower.
#' @param wparam The corresponding nondimensional parameter for the wavelet function
#' (Morlet, DoG or Paul).
#' @param waverad Numeric. The radius of the wavelet used in the computations for the cone
#' of influence. If it is not specified, it is asumed to be \eqn{\sqrt{2}} for Morlet and DoG,
#' \eqn{1/\sqrt{2}} for Paul and 0.5 for Haar.
#' @param border_effects String, equal to "BE", "INNER", "PER" or "SYM", which indicates
#' how to manage the border effects which arise usually when a convolution is performed on
#' finite-lenght signals.
#' \itemize{
#' \item "BE": With border effects, padding time series with zeroes.
#' \item "INNER": Normalized inner scalogram with security margin adapted for each
#'     different scale. Although there are no border effects, it is shown as a regular COI
#'     the zone in which the length of \eqn{J(s)} (see Benítez et al. 2010) is smaller and
#'     it has to be normalized.
#' \item "PER": With border effects, using boundary wavelets (periodization of the
#' original time series).
#' \item "SYM": With border effects, using a symmetric catenation of the original time
#' series.
#' }
#' @param energy_density Logical. If TRUE (default), divide the scalograms by the square
#' root of the scales for convert them into energy density.
#' @param makefigure Logical. If TRUE (default), a figure with the scalograms is plotted.
#' @param time_values A numerical vector of length \code{length(signal)} containing custom
#' time values for the figure. If NULL (default), it will be computed starting at 0.
#' @param figureperiod Logical. If TRUE (default), periods are used in the figure instead
#' of scales.
#' @param xlab A string giving a custom X axis label.
#' @param ylab A string giving a custom Y axis label. If NULL (default) the Y label is
#' either "Scale" or "Period" depending on the value of \code{figureperiod} if
#' \code{length(scales) > 1}, or "Windowed Scalogram" if \code{length(scales) == 1}.
#' @param main A string giving a custom main title for the figure.
#' @param zlim A vector of length 2 with the limits for the z-axis (the color bar).
#'
#' @return A list with the following fields:
#' \itemize{
#' \item \code{wsc}: A matrix of size \code{length(tcentral)} x \code{length(scales)}
#' containing the values of the windowed scalograms at each scale and at each time window.
#' \item \code{tcentral}: The vector of central times at which the windows are centered.
#' \item \code{scales}: The vector of the scales.
#' \item \code{fourierfactor}: A factor for converting scales into periods.
#' \item \code{coi_maxscale}: A vector of length \code{length(tcentral)} containing the
#' values of the maximum scale from which there are border effects for the respective
#' central time.
#' }
#'
#' @importFrom graphics abline
#'
#' @examples
#' dt <- 0.1
#' time <- seq(0, 50, dt)
#' signal <- c(sin(pi * time), sin(pi * time / 2))
#' wscalog <- windowed_scalogram(signal = signal, dt = dt)
#'
#'
#' @section References:
#'
#' C. Torrence, G. P. Compo. A practical guide to wavelet analysis. B. Am. Meteorol. Soc.
#' 79 (1998), 61–78.
#'
#' V. J. Bolós, R. Benítez, R. Ferrer, R. Jammazi. The windowed scalogram difference: a
#' novel wavelet tool for comparing time series. Appl. Math. Comput., 312 (2017), 49-65.
#'
#' R. Benítez, V. J. Bolós, M. E. Ramírez. A wavelet-based tool for studying
#' non-periodicity. Comput. Math. Appl. 60 (2010), no. 3, 634-641.
#'
#' @export
#'

windowed_scalogram <-
function(signal,
dt = 1,
scales = NULL,
powerscales = TRUE,
delta_t = NULL,
wname = c("MORLET", "DOG", "PAUL", "HAAR", "HAAR2"),
wparam = NULL,
border_effects = c("BE", "INNER", "PER", "SYM"),
energy_density = TRUE,
makefigure = TRUE,
time_values = NULL,
figureperiod = TRUE,
xlab = "Time",
ylab = NULL,
main = "Windowed Scalogram",
zlim = NULL) {

#  require(zoo)
#  require(Matrix)

wname <- toupper(wname)
wname <- match.arg(wname)

if ((wname == "MORLET") || (wname == "DOG")) {
} else if (wname == "PAUL") {
} else { # HAAR
}
}

border_effects <- toupper(border_effects)
border_effects <- match.arg(border_effects)
if (border_effects == "INNER") {
border_effects_cwt <- "BE"
} else {
border_effects_cwt <- border_effects
}

nt <- length(signal)

if (is.null(delta_t)) {
delta_t <- ceiling(nt / 256)
}

} else {
}

fourierfactor <- fourier_factor(wname = wname, wparam = wparam)

if (is.null(scales)) {
scmin <- 2 / fourierfactor
if (powerscales) {
scales <- pow2scales(c(scmin, scmax, ceiling(256 / log2(scmax / scmin))))
} else {
scales <- seq(scmin, scmax, by = (scmax - scmin) / 256)
}
scalesdt <- scales * dt
} else {
if (powerscales && length(scales) == 3) {
scales <- pow2scales(scales)
} else {
if (is.unsorted(scales)) {
warning("Scales were not sorted.")
scales <- sort(scales)
}
aux <- diff(log2(scales))
if (powerscales && ((max(aux) - min(aux)) / max(aux) > 0.05)) {
warning("Scales seem like they are not power 2 scales. Powerscales set to FALSE.")
powerscales <- FALSE
}
}
scalesdt <- scales
scales <- scales / dt
}

ns <- length(scales)

cwt <- cwt_wst(signal = signal,
dt = dt,
scales = scalesdt,
powerscales = FALSE,
wname = wname,
wparam = wparam,
border_effects = border_effects_cwt,
makefigure = FALSE)

coefs <- cwt\$coefs

if (border_effects == "INNER") {

if (tcentral_ini > tcentral_end) {
stop("We need a larger signal")
}
tcentral <- seq(from = tcentral_ini, to = tcentral_end, by = delta_t)
ntcentral <- length(tcentral)
wsc <- matrix(NA, nrow = ntcentral, ncol = ns)

abscoefs2 <- matrix(abs(coefs) ^ 2, nrow = nt, ncol = ns)

# Regular version
for (i in 1:ntcentral) {
for (j in 1:ns) {
t_ini <- max(tcentral[i] - windowrad, 1 + wrs[j])
t_end <- min(tcentral[i] + windowrad, nt - wrs[j])
if (t_ini <= t_end) {
wsc[i, j] <- sqrt(abs(sum(abscoefs2[t_ini:t_end, j]))) # abs: sometimes wsc is negative due to numerical errors
wsc[i, j] <- wsc[i, j] / sqrt(t_end - t_ini + 1) # Normalization
}
}
}

wsc <- as.matrix(wsc)

} else {

tcentral <- seq(from = tcentral_ini, to = tcentral_end, by = delta_t)
ntcentral <- length(tcentral)
wsc <- matrix(0, nrow = ntcentral, ncol = ns)

abscoefs2 <- matrix(abs(coefs) ^ 2, nrow = nt, ncol = ns)

if (delta_t < windowrad) { # Fast version
for (j in 1:ns) {
wsc[1, j] <- sum(abscoefs2[1:(1 + 2 * windowrad), j])
for (i in 2:ntcentral) {
wsc[i, j] <- wsc[i-1, j] - sum(abscoefs2[(tcentral[i] - windowrad - delta_t):(tcentral[i] - windowrad - 1), j]) + sum(abscoefs2[(tcentral[i] + windowrad - delta_t + 1):(tcentral[i] + windowrad), j])
}
}
} else { # Regular version
for (i in 1:ntcentral) {
for (j in 1:ns) {
}
}
}
wsc <- as.matrix(sqrt(abs(wsc))) # abs: sometimes wsc is negative due to numerical errors
wsc <- wsc / sqrt(2 * windowrad + 1) # Normalization

}

# COI
coi_maxscale <- numeric(ntcentral)
for (i in 1:ntcentral) {
}

tcentraldt <- tcentral * dt

# Energy density
if (energy_density) {
wsc <- t(t(wsc) / sqrt(scalesdt))
}

# Make figure
if (makefigure) {

if (figureperiod) {
Y <- fourierfactor * scalesdt
coi <- fourierfactor * coi_maxscale
if (is.null(ylab)) ylab <- "Period"
} else {
Y <- scalesdt
coi <- coi_maxscale
if (is.null(ylab)) ylab <- "Scale"
}

if (is.null(time_values)) {
X <- tcentraldt
} else {
if (length(time_values) != nt) {
warning("Invalid length of time_values vector. Changing to default.")
X <- tcentraldt
} else {
X <- time_values[tcentral]
}
}

if (length(Y) > 1) {
wavPlot(
Z = wsc,
X = X,
Y = Y,
Ylog = powerscales,
coi = coi,
Xname = xlab,
Yname = ylab,
Zname = main,
zlim = zlim
)
} else {
if (is.null(ylab)) ylab <- "Windowed Scalogram"
plot(X, wsc, type = "l", xlab = xlab, ylab = ylab, main = main, xaxt = "n")
axis(side = 1, at = X[1 + floor((0:8) * (ntcentral - 1) / 8)])
abline(v = range(X[(coi > Y)]), lty = 2)
}

}

return(list(
wsc = wsc,
tcentral = tcentraldt,
scales = scalesdt,