timeSpectrum

Time-averaged wavelet spectrum

Since R2020b

Syntax

tavgp = timeSpectrum(fb,x)

tavgp = timeSpectrum(fb,cfs)

[tavgp,f] = timeSpectrum(___)

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

timeSpectrum(___)

Description

tavgp = timeSpectrum(fb,x) returns the time-averaged wavelet power spectrum of the signal x using the continuous wavelet transform (CWT) filter bank fb. By default, tavgp is obtained by time-averaging the magnitude-squared scalogram over all times. The power of the time-averaged wavelet spectrum is normalized to equal the variance of x.

tavgp = timeSpectrum(fb,cfs) returns the time-averaged wavelet spectrum for the CWT coefficients cfs.

Note

When using this syntax, the power of the time-averaged wavelet spectrum is normalized to equal the variance of the last signal processed by the filter bank object function wt.

[tavgp,f] = timeSpectrum(___) returns the wavelet center frequencies or center periods for the time-averaged wavelet spectrum. f is a column vector or duration array depending on whether the sampling frequency or sampling period is specified in the CWT filter bank, fb.

[___] = timeSpectrum(___,Name,Value) specifies additional options using name-value pair arguments. These arguments can be added to any of the previous input syntaxes. For example, 'Normalization','none' specifies no normalization of the time-averaged wavelet spectrum.

timeSpectrum(___) with no output arguments plots the time-averaged wavelet power spectrum in the current figure.

Examples

Time-Averaged Wavelet Spectrum of Oceanic Eddy Data

Load the NPG2006 dataset [1]. The data is the trajectory of a subsurface float trapped in an eddy. Plot the eastward and northward displacement. The triangle marks the initial position.

load npg2006
plot(npg2006.cx)
hold on
plot(npg2006.cx(1),'^','markersize',11,'color','r', ...
    'markerfacecolor',[1 0 0 ])
hold off
grid on
xlabel('Eastward Displacement (km)')
ylabel('Northward Displacement (km)')

Create a CWT filter bank that can be applied to the data. Use the default Morse wavelet. The sampling period for the data is 4 hours.

fb = cwtfilterbank('SignalLength',length(npg2006.cx), ...
    'SamplingPeriod',hours(4));

Obtain the time-averaged wavelet power spectra and the center periods.

[tavgp,centerP] = timeSpectrum(fb,npg2006.cx);
size(tavgp)

ans = 1×3

    73     1     2

The first page is the time-averaged wavelet spectrum for the positive scales (analytic part or counterclockwise component), and the second page is the time-averaged wavelet spectrum for the negative scales (anti-analytic part or clockwise component). Plot both spectra.

tiledlayout(2,1)
nexttile
plot(centerP,tavgp(:,1,1))
title('Counterclockwise Component')
ylabel('Power')
xlabel('Period (hrs)')
nexttile
plot(centerP,tavgp(:,1,2))
title('Clockwise Component')
ylabel('Power')
xlabel('Period (hrs)')

If you omit the output arguments and execute timeSpectrum(fb,npg2006.cx) on the command line, the scalograms and time-averaged power spectra are plotted in the current figure. Note that the clockwise rotation of the float is captured in the clockwise rotary scalogram and the time-averaged spectrum.

Normalize Time-Averaged Wavelet Spectrum

Load a time series of solar magnetic field magnitudes recorded hourly over the south pole of the sun by the Ulysses spacecraft from 21:00 UT on December 4, 1993 to 12:00 UT on May 24, 1994. See [3] pp. 218–220 for a complete description of this data. Create a CWT filter bank that can be applied to the data. Plot the scalogram and the time-averaged wavelet spectrum.

load solarMFmagnitudes
fb = cwtfilterbank('SignalLength',length(sm),'SamplingPeriod',hours(1));
timeSpectrum(fb,sm)

Obtain the time-averaged wavelet spectrum of the signal using default values. By default, timeSpectrum normalizes the power of the time-averaged wavelet spectrum to equal the variance of the signal. Verify that the sum of the spectrum equals the variance of the signal.

tavg = timeSpectrum(fb,sm);
[var(sm) sum(tavg)]

ans = 1×2

    0.0448    0.0447

Obtain the time-averaged wavelet spectrum of the signal, but instead normalize the power as a probability density function. Verify that the sum is equal to 1.

tavg = timeSpectrum(fb,sm,'Normalization','pdf');
sum(tavg)

ans = 
1.0000

If you set SpectrumType to 'density', timeSpectrum normalizes the weighted integral of the wavelet spectrum according to the value of Normalization. The spectrum mimics a probability density function whose integral, numerically evaluated, equals the value specified by Normalization.

Plot the scalogram and the time-averaged wavelet spectrum with spectrum type 'density' and 'pdf' normalization.

figure
timeSpectrum(fb,sm,'SpectrumType','density', ...
    'Normalization','pdf')

To confirm the integral of the spectrum equals 1, first obtain the time-averaged wavelet spectrum with 'density' spectrum type and 'pdf' normalization.

tavg = timeSpectrum(fb,sm,'SpectrumType','density', ...
    'Normalization','pdf');

By default, the filter bank uses the analytic Morse (3,60) wavelet. Obtain the admissibility constant for the wavelet and numerically integrate the wavelet spectrum using the trapezoidal rule. Keep in mind that the CWT uses L1 normalization. Confirm that the integral equals 1.

