Test/WaveletTransform.R

In yifanzhang0842/CoFESWave: This is a Test for CoFESWave package

WaveletTransform <- function (x, dt = 1, dj = 1/20, lowerPeriod = 2 * dt, upperPeriod = floor(length(x) *
                                                                            dt/3))
{
  series.length = length(x)
  pot2 = trunc(log2(series.length) + 0.5)
  pad.length = 2^(pot2 + 1) - series.length
  omega0 = 6
  fourier.factor = (2 * pi)/omega0
  min.scale = lowerPeriod/fourier.factor
  max.scale = upperPeriod/fourier.factor
  J = as.integer(log2(max.scale/min.scale)/dj)
  scales = min.scale * 2^((0:J) * dj)
  scales.length = length(scales)
  periods = fourier.factor * scales
  N = series.length + pad.length
  omega.k = 1:floor(N/2)
  omega.k = omega.k * (2 * pi)/(N * dt)
  omega.k = c(0, omega.k, -omega.k[floor((N - 1)/2):1])
  morlet.wavelet.transform = function(x) {
    x = (x - mean(x))/sd(x)
    xpad = c(x, rep(0, pad.length))
    fft.xpad = fft(xpad)
    wave = matrix(0, nrow = scales.length, ncol = N)
    wave = wave + (0+1i) * wave
    for (ind.scale in (1:scales.length)) {
      my.scale = scales[ind.scale]
      norm.factor = pi^(1/4) * sqrt(2 * my.scale/dt)
      expnt = -((my.scale * omega.k - omega0)^2/2) * (omega.k >
                                                        0)
      daughter = norm.factor * exp(expnt)
      daughter = daughter * (omega.k > 0)
      wave[ind.scale, ] = fft(fft.xpad * daughter, inverse = TRUE)/N
    }
    wave = wave[, 1:series.length]
    return(wave)
  }
  Wave = morlet.wavelet.transform(x)
  Power = Mod(Wave)^2/matrix(rep(scales, series.length), nrow = scales.length)
  Phase = Arg(Wave)
  Ampl = Mod(Wave)/matrix(rep(sqrt(scales), series.length),
                          nrow = scales.length)
  output = list(Wave = Wave, Phase = Phase, Ampl = Ampl, Period = periods,
                Scale = scales, Power = Power, nc = series.length, nr = scales.length)
  return(invisible(output))
}