| Title: | Synchrosqueezed Wavelet Transform |
|---|---|
| Description: | The synchrosqueezed wavelet transform is implemented. The package is a translation of MATLAB Synchrosqueezing Toolbox, version 1.1 originally developed by Eugene Brevdo (2012). The C code for curve_ext was authored by Jianfeng Lu, and translated to Fortran by Dongik Jang. Synchrosqueezing is based on the papers: [1] Daubechies, I., Lu, J. and Wu, H. T. (2011) Synchrosqueezed wavelet transforms: An empirical mode decomposition-like tool. Applied and Computational Harmonic Analysis, 30. 243-261. [2] Thakur, G., Brevdo, E., Fukar, 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. |
| Authors: | Matlab original by Eugene Brevdo; R port by Dongik Jang, Hee-Seok Oh and Donghoh Kim |
| Maintainer: | Donghoh Kim <[email protected]> |
| License: | LGPL (>= 2) |
| Version: | 1.1.2 |
| Built: | 2025-02-07 03:11:03 UTC |
| Source: | https://github.com/cran/SynchWave |
This function extracts a maximum energy / minimum curvature curve from Synchrosqueezed Representation. This code is translated from MATLAB Synchrosqueezing Toolbox, version 1.1 developed by Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/).
curve_ext(Tx, fs, lambda=0)curve_ext(Tx, fs, lambda=0)
Tx |
synchrosqueezed output of x (columns associated with time t) |
fs |
frequencies associated with rows of |
lambda |
lambda should be greater than or equal to 0. Default: lambda=0 |
This function extracts a maximum energy, minimum curvature, curve from
Synchrosqueezed Representation. Note, energy is given as:
abs(Tx)^2.
This implements the solution to Eq. (8) of [1].
Original Author: Jianfeng Lu
C |
the curve locations (indices) |
E |
the (logarithmic) energy of the curve |
[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.
synsq_cwt_fw, curve_ext_multi, curve_ext_recon.
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/).
curve_ext_multi(Tx, fs, nc, lambda = 1e3, clwin = 4)curve_ext_multi(Tx, fs, nc, lambda = 1e3, clwin = 4)
Tx |
same as input to |
fs |
same as input to |
nc |
Number of curves to extract |
lambda |
same as input to |
clwin |
frequency clearing window; after curve extraction, a window
of frequencies ( |
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].
Cs |
|
Es |
|
[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.
synsq_cwt_fw, curve_ext, curve_ext_recon.
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)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)
This function reconstructs the curves given by curve_ext or curve_ext_multi.
This code is translated from MATLAB Synchrosqueezing Toolbox, version 1.1 developed by Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/).
curve_ext_recon(Tx, fs, Cs, opt=NULL, clwin=4)curve_ext_recon(Tx, fs, Cs, opt=NULL, clwin=4)
Tx |
same as input to |
fs |
same as input to |
opt |
same as input to |
Cs |
same as output to |
clwin |
same as input to |
This function reconstructs the curves given by curve_ext or curve_ext_multi.
xrs |
|
[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.
synsq_cwt_fw, curve_ext, curve_ext_multi, curve_ext_recon.
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)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)
This function performs forward continuous wavelet transform, discretized, as described in Sec. 4.3.3 of [1] and Sec. IIIA of [2]. This code is translated from MATLAB Synchrosqueezing Toolbox, version 1.1 developed by Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/).
cwt_fw(x, type, nv, dt, opt)cwt_fw(x, type, nv, dt, opt)
x |
input signal vector, length |
type |
wavelet type, string (see |
nv |
number of voices (suggest |
dt |
sampling period (default, |
opt |
list of options. |
This function performs forward continuous wavelet transform, discretized, as described in Sec. 4.3.3 of [1] and Sec. IIIA of [2]. This algorithm uses the FFT and samples the wavelet atoms in the fourier domain. Options such as padding of the original signal are allowed. Returns the vector of scales and, if requested, the analytic time-derivative of the wavelet transform (as described in Sec. IIIB of [2].
Wx |
|
asc |
|
dWx |
|
[1] Mallat, S (2009) A Wavelet Tour of Signal Processing, New York: Academic Press.
[2] 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.
cwt_iw, wfiltfn, est_riskshrink_thresh.
tt <- seq(0, 10, , 1024) 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 nv <- 32 opt <- list(type = "bump") cwtfit <- cwt_fw(f, opt$type, nv, tt[2]-tt[1], opt)tt <- seq(0, 10, , 1024) 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 nv <- 32 opt <- list(type = "bump") cwtfit <- cwt_fw(f, opt$type, nv, tt[2]-tt[1], opt)
This function performs
the inverse wavelet transform of signal Wx.
This code is translated from MATLAB Synchrosqueezing Toolbox, version 1.1 developed by Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/).
cwt_iw(Wx, type, opt=NULL)cwt_iw(Wx, type, opt=NULL)
Wx |
wavelet transform of a signal, see |
type |
wavelet used to take the wavelet transform,
see |
opt |
options list used for forward wavelet transform. |
This function performs
the inverse wavelet transform of signal Wx, and
implements Eq. (4.67) of [1].
x |
the signal, as reconstructed from |
[1] Mallat, S (2009) A Wavelet Tour of Signal Processing, New York: Academic Press.
cwt_fw, wfiltfn, est_riskshrink_thresh.
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 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)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 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)
This function estimates the RiskShrink hard thresholding level.
This code is translated from MATLAB Synchrosqueezing Toolbox, version 1.1 developed by Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/).
est_riskshrink_thresh(Wx, nv)est_riskshrink_thresh(Wx, nv)
Wx |
wavelet transform of a signal, see |
nv |
number of voices |
This function implements Defn. 1 of Sec. 2.4 in [1], using the suggested noise estimator from the discussion "Estimating the noise level" in that same section.
the RiskShrink hard threshold estimate
[1] Donoho, D. L. and Johnstone, I. M. (1994) Ideal spatial adaptation by wavelet shrinkage. Biometrika, 81, 425–455.
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.0tt <- 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
This function exchanges the left halves of a vector with the right halves.
fftshift(x)fftshift(x)
x |
a vector |
This function exchanges the left halves of a vector with the right halves. This function is adapted from Matlab.
shifted vector
x <- 1:4 fftshift(fftshift(x)) ifftshift(fftshift(x)) x <- 1:5 fftshift(fftshift(x)) ifftshift(fftshift(x))x <- 1:4 fftshift(fftshift(x)) ifftshift(fftshift(x)) x <- 1:5 fftshift(fftshift(x)) ifftshift(fftshift(x))
This function exchanges the left halves of a vector with the right halves.
ifftshift(x)ifftshift(x)
x |
a vector |
This function exchanges the left halves of a vector with the right halves. This function is adapted from Matlab.
shifted vector
x <- 1:4 fftshift(fftshift(x)) ifftshift(fftshift(x)) x <- 1:5 fftshift(fftshift(x)) ifftshift(fftshift(x))x <- 1:4 fftshift(fftshift(x)) ifftshift(fftshift(x)) x <- 1:5 fftshift(fftshift(x)) ifftshift(fftshift(x))
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/).
synsq_cwt_fw(tt, x, nv=16, opt=NULL)synsq_cwt_fw(tt, x, nv=16, opt=NULL)
tt |
vector of times samples are taken (e.g. |
x |
vector of signal samples (e.g. |
nv |
number of voices (default: 16, recommended: 32 or 64) |
opt |
list of options. |
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.
Tx |
synchrosqueezed output of |
fs |
frequencies associated with rows of |
Wx |
wavelet transform of |
asc |
scales associated with rows of |
w |
instantaneous frequency estimates for each element of |
[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.
cwt_fw, est_riskshrink_thresh, wfiltfn.
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)))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)))
This function performs inverse synchrosqueezing transform.
This code is translated from MATLAB Synchrosqueezing Toolbox, version 1.1 developed by Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/).
synsq_cwt_iw(Tx, fs, opt=NULL)synsq_cwt_iw(Tx, fs, opt=NULL)
Tx |
synchrosqueezed output of |
fs |
frequencies associated with rows of |
opt |
list of options. See |
This function performs
inverse synchrosqueezing transform of Tx with associated
frequencies in fs. This implements Eq. 5 of [1].
x |
reconstructed signal |
[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.
synsq_cwt_fw, synsq_filter_pass.
set.seed(7) n <- 2048 tu <- seq(0,10,, n) dt <- tu[2]-tu[1] feq1 <- function(t) (1+0.2*cos(t))*cos(2*pi*(2*t+0.3*cos(t))) feq2 <- function(t) (1+0.3*cos(2*t))*exp(-t/15)*cos(2*pi*(2.4*t+0.5*t^(1.2)+0.3*sin(t))) feq3 <- function(t) cos(2*pi*(5.3*t-0.2*t^(1.3))) feq <- function(t) feq1(t) + feq2(t) + feq3(t) s2 <- 2.4 noise <- sqrt(s2)*rnorm(length(tu)) fu0 <- feq(tu); fu <- fu0 + noise; fus <- cbind(feq1(tu), feq2(tu), feq3(tu)) # Continuous wavelet transform nv <- 32 opt <- list(type = "bump") cwtfit <- cwt_fw(fu, opt$type, nv, dt, opt) thresh <- est_riskshrink_thresh(cwtfit$Wx, nv) # Hard thresholding and Reconstruction cwtfit$Wx[which(abs(cwtfit$Wx) < thresh)] <- 0.0 fur <- cwt_iw(cwtfit$Wx, opt$type, opt) # Synchrosqueezed wavelet transform using denoised signal sstfit <- synsq_cwt_fw(tu, fur, nv, opt) #par(mfrow=c(1,2)) #image.plot(list(x=tu, 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(1, 8)) # Extracting the second component by filtering of synchrosqueezed wavelet transform fm <- fM <- (2.4+0.5*1.2*tu^0.2+0.3*cos(tu)) #lines(tu, 0.88*fm, col="red", lty=3, lwd=2) #lines(tu, 1.22*fM, col="red", lty=3, lwd=2) tmp <- synsq_filter_pass(sstfit$Tx, sstfit$fs, 0.88*fm, 1.12*fM); fursst <- synsq_cwt_iw(tmp$Txf, w, opt); #plot(tu, fursst, type="l", main="SST", xlab="time", ylab="f", col="red", # xlim=c(1.5,8.5), ylim=c(-1,1)) #lines(tu, feq2(tu), col="blue") # Example: # tmp <- synsq_cwt_fw(t, x, 32) # Synchrosqueezing # Txf <- synsq_filter_pass(tmp$Tx, tmp$fs, -Inf, 1) # Pass band filter # xf <- synsq_cwt_iw(Txf, fs) # Filtered signal reconstructionset.seed(7) n <- 2048 tu <- seq(0,10,, n) dt <- tu[2]-tu[1] feq1 <- function(t) (1+0.2*cos(t))*cos(2*pi*(2*t+0.3*cos(t))) feq2 <- function(t) (1+0.3*cos(2*t))*exp(-t/15)*cos(2*pi*(2.4*t+0.5*t^(1.2)+0.3*sin(t))) feq3 <- function(t) cos(2*pi*(5.3*t-0.2*t^(1.3))) feq <- function(t) feq1(t) + feq2(t) + feq3(t) s2 <- 2.4 noise <- sqrt(s2)*rnorm(length(tu)) fu0 <- feq(tu); fu <- fu0 + noise; fus <- cbind(feq1(tu), feq2(tu), feq3(tu)) # Continuous wavelet transform nv <- 32 opt <- list(type = "bump") cwtfit <- cwt_fw(fu, opt$type, nv, dt, opt) thresh <- est_riskshrink_thresh(cwtfit$Wx, nv) # Hard thresholding and Reconstruction cwtfit$Wx[which(abs(cwtfit$Wx) < thresh)] <- 0.0 fur <- cwt_iw(cwtfit$Wx, opt$type, opt) # Synchrosqueezed wavelet transform using denoised signal sstfit <- synsq_cwt_fw(tu, fur, nv, opt) #par(mfrow=c(1,2)) #image.plot(list(x=tu, 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(1, 8)) # Extracting the second component by filtering of synchrosqueezed wavelet transform fm <- fM <- (2.4+0.5*1.2*tu^0.2+0.3*cos(tu)) #lines(tu, 0.88*fm, col="red", lty=3, lwd=2) #lines(tu, 1.22*fM, col="red", lty=3, lwd=2) tmp <- synsq_filter_pass(sstfit$Tx, sstfit$fs, 0.88*fm, 1.12*fM); fursst <- synsq_cwt_iw(tmp$Txf, w, opt); #plot(tu, fursst, type="l", main="SST", xlab="time", ylab="f", col="red", # xlim=c(1.5,8.5), ylim=c(-1,1)) #lines(tu, feq2(tu), col="blue") # Example: # tmp <- synsq_cwt_fw(t, x, 32) # Synchrosqueezing # Txf <- synsq_filter_pass(tmp$Tx, tmp$fs, -Inf, 1) # Pass band filter # xf <- synsq_cwt_iw(Txf, fs) # Filtered signal reconstruction
This function filters the Synchrosqueezing representation Tx, having associated
frequencies fs (see synsq_cwt_fw).
This band-pass filter keeps frequencies in the
range [fm, fM].
This code is translated from MATLAB Synchrosqueezing Toolbox, version 1.1 developed by Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/).
synsq_filter_pass(Tx, fs, fm, fM)synsq_filter_pass(Tx, fs, fm, fM)
Tx |
synchrosqueezed output of |
fs |
frequencies associated with rows of |
fm |
Minimum band pass values. scalars, or vectors with each value associated
with a time index ( |
fM |
Maximum band pass values. scalars, or vectors with each value associated
with a time index ( |
This function filters the Synchrosqueezing representation Tx, having associated
frequencies fs (see synsq_cwt_fw).
This band-pass filter keeps frequencies in the
range [fm, fM].
Txf |
Filtered version of |
fmi |
time-length vector of min-frequency row indices |
fMi |
time-length vector of max-frequency row indices |
set.seed(7) n <- 2048 tu <- seq(0,10,, n) dt <- tu[2]-tu[1] feq1 <- function(t) (1+0.2*cos(t))*cos(2*pi*(2*t+0.3*cos(t))) feq2 <- function(t) (1+0.3*cos(2*t))*exp(-t/15)*cos(2*pi*(2.4*t+0.5*t^(1.2)+0.3*sin(t))) feq3 <- function(t) cos(2*pi*(5.3*t-0.2*t^(1.3))) feq <- function(t) feq1(t) + feq2(t) + feq3(t) s2 <- 2.4 noise <- sqrt(s2)*rnorm(length(tu)) fu0 <- feq(tu); fu <- fu0 + noise; fus <- cbind(feq1(tu), feq2(tu), feq3(tu)) # Continuous wavelet transform nv <- 32 opt <- list(type = "bump") cwtfit <- cwt_fw(fu, opt$type, nv, dt, opt) thresh <- est_riskshrink_thresh(cwtfit$Wx, nv) # Hard thresholding and Reconstruction cwtfit$Wx[which(abs(cwtfit$Wx) < thresh)] <- 0.0 fur <- cwt_iw(cwtfit$Wx, opt$type, opt) # Synchrosqueezed wavelet transform using denoised signal sstfit <- synsq_cwt_fw(tu, fur, nv, opt) #par(mfrow=c(2,2)) #image.plot(list(x=tu, 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=c(1, 0.0625)) # Extracting the second component by filtering of continuous wavelet transform am <- 0.2 * rep(1, length(tu)) aM <- 0.3 * rep(1, length(tu)) #lines(tu, am, col="red", lty=3, lwd=2) #lines(tu, aM, col="red", lty=3, lwd=2) tmp <- synsq_filter_pass(sstfit$Wx, sstfit$asc, am, aM); furcwt <- cwt_iw(tmp$Txf, opt$type, opt); #image.plot(list(x=tu, 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(1, 8)) # Extracting the second component by filtering of synchrosqueezed wavelet transform fm <- fM <- (2.4+0.5*1.2*tu^0.2+0.3*cos(tu)) #lines(tu, 0.88*fm, col="red", lty=3, lwd=2) #lines(tu, 1.22*fM, col="red", lty=3, lwd=2) tmp <- synsq_filter_pass(sstfit$Tx, sstfit$fs, 0.88*fm, 1.12*fM); fursst <- synsq_cwt_iw(tmp$Txf, w, opt); #plot(tu, fursst, type="l", main="SST", xlab="time", ylab="f", col="red", # xlim=c(1.5,8.5), ylim=c(-1,1)) #lines(tu, feq2(tu), col="blue") #plot(tu, furcwt, type="l", main="CWT", xlab="time", ylab="f", col="red", # xlim=c(1.5,8.5), ylim=c(-1,1)) #lines(tu, feq2(tu), col="blue") # Remove all energy for normalized frequencies above 1. # synsq_filter_pass(Tx, fs, -Inf, 1)set.seed(7) n <- 2048 tu <- seq(0,10,, n) dt <- tu[2]-tu[1] feq1 <- function(t) (1+0.2*cos(t))*cos(2*pi*(2*t+0.3*cos(t))) feq2 <- function(t) (1+0.3*cos(2*t))*exp(-t/15)*cos(2*pi*(2.4*t+0.5*t^(1.2)+0.3*sin(t))) feq3 <- function(t) cos(2*pi*(5.3*t-0.2*t^(1.3))) feq <- function(t) feq1(t) + feq2(t) + feq3(t) s2 <- 2.4 noise <- sqrt(s2)*rnorm(length(tu)) fu0 <- feq(tu); fu <- fu0 + noise; fus <- cbind(feq1(tu), feq2(tu), feq3(tu)) # Continuous wavelet transform nv <- 32 opt <- list(type = "bump") cwtfit <- cwt_fw(fu, opt$type, nv, dt, opt) thresh <- est_riskshrink_thresh(cwtfit$Wx, nv) # Hard thresholding and Reconstruction cwtfit$Wx[which(abs(cwtfit$Wx) < thresh)] <- 0.0 fur <- cwt_iw(cwtfit$Wx, opt$type, opt) # Synchrosqueezed wavelet transform using denoised signal sstfit <- synsq_cwt_fw(tu, fur, nv, opt) #par(mfrow=c(2,2)) #image.plot(list(x=tu, 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=c(1, 0.0625)) # Extracting the second component by filtering of continuous wavelet transform am <- 0.2 * rep(1, length(tu)) aM <- 0.3 * rep(1, length(tu)) #lines(tu, am, col="red", lty=3, lwd=2) #lines(tu, aM, col="red", lty=3, lwd=2) tmp <- synsq_filter_pass(sstfit$Wx, sstfit$asc, am, aM); furcwt <- cwt_iw(tmp$Txf, opt$type, opt); #image.plot(list(x=tu, 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(1, 8)) # Extracting the second component by filtering of synchrosqueezed wavelet transform fm <- fM <- (2.4+0.5*1.2*tu^0.2+0.3*cos(tu)) #lines(tu, 0.88*fm, col="red", lty=3, lwd=2) #lines(tu, 1.22*fM, col="red", lty=3, lwd=2) tmp <- synsq_filter_pass(sstfit$Tx, sstfit$fs, 0.88*fm, 1.12*fM); fursst <- synsq_cwt_iw(tmp$Txf, w, opt); #plot(tu, fursst, type="l", main="SST", xlab="time", ylab="f", col="red", # xlim=c(1.5,8.5), ylim=c(-1,1)) #lines(tu, feq2(tu), col="blue") #plot(tu, furcwt, type="l", main="CWT", xlab="time", ylab="f", col="red", # xlim=c(1.5,8.5), ylim=c(-1,1)) #lines(tu, feq2(tu), col="blue") # Remove all energy for normalized frequencies above 1. # synsq_filter_pass(Tx, fs, -Inf, 1)
This function provides wavelet transform function of the wavelet filter in question, fourier domain. This code is translated from MATLAB Synchrosqueezing Toolbox, version 1.1 developed by Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/).
wfiltfn(type, opt)wfiltfn(type, opt)
type |
wavelet type. ‘gauss’, ‘mhat’, ‘cmhat’, ‘morlet’, ‘shannon’, ‘hshannon’, ‘hhat’, ‘hhhat’, ‘bump’ |
opt |
list of options for wavelet type.
For ‘gauss’, |
This function provides wavelet transform function of the wavelet filter in question, fourier domain.
wavelet transform function
psihfn <- wfiltfn("bump", list(mu=1, s=.5)) plot(psihfn(seq(-5, 5, by=0.01)), type="l", xlab="", ylab="")psihfn <- wfiltfn("bump", list(mu=1, s=.5)) plot(psihfn(seq(-5, 5, by=0.01)), type="l", xlab="", ylab="")
This function outputs the FFT of the wavelet. This code is translated from MATLAB Synchrosqueezing Toolbox, version 1.1 developed by Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/).
wfilth(type, N, a, opt)wfilth(type, N, a, opt)
type |
wavelet type. See |
N |
number of samples to calculate |
a |
wavelet scale parameter (default = 1) |
opt |
list of options for wavelet type. See |
This function outputs the FFT of the wavelet of family 'type' with parameters in 'opt', of length N at scale a: (psi(-t/a)).
Note that the output is made so that the inverse fft of the
result is zero-centered in time. This is important for
convolving with the derivative(dpsih). To get the correct
output, perform an ifftshift. That is,
psi = ifftshift(fft(psih, inverse=TRUE) / length(psih)),
xfilt = ifftshift(fft(fft(x) * psih, inverse=TRUE) / length(fft(x) * psih))
psih |
wavelet sampling in frequency domain (for use in |
dpsih |
derivative of same wavelet, sampled in frequency domain (for |
xi |
associated fourier domain frequencies of the samples. |
tmp <- wfilth("morlet", 1024, 4) plot(fftshift(tmp$xi/(2*pi)), fftshift(abs(tmp$psih)), type="l", col="blue", xlab="", ylab="")tmp <- wfilth("morlet", 1024, 4) plot(fftshift(tmp$xi/(2*pi)), fftshift(abs(tmp$psih)), type="l", col="blue", xlab="", ylab="")