ga = 3;
tbw = 60;

be = tbw/ga;
anorm = 2*exp(be/ga*(1+(log(ga)-log(be))));
cPsi = anorm^2/(2*ga).*(1/2)^(2*(be/ga)-1)*gamma(2*be/ga);

rawScales = scales(fb);
numInt = 2/cPsi*1/length(sm)*trapz(rawScales(:),tavg./rawScales(:))

numInt = 
1.0000

Input Arguments

fb — Continuous wavelet transform filter bank

cwtfilterbank object

Continuous wavelet transform (CWT) filter bank, specified as a cwtfilterbank object.

x — Input data

vector

Input data, specified as a real- or complex-valued vector. The input data x must have at least four samples.

Data Types: single | double
Complex Number Support: Yes

cfs — CWT coefficients

matrix | 3-D array

CWT coefficients, specified as a 2-D matrix or as an M-by-N-by-2 array. cfs should be the output of the wt object function of the CWT filter bank fb.

Data Types: single | double
Complex Number Support: Yes

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.

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

Example: timeSpectrum(fb,x,'TimeLimits',[100 500],'Normalization','none') returns the time-averaged wavelet spectrum averaged over the time limits specified in samples without normalizing the spectrum.

Normalization — Normalization

'var' (default) | 'pdf' | 'none'

Normalization of the time-averaged wavelet spectrum, specified as a comma-separated pair consisting of 'Normalization' and one of the following:

• 'var' — Normalize to equal the variance of the time series x. If you provide the cfs input, the timeSpectrum function uses the variance of the last time series processed by the filter bank object function wt.

• 'pdf' — Normalize to equal 1.

• 'none' — No normalization is applied.

SpectrumType — Type of wavelet spectrum

'power' (default) | 'density'

Type of wavelet spectrum to return, specified as a comma-separated pair consisting of 'SpectrumType' and either 'power' or 'density'. If specified as 'power', the averaged sum of the time-averaged wavelet spectrum over all times is normalized according to the value specified in 'Normalization'. If specified as 'density', the weighted integral of the wavelet spectrum over all times is normalized according to the value specified in 'Normalization'.

Note

With regards to the numerical implementation of the continuous wavelet transform, the integral over scale is performed using L1 normalization. With L1 normalization, if you have equal amplitude oscillatory components in your data at different scales, they will have equal magnitude in the CWT. Using L1 normalization provides a more accurate representation of the signal. For more information, see L1 Norm for CWT.

TimeLimits — Time limits

[1 length(x)] or [1 size(cfs,2)] (default) | two-element vector

Time limits over which to average the wavelet spectrum, specified in samples. TimeLimits is specified as a comma-separated pair consisting of 'TimeLimits' and a two-element vector with nondecreasing elements. When you specify the input data as a signal, the elements are between 1 and the length of x. When you specify the input data as CWT coefficients, the elements are between 1 and size(cfs,2).

Output Arguments

tavgp — Time-averaged wavelet power spectrum

real-valued vector | real-valued 3-D array

Time-averaged wavelet power spectrum, returned as a real-valued vector or real-valued 3-D array. If x is real-valued, tavgp is an F-by-1 vector, where F is the number of wavelet center frequencies or center periods in the CWT filter bank fb. If x is complex-valued, tavgp is an F-by-1-by-2 array, where the first page is the time-averaged wavelet spectrum for the positive scales (analytic part or counterclockwise component), and the second page is the time-averaged wavelet spectrum for the negative scales (anti-analytic part or clockwise component).

f — Center frequencies or center periods

column vector | duration array

Center frequencies or center periods for the time-averaged wavelet spectrum, returned as a column vector or duration array, respectively. If the sampling frequency is specified in fb, then the elements of f are the center frequencies ordered from high to low. If the sampling period is specified in fb, then the elements of f are the center periods.

References

[1] Lilly, J. M., and J.-C. Gascard. “Wavelet Ridge Diagnosis of Time-Varying Elliptical Signals with Application to an Oceanic Eddy.” Nonlinear Processes in Geophysics 13, no. 5 (September 14, 2006): 467–83. https://doi.org/10.5194/npg-13-467-2006.

[2] Torrence, Christopher, and Gilbert P. Compo. “A Practical Guide to Wavelet Analysis.” Bulletin of the American Meteorological Society 79, no. 1 (January 1, 1998): 61–78. https://doi.org/10.1175/1520-0477(1998)079<0061:APGTWA>2.0.CO;2.

[3] Percival, Donald B., and Andrew T. Walden. Wavelet Methods for Time Series Analysis. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge ; New York: Cambridge University Press, 2000.

[4] Lilly, J.M., and S.C. Olhede. “Higher-Order Properties of Analytic Wavelets.” IEEE Transactions on Signal Processing 57, no. 1 (January 2009): 146–60. https://doi.org/10.1109/TSP.2008.2007607.

Extended Capabilities

C/C++ Code Generation

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

Usage notes and limitations:

• Plotting is not supported.

GPU Arrays

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

The timeSpectrum function fully supports GPU arrays. To run the function on a GPU, specify the input data as a gpuArray (Parallel Computing Toolbox). For more information, see Run MATLAB Functions on a GPU (Parallel Computing Toolbox).

Version History

Introduced in R2020b