wsstridge

R2026b

Time-frequency ridges from wavelet synchrosqueezing

Syntax

fridge = wsstridge(sst)

[fridge,iridge] = wsstridge(sst)

[___] = wsstridge(sst,penalty)

[___] = wsstridge(___,f)

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

Description

fridge = wsstridge(sst) extracts the maximum energy time-frequency ridge in cycles per sample from the wavelet synchrosqueezed transform, sst. The sst input is the output of wsst using the analytic Morlet wavelet. Each ridge is a separate signal mode.

If you used a wavelet other than the default Morlet wavelet ("amor") in wsst, or a nondefault VoicesPerOctave value in wsst, you can use the optional continuous wavelet transform (CWT) filter bank parameters output from wsst and input that to wsstridge using the FilterBankParameters name-value argument. Alternatively, you can use the frequencies output of wsst directly in wsstridge.

[fridge,iridge] = wsstridge(sst) returns in iridge the row indices of sst. The row indices are the maximum time-frequency ridge at each sample. Use iridge to reconstruct the signal mode along a time-frequency ridge using iwsst.

[___] = wsstridge(sst,penalty) multiplies the squared distance between frequency bins by the penalty value. You can include any of the output arguments from previous syntaxes.

[___] = wsstridge(___,f) returns the maximum energy time-frequency ridge in cycles per unit time based on the f input frequency vector. f is the frequency output of wsst. The f input and fridge output have the same units.

[___]= wsstridge(___,Name=Value) returns the time-frequency ridge with additional options specified by one or more Name=Value arguments.

Examples

Extract Time-Frequency Ridge from Chirp Signal

Obtain the wavelet synchrosqueezed transform of a quadratic chirp and extract the maximum time-frequency ridge, in fridge, and the associated row indices, in iridge.

Load the chirp signal and obtain its synchrosqueezed transform.

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

Extract the maximum time-frequency ridge.

[fridge,iridge] = wsstridge(sst);

Plot the synchrosqueezed transform.

pcolor(tquad,f,abs(sst))
shading interp
title("Synchrosqueezed Transform")

Overlay the plot of the maximum energy frequency ridge.

hold on
plot(tquad,fridge,"w")
hold off
title("Synchrosqueezed Transform with Overlaid Ridge")

Extract Time-Frequency Ridge from Multicomponent Signal

Extract the two highest energy modes from a multicomponent signal using two methods:

• Use the frequency vector output of the wsst function.

• Use the CWT filter bank parameters output of wsst.

Load a multicomponent signal. Obtain the CWT of the signal using the analytic Morlet wavelet.

load multicompsig
sig = sig1+sig2;

Method 1 — Frequency Vector

Obtain the wavelet synchrosqueezed transform of the signal. Also obtain the frequency vector output. Plot the output.

[sst,F] = wsst(sig,sampfreq);
contour(t,F,abs(sst))
xlabel("Time")
ylabel("Hz")
grid on
title("Synchrosqueezed Transform of Two-Component Signal")

Using a penalty of 10, extract the two highest energy modes and plot the result.

[fridge,iridge] = wsstridge(sst,10,F,NumRidges=2);
hold on
plot(t,fridge,"k",LineWidth=2)
hold off

Method 2 — CWT Filter Bank Parameters

Obtain the wavelet synchrosqueezed transform of the signal. Also obtain the filter bank parameters. Plot the output.

[sst,~,fbparams] = wsst(sig,sampfreq);
contour(t,F,abs(sst))
xlabel("Time")
ylabel("Hz")
grid on
title("Synchrosqueezed Transform of Two-Component Signal")

Using a penalty of 10, extract the two highest energy modes using the filter bank parameters. Plot the result.

[fridge,iridge] = wsstridge(sst,10,NumRidges=2, ...
    FilterbankParameters=fbparams);
contour(t,F,abs(sst))
xlabel("Time")
ylabel("Hz")
grid on
title("Using Filter Bank Parameters")
hold on
plot(t,fridge,"k",LineWidth=2)
hold off

Input Arguments

sst — Synchrosqueezed transform

matrix

Synchrosqueezed transform, specified as a matrix. sst is a time-frequency matrix and is the output of wsst.

Data Types: single | double

penalty — Frequency bins scaling penalty

0 (default) | nonnegative scalar

Frequency bins scaling penalty, specified as a nonnegative scalar. This input penalizes changes in frequency by multiplying the penalty value by the squared distance between frequency bins. Use a penalty term when you extract multiple ridges, or when you have a single modulated component in additive noise. The penalty term prevents jumps in frequency that occur when the region of highest energy in the time-frequency plane changes abruptly.

