Documentation

cwtft

Continuous wavelet transform using FFT algorithm

Syntax

cwtstruct = cwtft(sig)
cwtstruct = cwtft(sig,Name,Value)
cwtstruct = cwtft(...,'plot')

Description

cwtstruct = cwtft(sig) returns the continuous wavelet transform (CWT) of the 1–D input signal sig. cwtft uses an FFT algorithm to compute the CWT. sig can be a vector, a structure array, or a cell array. If the sampling interval of your signal is not equal to 1, you must input the sampling period with sig in a cell array or a structure array to obtain correct results. If sig is a cell array, sig{1} is equal to your signal and sig{2} is equal to the sampling interval. If sig is a structure array, the field sig.val contains your signal and sig.period contains the sampling interval.

By default, cwtft uses the analytic Morlet wavelet. See More About for descriptions of valid analyzing wavelets.

For additional default values, see scales in Name-Value Pair Arguments.

cwtstruct = cwtft(sig,Name,Value) returns the continuous wavelet transform (CWT) of the 1–D input signal sig with additional options specified by one or more Name,Value pair arguments. See Name-Value Pair Arguments for a comprehensive list.

cwtstruct = cwtft(...,'plot') plots the continuous wavelet transform. If the analyzing wavelet is real-valued, the original signal along with the CWT coefficient magnitudes and signed CWT coefficients are plotted. If the analyzing wavelet is complex-valued, the original signal is plotted along with the moduli, real parts, imaginary parts, and angles of the CWT coefficients. You can select the radio button in the bottom left of the plot to superimpose the signal's reconstruction using icwtft.

Input Arguments

sig

The 1–D input signal. sig can be a vector, a structure array, or a cell array. If sig is a structure array, sig contains two fields: val and period. sig.val is the signal vector and sig.period is the sampling period. If sig is a cell array, sig{1} is the signal vector and sig{2} is the sampling period.

If sig is a vector, the sampling period defaults to 1.

    Note:   If the sampling interval of your input signal is not 1, you must input the sampling interval with sig in a cell array or structure array to obtain correct results. If sig is a cell array, sig{1} is the 1–D input signal and sig{2} is the sampling period. If sig is a structure array, the field sig.val is the 1–D input signal and sig.period is the sampling interval.

Name-Value Pair Arguments

'scales'

Scales over which to compute the CWT. The value of scales can be a vector, a structure array, or a cell array. If scales is a structure array, it contains at most five fields. The first three fields are mandatory. The last two fields are optional.

  1. s0 — The smallest scale. The default s0 depends on the wavelet. See the More About for descriptions of the default for each wavelet.

  2. ds — Spacing between scales. The default ds depends on the wavelet. See the More About for descriptions of the default for eacj wavelet. You can construct a linear or logarithmic scale vector using ds. See type for a description of the type of spacing.

  3. nb — Number of scales. The default nb depends on the wavelet. See the More About for descriptions of the default for each wavelet.

  4. type — Type of spacing between scales. type can be one of 'pow' or 'lin'. The default is 'pow'. If type is equal to 'pow', the CWT scales are s0*pow.^((0:nb-1)*ds). This results in a constant spacing of ds if you take the logarithm to the base power of the scales vector. If type is equal to 'lin', the CWT scales are linearly spaced by s0 + (0:nb-1)*ds.

    Use the default power of two spacing to ensure an accurate approximation to the original signal based only on select scales. See the second example in Examples for a signal approximation based on select scales.

  5. pow — The base for 'pow' spacing. The default is 2. This input is valid only if the type argument is 'pow'.

If scales is a cell array, the first three elements of the cell array are identical to the first three elements of the structure array described in the preceding list. The last two elements of the cell array are optional and match the two optional inputs in the structure array described in the preceding list.

'wavelet'

Analyzing wavelet. The supported analyzing wavelets are:

  • 'dog'm-th order derivative of a Gaussian wavelet where m is a positive even integer. The default value of m is 2.

  • 'morl' — Morlet wavelet. Results in an analytic Morlet wavelet. The Fourier transform of an analytic wavelet is zero for negative frequencies.

  • 'morlex' — non-analytic Morlet wavelet

  • 'morl0' — non-analytic Morlet wavelet with zero mean

  • 'mexh' — Mexican hat wavelet. The Mexican hat wavelet is a special case of the m-th order derivative of a Gaussian wavelet with m=2.

  • 'paul' — Paul wavelet

  • 'bump' — Bump wavelet

