Continuous 1-D wavelet transform
cwt for information on the older
version of the
cwt. The older version is no longer
wt = cwt(x)
wt = cwt(x,wname)
[wt,f] = cwt(___,fs)
[wt,period] = cwt(___,ts)
[wt,f,coi] = cwt(___,fs)
[wt,period,coi] = cwt(___,ts)
[___] = cwt(___,Name,Value)
[___,coi,fb] = cwt(___)
[___,fb,scalingcfs] = cwt(___)
returns the continuous wavelet transform (CWT) of
wt = cwt(
x, is a double-precision real- or complex-valued
vector, or a single-variable regularly sampled timetable and must have at least
four samples. The CWT is obtained using the analytic Morse wavelet with the
symmetry parameter (gamma) equal to 3 and the time-bandwidth product equal to
cwt uses 10 voices per octave. The minimum and maximum
scales are determined automatically based on the energy spread of the wavelet in
frequency and time. If
x is real-valued,
wt is a 2-D matrix where each row corresponds to one
scale. The column size of
wt is equal to the length of
x is complex-valued,
wt is a 3-D matrix, where the first page is the CWT for
the positive scales (analytic part or counterclockwise component) and the second
page is the CWT for the negative scales (anti-analytic part or clockwise
cwt function uses 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 shows a more accurate representation of the signal. See L1 Norm for CWT and Sinusoid and Wavelet Coefficient Amplitudes.
cwt function supports GPU code generation. See GPU Code Generation for more
specifies the sampling frequency,
fs, in Hz as a positive
fs to determine the
scale-to-frequency conversions and returns the frequencies
f in Hz. If you do not specify a sampling frequency,
f in cycles per
sample. If the input
x is complex, the scale-to-frequency
conversions apply to both pages of
x is a timetable, you cannot specify
fs is determined from the
RowTimes of the timetable.
specifies the sampling period,
ts, as a positive
duration scalar. The
duration can be in years, days, hours, minutes, or
ts to compute the
scale-to-period conversion and returns the time periods in
period. The array of durations in
period has the same format property as
ts. If the input
x is complex, the
scale-to-period conversions apply to both pages of
x is a timetable, you cannot specify
ts is determined from the
RowTimes of the timetable when you set the
[___] = cwt(___,
returns the CWT with additional options specified by one or more
Name,Value pair arguments.
cwt(___) with no output arguments plots the
CWT scalogram. The scalogram is the absolute value of the CWT as a function of
time and frequency. Frequency is plotted on a logarithmic scale. The cone of
influence showing where edge effects become significant is also plotted. Gray
regions outside the dashed white line delineate regions where edge effects are
significant. If the input signal is complex-valued, the positive
(counterclockwise) and negative (clockwise) components are plotted in separate
If you do not specify a sampling frequency or sampling period, the frequencies are plotted in cycles per sample. If you specify a sampling frequency, the frequencies are in Hz. If you specify a sampling period, the scalogram is plotted as a function of time and periods. If the input signal is a timetable, the scalogram is plotted as a function of time and frequency in hertz and uses the RowTimes as the basis for the time axis.
To see the time, frequency, and magnitude of a scalogram point, enable data tips in the figure axes toolbar and click the desired point in the scalogram.
cwt clears (
clf) the current
figure. To plot the scalogram in a subplot, use a plotting
Obtain the continuous wavelet transform of a speech sample using default values.
load mtlb; w = cwt(mtlb);
Obtain the continuous wavelet transform of a speech sample using the bump wavelet instead of the default Morse wavelet.
load mtlb; cwt(mtlb,'bump',Fs);
Compare the result obtained from the CWT using the default Morse wavelet.
Create two complex exponentials, of different amplitudes, with frequencies of 32 and 64 Hz. The data is sampled at 1000 Hz. The two complex exponentials have disjoint support in time.
Fs = 1e3; t = 0:1/Fs:1; z = exp(1i*2*pi*32*t).*(t>=0.1 & t<0.3)+2*exp(-1i*2*pi*64*t).*(t>0.7);
Add complex white Gaussian noise with a standard deviation of 0.05.
wgnNoise = 0.05/sqrt(2)*randn(size(t))+1i*0.05/sqrt(2)*randn(size(t)); z = z+wgnNoise;
Obtain and plot the cwt using a Morse wavelet.
Note the magnitudes of the complex exponential components in the colorbar are essentially their amplitudes even though they are at different scales. This is a direct result of the L1 normalization. You can verify this by executing this script and exploring each subplot with a data cursor.
Obtain the CWT of Kobe earthquake data. The sampling frequency is 1 Hz.
Plot the earthquake data.
plot((1:numel(kobe))./60,kobe) xlabel('mins') ylabel('nm/s^2') grid on title('Kobe Earthquake Data')
Obtain the CWT, frequencies, and cone of influence.
[wt,f,coi] = cwt(kobe,1);
Plot the data, including the cone of influence.
Obtain the CWT, time periods, and cone of influence by specifying a sampling period instead of a sampling frequency.
[wt,periods,coi] = cwt(kobe,minutes(1/60));
View the same data by specifying a sampling period input instead of a frequency.
This example shows that the amplitudes of oscillatory components in a signal agree with the amplitudes of the corresponding wavelet coefficients.
Create a signal composed of two sinusoids with disjoint support in time. One sinusoid has a frequency of 32 Hz and amplitude equal to 1. The other sinusoid has a frequency of 64 Hz and amplitude equal to 2. The signal is sampled for one second at 1000 Hz. Plot the signal.
frq1 = 32; amp1 = 1; frq2 = 64; amp2 = 2; Fs = 1e3; t = 0:1/Fs:1; x = amp1*sin(2*pi*frq1*t).*(t>=0.1 & t<0.3)+amp2*sin(2*pi*frq2*t).*(t>0.6 & t<0.9); plot(t,x) grid on xlabel('Time (sec)') ylabel('Amplitude') title('Signal')
Create a CWT filter bank that can be applied to the signal. Since the signal component frequencies are known, set the frequency limits of the filter bank to a narrow range that includes the known frequencies. To confirm the range, plot the magnitude frequency responses for the filter bank.
fb = cwtfilterbank('SignalLength',numel(x),'SamplingFrequency',Fs,... 'FrequencyLimits',[20 100]); figure freqz(fb)
cwt and the filter bank to plot the scalogram of the signal.
Execute this script and use a data cursor to confirm that the amplitudes of the wavelet coefficients are essentially equal to the amplitudes of the sinusoidal components.
This example shows how using a CWT filter bank improves computational efficiency when taking the CWT of multiple time series.
Load the seismograph data recorded during the 1995 Kobe earthquake. The data are seismograph (vertical acceleration, nm/sq.sec) measurements recorded at Tasmania University, Hobart, Australia on 16 January 1995 beginning at 20:56:51 (GMT) and continuing for 51 minutes at 1 second intervals. Create a CWT filter bank that can be applied to the data.
load kobe fb = cwtfilterbank('SignalLength',numel(kobe),'SamplingFrequency',1);
cwt function and take the CWT of the data 250 times. Display the elapsed time used.
num = 250; tic; for k=1:num cfs = cwt(kobe); end toc
Elapsed time is 6.551628 seconds.
Now use the
wt object function of the filter bank to take the CWT of the data. Confirm using the filter bank is faster.
tic; for k=1:num cfs = wt(fb,kobe); end toc
Elapsed time is 3.782376 seconds.
This example shows how to change the default frequency axis labels for the CWT when you obtain a plot with no output arguments.
Create two sine waves with frequencies of 32 and 64 Hz. The data is sampled at 1000 Hz. The two sine waves have disjoint support in time. Add white Gaussian noise with a standard deviation of 0.05. Obtain and plot the CWT using the default Morse wavelet.
Fs = 1e3; t = 0:1/Fs:1; x = cos(2*pi*32*t).*(t>=0.1 & t<0.3)+sin(2*pi*64*t).*(t>0.7); wgnNoise = 0.05*randn(size(t)); x = x+wgnNoise; cwt(x,1000)
The plot uses a logarithmic frequency axis because frequencies in the CWT are logarithmic. In MATLAB, logarithmic axes are in powers of 10 (decades). You can use
cwtfreqbounds to determine what the minimum and maximum wavelet bandpass frequencies are for a given signal length, sampling frequency, and wavelet.
[minf,maxf] = cwtfreqbounds(numel(x),1000);
You see that by default MATLAB has placed frequency ticks at 10 and 100 because those are the powers of 10 between the minimum and maximum frequencies. If you wish to add more frequency axis ticks, you can obtain a logarithmically spaced set of frequencies between the minimum and maximum frequencies using the following.
numfreq = 10; freq = logspace(log10(minf),log10(maxf),numfreq);
Next, get the handle to the current axes and replace the frequency axis ticks and labels with the following.
AX = gca; AX.YTickLabelMode = 'auto'; AX.YTick = freq;
In the CWT, frequencies are computed in powers of two. To create the frequency ticks and tick labels in powers of two, you can do the following.
newplot; cwt(x,1000); AX = gca; freq = 2.^(round(log2(minf)):round(log2(maxf))); AX.YTickLabelMode = 'auto'; AX.YTick = freq;
This example shows how to scale scalogram values by maximum absolute value at each level for plotting.
Load in a signal and display the default scalogram. Change the colormap to
load noisdopp cwt(noisdopp) colormap(pink(240))
Take the CWT of the signal and obtain the wavelet coefficients and frequencies.
[cfs,frq] = cwt(noisdopp);
To efficiently find the maximum value of the coefficients at each frequency (level), first transpose the absolute value of the coefficients. Find the minimum value at every level. At each level, subtract the level's minimum value.
tmp1 = abs(cfs); t1 = size(tmp1,2); tmp1 = tmp1'; minv = min(tmp1); tmp1 = (tmp1-minv(ones(1,t1),:));
Find the maximum value at every level of
tmp1. For each level, divide every value by the maximum value at that level. Multiply the result by the number of colors in the colormap. Set equal to 1 all zero entries. Transpose the result.
maxv = max(tmp1); maxvArray = maxv(ones(1,t1),:); indx = maxvArray<eps; tmp1 = 240*(tmp1./maxvArray); tmp2 = 1+fix(tmp1); tmp2(indx) = 1; tmp2 = tmp2';
Display the result. The scalogram values are now scaled by the maximum absolute value at each level. Frequencies are displayed on a linear scale.
t = 0:length(noisdopp)-1; pcolor(t,frq,tmp2); shading interp ylabel('Frequency') title('Scalogram Scaled By Level') colormap(pink(240))
This example shows that increasing the time-bandwidth product of the Morse wavelet creates a wavelet with more oscillations under its envelope. All other things being equal, increasing makes the wavelet more frequency-localized as a result.
Create two filter banks. One filter bank has the default
TimeBandwidth value of 60. The second filter bank has a
TimeBandwidth value of 10. The
SignalLength for both filter banks is 4096 samples.
sigLen = 4096; fb60 = cwtfilterbank('SignalLength',sigLen); fb10 = cwtfilterbank('SignalLength',sigLen,'TimeBandwidth',10);
Obtain the time-domain wavelets for the filter banks.
[psi60,t] = wavelets(fb60); [psi10,~] = wavelets(fb10);
scales function to find the mother wavelet for each filter bank.
sca60 = scales(fb60); sca10 = scales(fb10); [~,idx60] = min(abs(sca60-1)); [~,idx10] = min(abs(sca10-1)); m60 = psi60(idx60,:); m10 = psi10(idx10,:);
Since the time-bandwidth product is larger for the
fb60 filter bank, verify the
m60 wavelet has more oscillations under its envelope than the
subplot(2,1,1) plot(t,abs(m60)) grid on hold on plot(t,real(m60)) plot(t,imag(m60)) xlim([-30 30]) legend('abs(m60)','real(m60)','imag(m60)') title('TimeBandwidth = 60') subplot(2,1,2) plot(t,abs(m10)) grid on hold on plot(t,real(m10)) plot(t,imag(m10)) xlim([-30 30]) legend('abs(m10)','real(m10)','imag(m10)') title('TimeBandwidth = 10')
Align the peaks of the
m10 magnitude frequency responses. Verify the frequency response of the
m60 wavelet is more narrow than the frequency response for the
cf60 = centerFrequencies(fb60); cf10 = centerFrequencies(fb10); m60cFreq = cf60(idx60); m10cFreq = cf10(idx10); freqShift = 2*pi*(m60cFreq-m10cFreq); x10 = m10.*exp(1j*freqShift*(-sigLen/2:sigLen/2-1)); figure plot([abs(fft(m60)).' abs(fft(x10)).']) grid on legend('Time-bandwidth = 60','Time-bandwidth = 10') title('Magnitude Frequency Responses')
This example shows how to plot the CWT scalogram in a figure subplot.
Load the speech sample. The data is sampled at 7418 Hz. Plot the default CWT scalogram.
load mtlb cwt(mtlb,Fs)
Obtain the continuous wavelet transform of the signal, and the frequencies of the CWT.
[cfs,frq] = cwt(mtlb,Fs);
cwt function sets the time and frequency axes in the scalogram. Create a vector representing the sample times.
tms = (0:numel(mtlb)-1)/Fs;
In a new figure, plot the original signal in the upper subplot and the scalogram in the lower subplot. Plot the frequencies on a logarithmic scale.
figure subplot(2,1,1) plot(tms,mtlb) axis tight title('Signal and Scalogram') xlabel('Time (s)') ylabel('Amplitude') subplot(2,1,2) surface(tms,frq,abs(cfs)) axis tight shading flat xlabel('Time (s)') ylabel('Frequency (Hz)') set(gca,'yscale','log')
x— Input signal
Input signal, specified as a double-precision real- or complex-valued
vector or a single-variable regularly sampled timetable.
x must have at least four samples.
wname— Analytic wavelet
Analytic wavelet used to compute the CWT, specified as
'bump'. These character
vectors specify the analytic Morse, Morlet (Gabor), and bump wavelet,
The default Morse wavelet has symmetry parameter () equal to 3 and time-bandwidth product equal to 60.
fs— Sampling frequency
Sampling frequency, in Hz, specified as a positive scalar. If
fs, then you cannot specify
ts— Sampling period
Sampling period, also known as the time duration, specified as a scalar
duration. Valid durations are
seconds. You cannot
use calendar durations. If you specify
ts, then you
wt = cwt(x,hours(12))
comma-separated pairs of
the argument name and
Value is the corresponding value.
Name must appear inside quotes. You can specify several name and value
pair arguments in any order as
'ExtendSignal',falseindicates that the signal is not extended.
'ExtendSignal'— Extend input signal symmetrically
Option to extend the input signal symmetrically by reflection,
specified as the comma-separated pair consisting of
the signal symmetrically can mitigate boundary effects. If you specify
then the signal is extended. If you specify
then the signal is not extended.
'FrequencyLimits'— Frequency limits
Frequency limits to use in the CWT, specified as a two-element vector with positive strictly increasing entries. The first element specifies the lowest peak passband frequency and must be greater than or equal to the product of the wavelet peak frequency in hertz and two time standard deviations divided by the signal length. The second element specifies the highest peak passband frequency and must be less than or equal to the Nyquist frequency. The base 2 logarithm of the ratio of the maximum frequency to the minimum frequency must be greater than or equal to 1/NV where NV is the number of voices per octave.
If you specify frequency limits outside the permissible range,
cwt truncates the limits to the minimum and
maximum valid values. Use
cwtfreqbounds to determine frequency limits for
different parameterizations of the CWT. For complex-valued signals,
flimits is used for the anti-analytic part,
flimits is the vector specified by
'PeriodLimits'— Period limits
Period limits to use in the CWT, specified as a two-element duration
array with strictly increasing positive entries.
The first element must be greater than or equal to
ts is the
sampling period. The base 2 logarithm of the ratio of the minimum period
to the maximum period must be less than or equal to -1/NV where NV is
the number of voices per octave. The maximum period cannot exceed the
signal length divided by the product of two time standard deviations of
the wavelet and the wavelet peak frequency.
If you specify period limits outside the permissible range,
cwt truncates the limits to the minimum and
maximum valid values. Use
cwtfreqbounds to determine period limits for different
parameterizations of the wavelet transform. For complex-valued signals,
plimits is used for the anti-analytic part,
plimits is the vector specified by
'VoicesPerOctave'— Number of voices per octave
10(default) | even integer from 4 to 48
Number of voices per octave to use for the CWT, specified as
the comma-separated pair consisting of
an even integer from 4 to 48. The CWT scales are discretized using
the specified number of voices per octave. The energy spread of the
wavelet in frequency and time automatically determines the minimum
and maximum scales.
'TimeBandwidth'— Time-bandwidth product of the Morse wavelet
Time-bandwidth product of the Morse wavelet, specified as the
comma-separated pair consisting of
and a scalar greater than 3 and less than or equal to 120. The symmetry
parameter, gamma (), is fixed at 3. Wavelets with larger time-bandwidth
products have larger spreads in time and narrower spreads in frequency.
The standard deviation of the Morse wavelet in time is approximately
sqrt(TimeBandwidth/2). The standard deviation of
the Morse wavelet in frequency is approximately
If you specify
'TimeBandwidth', you cannot specify
'WaveletParameters'. To specify both the symmetry
and time-bandwidth product, use
In the notation of Morse Wavelets,
'WaveletParameters'— Symmetry and time-bandwidth product of the Morse wavelet
(3,60)(default) | two-element vector of scalars
Symmetry and time-bandwidth product of the Morse wavelet, specified as
the comma-separated pair consisting of
'WaveletParameters' and a two-element vector of
The first element is the symmetry, , which must be greater than or equal to 1. The second element is the time-bandwidth product, which must be strictly greater than . The ratio of the time-bandwidth product to cannot exceed 40.
When is equal to 3, the Morse wavelet is perfectly symmetric in the frequency domain and the skewness is 0. When is greater than 3, the skewness is positive. When is less than 3, the skewness is negative.
If you specify
'WaveletParameters', you cannot
'NumOctaves'— Number of octaves
Number of octaves, specified as the comma-separated pair consisting of
'NumOctaves' and a positive integer. The number
of octaves cannot exceed
fmin are the maximum
and minimum CWT frequencies (or periods) as determined by the signal
length, sampling frequency, and wavelet. See
cwtfreqbounds for details.
'NumOctaves' name-value pair is not recommended
and will be removed in a future release. The recommended way to modify
the frequency or period range of the CWT is with the
'PeriodLimits' name-value pairs. You cannot
specify both the
'PeriodLimits' name-value pairs.
'FilterBank'— CWT filter bank
CWT filter bank to use to compute the CWT, specified as a CWT filter
bank object. If you use the
pair, you cannot specify any other options. All options for the
computation of the CWT are defined as properties of
x is a timetable, the sampling frequency or
sampling period in
fb must agree with the sampling
frequency or sampling period determined by the
RowTimes of the timetable.
wt = cwt(x,'FilterBank',fb)
wt— Continuous wavelet transform
Continuous wavelet transform, returned as a matrix of complex values. By
cwt uses the analytic Morse (3,60) wavelet,
where 3 is the symmetry and 60 is the time-bandwidth product.
cwt uses 10 voices per octave. If
x is real-valued,
wt is an
Na-by-N matrix, where
Na is the number of scales, and N
is the number of samples in
x is complex-valued,
wt is a
3-D matrix, where the first page is the CWT for the positive scales
(analytic part or counterclockwise component) and the second page is the CWT
for the negative scales (anti-analytic part or clockwise component). The
minimum and maximum scales are determined automatically based on the energy
spread of the wavelet in frequency and time. See Algorithms for information on how the scales are
Frequencies of the CWT, returned as a vector. If you specify
a sampling frequency,
in Hz. If you do not specify
cycles per sample.
period— Time periods
Time periods, returned as an array of durations. The durations
are in the same format as
ts. Each row corresponds
to a period.
coi— Cone of influence
Cone of influence for the CWT, returned as either an array of doubles or array of durations.
The cone of influence indicates where edge effects occur in the CWT. If you
specify a sampling frequency,
fs, the cone of influence
is in Hz. If you specify a scalar duration,
cone of influence is in periods. Due to the edge effects, give less credence
to areas that are outside or overlap the cone of influence.
For additional information, see Boundary Effects and the Cone of Influence.
fb— CWT filter bank
CWT filter bank used in the CWT, returned as a CWT filter bank object. See
scalingcfs— Scaling coefficients
Scaling coefficients for the CWT if the analyzing wavelet is
'amor', returned as a
real- or complex-valued vector. The length of
is equal to the length of the input
Scaling coefficients are not supported for the bump wavelet.
Analytic wavelets are complex-valued wavelets whose Fourier transform vanish for negative frequencies. Analytic wavelets are a good choice when doing time-frequency analysis with the CWT. Because the wavelet coefficients are complex-valued, the coefficients provide phase and amplitude information of the signal being analyzed. Analytic wavelets are well suited for studying how the frequency content in real world nonstationary signals evolves as a function of time.
Analytic wavelets are almost exclusively based on rapidly decreasing functions. If is an analytic rapidly decreasing function in time, then its Fourier transform is a rapidly decreasing function in frequency and is small outside of some interval where . Orthogonal and biorthogonal wavelets are typically designed to have compact support in time. Wavelets with compact support in time have relatively poorer energy concentration in frequency than wavelets which rapidly decrease in time. Most orthogonal and biorthogonal wavelets are not symmetric in the Fourier domain.
Morse Wavelet Family (default)
Analytic Morlet (Gabor) Wavelet
If you want to do time-frequency analysis using orthogonal or
biorthogonal wavelets, we recommend
When using wavelets for time-frequency analysis, you usually convert scales to
frequencies or periods to interpret results.
cwtfilterbank do the conversion. You can obtain the corresponding
scales associated by using
on the optional
cwt output argument
cwt supports CUDA®code generation. You must have MATLAB®
Coder™ and GPU Coder™ to generate GPU code.
Usage notes and limitations:
Timetable input signal is not supported.
All inputs must be constant and specified at compilation time.
Only analytic Morse (
'morse') and Morlet
'amor') wavelets are supported.
The following input arguments are not supported:
PeriodLimits name-value pair,
NumOctave name-value pair, and
FilterBank name-value pair.
Scaling coefficient output and filter bank output are not supported.
Plotting is not supported.
If you are taking the CWT of multiple time series, it is more efficient to
precompute and use a CWT filter bank than repeatedly apply the
cwt function. See Using CWT Filter Bank on Multiple Time Series.
To determine the minimum scale, find the peak frequency of the base wavelet. For Morse wavelets, dilate the wavelet so that the Fourier transform of the wavelet at radians is equal to 10% of the peak frequency. The smallest scale occurs at the largest frequency:
As a result, the smallest scale is the minimum of (2, ). For Morse wavelets, the smallest scale is usually . For the Morlet wavelet, the smallest scale is usually 2.
Both the minimum and maximum scales of the CWT are determined automatically based on the energy spread of the wavelet in frequency and time. To determine the maximum scale, CWT uses the following algorithm.
The standard deviation of the Morse wavelet in time, , is approximately , where is the time-bandwidth product. The standard deviation in frequency, , is approximately . If you scale the wavelet by some , the time duration changes to , which is the wavelet stretched to equal the full length (N samples) of the input. You cannot translate this wavelet or stretch it further without causing it to wrap, so the largest scale is .
Wavelet transform scales are powers of 2 and are denoted by . NV is the number of voices per octave, and j ranges from 0 to the largest scale. For a specific small scale, :
Converting to log2:
Therefore, the maximum scale is
In integral form, the CWT preserves energy. However, when you implement the CWT
numerically, energy is not preserved. In this case, regardless of the normalization
you use, the CWT is not an orthonormal transform. The
function uses L1 normalization.
Wavelet transforms commonly use L2 normalization of the wavelet. For the L2 norm, dilating a signal by 1/s, where s is greater than 0, is defined as follows:
The energy is now s times the original energy. When included in the Fourier transform, multiplying by produces different weights being applied to different scales, so that the peaks at higher frequencies are reduced more than the peaks at lower frequencies.
In many applications, L1 normalization is better. The L1 norm definition does not include squaring the value, so the preserving factor is 1/s instead of . Instead of high-frequency amplitudes being reduced as in the L2 norm, for L1 normalization, all frequency amplitudes are normalized to the same value. Therefore, using the L1 norm shows a more accurate representation of the signal. See example Continuous Wavelet Transform of Two Complex Exponentials.
Not recommended starting in R2018a
'NumOctaves' name-value pair argument will be removed in a
future release. Use either:
Name-value pair argument
modify the frequency range of the CWT.
Name-value pair argument
'PeriodLimits' to modify
the period range of the CWT.
cwtfreqbounds for additional information.
 Lilly, J. M., and S. C. Olhede. “Generalized Morse Wavelets as a Superfamily of Analytic Wavelets.” IEEE Transactions on Signal Processing. Vol. 60, No. 11, 2012, pp. 6036–6041.
 Lilly, J. M., and S. C. Olhede. “Higher-Order Properties of Analytic Wavelets.” IEEE Transactions on Signal Processing. Vol. 57, No. 1, 2009, pp. 146–160.
 Lilly, J. M. jLab: A data analysis package for Matlab, version 1.6.2. 2016. http://www.jmlilly.net/jmlsoft.html.
 Lilly, J. M. “Element analysis: a wavelet-based method for analysing time-localized events in noisy time series.” Proceedings of the Royal Society A. Volume 473: 20160776, 2017, pp. 1–28. dx.doi.org/10.1098/rspa.2016.0776.