The Synchrosqueezing algorithm for time-varying spectral analysis: robustness properties and new paleoclimate applications

Gaurav Thakur; Eugene Brevdo; Neven S. Fučkar; Hau-Tieng Wu

Source-linked AI summary

The Synchrosqueezing algorithm for time-varying spectral analysis: robustness properties and new paleoclimate applications

Gaurav Thakur, Eugene Brevdo, Neven S. Fučkar, Hau-Tieng Wu

arXiv:1105.0010v3Apr 29, 2011math.CAcs.CEmath.NAphysics.data-an

TL;DR

The paper addresses whether Synchrosqueezing can remain reliable for perturbed, noisy, and nonuniformly sampled signals while extracting time-varying oscillatory components. It develops stability results and a practical discretized implementation, then applies the method to paleoclimate data. The method delineates orbital cycles and improves the representation of orbital signals in δ18O records while supporting analysis of climate evolution over the past 2.5 million years.

Problem
Reliable time-frequency analysis and component extraction are needed for signals with time-varying oscillatory structure in the presence of noise, perturbations, and nonuniform sampling.
Method
The paper proves Synchrosqueezing stability, reformulates it in a numerically viable discretized form with parameter guidance, and applies it to paleoclimate records.
Results
Synchrosqueezing clearly delineates orbital cycles and provides a greatly improved representation of orbital signals in δ18O records over the past 2.5 million years.
Takeaways & Limitations
The method supports high-precision spectral analysis and decomposition of signals and provides insight into climate evolution over the past 2.5 million years.

Abstract

from arXiv · show

We analyze the stability properties of the Synchrosqueezing transform, a time-frequency signal analysis method that can identify and extract oscillatory components with time-varying frequency and amplitude. We show that Synchrosqueezing is robust to bounded perturbations of the signal and to Gaussian white noise. These results justify its applicability to noisy or nonuniformly sampled data that is ubiquitous in engineering and the natural sciences. We also describe a practical implementation of Synchrosqueezing and provide guidance on tuning its main parameters. As a case study in the geosciences, we examine characteristics of a key paleoclimate change in the last 2.5 million years, where Synchrosqueezing provides significantly improved insights.

I. INTRODUCTION

The paper develops Synchrosqueezing as a reconstructive time-frequency method for signals with time-varying oscillatory components, establishes robustness to several errors, and applies it to paleoclimate records. Its analysis and implementation support component extraction, discrete approximations, and more precise time-frequency representations of orbital climate signals.

Signal model: Synchrosqueezing decomposes signals into Fourier-like oscillatory components with slowly varying amplitudes and sufficiently smooth, time-varying frequencies.The components are modeled as f_k(t) = A_k(t) cos(2πφ_k(t)).
Method motivation: Unlike classical time-frequency reassignment methods, Synchrosqueezing can reconstruct the extracted components.Standard time-frequency reassignment sharpens spectrograms but does not allow synthesis of the components.
Stability: The paper proves stability under bounded deterministic perturbations and Gaussian white noise, supporting use with acquisition, quantization, interpolation, and sampling errors.The stability theory also implies that discrete approximations, including nonuniformly sampled data handled with splines, can replace the continuous signal under stated conditions.
Implementation and evaluation: The implementation reformulates the approach in discretized form, supplies parameter-selection guidance, and provides a freely available MATLAB toolbox.The paper also evaluates the algorithm on numerical test cases against established time-frequency and reassignment techniques.
Component extraction: Synchrosqueezing identifies and extracts components while enabling reconstruction through a partial, then completed, inversion of the continuous wavelet transform.The transform concentrates representation around instantaneous-frequency curves rather than inverting across the entire time-scale plane.

III. IMPLEMENTATION OVERVIEW

The implementation proceeds from a sampled signal through a discretized wavelet transform and phase transform to a reassigned time-frequency representation. Thresholding unreliable low-magnitude coefficients and using FFT-based calculations supports practical computation.