See the More About for formal definitions of the supported analyzing wavelets and associated defaults.

Default: 'morl'

'padmode'

Signal extension mode. See dwtmode for supported extension modes. By default, cwtft does not extend the signal prior to computing the CWT. In a Fourier-transform-based CWT algorithm, extending a signal can mitigate wrap-around effects. The number of CWT coefficients in each row of the output matrix cwtstruct.cfs is truncated to match the length of the input signal.

Output Arguments

cwtstruct

A structure array with six fields. The fields of the structure array are:

  • cfs — The CWT coefficient matrix. cwtstruct.cfs is an nb-by-N matrix where nb is the number of scales and N is the length of the input signal.

  • scales — Vector of scales at which the CWT is computed. The length of cwtstruct.scales is equal to the row dimension of cwtstruct.cfs.

  • frequencies — Frequencies in cycles per unit time (or space) corresponding to the scales. If the sampling period units are seconds, the frequencies are in hertz. The elements of frequencies are in decreasing order to correspond to the elements in the scales vector. Use this field to examine the CWT in the time-frequency plane.

  • omega — Vector of angular frequencies used in the Fourier transform of the wavelet. This field is used in icwtft and icwtlin for the inversion of the CWT for all wavelets except the bump wavelet.

  • meanSIG — Mean of the analyzed signal

  • dt — The sampling interval of the 1–D input signal

  • wav — Analyzing wavelet

Examples

Compute and display the CWT of sine waves with disjoint support. The sampling interval is 1/1023.

N = 1024;
% Sampling interval is 1/1023
t = linspace(0,1,N);
y = sin(2*pi*4*t).*(t<=0.5)+sin(2*pi*8*t).*(t>0.5);
% Because the sampling interval differs from the default
% you must input it along with the signal
% Using cell array input
sig = {y,1/1023};
cwtS1 = cwtft(sig,'plot');

You can display or hide the reconstructed signal using the radio button at the bottom left of the figure. When you select the radio button, the maximum and quadratic relative errors are computed and displayed along with the reconstructed signal.

Reconstruct an approximation to a sum of disjoint sine waves in noise using cwtft to decompose the signal and icwtft to reconstruct the approximation. Use the CWT coefficients to identify the scales isolating the sinusoidal components. Reconstruct an approximation to the signal based on those scales using the inverse CWT. To ensure an accurate approximation to the based on select scales, use the default power of two spacing in the CWT.

rng default % Reset random number generator for reproducible results
N = 1024;
% Sampling interval is 1/1023
t = linspace(0,1,N);
y = sin(2*pi*4*t).*(t<=0.5)+sin(2*pi*8*t).*(t>0.5);
ynoise = y+randn(size(t));
% Because the sampling interval differs from the default
% you must input it along with the signal
% Using structure array input
sig = struct('val',ynoise,'period',1/1023);
cwtS1 = cwtft(sig);
scales = cwtS1.scales;
MorletFourierFactor = 4*pi/(6+sqrt(2+6^2));
freq = 1./(scales.*MorletFourierFactor);
contour(t,freq,real(cwtS1.cfs));
xlabel('Seconds'); ylabel('Pseudo-frequency');
axis([0 t(end) 0 15]);

Extract the scales dominated by energy from the two sine waves and reconstruct a signal approximation using the inverse CWT.

cwtS2 = cwtS1;
cwtS2.cfs = zeros(size(cwtS1.cfs));
cwtS2.cfs(13:15,:) = cwtS1.cfs(13:15,:);
xrec = icwtft(cwtS2);
subplot(2,1,1);
plot(t,ynoise);
title('Sum of Disjoint Sinusoids in Noise');
subplot(2,1,2);
plot(t,xrec,'b'); hold on; axis([0 1 -4 4]);
plot(t,y,'r');
legend('Reconstructed Signal','Original Signal',...
    'Location','NorthWest');
xlabel('Seconds'); ylabel('Amplitude');

Alternatives

  • cwt — Computes the CWT using convolutions. cwt supports a wider choice of analyzing wavelets than cwtft, but may be more computationally expensive. The output of cwt is not compatible with the inverse CWT implemented with icwtft. To use icwtft, obtain the CWT with cwtft.

More About

expand all

Morlet Wavelet

Both non-analytic and analytic Morlet wavelets are supported. The analytic Morlet wavelet, 'morl', is defined in the Fourier domain by:

