wsst

Wavelet synchrosqueezed transform

Syntax

sst = wsst(x)

[sst,f] = wsst(x)

[sst,f,fbparam] = wsst(x)

[___] = wsst(x,fs)

[___] = wsst(x,ts)

[___] = wsst(___,wav)

wsst(___)

[___] = wsst(___,Name=Value)

Description

sst = wsst(x) returns the wavelet synchrosqueezed transform of x. The wsst function computes the synchrosqueezed transform using the analytic Morlet wavelet.

The synchrosqueezed transform has reduced energy smearing when compared to the continuous wavelet transform (CWT). The wsst function normalizes the analyzing wavelets to preserve the L1 norm. For more information, see Algorithms.

[sst,f] = wsst(x) returns a vector of frequencies, f, in cycles per sample.

[sst,f,fbparam] = wsst(x) returns the CWT filter bank parameters used in the synchrosqueezed transform.

[___] = wsst(x,fs) returns the synchrosqueezed transform using the specified sampling frequency, fs.

[___] = wsst(x,ts) returns the synchrosqueezed transform using the specified sampling interval, ts.

[___] = wsst(___,wav) uses the analytic wavelet specified by wav to compute the synchrosqueezed transform.

wsst(___) with no output arguments plots the magnitude of the synchrosqueezed transform as a function of time and frequency. If you do not specify a sampling frequency or interval, the synchrosqueezed transform is plotted in cycles per sample. If you specify a sampling frequency, the synchrosqueezed transform is plotted in hertz. If you specify a sampling interval using a duration, the plot is in cycles per unit time. The time units are derived from the duration.

[___] = wsst(___,Name=Value) returns the synchrosqueezed transform with additional options specified by one or more name-value arguments.

Examples

Synchrosqueezed Transform of Speech Signal

Obtain the wavelet synchrosqueezed transform of a speech sample using default values.

load mtlb
sst = wsst(mtlb);

Synchrosqueezed Transform and Reconstruction of Speech Signal

Load a speech signal. The signal is sampled at $$F_s=7418$$ Hz. Obtain the synchrosqueezed transform.

load mtlb
% To hear, type sound(mtlb,Fs)
dt = 1/Fs;
t = 0:dt:numel(mtlb)*dt-dt;
[sst,f] = wsst(mtlb,Fs);

Plot the magnitude of the synchrosqueezed transform.

pcolor(t,f,abs(sst))
shading interp
xlabel("Seconds")
ylabel("Frequency (Hz)")
title("Synchrosqueezed Transform")

Obtain the inverse synchrosqueezed transform. Plot the original signal and the reconstruction.

xrec = iwsst(sst);
% To hear, type soundsc(xrec,Fs)
plot(t,mtlb)
hold on
plot(t,xrec)
hold off
xlabel("Seconds")
ylabel("Amplitude")
legend("Original","Reconstruction")

Compare Synchrosqueezed Transform with CWT

Obtain the synchrosqueezed transform of a two-component signal and compare with the CWT of the same signal. The input is a combination of an amplitude-modulated and a frequency-modulated signal.

Create a two-component signal. The signal is a combination of an amplitude-modulated and a frequency-modulated signal. Sample the signal for one-tenth of a second at a sample rate of 1000 Hz.

Fs = 1e3;
t = 0:1/Fs:0.1;
x1 = (2+0.5*cos(2*pi*10*t)).*cos(2*pi*200*t+10*sin(2*pi*5*t));
x2 = cos(2*pi*75*t);
sig = x1+x2;

Obtain the synchrosqueezed transform of the signal. To minimize boundary effects, extend the signal symmetrically. Also obtain the parameters of the underlying CWT filter bank. Plot the magnitude of the synchrosqueezed transform.

[sst,sstF,fbparams] = wsst(sig,Fs,ExtendSignal=true);
pcolor(t,sstF,abs(sst))
shading interp
colorbar
title("WSST Magnitude")
xlabel("Time (s)")
ylabel("Frequency (Hz)")

