curve_ext_multi: Extract a Maximum Energy / Minimum Curvature Curves

In SynchWave: Synchrosqueezed Wavelet Transform

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)

SynchWave documentation built on May 7, 2022, 5:05 p.m.