single-gaussian.R

require(wmtsa)

Gaussian <- function (x, nu, sigma, H) {
  return(H * exp(-(x - nu)^2/(2*sigma^2)))
}

PLOT_P <- TRUE

# Artificial data
set.seed(1)
N <- 10000
noise <- 0.1
freq <- seq(from=0.01, to = 250, length.out = N)
numax <- 75
senv <- 5
Hosc <- 0.05
psd <- Gaussian(freq, numax, senv, Hosc) + rnorm(length(freq), 0, noise)
#plot(freq,psd, type="l", xlab=expression(nu), 
#     ylab=expression(paste("f(",nu,")", sep="")))

# Make the wavelet transforms
psd.series <- create.signalSeries(psd,
                                  position = list(from=min(freq),
                                                  by=(max(freq)-min(freq))/N, 
                                                  units="frequency"),
                                  units = "frequency")
psd.cwt <- wavCWT(psd.series, wavelet = "gaussian2")
if (PLOT_P) plot(psd.cwt, series = TRUE)

# Find the scale with the highest wavelet transform value
a_ind <- floor(which.max(psd.cwt) / attr(psd.cwt, "n.sample"))
a_max <- attr(psd.cwt, "scale")[a_ind]
numax_predicted <- freq[which.max(psd.cwt[,a_ind])]

# Make a "shrinkage" to reduce noise.
psd.shrink <- wavShrink(psd)
if (PLOT_P) {
  par(mfrow = c(3,1), cex = 1.1)
  plot(freq, Gaussian(freq, numax, senv, Hosc), type="l",
       ylim = c(-0.05,0.05),
       xlab = expression(nu),
       ylab = expression(paste("f(", nu, ")", sep="")))
  plot(freq, psd, type="l", ylim = c(-0.3,0.3),
       xlab = expression(nu),
       ylab = expression(paste("f(", nu, ")", sep="")))
  plot(freq, psd.shrink, type="l",ylim = c(-0.05,0.05), 
       xlab = expression(nu),
       ylab = expression(paste("f(", nu, ")", sep="")))
  par(mfrow = c(1,1))
}