SynchWave (version 1.1.2)

synsq_cwt_fw: Synchrosqueezing Transform

Description

This function calculates the synchrosqueezing transform.

This code is translated from MATLAB Synchrosqueezing Toolbox, version 1.1 developed by Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/).

Usage

synsq_cwt_fw(tt, x, nv=16, opt=NULL)

Arguments

tt

vector of times samples are taken (e.g. seq(0, 1,length=2))

x

vector of signal samples (e.g. x = cos(20*pi*tt))

nv

number of voices (default: 16, recommended: 32 or 64)

opt

list of options. opt.type: type of wavelet (see wfiltfn), opt$s, opt$mu, etc (wavelet parameters: see opt from wfiltfn), opt.disp: display debug information?, opt.gamma: wavelet hard thresholding value (see cwt_freq_direct, default: sqrt(machine epsilon)) opt.dtype: direct or numerical differentiation (default: 1, uses direct). if dtype=0, uses differentiation algorithm instead (see cwt_freq), which is faster and uses less memory, but may be less accurate.

Value

Tx

synchrosqueezed output of x (columns associated with time tt)

fs

frequencies associated with rows of Tx

Wx

wavelet transform of x (see cwt_fw)

asc

scales associated with rows of Wx

w

instantaneous frequency estimates for each element of Wx (see cwt_freq_direct)

Details

This function calculates the synchrosqueezing transform of vector x, with samples taken at times given in vector tt. Uses nv voices. opt is a struct of options for the way synchrosqueezing is done, the type of wavelet used, and the wavelet parameters. This implements the algorithm described in Sec. III of [1].

Note that Wx itself is not thresholded by opt$gamma. The instantaneoues frequency w is calculated using Wx by cwt_freq_direct and cwt_freq. After the calculation, instantaneoues frequency w is treated as NA where abs(Wx) < opt$gamma.

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

cwt_fw, est_riskshrink_thresh, wfiltfn.

Examples

# NOT RUN {
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

# Reconstruction
opt$gamma <- thresh
#[1] 0.0593984
#opt$gamma <- 10^-5
cwtrec <- cwt_iw(cwtfit$Wx, opt$type, opt)

par(mfrow=c(1,1))
plot(tt, f, type="p", lty=2, xlab="time", ylab="f", col="red", cex=0.1)
lines(tt, f0, col="blue")
lines(tt, cwtrec)

# Synchrosqueezed wavelet transform
sstfit <- synsq_cwt_fw(tt, f, nv, opt)

#par(mfrow=c(2,2))
#plot(tt, f, type="p", lty=2, xlab="time", ylab="f", col="red", cex=0.1)
#lines(tt, f0, col="blue")

#image.plot(list(x=tt, y=sstfit$asc, z=t(abs(sstfit$Wx))), log="y", 
#    xlab="Time", ylab="Scale", main="Time-Scale Representation by CWT",  
#    col=designer.colors(64, c("azure", "cyan", "blue", "darkblue")), ylim=rev(range(sstfit$asc)))
#image.plot(list(x=tt, y=sstfit$fs, z=t(abs(sstfit$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))
#image.plot(list(x=tt, y=sstfit$asc, z=t(sstfit$w)), log="y", 
#    xlab="Time", ylab="Scale", main="Instantaneous Frequency",  
#    col=designer.colors(64, c("azure", "cyan", "blue", "darkblue")), ylim=rev(range(sstfit$asc)))
# }