Time-frequency Analysis Tutorial

This page demonstrates the use of our time-frequency analysis add-on package for Mathematica.

Tutorial

This tutorial notebook demonstrates the capabilities of the Continuous Wavelet Transform (CWT) for the time-frequency analysis of signals.

Loading

Once purchased and installed, our CWT Mathematica package can be loaded with a single line of code:

Needs["CWT`"] ;

SetOptions[ContourPlot, RotateLabelFalse] ;

Executing this notebook will then reproduce the results shown below.

Introduction

Techniques which decompose a function into a concentration (such as instantaneous amplitude or energy density) over the time-frequency plane are collectively known as time-frequency methods. We refer to such a concentration as a time-frequency representation (TFR) of a function.

For example, given a signal f(t)=sin(2π t^2) with linearly increasing frequency:

ϕ[t_] := π t^2

f[t_] := Sin[ϕ[t]]

Plot[f[t], {t, 0, 4}, AxesLabel {"t", "f(t)"}] ;

[Graphics:HTMLFiles/tutorial_7.gif]

The instantaneous frequency ν(t) of f(t) is the derivative of the phase: ν(t)=1/(2π)|(d ϕ)/(d t)|=|t|.

ν[t_] = 1/(2π) Abs[ϕ '[t]]

Abs[t]

This is most easily illustrated by expanding the signal about a point as a sine wave, say about t=2:

Plot[{f[t], Sin[ϕ '[2] t]}, {t, 0, 4}, AxesLabel {"t", "f(t)"}] ;

[Graphics:HTMLFiles/tutorial_13.gif]

On the basis of the instantaneous frequency, the time-frequency representation of this signal is expected to be a single line ν(t)=|t| in the t-ν plane:

Plot[ν[t], {t, 0, 4}, AxesLabel {"t", ν  (t)}, AspectRatio1] ;                                                                    f

[Graphics:HTMLFiles/tutorial_15.gif]

However, obtaining the time-frequency representation of an arbitrary signal is decidedly tricky. There are many approaches to solving this problem: the Hilbert transform, the Wigner-Ville distribution, the Choi-Williams distribution, the short-time Fourier transform, the windowed Fourier transform and the wavelet transform. Due to its properties, the wavelet transform is the most generally applicable of these time-frequency methods. In particular, the continuous wavelet transform with a suitable wavelet is a very powerful tool for analysing the time-frequency content of arbitrary signals.

Artificial Example Analyses

The following examples illustrate the use of the CWT package for time-frequency analysis.

Sine wave

The simplest function to analyse is a sine wave. The time-frequency representation of a sine wave should simply be a line of constant frequency. Consider the function:

f[t_] := Cos[t]

Plot[f[t], {t, 0, 6π}, AxesLabel {"t", "f(t)"}] ;

[Graphics:HTMLFiles/tutorial_18.gif]

The CWT of f(t) is easily computed using the FunctionCWT function provided by the CWT package:

? FunctionCWT

FunctionCWT[f, {t, bmin, bmax, n}, options...] computes an approximation to the continuous w ... Parameter option can be used to control the trade-off between temporal and spectral resolutions.

sinecwt = FunctionCWT[f[t], {t, 0, 2π, 128}] ;//Timing

{0.897863 Second, Null}

Note that this computation was performed in around 1 second.

The resulting expression encapsulates the information used to compute the CWT, the extent of the validity of the result and the result itself:

sinecwt

CWT[{{-4.15888, 0.693147}, {0., 6.28319}}, 0.853553, 128, InterpolatingFunction[{{-4.15888, 0.693147}, {0., 6.28319}}, <>]]

This result contains thousands of carefully computed accurate numerical approximations to the true values of the CWT. These numerical approximations are interpolated to provide the illusion of a truly continuous transform. For example, the complex value of the CWT at a=1, b=1 may be obtained by suppling ln(a) and b as arguments to the CWT:

sinecwt[Log[a], b]/.{a1, b1}

0.583015 + 0.937239 

The range over ln(a) and b for which the CWT is valid can be extracted using GetCWTlogaRange and GetCWTbRange. In this case -4.16≤ln(a)≤0.69 and 0≤b≤6.28:

{logamin, logamax} = GetCWTlogaRange[sinecwt]

{-4.15888, 0.693147}

{bmin, bmax} = GetCWTbRange[sinecwt]

{0., 6.28319}

The CWT is most easily visualised as a 3D surface plot of the modulus |Overscript[f, ~] _Ψ(a,b)|:

Plot3D[Abs[sinecwt[Log[a], b]], {b, bmin, bmax}, {a, ^logamin, 2}, AxesLabel ...                                                                                           Ψ

[Graphics:HTMLFiles/tutorial_33.gif]

Some forms of analysis act directly upon the CWT. However, the TFR of the function is typically much more useful and can either be obtained directly (computed using FunctionTFR) or by wrapping the CWT in TFR:

sinetfr = TFR[sinecwt] ;

As for the CWT, the complex value of the TFR may be obtained by suppling a time t and angular frequency ω as arguments to the TFR:

sinetfr[t, ω]/.{t1, ω1}

0.527199 + 0.84751 

The range over t and ω for which the TFR is valid can be extracted using GetTFRωRange and GetTFRtRange. In this case 0.5≤ω≤64 and 0≤t≤6.28:

{ωmin, ωmax} = GetTFRωRange[sinetfr]

{0.5, 64.}

{tmin, tmax} = GetTFRtRange[sinetfr]

{0., 6.28319}

The TFR is most easily visualised as a 3D surface or contour plot of the modulus of the TFR, which is the instantaneous amplitude A_f of the signal f as a function of time t and angular frequency ω:

Plot3D[Abs[sinetfr[t, ω]], {t, tmin, tmax}, {ω, ωmin, 2}, AxesLabel { ...                                                                                               f\

[Graphics:HTMLFiles/tutorial_43.gif]

In this case, the TFR exhibits a ridge at constant frequency ω=1, indicating the presence of a single, constant-frequency component as expected.

Chirp

Consider a more interesting example signal f(t)=Sin(t^2):

ϕ[t_] := t^2

f[t_] := Sin[ϕ[t]]

A direct plot of this function shows that its instantaneous frequency ω_f(t) increases away from t=0 but the functional form of ω_f(t) is not easily determined from this plot:

Plot[f[t], {t, -3π, 3π}, AxesLabel {"t", "f(t)"}] ;

[Graphics:HTMLFiles/tutorial_50.gif]

Let's skip the CWT this time and compute the time-frequency representation directly, using more samples than before to capture the added complexity of this signal:

chirptfr = FunctionTFR[f[t], {t, -4π, 4π, 512}] ;//Timing

{5.67214 Second, Null}

The ranges of time t and angular frequency ω for which this TFR is valid are:

{{ωmin, ωmax}, {tmin, tmax}} = {GetTFRωRange[chirptfr], GetTFRtRange[chirptfr]}

{{0.125, 64.}, {-12.5664, 12.5664}}

This TFR is substantially more interesting that the previous one, and clearly conveys useful information about the instantaneous frequency and amplitude of the signal as a function of time:

Plot3D[Abs[chirptfr[t, ω]], {t, -3π, 3π}, {ω, ωmin, 6π}, AxesL ...                                                                                                f

[Graphics:HTMLFiles/tutorial_56.gif]

Specifically, the "V"-shaped ridge indicates an instantaneous frequency which can be seen to increase linearly away from t=0. This is even clearer on a contour plot:

ContourPlot[Abs[chirptfr[t, ω]], {t, -3π, 3π}, {ω, ωmin, 6π}, FrameLabel {"t", "ω"}, PlotRangeAll] ;

[Graphics:HTMLFiles/tutorial_58.gif]

The instantaneous angular frequency of a component in the signal is the position of its ridge in its time-frequency representation. As expected, there in a single dominant ridge which is seen to lie mostly along the line ω(t)=2|t|.

The trajectory of the ridge is easily extracted numerically using Mathematica's built-in FindMaximum function:

ω[t_] := ω/.Last[FindMaximum[Abs[chirptfr[t, ω]], {ω, Max[1, 2Abs[t]] - 1/2, Max[1, 2Abs[t]] + 1/2}]]

Note: supplying a bound of two initial values instead of one value makes Mathematica's FindMaximum function considerably more efficient.

Comparing the numerical approximation (blue) with the expected result (red) shows that the instantaneous frequency has been determined very accurately, except for the region t=0 where the expected value (oscillations with zero frequency) is dubious:

Plot[{2Abs[t], ω[t]}, {t, -3π, 3π}, PlotRange {0, 6π}, AxesLabel ...                                                                                                f

[Graphics:HTMLFiles/tutorial_61.gif]

Due to the feature at t=0, the instantaneous amplitude A_f of this signal is also interesting. The instantaneous amplitude of the signal is simply the modulus of the TFR on the ridge:

Plot[{1, Abs[chirptfr[t, ω[t]]]}, {t, -3π, 3π}, PlotRange {0, 1.1}, A ...                                                                                                f

[Graphics:HTMLFiles/tutorial_64.gif]

In addition to the expected dip in instantaneous amplitude at t=0, this function also exhibits significant unwanted oscillation and systematic deviation from the desired value of 1 for |t|>>0. These artefacts can be reduced by increasing the spectral resolution of the transform, at the cost of lower temporal resolution. This is done by specifying a larger value for the parameter σ of the wavelet. The previous transform used the default parameter value:

GetCWTσ[chirptfr〚1〛]

1/4 (2 + 2^(1/2))

N[%]

0.853553

So, let's recalculate the transform with the parameter σ=2:

chirptfr = FunctionTFR[f[t], {t, -4π, 4π, 512}, Parameter2] ;//Timing

{6.18206 Second, Null}

Plot3D[Abs[chirptfr[t, ω]], {t, -3π, 3π}, {ω, ωmin, 6π}, AxesL ...                                                                                                f

[Graphics:HTMLFiles/tutorial_72.gif]

Export["tfr.eps", %]

tfr.eps

The instantaneous frequency is comparatively unchanged in this case:

Plot[{2Abs[t], ω[t]}, {t, -3π, 3π}, PlotRange {0, 6π}, AxesLabel ...                                                                                                f

[Graphics:HTMLFiles/tutorial_76.gif]

In contrast, the oscillations in the instantaneous amplitude are eliminated:

Plot[{1, Abs[chirptfr[t, ω[t]]]}, {t, -3π, 3π}, PlotRange {0, 1.1}, A ...                                                                                                f

[Graphics:HTMLFiles/tutorial_78.gif]

Superimposed Chirps

The CWT-based approach to time-frequency analysis is also able to handle polychromatic signals, those which contain multiple frequency components at the same time.

Consider a polychromatic signal composed of two superimposed chirps:

f[t_] := Sin[t^2] + Sin[11^(1/2) t^2]

Plot[f[t], {t, π, 3π}, AxesLabel {"t", "f(t)"}] ;

[Graphics:HTMLFiles/tutorial_81.gif]

chirp2tfr = FunctionTFR[f[t], {t, 0, 4π, 512}, Parameter7] ;//Timing

{6.16306 Second, Null}

{{ωmin, ωmax}, {tmin, tmax}} = {GetTFRωRange[chirp2tfr], GetTFRtRange[chirp2tfr]}

{{0.25, 128.}, {0., 12.5664}}

Plot3D[Abs[chirp2tfr[t, ω]], {t, π, 3π}, {ω, ωmin, 60}, AxesLabel&# ...                                                                                                f

[Graphics:HTMLFiles/tutorial_87.gif]

ContourPlot[Abs[chirp2tfr[t, ω]], {t, π, 3π}, {ω, ωmin, 211^(1/2) 3π}, FrameLabel {"t", "ω"}, PlotRangeAll] ;

[Graphics:HTMLFiles/tutorial_89.gif]

Two ridges are clearly visible, although the results are not as accurate as those in the preceding examples due to interference between the two signal components.

The instantaneous frequencies can be extracted as before:

ω1[t_] := ω/.Last[FindMaximum[Abs[chirp2tfr[t, ω]], {ω, 0.5 2Abs[t], 2 2Abs[t]}]]

ω2[t_] := ω/.Last[FindMaximum[Abs[chirp2tfr[t, ω]], {ω, 0.8 211^(1/2) Abs[t], 1.2 211^(1/2) Abs[t]}]]

Plot[{2Abs[t], ω1[t], 211^(1/2) Abs[t], ω2[t]}, {t, π, 3π}, AxesLabel {"t", "ω(t)"}] ;

[Graphics:HTMLFiles/tutorial_93.gif]

Despite the interference, the extracted instantaneous frequencies are clearly accurate enough for a wide variety of applications.

Noise

When analyzing measured signals, knowledge of the CWT of random noise can be useful in order to distinguish between noise and signal features.

As a trivial example, consider the following noise:

f = Interpolation[Re[Fourier[Table[Random[] ^( 2π Random[]), {t, 256}]]]] ;

Plot[f[t], {t, 1, 256}, PlotRangeAll, AxesLabel {"t", "f(t)"}] ;

[Graphics:HTMLFiles/tutorial_96.gif]

noisecwt = FunctionCWT[f[t], {t, 1, 256, 256}] ;//Timing

{1.46078 Second, Null}

{{logamin, logamax}, {bmin, bmax}} = {GetCWTlogaRange[noisecwt], GetCWTbRange[noisecwt]}

{{-1.14864, 4.39653}, {1., 256.}}

The CWT of this noise possesses low coefficients at large scales and much higher coefficients at smaller scales:

Plot3D[Abs[noisecwt[loga, b]], {b, bmin, bmax}, {loga, logamin, logamax}, AxesLabel  ...                                                                                           Ψ

[Graphics:HTMLFiles/tutorial_102.gif]

Scaling by a^(1/2) to recover the CWT as an "energy density" (after Morlet et al.) shows a more constant distribution of energies over different scales and translations:

Plot3D[^loga^(1/2) Abs[noisecwt[loga, b]], {b, bmin, bmax}, {loga, logamin, logamax} ...                                                                                           Ψ

[Graphics:HTMLFiles/tutorial_105.gif]

Create a TFR from the CWT:

noisetfr = TFR[noisecwt] ;

{{ωmin, ωmax}, {tmin, tmax}} = {GetTFRωRange[noisetfr], GetTFRtRange[noisetfr]}

{{0.01232, 3.15391}, {1., 256.}}

As expected, the TFR exhibits high amplitudes at high frequencies, due to the sharp features in the noisy signal:

Plot3D[Abs[noisetfr[t, ω]], {t, tmin, tmax}, {ω, ωmin, ωmax}, AxesLabel& ...                                                                                                f

[Graphics:HTMLFiles/tutorial_110.gif]

Advanced Options

Several facilities are provided for more advanced users:

? CWT`*

CWT`
calcFdata FunctionCWT GetCWTF GetCWTn GetTFRωRange SpectralRange TFRCombinedPlot
CWT FunctionTFR GetCWTlogaRange GetCWTσ LogScaleRange TemporalRange TranslationRange
CWTContourPlot GetCWTbRange GetCWTMatrix GetTFRtRange Parameter TFR
Restricted Scale and Translation

By default, the CWT is computed over ranges of scale equivalent to the Nyquist frequency in Fourier analysis, and translation over the whole signal:

ω_Ψ/(π n)(t_max-t_min)≤a<ω_Ψ/π(t_max-t_min)
t_min≤b<t_max
where ω_Ψ is the central frequency of the mother wavelet.

The TFR is computed over equivalent ranges:

(t_max - t_min)/(π n)≤ω<(t_max - t_min)/π
t_min≤t<t_max

Either of these ranges may be arbitrarily restricted to smaller regions. Restricting the range of scales reduces both the amount of computation and the amount of storage required for the transform. As the CWT implementation is based upon cyclic convolutions, restricting the range of translations only reduces the overall storage requirements.

The ranges that a CWT is computed over can be restricted by specifying either or both of the LogScaleRange and TranslationRange optional arguments to FunctionCWT. For example:

FunctionCWT[Sin, {0, 8π, 256}, 3, LogScaleRange {-2, 0}, TranslationRange {0, 2π}]

FunctionCWT[Sin, {0, 8 π, 256}, 3, LogScaleRange {-2, 0}, TranslationRange {0, 2 π}]

The ranges that a TFR is computed over can be restricted by specifying either or both of the SpectralRange and TemporalRange optional arguments to FunctionCWT. For example:

ByteCount[FunctionTFR[Sin[t], {t, 0, 8π, 256}]]//Timing

{1.46578 Second, 1316324}

ByteCount[FunctionTFR[Sin[t], {t, 0, 8π, 256}, SpectralRange {8, 32}]]//Timing

{0.439933 Second, 336876}

ByteCount[FunctionTFR[Sin[t], {t, 0, 8π, 256}, TemporalRange {0, 2π}]]//Timing

{1.48677 Second, 336876}

Thus, restricting the ranges of scales, translations, frequencies and times is the simplest way to make an analysis tractable.

Property Extraction

When performing wavelet analysis on many data sets, the ability to extract the properties of a transform can be useful. For example, to determine the number of samples used to perform the transform.

The CWT package provides several functions for extracting transform properties:

? GetCWTn

GetCWTn[cwt] extracts the number of samples contained in the original sampling of the input function which was used to generate the given continuous wavelet transform object.

? GetCWTσ

GetCWTσ[cwt] extracts the wavelet parameter used to generate the given continuous wavelet transform function.

The most advanced users may even wish to extract the innards of the transform, in discrete form:

? GetCWTF

GetCWTF[cwt] extracts the interpolating function which represents the continuous wavelet transform over log scale and translation.

As these functions can be used to extract the properties of a transform from the result, sets of transformed data can be marshalled to and from disc using the Mathematica functions DumpSave and Get.

Examples

We have analysed a suite of example signals from various disciplines: