# Synchrosqueezing Transform 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/). ```r 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`. ## Returns - **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 ```r 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))) ```

synsq_cwt_fw function