f — Synchrosqueezed transform frequencies

vector

Synchrosqueezed transform frequencies corresponding to the rows of the synchrosqueezed transform, which is the vector output of wsst. The number of elements in the frequency vector is equal to the number of rows in the sst input.

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: NumRidges=3

FilterBankParameters — CWT filter bank parameters

structure array

Since R2024a

CWT filter bank parameters, specified as a structure array. The filter bank parameters are the optional structure array output fbparam of wsst. The wsstridge function uses the parameters to determine the frequencies used in the construction of the time-frequency ridges. Specifying FilterBankParameters is only valid if you do not use the input frequency vector, f.

NumRidges — Number of highest energy time-frequency ridges

1 (default) | positive integer

Number of highest energy time-frequency ridges to extract, specified as a positive integer. If this integer is greater than 1, wsstridge iteratively determines the maximum energy time-frequency ridge by removing the previously computed ridges and the default or specified NumFrequencyBins on either side of each ridge bin.

NumFrequencyBins — Number of frequency bins to remove

4 (default) | positive integer

Number of frequency bins to remove from synchrosqueezed transform sst when extracting multiple ridges, specified as a positive integer. This integer must be less than or equal to round(size(sst,1)/4). You can specify the number of frequency bins to remove only if you extract more than one ridge. After extracting the highest energy time-frequency ridge, wsstridge removes the sst values corresponding to the iridge indices at each time step. The energy is removed along the time-frequency ridge extended on both sides of the iridge index by the specified number of frequency bins. If the index of the extended time-frequency ridge exceeds the number of frequency bins at any time step, wsstridge truncates the removal region at the first or last frequency bin. To specify 'NumFrequencyBins', you must specify 'NumRidges'.

Output Arguments

fridge — Time-frequency ridge frequencies

vector or matrix

Time-frequency ridge frequencies, returned as a vector or matrix. The frequencies correspond to the time-frequency ridge at each time step. fridge is an N-by-nr matrix where N is the number of time samples (columns) in sst and nr is the number of ridges. The first column of the matrix contains the frequencies for the maximum energy time-frequency ridge in sst. Subsequent columns contain the frequencies for the time-frequency ridges in decreasing energy order. By default, fridge contains frequencies in cycles per sample.

Data Types: single | double

iridge — Time-frequency ridge indices

vector or matrix

Time-frequency ridge row indices of sst, returned as a vector or matrix. The row indices in iridge correspond to the row index of the maximum time-frequency ridge for each sst column. iridge is an N-by-nr matrix where N is the number of time samples (columns) in sst, and nr is the number of ridges. The first column of the matrix contains the indices for the maximum energy time-frequency ridge in sst. Subsequent columns contain the indices for the time-frequency ridges in decreasing energy order.

Data Types: single | double

Algorithms

The function uses a penalized forward-backward greedy algorithm to extract the maximum-energy ridges from a time-frequency matrix. For unpenalized ridge extraction —where the penalty for changing frequency is 0— the function instead selects the maximum energy at each time point.

The penalized forward-backward greedy algorithm finds the maximum time-frequency ridge by minimizing –ln S at each time point. Minimizing –ln S is equivalent to maximizing the value of S. The algorithm optionally constrains jumps in frequency with a penalty that is proportional to the distance between frequency bins.

Assume a real-valued matrix E(j,i) = – ln(S(j,i)), where j = 1, 2, …, M indexes frequency bins and i = 1, 2, …, N indexes time bins. S can be, for instance, the absolute value of the Fourier or wavelet synchrosqueezed transform of a signal. The algorithm seeks a curve γ(i) that minimizes the quadratic regularized cost

$$J(\gamma )={\displaystyle \sum _{i=1}^{N}E}(\gamma (i),i)+\lambda {\displaystyle \sum _{i=2}^N{(\gamma (i)-\gamma (i-1))}^2}$$,

where the first term enforces the curve to follow low-energy bins and the second term penalizes large frequency jumps. The parameter λ > 0 controls smoothness.

Let the state cost be F(j,i) = – ln(S(j,i)), which is the minimum total cost of all paths ending at bin j at time i. The algorithm performs these steps:

1. Initialization: $$F(j,1)=E(j,1)$$

2. Recursion: $$F(j,i)=E(j,i)+\underset{k=1,\dots ,M}{\min }F(k,i-1)+\lambda {(j-k)}^2$$

