Wavelet Scattering Spectra

Wavelet Sparsity Factors

The wavelet sparsity factors are the ratios of the first and second moments of the first-order wavelet scalogram coefficients. The wavelet sparsity factors are defined as:

where N is the number of wavelet coefficients, k is the wavelet filter index in the first filter bank, and $${\sigma }_{1,k}$$ is the biased standard deviation of the wavelet coefficients.

By default, the scatteringSpectra function returns sparsity factors that include a term for the lowpass filter φ(t):

where σ_φ is the root mean square (RMS) instead of the standard deviation. Because the scaling function is real-valued, the absolute value is not necessary.

By using the RMS, the scatteringSpectra function protects against edge cases where the sparsity coefficient corresponding to the lowpass filter would otherwise blow up in value. Consider the following equality:

where μ is the mean and $${\sigma }_z^2$$ is the biased estimate of the standard deviation:

Because the mean of the wavelet coefficients is zero, the RMS and the biased standard deviation are equal. However, for φ(t), the scaling function, the situation is different. For the constant signal, the standard deviation is zero but the mean is equal to the value of each element in the signal. Therefore, the denominator σ_φ is not zero.

Variance (Wavelet Spectrum)

The variance, or wavelet spectrum, is defined as:

These features tend to differentiate signals with differing long-range dependencies. By default, there is an additional term corresponding to the lowpass filter:

Phase-Modulus Cross Spectrum

The phase-modulus cross spectrum is the correlation between the wavelet coefficients and their modulus:

where the filter indices satisfy l > k. Morel et al. [7] show that taking the modulus enables you to obtain some correlation with wavelet coefficients at a different scale. The power spectrum of $$\left|x\ast {\psi }_{j-1}\right|$$ overlaps the power spectrum of $$x\ast {\psi }_j$$.

As implemented, for a given $${\psi }_{1,k}$$, the scatteringSpectra function chooses those wavelets $${\psi }_{1,l}$$ such that 1/2 the 3-dB bandwidth of $${\psi }_{1,k}$$ overlaps the center frequency of $${\psi }_{1,l}$$ minus 1/2 its 3-dB bandwidth. For more information, see Visualize Wavelet Filter Bandwidth Criterion.

The phase-modulus coefficients tend to differentiate signals based on their skewness. In other words, the phase modulus is sensitive to time asymmetries in the signal. Note that the coefficients involve only wavelets from the first filter bank, which can potentially contain more wavelets than the second filter bank.

Modulus-Modulus Correlation (Scattering Cross-Spectrum)

The modulus-modulus correlations, or scattering cross-spectrum, are defined as

where the filter indices satisfy k ≥ k', l − k ≠ 0, and l − k' ≠ 0. The scattering cross-spectrum depends on the wavelet filters in the first and second filter banks. These coefficients quantify envelope dependencies across scales because they capture correlations between the envelopes of the wavelet coefficients, $$\left|x\ast {\tilde{\psi }}_k\right|$$.

The scatteringSpectra function applies a bandwidth criterion similar to the one used for the phase-modulus cross spectrum. However, in this case, the criterion involves filters in the first and second filter banks. For a filter $${\psi }_{1,k}$$, the function chooses those filters $${\psi }_{2,l}$$ such that 1/2 the 3-dB bandwidth of $${\psi }_{1,k}$$ overlaps the center frequency of $${\psi }_{2,l}$$ minus 1/2 its 3-dB bandwidth.

Similar to the situation with the phase-modulus cross spectrum, wavelet filter indices must satisfy l > k and l > k' when adjusted for the possibility of different quality factors. For example, suppose both filter banks in the wavelet time scattering network have the same Q factor. Then

is valid but scatteringSpectra does not return its conjugate (Hermitian adjoint term). The function also only pairs lowpass filters with lowpass filters: $$\mathbb{E}\left\{\frac{\left|x\ast {\tilde{\psi }}_{1,k}\right|}{{\sigma }_{1,k}}\ast \varphi ,\frac{\left|x\ast {\tilde{\psi }}_{1,k^{\prime }}\right|}{{\sigma }_{1,k^{\prime }}}\ast \varphi \right\}.$$ For more information, see Network Quality Factors and Scattering Spectra.

Wavelet Phase Harmonics

Other strategies for correlating wavelet coefficients at different scales involve using nonlinear transformations to overcome the issue whereby wavelets at different scales might have non-overlapping spectra. One such approach is to compute the wavelet phase harmonics. For a complex number $$z$$, the phase harmonics are computed as $${[z]}^k\to |z|e^{ik\varphi (z)}$$, where $$\varphi (z)$$ is the phase angle of $$z$$ and $$k$$ is a positive integer. Note that the phase harmonics operator does not alter the modulus.