Create a CWT filter bank using the underlying filter bank parameters. Then use the filter bank to obtain the CWT of the signal. Plot the scalogram of the CWT and compare with the synchrosqueezed transform results.

cfb = cwtfilterbank(fbparams);
[cfs,cfsF] = cwt(sig,FilterBank=cfb);
pcolor(t,cfsF,abs(cfs))
colorbar
shading interp
title("Scalogram")
xlabel("Time (s)")
ylabel("Frequency (Hz)")

Synchrosqueezed Transform of Quadratic Chirp

Obtain the wavelet synchrosqueezed transform of a quadratic chirp. The chirp is sampled at 1000 Hz.

load quadchirp
[sst,f] = wsst(quadchirp,1000);

Plot the magnitude of the transform.

hp = pcolor(tquad,f,abs(sst));
hp.EdgeColor = "none";
title("Wavelet Synchrosqueezed Transform")
xlabel("Time (s)")
ylabel("Frequency (Hz)")

Synchrosqueezed Transform of Sunspot Data

Plot the wavelet synchrosqueezed transform of sunspot data using the analytic Morse (3,30) wavelet. Specify the sampling interval to be one year.

load sunspot
wsst(sunspot(:,2),years(1),"Morse",WaveletParameters=[3 30])

Synchrosqueezed Transform of Sunspot Data Using Bump Wavelet

Plot the wavelet synchrosqueezed transform of sunspot data using the bump wavelet. Specify the sampling interval to be one year.

load sunspot
wsst(sunspot(:,2),years(1),"bump")

Input Arguments

x — Input signal

vector

Input signal, specified as a vector. x must have at least four samples.

Data Types: single | double

fs — Sampling frequency

positive scalar

Sampling frequency in hertz, specified as a positive scalar.

ts — Sampling interval

duration scalar

Sampling interval, also known as the sampling period, specified as a duration scalar. ts is the time between samples of x. Valid durations are years, days, hours, minutes and seconds. You cannot specify both ts and fs.

Example: seconds(1) is a duration scalar representing a 1-second time difference between consecutive signal samples.

wav — Analytic wavelet

"amor" (default) | "Morse" | "bump"

Analytic wavelet used to compute the synchrosqueezed transform, specified as one of the following:

• "amor" — Analytic Morlet wavelet

• "Morse" — Analytic Morse wavelet (since R2024a)

• "bump" — Bump wavelet

For use in the inverse synchrosqueezed transform, the analytic Morlet or Morse (3,30) wavelet generally provide the best results.

Name-Value Arguments

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Example: VoicesPerOctave=26

Before R2021a, use commas to separate each name and value, and enclose Name in quotes.

Example: 'VoicesPerOctave',26

VoicesPerOctave — Number of voices per octave

32 (default) | even integer from 10 to 48

Number of voices per octave to use in the synchrosqueezed transform, specified as an even integer between 10 and 48. The number of scales is the product of the number of voices per octave and the number of octaves.

WaveletParameters — Morse wavelet parameters

[3,60] (default) | two-element vector of scalars

Since R2024a

Morse wavelet parameters, specified as a two-element vector. The first element is the symmetry parameter (gamma), which must be greater than or equal to 1. The second element is the time-bandwidth product, which must be greater than or equal to gamma. The ratio of the time-bandwidth product to gamma cannot exceed 40. Note the ratio of the time-bandwidth parameter to gamma is often referred to as the beta parameter in the literature on Morse wavelets.

When gamma is equal to 3, the Morse wavelet is perfectly symmetric in the frequency domain. The skewness is equal to 0. Values of gamma greater than 3 result in positive skewness, while values of gamma less than 3 result in negative skewness.

If you are using synchrosqueezing for mode reconstruction, the Morse (3,30) wavelet generally provides superior inverse synchrosqueezing results.

WaveletParameters is only valid if you specify the analytic Morse wavelet.

For more information, see Morse Wavelets.

ExtendSignal — Extend signal symmetrically

false or 0 (default) | true or 1

Extend signal symmetrically logical, specified as a numeric or logical 0 (false) or 1 (true). If specified as true, the wsst function extends the signal symmetrically. Extending the signal symmetrically can mitigate boundary effects.

