synsq_cwt_fw {SynchWave}

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.

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.

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)

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

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)))