Computing wavelet phase harmonics amounts to a nonlinear dilation in the Fourier domain without affecting the temporal (or spatial) localization of the wavelet. See [5] for the theoretical development of wavelet phase harmonics and [1], [9], and [12] to see wavelet phase harmonics applied to data.

Demonstrate how the phase harmonic operator can align phases for wavelet coefficients at different scales. For reproducibility, initialize the random seed to the default value. Create a random signal with 1024 samples. Create a wavelet scattering network appropriate for the signal. Each network filter bank has a quality factor of 1. In order to obtain the scattering spectra, set the filter downsampling criterion to "bandlimited".

rng("default")
sig = randn(1024,1);
tsn = waveletScattering(QualityFactors=[1 1], ...
    FilterDownsampling="bandlimited");

Use the scatteringSpectra object function to obtain the metadata associated with the scattering spectra of the signal. To obtain metadata associated with the phase harmonics, set the scattering spectra method to "phaseharmonic".

[~,metadataTable] = scatteringSpectra(tsn,sig,Method="phaseharmonic");

Use the filterbank object function to obtain the network filter banks. Extract the wavelet filters from the first filter bank.

psif = filterbank(tsn);
psift = psif{2}.psift;

The wavelet filters are defined in the Fourier domain. Use the filters to obtain the wavelet transform of the signal.

wavcfs = ifft(fft(sig).*psift);

In the scattering spectra, the wavelet phase harmonic coefficients are defined as $$\mathbb{E}\left\{\frac{{[x\ast {\psi }_{1,k}]}^p}{{\sigma }_{1,k}},\frac{x\ast {\psi }_{1,l}}{{\sigma }_{1,l}}\right\}$$, where $$x$$ is the signal, $$p$$ is the order of the harmonic, and $${\sigma }_{1,k}$$ and $${\sigma }_{1,l}$$ are the RMS of the wavelet coefficients at scales $$k$$ and $$l$$, respectively.

From the scattering spectra metadata table, extract the metadata associated with the phase harmonic coefficients. Extract the table variables jl1, jr1, and harmonic. The variables correspond to $$k$$, $$l$$, and $$p$$, respectively, in the inner product.

ind = ismember(metadataTable.type,"ph");
phTable = metadataTable(ind,["jl1" "jr1" "harmonic"])

phTable = 25×3 table
    jl1    jr1    harmonic
    ___    ___    ________

     2      1        2    
     3      2        2    
     4      3        2    
     5      4        2    
     6      5        2    
     7      6        2    
     3      1        3    
     4      2        3    
     5      3        3    
     6      4        3    
     7      5        3    
     3      1        4    
     4      2        4    
     5      3        4    
     6      4        4    
     7      5        4    
      ⋮

Choose a row. Obtain the wavelet coefficients $$x\ast {\psi }_{1,k}$$ and $$x\ast {\psi }_{1,l}$$, where $$k$$ equals jl1 and $$l$$ equals jr1.

row = 1;
k = phTable.jl1(row);
l = phTable.jr1(row);

wavA = wavcfs(:,k);
wavB = wavcfs(:,l);

Let p equal the harmonic value of the chosen row. Determine the phase angles of the coefficients $$x\ast {\psi }_{1,k}$$. Then apply the phase harmonic operator of order $$p$$ to those coefficients.

p = phTable.harmonic(row);
theta = atan2(imag(wavA),real(wavA));
wavAph = abs(wavA).*exp(1i*theta*p);

Show how applying the phase harmonic operator affects the scale $$k$$ wavelet coefficients in the Fourier domain. Compare the magnitudes of the Fourier transforms of $$x\ast {\psi }_{1,k}$$ and $$x\ast {\psi }_{1,l}$$. Then compare the magnitudes of the Fourier transforms of $${[x\ast {\psi }_{1,k}]}^p$$ and $$x\ast {\psi }_{1,l}$$. The phase harmonic operator increases the overlap of the spectra of the level $$k$$ and $$l$$ coefficients without affecting the temporal localization of the wavelet.

t=tiledlayout(2,1);
nexttile
plot([abs(fft(wavA)) abs(fft(wavB))])
axis tight
legend("Scale k","Scale l")
nexttile
plot([abs(fft(wavAph)) abs(fft(wavB))])
axis tight
legend("Phase Harmonic(Scale k)","Scale l")
title(t,"Effect of Phase Harmonic Operator in Fourier Domain")