A. DWT of sampled signal: The method starts with a uniformly sampled, padded signal and computes its discretized continuous wavelet transform using a selected mother wavelet.The signal has n = 2L+1 samples, is padded to reduce boundary effects, and uses a wavelet concentrated around a positive frequency.
B. Phase transform: The phase transform is computed from the wavelet transform and its time derivative, with Fourier methods providing a practical derivative estimate.The resulting normalized frequency information supports reassignment from the time-scale plane to the time-frequency plane.
B. Phase transform: A hard threshold γ excludes points where the wavelet magnitude is too small for a trustworthy phase estimate.This addresses phase instability when |Wf| is approximately zero; γ = 10^-5 is used in the example figure.
B. Phase transform: The reassignment map sends each time-scale location (a, b) to a frequency-time location (ω(a, b), b), where ω supplies frequency information for recovery.The phase transform is not generally the instantaneous frequency itself, but it is used to recover actual frequencies in the next step.

C. Synchrosqueezing in the time-frequency plane: Tf(ω, b)

Synchrosqueezing constructs a discrete frequency-plane representation by assigning phase-transform values to frequency bins. Logarithmic frequency divisions respect sampling bounds and enable an efficient implementation.

C. Synchrosqueezing in the time-frequency plane: The discretized Synchrosqueezing transform assigns each phase-transform value to the frequency bin containing its reassigned frequency.The discrete definition replaces the continuous frequency variable with intervals W_l and corresponds to the limiting case as δ → 0.
C. Synchrosqueezing in the time-frequency plane: The frequency range is bounded by the sampling interval and signal duration, with maximum frequency 1/(2∆t) and minimum frequency 1/(n∆t).These bounds determine the divisions used for the frequency plane.
C. Synchrosqueezing in the time-frequency plane: The implementation uses logarithmically spaced frequency divisions defined by w_l = 2^(l∆w)w_0.The spacing is chosen consistently with the logarithmic scale discretization of the wavelet transform.
C. Synchrosqueezing in the time-frequency plane: The complete discretized Synchrosqueezed plane requires O(n_an) computations after the reassigned values and frequency bins are available.The algorithm finds the associated bin for each scale-time location and adds the value to the corresponding sum.
C. Synchrosqueezing in the time-frequency plane: Linear frequency spacing produces no significant practical differences from the logarithmic spacing used in the implementation.The logarithmic approach is preferred because it matches the CWT discretization and is computationally efficient.

D. Component reconstruction

Component reconstruction integrates the Synchrosqueezed transform over a narrow frequency band surrounding each component’s curve. The band may be selected manually or estimated by ridge extraction.

D. Component reconstruction: Each component f_k is recovered by integrating the Synchrosqueezed transform over frequencies associated with that component.This parallels filtering within a conventional time-frequency plot.
D. Component reconstruction: The relevant frequency indices L_k(t_m) form a small band around the kth component’s curve in phase-transform space.Theoretical results guide the choice of this band.
D. Component reconstruction: The component band can be selected manually or estimated with a standard least-squares ridge extraction method.The Synchrosqueezing Toolbox uses the ridge-extraction approach.
D. Component reconstruction: For real-valued components, the reconstruction uses the CWT normalization constant R_ψ.The reconstruction formula exploits the fact that f_k is real.

E. Selecting the threshold γ

The threshold γ determines which wavelet coefficients are trusted for phase estimation. The paper recommends a noise-adaptive choice based on the finest-scale coefficient magnitudes, while noiseless data can use machine epsilon.

