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:

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π ) with linearly increasing frequency:

The instantaneous frequency ν(t) of f(t) is the derivative of the phase: ν(t)=||=|t|.

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

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:

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:

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

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:

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:

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:

The CWT is most easily visualised as a 3D surface plot of the modulus |(a,b)|:

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:

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:

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:

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

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():

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

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:

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

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:

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:

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:

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:

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

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:

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

The instantaneous frequency is comparatively unchanged in this case:

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

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:

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:

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:

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

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

Create a TFR from the CWT:

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

Advanced Options

Several facilities are provided for more advanced users:

CWT`

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:

(-)≤a<(-)
≤b<
where is the central frequency of the mother wavelet.

The TFR is computed over equivalent ranges:

≤ω<
≤t<

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:

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:

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:

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

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:

• Physics
• Finance
• Biology
• Bioinformatics
• Engineering