Ψ^(sω)=π1/4e(sωω0)2/2U(sω)

where U(ω) is the Heaviside step function [5].

The non-analytic Morlet wavelet, 'morlex', is defined in the Fourier domain by:

Ψ^(sω)=π1/4e(sωω0)2/2

'morl0' defines a non-analytic Morlet wavelet in the Fourier domain with exact zero mean:

Ψ^(sω)=π1/4{e(sωω0)2/2eω02/2}

The default value of ω0 is 6.

The default smallest scale for the Morlet wavelets is s0 = 2*dt where dt is the sampling period.

The default spacing between scales for the Morlet wavelets is ds=0.4875.

The default number of scales for the Morlet wavelets is NbSc = fix(log2(length(sig)*dt/s0)/ds).

The default scales for the Morlet wavelet are s0*2.^((0:NbSc-1)*ds).

m-th Order Derivative of Gaussian Wavelets

In the Fourier domain, the m-th order derivative of Gaussian wavelets, 'dog', are defined by:

Ψ^(sω)=1Γ(m+1/2)(jsω)me(sω)2/2

where Г( ) denotes the gamma function [5].

The derivative must be an even order. The default order of the derivative is 2, which is also known as the Mexican hat wavelet .

The default smallest scale for the DOG wavelet is s0 = 2*dt where dt is the sampling period.

The default spacing between scales for the DOG wavelet is ds=0.4875.

The default number of scales for the DOG wavelet is NbSc = fix(log2(length(sig)*dt/s0)/ds).

The default scales for the DOG wavelet are s0*2.^((0:NbSc-1)*ds).

Paul Wavelet

The Fourier transform of the analytic Paul wavelet, 'paul', of order m is:

Ψ^(sω)=2mm(2m1)!(sω)mesωU(sw)

where U(ω) is the Heaviside step function [5].

The default order of the Paul wavelet is 4.

The default smallest scale for the Paul wavelet is s0 = 2*dt where dt is the sampling period.

The default spacing between scales for the Paul wavelet is ds=0.4875.

The default number of scales for the Paul wavelet is NbSc = fix(log2(length(sig)*dt/s0)/ds).

The default scales for the Paul wavelet are s0*2.^((0:NbSc-1)*ds).

Bump wavelet

The Fourier transform of the analytic bump wavelet, 'bump', with parameters μ and σ is

ψ^(sω)=e(111(sωμ)2/σ2)1[(μσ)/s,(μ+σ)/s]

where 1[(μσ)/s,(μ+σ)/s] is the indicator function for the interval (μσ)/sω(μ+σ)/s.

Valid values for μ are [3,6]. Valid values for σ are [0.1, 1.2]. Smaller values of σ result in a wavelet with superior frequency localization but poorer time localization. Larger values of σ produce a wavelet with better time localization and poorer frequency localization.

The default values for μ and σ are 5 and 0.6 respectively.

The default smallest scale for the bump wavelet is s0 = 2*dt where dt is the sampling period.

The default spacing between scales for the bump wavelet is ds=1/10.

The default number of scales for the bump wavelet is NbSc = fix(log2(length(sig)*dt/s0)/ds).

The default scales for the bump wavelet are s0*2.^((0:NbSc-1)*ds).

Algorithms

cwtft implements the following algorithm:

  • Obtain the discrete Fourier transform (DFT) of the signal using fft.

  • Obtain the DFT of the analyzing wavelet at the appropriate angular frequencies. Scale the DFT of the analyzing wavelet at different scales to ensure different scales are directly comparable.

  • Take the product of the signal DFT and the wavelet DFT over all the scales. Invert the DFT to obtain the CWT coefficients.

For a mathematical motivation for the FFT-based algorithm see Continuous Wavelet Transform and Scale-Based Analysis.

References

[1] Daubechies, I. Ten Lectures on Wavelets, Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM), 1992.

[2] Farge, M. "Wavelet Transforms and Their Application to Turbulence", Ann. Rev. Fluid. Mech., 1992, 24, 395–457.

[3] Mallat, S. A Wavelet Tour of Signal Processing, San Diego, CA: Academic Press, 1998.

[4] Sun,W. "Convergence of Morlet's Reconstruction Formula", preprint, 2010.

[5] Torrence, C. and G.P. Compo. "A Practical Guide to Wavelet Analysis", Bull. Am. Meteorol. Soc., 79, 61–78, 1998.

See Also

| |

Introduced in R2011a

Was this topic helpful?