3. Backtracking: $$\{\begin{array}{l}I(j,i)=\min (\arg \underset{k=1,\dots ,M}{\min }(F(k,i-1)+\lambda {(}^j)\\ \gamma (i-1)=I(\gamma (i),i),i=N,N-1,\dots ,2\end{array}$$

The following example illustrates the time-frequency ridge algorithm using a penalty that is two times the distance between frequency bins. Specifically, the distance between the elements (j,k) and (m,n) is defined as (j-m)2. The time-frequency matrix has three frequency bins and three time steps. The matrix columns correspond to time steps, and the matrix rows correspond to frequency bins. The values in the second row represent a sine wave.

1. Suppose you have the matrix:

   1   4   4
   2   2   2
   5   5   4

2. Update the value for the (1,2) element as follows.

   1. Leave the values at the first time point unaltered. Begin the algorithm with the (1,2) element of the matrix, which presents the first frequency bin at the second time point. The bin value is 4. Penalize the values in the first column based on their distance from the (1,2) element. Applying the penalty to the first column produces

      originalValues + penalty × distance = penalizedValues
              1      +       2 × 0        =  1
              2      +       2 × 1        =  4
              5      +       2 × 4        = 13
      

      The minimum value of the penalizedValues vector is 1, which is from the frequency bin (row) 1.

      Note

      If there are two or more rows containing the minimum value of the penalizedValues vector, the algorithm chooses the frequency bin with the highest row index that satisfies the minimum value criterion because the algorithm loops down the rows and updates the minimum as it goes. While ties are possible to occur in synthetic data, such ties are unlikely for data coming from real-world time-frequency data.

   2. Add the minimum value in column 1 to the current bin value, 4. The updated value for (1,2) becomes 5, which came from bin 1.

3. Update the values for the remaining elements in column 2 as follows.

   Recompute the original column 1 values with the penalty factor using the same process as in Step 2a. Obtain the remaining second column values using the same process as in Step 2b. For example, when updating the (2,2) element, which has bin value 2, applying the penalty to the column yields

   originalValues + penalty × distance = penalizedValues
           1      +       2 × 1        =  3
           2      +       2 × 0        =  2
           5      +       2 × 1        =  7
   

   Add the minimum value, 2, to the current bin value. The updated value for (2,2) becomes 4. After updating the (3,2) element, the matrix is

   1   5(1)  4
   2   4(2)  2
   5   9(2)  4

   Only the second column has been updated. The subscripts indicate the index of the bin in the previous column from which a value came.

4. Repeat Step 2 for the third column. But now the penalty is applied to the updated second column. For example, when updating the (1,3) element, the penalized values are

   5 + 2 × 0 =  5
   4 + 2 × 1 =  6
   9 + 2 × 4 = 17
   

   The minimum value, 5, which is in the first bin, is added to the (1,3) bin value. After updating all the values in the third column, the final matrix is

   1   5(1)   9(1)
   2   4(2)   6(2)
   5   9(2)  10(2)

5. Starting at the last column of the matrix, find the minimum value. Walk back in time through the matrix by going from the current bin to the origin of that bin at the previous time point. Keep track of the bin indices, which form the path composing the ridge. The algorithm smooths the transition by using the origin bin instead of the bin with the minimum value. For this example, the ridge indices are 2, 2, 2, which match the energy path of the sine wave in row 2 of the matrix shown in Step 1.

If you are extracting multiple ridges, the algorithm removes the first ridge from the time-frequency matrix and repeats the process.

References

[1] Daubechies, I., J. Lu, and H.-T. Wu. "Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool." Applied and Computational Harmonic Analysis. Vol. 30, Number 2, 2011, pp. 243–261.

[2] Thakur, G., E. Brevdo, N. S. Fučkar, and H.-T. Wu. "The Synchrosqueezing algorithm for time-varying spectral analysis: Robustness properties and new paleoclimate applications." Signal Processing. Vol. 93, Number 4, 2013, pp. 1079–1094.

Extended Capabilities

Thread-Based Environment

Run code in the background using MATLAB® backgroundPool or accelerate code with Parallel Computing Toolbox™ ThreadPool.

The wsstridge function fully supports thread-based environments. For more information, see Run MATLAB Functions in Thread-Based Environment.

Version History

Introduced in R2016a

R2026b: Use thread-based environments

You can use wsstridge in thread-based environments. For more information, see Run MATLAB Functions in Thread-Based Environment.

R2024a: wsstridge supports analytic Morse wavelets

The wsstridge function supports analytic Morse wavelets. To use a Morse wavelet, set the FilterBankParameters name-value argument.