E. Selecting the threshold γ: The hard threshold γ sets the lowest CWT magnitude at which the phase estimate ω is considered trustworthy.It also acts as a hard threshold on the wavelet representation by shrinking small coefficients to zero.
E. Selecting the threshold γ: For an uncorrupted signal, γ can be set to machine epsilon, with 10^-8 suggested for double-precision floating-point systems.This setting applies in the ideal noiseless case.
E. Selecting the threshold γ: Wavelet thresholding is related to established denoising procedures that threshold corrupted-signal coefficients before inversion.Prior work described both soft- and hard-thresholding and reported near-optimal root mean square error for its estimator.
E. Selecting the threshold γ: The noise power needed for threshold selection can be estimated from the wavelet representation.The paper follows prior work in using this representation for noise-power estimation.
E. Selecting the threshold γ: For noisy signals, the suggested threshold uses the Median Absolute Deviation of coefficient magnitudes at the nv finest scales.The factor 1.4826 relates Gaussian MAD to standard deviation, and the first octave supplies the finest-scale coefficients.

IV. NUMERICAL SIMULATIONS

The numerical simulations illustrate Synchrosqueezing’s analysis, extraction, and comparison with established time-frequency transforms, using scripts available for reproducing the figures.

The simulations illustrate Synchrosqueezing’s theoretical ideas and compare it with several time-frequency transforms.
MATLAB scripts used to generate the numerical-example figures are publicly available.

A. Comparison of Synchrosqueezing with the CWT, STFT and EEMD

On noisy multicomponent signals, Synchrosqueezing more clearly resolves time-varying frequencies than the CWT and STFT, while enabling component reconstruction comparable to those methods and better than EEMD in the reported example.

The synthetic signal contains three components with time-varying amplitude or frequency and additive Gaussian white noise with power σ2 = 2.4, producing an SNR of −2.6 dB.
Synchrosqueezing’s selected parameters balance time and frequency resolution, while shorter windows or wavelets increase temporal resolution at the cost of frequency resolution and greater component smearing.
Synchrosqueezing clearly distinguishes time-varying frequencies, whereas the STFT and CWT show more smearing and distortion.The nearly harmonic third component is identified by both Synchrosqueezing and STFT, but its frequency estimate is more precise in Synchrosqueezing.
All three time-frequency methods reconstruct component s2 reasonably accurately, although noise slightly smooths its amplitude modulation around t ∈[5, 7].The reconstruction is obtained by inverting each transform in a small band around s2’s curve.
EEMD shows poor amplitude recovery and drifting phase even without noise, while Synchrosqueezing and the other time-frequency methods produce good results under suitable bandwidth conditions.The cited comparison attributes EMD/EEMD sensitivity to amplitude changes and strong frequency-separation requirements.

B. Comparison of Synchrosqueezing with Reassignment Techniques

The reassignment comparison evaluates Synchrosqueezing against reassigned spectrogram and smoothed pseudo-Wigner–Ville methods under increased noise, highlighting differences in frequency reassignment and reconstruction.

At noise power σ2 = 5, corresponding to −5.8 dB SNR, Synchrosqueezing’s time-frequency plot contains fewer spurious components than the other reassignment plots.
Unlike standard reassignment, Synchrosqueezing reassigns only in frequency, preserving the time coordinate and allowing signal reconstruction.
The reassigned Wigner–Ville plot contains an extra curve between the second and third components from quadratic cross-terms.

C. Nonuniform Samples and Spline Fitting

For nonuniform samples, the simulations interpolate perturbed sampling times with a cubic spline before applying Synchrosqueezing; the resulting errors remain localized and mainly affect the highest-frequency component.

Nonuniform sample spacing is handled by fitting a cubic spline and discretizing the interpolant on a finer grid.The discretized signal equals the original signal plus a spline error term e(t).
Approximately 165 samples on [2, 8] yield an average sampling rate of 27.2, about three times the maximum instantaneous frequency of 9.85.
The spline interpolant approximates the original signal closely, with deviations concentrated in a few oscillations for t > 7.3 at the highest frequencies.
Synchrosqueezing is largely unaffected by interpolation errors and introduces no spurious spectral information in the time-frequency plot.
The localized errors affect amplitude-modulation recovery mainly for the third, highest-frequency component, which is the most difficult to recover.
Components near the Nyquist frequency are picked up fairly accurately when the mother wavelet follows Theorem II.1 and components are sufficiently separated.

