cmorwavf: Complex Morlet Wavelet

In gjmvanboxtel/gsignal: Signal Processing

Complex Morlet Wavelet

Description

Compute the complex Morlet wavelet on a grid.

Usage

cmorwavf(lb = -8, ub = 8, n = 1000, fb = 5, fc = 1)

Arguments

lb, ub	Lower and upper bounds of the interval to evaluate the complex Morlet waveform on. Default: -8 to 8.
n	Number of points on the grid between lb and ub (length of the wavelet). Default: 1000.
fb	Time-decay parameter of the wavelet (bandwidth in the frequency domain). Must be a positive scalar. Default: 5.
fc	Center frequency of the wavelet. Must be a positive scalar. Default: 1.

Details

The Morlet (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 fb, and center frequency fc. The general expression for the complex Morlet wavelet is

 Psi(x) = ((pi * fb)^-0.5) *  exp(2 * pi * i * fc * x) * exp(-(x^2) / fb)

x is evaluated on an n-point regular grid in the interval (lb, ub).

fb controls the decay in the time domain and the corresponding energy spread (bandwidth) in the frequency domain. fb is the inverse of the variance in the frequency domain. Increasing fb makes the wavelet energy more concentrated around the center frequency and results in slower decay of the wavelet in the time domain. Decreasing fb results in faster decay of the wavelet in the time domain and less energy spread in the frequency domain. The value of 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.

Value

A list containing 2 variables; x, the grid on which the complex Morlet wavelet was evaluated, and psi (\Psi), the evaluated wavelet on the grid x.

Author(s)

Sylvain Pelissier, [e-mail].
Conversion to R by Geert van Boxtel, [e-mail].

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))

gjmvanboxtel/gsignal documentation built on Nov. 22, 2023, 8:19 p.m.