Accelerating the pace of engineering and science

# 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 Wavelet Definitions for a description 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. s0 — The smallest scale. The default s0 depends on the wavelet. See Wavelet Definitions for the wavelet-dependent default.ds — Spacing between scales. The default ds depends on the wavelet. See Wavelet Definitions for the wavelet-dependent default. You can construct a linear or logarithmic scale vector using ds. See type for a description of the type of spacing.nb — Number of scales. The default nb depends on the wavelet. See Wavelet Definitions for the wavelet-dependent default.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.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 See Wavelet Definitions for formal definitions of the supported analyzing wavelets and associated defaults. Default: 'morl'

## Output Arguments

 cwtstruct A structure array with six fields. The fields of the structure array are:dt — The sampling interval of the 1–D input signalcfs — 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.meanSIG — Mean of the analyzed signalomega — Vector of angular frequenciesscales — Vector of scales at which the CWT is computed. The length of cwtstruct.scales is equal to the row dimension of cwtstruct.cfs.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.

expand all

### Wavelet Definitions

#### Morlet Wavelet

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

$\stackrel{^}{\Psi }\left(s\omega \right)={\pi }^{-1/4}{e}^{{\left(s\omega -{\omega }_{0}\right)}^{2}/2}U\left(s\omega \right)$

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

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

$\stackrel{^}{\Psi }\left(s\omega \right)={\pi }^{-1/4}{e}^{{\left(s\omega -{\omega }_{0}\right)}^{2}/2}$

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

$\stackrel{^}{\Psi }\left(s\omega \right)={\pi }^{-1/4}\left\{{e}^{{\left(s\omega -{\omega }_{0}\right)}^{2}/2}-{e}^{{\omega }_{0}^{2}/2}\right\}$

The default value of ω0 is 6.

The scale-to-frequency Fourier factor for the Morlet wavelet is:

$\frac{4\pi s}{{\omega }_{0}+\sqrt{2+{\omega }_{0}^{2}}}$

The default smallest scale for the Morlet wavelets is 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 fix(log2(length(sig))/ds)+1.

#### m-th Order Derivative of Gaussian Wavelets

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

$\stackrel{^}{\Psi }\left(s\omega \right)=-\frac{1}{\sqrt{\Gamma \left(m+1/2\right)}}{\left(js\omega \right)}^{m}{e}^{-{\left(s\omega \right)}^{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 scale-to-frequency Fourier factor for the DOG wavelet is:

$\frac{2\pi s}{\sqrt{m+\frac{1}{2}}}$

The default smallest scale for the DOG wavelet is 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 max([fix(log2(length(sig))/ds),1]).

#### Paul Wavelet

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

$\stackrel{^}{\Psi }\left(s\omega \right)={2}^{m}\sqrt{m\left(2m-1\right)!}\text{ }\text{ }{\left(s\omega \right)}^{m}{e}^{-s\omega }U\left(sw\right)$

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

The default order of the Paul wavelet is 4.

The scale-to-frequency Fourier factor for the Paul wavelet is:

$\frac{4\pi s}{2m+1}$

The default smallest scale for the Paul wavelet is 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 fix(log2(length(sig))/ds)+1.

### 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 DFT-Based Continuous Wavelet Transform.

## 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.