D. Invariance to the underlying transform

Synchrosqueezing is largely invariant to the mother wavelet, with practical differences driven mainly by each wavelet’s time-frequency concentration. Among the tested wavelets, the bump wavelet gives the most accurate representation.

Synchrosqueezing is largely invariant to the mother wavelet’s precise shape.Observed differences instead reflect relative concentration in time and frequency, including distance from zero frequency.
The bump wavelet ψc provides the most accurate representation among the tested choices.Its frequency support is the smallest and exactly positive and finite rather than approximately so.
Figure 6 compares the wavelet spectrum, CWT, and Synchrosqueezing transform across the three mother-wavelet choices.Columns correspond to ψa, ψb, and ψc; rows show |2 bψ(4ξ)|, |Wfs|, and |Tfs|.

V. ASPECTS OF THE MID-PLEISTOCENE TRANSITION

Synchrosqueezing gives sharper, more precise time-periodicity representations of solar forcing and benthic δ18O records than CWT or STFT. This reveals changing climate periodicities across the mid-Pleistocene transition, including the emergence of a dominant late-Pleistocene 100 kyr component.

Solar forcing and climate response: Synchrosqueezing confirms solar-forcing precession cycles at 19 kyr and 23 kyr, obliquity cycles at 41 kyr and 54 kyr, and weak eccentricity cycles at 95 kyr, 124 kyr, and 400 kyr.Its spectral structure is more concentrated along the periodicity direction than the CWT.
Solar forcing and climate response: Synchrosqueezing provides a far more precise decomposition of DSDP607 and stacked benthic δ18O records than CWT or STFT.The stacks also reduce noise relative to individual cores, especially at finer scales and higher frequencies.
Solar forcing and climate response: 23 kyr precession becomes mostly coherent over the last 1 Myr, particularly in the stacked records compared with CWT decompositions.The stacks show less stochasticity above the obliquity band and less noisy time-periodicity evolution overall.
Mid-Pleistocene transition: The 41 kyr obliquity cycle persists through nearly the entire Pleistocene in the global climate response.Synchrosqueezing precisely reveals modulated signals that rise above stochastic variability.
Mid-Pleistocene transition: A lower-frequency signal near 70 kyr begins about 1.2 Myr ago and gradually evolves into the dominant 100 kyr component in the late Pleistocene.The analysis indicates this is not a direct cause-and-effect response to eccentricity variability.
Mid-Pleistocene transition: Synchrosqueezing shows increasing climate-response variability from the precession toward eccentricity bands after approximately 0.6 Myr.Internal feedback mechanisms, likely associated with long-term global cooling, amplify the eccentricity-band response after approximately 1.2 Myr.
Mid-Pleistocene transition: Cross-band differences indicate that superposed precession cycles can modulate lower-frequency climate responses as the climate becomes colder and potentially more nonlinear.These results address unresolved questions about the emergence of the strong 100 kyr cycle and earlier precession-band incoherence.
Implications: The analysis improves spectral tracing, noise characterization, and reconstruction of climate-signal components within selected frequency bands.These contributions support more accurate inversion across insolation and climate evolution.

VI. CONCLUSIONS AND FUTURE DIRECTIONS

The paper concludes that Synchrosqueezing analyzes and decomposes diverse signals with high time-frequency precision while remaining stable against signal errors. Its paleoclimate application provides further insight into climate evolution over the past 2.5 million years.

Synchrosqueezing can spectrally analyze and decompose a wide variety of signals with high precision in time and frequency.
An efficient implementation runs in O(nvn log2 n) time and is stable against errors in the signals.
The method provides further insight into climate evolution over the past 2.5 million years.
The authors report applications in climate dynamics, meteorology, oceanography, ECG analysis, mechanical-transmission analysis, and automated trading systems.