R/cmorwavf.R

# cmorwavf.R
# Copyright (C) 2019 Geert van Boxtel <[email protected]>
# Matlab/Octave signal package:
# Copyright (C) 2007 Sylvain Pelissier <[email protected]>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; see the file COPYING. If not, see
# <https://www.gnu.org/licenses/>.
#
# 20191122 Geert van Boxtel          First version for v0.1.0
#---------------------------------------------------------------------------------------------------------------------------------

#' Complex Morlet Wavelet
#' 
#' the Morlet wavelet (or Gabor wavelet) is a wavelet composed of a complex exponential (carrier) multiplied by a Gaussian window
#' (envelope). The wavelet exists as a complex version or a purely real-valued version. Some distinguish between the "real Morlet"
#' versus the "complex Morlet". Others consider the complex version to be the "Gabor wavelet", while the real-valued version is the
#' "Morlet wavelet". This function returns the complex Morlet wavelet, with time-decay parameter \code{fb}, and center frequency 
#' \code{fc}. The general expression for the complex Morlet wavelet is
#' \deqn{\Psi(x) = ((\pi * fb)^{-0.5}) * \exp{2 * \pi * i * fc * x} * exp{-(x^2) / fb}}
#' \code{x} is evaluated on an \code{n}-point regular grid in the interval (lb, ub).
#' 
#' \code{fb} controls the decay in the time domain and the corresponding energy spread (bandwidth) in the frequency domain.
#' \code{fb} is the inverse of the variance in the frequency domain. Increasing \code{fb} makes the wavelet energy more concentrated
#' around the center frequency and results in slower decay of the wavelet in the time domain. Decreasing \code{fb} results in faster
#' decay of the wavelet in the time domain and less energy spread in the frequency domain. The value of \code{fb} does not affect
#' the center frequency. When converting from scale to frequency, only the center frequency affects the frequency values.
#' The energy spread or bandwidth parameter affects how localized the wavelet is in the frequency domain. See the examples.
#' 
#' @param lb,ub Lower and upper bounds of the interval to evaluate the complex Morlet waveform on. Default: -8 to 8.
#' @param n Number of points on the grid between \code{lb} and \code{ub} (length of the wavelet). Default: 1000.
#' @param fb Time-decay parameter of the wavelet (bandwidth in the frequency domain). Must be a positive scalar. Default: 5.
#' @param fc Center frequency of the wavelet. Must be a positive scalar. Default: 1.
#' 
#' @return A list containing 2 variables; \code{x}, the grid on which the complex Morlet wavelet was evaluated, and \code{psi}
#' (\eqn{\Psi}), the evaluated wavelet on the grid \code{x}.
#' 
#' @examples
#' 
#' ## Construct a complex-valued Morlet wavelet with a bandwidth parameter of 1.5
#' ## and a center frequency of 1. Set the effective support to [-8,8] and the
#' ## length of the wavelet to 1000.
#' cmw <- cmorwavf(-8, 8, 1000, 1.5, 1)
#' 
#' # Plot the real and imaginary parts of the wavelet.
#' op <- par(mfrow = c(2, 1))
#' plot(cmw$x, Re(cmw$psi), type = "l", main = "Real Part")
#' plot(cmw$x, Im(cmw$psi), type = "l", main = "Imaginary Part")
#' par(op)
#' 
#' ## This example shows how the complex Morlet wavelet shape in the frequency domain is affected
#' ## by the value of the bandwidth parameter (fb). Both wavelets have a center frequency of 1.
#' ## One wavelet has an fb value of 0.5 and the other wavelet has a value of 8.
#' 
#' op <- par(mfrow = c(2,1))
#' cmw1 <- cmorwavf(fb = 0.5)
#' cmw2 <- cmorwavf(fb = 8)
#' 
#' # time domain plot
#' plot(cmw1$x, Re(cmw1$psi), type = "l", xlab = "Time", ylab = "", main = "Time domain, Real part")
#' lines(cmw2$x, Re(cmw2$psi), col = "red")
#' legend("topright", legend = c("fb = 0.5", "fb = 8"), lty= 1, col = c(1,2))
#' 
#' # frequency domain plot
#' f <- seq(-5, 5, .01)
#' Fc <- 1  
#' Fb1 <- 0.5  
#' Fb2 <- 8
#' PSI1 <- exp(-pi^2 * Fb1 * (f-Fc)^2)  
#' PSI2 <- exp(-pi^2 * Fb2 * (f-Fc)^2)  
#' plot(f, PSI1, type="l", xlab = "Frequency", ylab = "", main = "Frequency domain")
#' lines(f, PSI2, col = "red")
#' legend("topright", legend = c("fb = 0.5", "fb = 8"), lty= 1, col = c(1,2))
#' par(op)
#' 
#' ## The fb bandwidth parameter for the complex Morlet wavelet is the inverse of the
#' ## variance in frequency. Therefore, increasing Fb results in a narrower concentration
#' ## of energy around the center frequency.
#' 
#' ## alternative to the above frequency plot:
#' fs <- length(cmw1$x) / sum(abs(range(cmw1$x)))
#' hz  <- seq(0, fs/2, len=floor(length(cmw1$psi)/2)+1)
#' PSI1 <- fft(cmw1$psi) / length(cmw1$psi)
#' PSI2 <- fft(cmw2$psi) / length(cmw2$psi)
#' plot(hz, 2 * abs(PSI1)[1:length(hz)], type="l", xlab = "Frequency", 
#'      ylab = "", main = "Frequency domain", xlim=c(0,5))
#' lines(hz, 2 * abs(PSI2)[1:length(hz)], col = 2)
#' legend("topright", legend = c("fb = 0.5", "fb = 8"), lty= 1, col = c(1,2))
#' 
#' @author Original Matlab/Octave code Copyright (C) 2007 Sylvain Pelissier \email{[email protected]@gmail.com}.
#' Port to R by Geert van Boxtel \email{[email protected]@gmail.com}.
#
#' 
#' @export

cmorwavf <- function (lb = -8, ub = 8, n = 1000, fb = 5, fc = 1) {
  
  if (!isPosscal(n) || !isWhole(n) || n <= 0) stop('n must be an integer strictly positive')
  if (!isPosscal(fb) || fb <= 0) stop('fb must be a positive scalar > 0')
  if (!isPosscal(fc) || fc <= 0) stop('fc must be a positive scalar > 0')
  
  x <- seq(lb, ub, length.out = n)
  psi <- ((pi * fb)^(-0.5)) * exp(2 * 1i * pi * fc * x) * exp(-x^2 / fb)
  
  list(x = x, psi = psi)
}
gjmvanboxtel/gsignal documentation built on Dec. 9, 2019, 6:43 p.m.