Output Arguments

sst — Synchrosqueezed transform

matrix

Synchrosqueezed transform, returned as an Na-by-N matrix, where Na is the number of scales and N is the number of samples in x.

The number of scales is determined by the underlying CWT filter bank object using frequency limits from the minimum frequency to the Nyquist frequency (half the sample rate). You can use the frequency vector f and filter bank parameters fbparam to obtain the scales used in the synchrosqueezed transform.

Data Types: single | double

f — Frequencies

vector

Frequencies of the synchrosqueezed transform, returned as a vector. The frequencies correspond to the rows of the sst. The length of f and the number of rows in sst are equal.

• If you do not specify sampling frequency fs or sampling interval ts, the frequencies are in cycles per sample.

• If you specify a sampling frequency, the frequencies are in hertz.

• If you specify a sampling interval, the frequencies are in cycles per unit time, where the time unit is the same time unit as ts. The function uses ts to compute the scale-to-frequency conversion for f.

fbparam — Parameters

structure array

Since R2024a

Parameters of the underlying CWT filter bank used in the synchrosqueezed transform, returned as a structure array. The field names denote the properties of the CWT filter bank used in the synchrosqueezed transform. The values for each field are the property values.

You can use f and fbparam to obtain the scales used in the synchrosqueezed transform.

Algorithms

The wsst function normalizes the analyzing wavelets to preserve the L1 norm. An equivalent way to state this is that wsst does not multiply the Fourier transforms of the wavelet bandpass filters by the square root of the scale. Multiplying by the square root of the scale would unequally weight different bandpass contributions.

With L1 normalization, if you have equal amplitude oscillatory components in your data at different scales, they will have equal magnitude in the CWT. The cwt function also uses L1 normalization. For more information, see L1 Norm for CWT.

References

[1] Daubechies, Ingrid, Jianfeng Lu, and Hau-Tieng Wu. “Synchrosqueezed Wavelet Transforms: An Empirical Mode Decomposition-like Tool.” Applied and Computational Harmonic Analysis 30, no. 2 (March 2011): 243–61. https://doi.org/10.1016/j.acha.2010.08.002.

[2] Thakur, Gaurav, Eugene Brevdo, Neven S. Fučkar, and Hau-Tieng Wu. “The Synchrosqueezing Algorithm for Time-Varying Spectral Analysis: Robustness Properties and New Paleoclimate Applications.” Signal Processing 93, no. 5 (May 2013): 1079–94. https://doi.org/10.1016/j.sigpro.2012.11.029.

Extended Capabilities

C/C++ Code Generation

Generate C and C++ code using MATLAB® Coder™.

Usage notes and limitations:

• The input wavelet name wav must be constant at compile time. Use coder.Constant (MATLAB Coder).

• The wavelet synchrosqueezed transform algorithm is sensitive to single-precision data. As a result, there may be a mismatch between results obtained using MATLAB^® code and generated C/C++ code. The mismatches are due to off-by-one differences in frequency reassignments.

GPU Arrays

Accelerate code by running on a graphics processing unit (GPU) using Parallel Computing Toolbox™.

The wsst function supports GPU array input with these usage notes and limitations:

• The wavelet synchrosqueezed transform algorithm is sensitive to single-precision data. As a result, there may be a mismatch between results obtained on a CPU and GPU. The mismatches are due to off-by-one differences in frequency reassignments.

For more information, see Run MATLAB Functions on a GPU (Parallel Computing Toolbox).

Version History

Introduced in R2016a

R2024a: Obtain CWT filter bank parameters

You can obtain the parameters of the underlying CWT filter bank the wsst function uses in the synchrosqueezing transform.

R2024a: wsst supports analytic Morse wavelets

The wsst function supports analytic Morse wavelets.

R2024a: wsst supports C/C++ code generation and gpuArray objects

The wsst function supports:

• C/C++ code generation. You must have MATLAB Coder™ to generate C/C++ code.

• gpuArray objects. You must have Parallel Computing Toolbox™ to use this functionality.