curve_ext_multi {SynchWave}

Extract a Maximum Energy / Minimum Curvature Curves

Description

This function consecutively extracts maximum energy / minimum curvature curves from a synchrosqueezing representation. This code is translated from MATLAB Synchrosqueezing Toolbox, version 1.1 developed by Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/).

Usage

curve_ext_multi(Tx, fs, nc, lambda = 1e3, clwin = 4)

Arguments

Tx

same as input to curve_ext. synchrosqueezed output of x (columns associated with time t)

fs

same as input to curve_ext. frequencies associated with rows of Tx

nc

Number of curves to extract

lambda

same as input to curve_ext. lambda default = 1e3.

clwin

frequency clearing window; after curve extraction, a window of frequencies (Cs[,i]-clwin) : (Cs[,i]+clwin) is removed from the original representation. (default = 2)

Details

This function consecutively extracts maximum energy / minimum curvature curves from a synchrosqueezing representation. As curves are extracted, their associated energies are zeroed out in the synsq representation. Curve extraction is then again performed on the remaining representaton.

For more details, see help curve_ext and Sec. IIID of [1].

Value

Cs

[N x nc] matrix of curve indices (N=ncol(Tx))

Es

[nc x 1] vector of associated (logarithmic) energies

References

[1] Thakur, G., Brevdo, E., Fuckar, N. S. and Wu, H-T. (2013) The Synchrosqueezing algorithm for time-varying spectral analysis: Robustness properties and new paleoclimate applications. Signal Processing, 93, 1079–1094.

See Also

synsq_cwt_fw, curve_ext, curve_ext_recon.

Examples

set.seed(7)
tt <- seq(0, 10, , 1024)
nv <- 32
f0 <- (1+0.6*cos(2*tt))*cos(4*pi*tt+1.2*tt^2)
sigma <- 0.5
f <- f0 + sigma*rnorm(length(tt))

# Continuous wavelet transform
opt <- list(type = "bump")
cwtfit <- cwt_fw(f, opt$type, nv, tt[2]-tt[1], opt)

# Hard thresholing
thresh <- est_riskshrink_thresh(cwtfit$Wx, nv)
cwtfit$Wx[which(abs(cwtfit$Wx) < thresh)] <- 0.0

# Denoised signal
opt$gamma <- thresh
fr <- cwt_iw(cwtfit$Wx, opt$type, opt)

# Synchrosqueezed wavelet transform using denoised signal
sstfit2 <- synsq_cwt_fw(tt, fr, nv, opt)

# Ridge extraction
lambda <- 1e+04
nw <- 16
imtfit <- curve_ext_multi(sstfit2$Tx, log2(sstfit2$fs), 1, lambda, nw)

# Reconstruction
curvefit <- curve_ext_recon(sstfit2$Tx, sstfit2$fs, imtfit$Cs, opt, nw)

par(mfrow=c(2,1))
image.plot(list(x=tt, y=sstfit2$fs, z=t(abs(sstfit2$Tx))), log="y", 
    xlab="Time", ylab="Frequency", main="Time-Frequency Representation by SST", 
    col=designer.colors(64, c("azure", "cyan", "blue", "darkblue")), ylim=c(0.5, 25))
lines(tt, sstfit2$fs[imtfit$Cs[,1]], col="red", lty=3, lwd=2)

plot(tt, f0, type="l")
lines(tt, curvefit, lty=2)