Theory

This page describes the mathematical background of the Stockwell transform (S-transform) and its implementation in this package.

Continuous S-Transform

The S-transform of a continuous signal \(x(t)\) is defined by Stockwell et al. [1996]:

\[S_x(t, f) = \int_{-\infty}^{\infty} x(\tau) \, w(t - \tau, f) \, e^{-j 2\pi f \tau} \, d\tau\]

where \(w(t, f)\) is a frequency-dependent window function. The choice of window determines the time-frequency resolution trade-off; two window types are available (see Window functions below).

The inverse S-transform recovers the original signal by integrating over time and frequency:

\[x(\tau) = \int_{-\infty}^{\infty} \left[ \int_{-\infty}^{\infty} S_x(t, f) \, dt \right] e^{j 2\pi f \tau} \, df\]

Discrete Implementation

This package implements the discrete S-transform in the frequency domain, following the fast algorithm of Brown et al. [2010] (see also Stockwell [1999]). Given a discrete signal \(x[k]\) of length \(N\), the algorithm proceeds as follows:

  1. Compute the discrete Fourier transform (DFT) \(X[m]\) of \(x[k]\), where \(m = 0, \dots, N-1\) is the frequency bin index.

  2. Apply the Hilbert transform in the frequency domain: double the positive frequencies and zero the negative frequencies.

  3. For each frequency voice \(n\) (from lo to hi), where \(n\) is the frequency row index in the output S-transform:

    1. Compute the frequency-domain window \(W[n, m]\).

    2. Shift the spectrum: multiply \(X[m]\) by \(W[n, (m-n) \bmod N]\).

    3. Inverse DFT to obtain the \(n\)-th row of the S-transform.

The result is a complex matrix of shape \((\text{nfreqs}, N)\), where the first dimension runs over frequency and the second over time.

Window functions

The S-transform is implemented in the frequency domain: for each frequency voice \(n\), the spectrum \(X[m]\) is multiplied by a shifted frequency-domain window \(W[n, m]\) and then inverse-transformed.

Two window functions are available. Both depend on the gamma parameter described below.

Gaussian window (default, win_type='gauss')

In the time domain the window is a scaled Gaussian, following the original definition of Stockwell et al. [1996]:

\[w(t, f) = \frac{|f|}{\gamma\sqrt{2\pi}} \, e^{-t^2 f^2 / (2\gamma^2)}\]

Its Fourier transform gives the frequency-domain form used in the implementation:

\[G(n, m, \gamma) = \exp\left( -\frac{2\pi^2 m^2 \gamma^2}{n^2} \right)\]

where \(n\) is the frequency voice index and \(m\) the frequency sample index.

Note

The Wikipedia article on the S-transform uses an alternative Gaussian convention: \(|f| e^{-\pi t^2 f^2}\), whose frequency-domain counterpart is \(e^{-\pi m^2 / n^2}\). Both forms are equivalent up to a different choice of the \(\gamma\) scaling. This package follows the original convention of Stockwell et al. [1996].

Kazemi window (win_type='kazemi')

Introduced by Kazemi et al. [2014], this window provides sharper frequency localisation. In the frequency domain:

\[K(n, m, \gamma) = \frac{1}{1 + \left( \dfrac{m^2 \gamma}{n} \right)^2}\]

Gamma parameter

The parameter \(\gamma > 0\) (default 1) controls the number of Fourier sinusoidal periods within one standard deviation of the window: raising \(\gamma\) increases frequency resolution at the expense of time resolution, and vice versa. Typical values are around 1; very large values produce an extremely narrow window in time.

Hilbert Transform

The Hilbert transform \(\mathcal{H}\{x(t)\}\) of a real signal \(x(t)\) is computed in the frequency domain:

\[\begin{split}\mathcal{F}\{\mathcal{H}\{x\}\}[m] = \begin{cases} 2 \, X[m] & 0 < m < \lceil N/2 \rceil \\ X[m] & m = 0 \;\text{or}\; m = N/2 \;(\text{if } N \text{ even}) \\ 0 & m \geq \lfloor N/2 \rfloor + 1 \end{cases}\end{split}\]

where \(\mathcal{F}\) denotes the Fourier transform and \(X[m] = \mathcal{F}\{x\}[m]\). In words: positive frequencies (excluding DC and Nyquist) are doubled, negative frequencies are zeroed, and DC (and Nyquist for even \(N\)) are left unchanged.

Sine Tapers

The \(K\)-th sine taper of length \(N\) is defined as Riedel and Sidorenko [1979]:

\[w_K[n] = \sqrt{\frac{2}{N+1}} \, \sin\left( \frac{\pi (K+1) (n+1)}{N+1} \right), \quad n = 0, \dots, N-1\]

Acknowledgments

The C implementation of the Stockwell transform is based on original code from the NIMH MEG Core Facility, released into the public domain.

The following modifications have been made to the original code:

  • Windows compatibility: added portable defines and MSVC export declarations.

  • ctypes interface (v1.1): replaced Python C extension modules with ctypes-based loading for easier cross-platform packaging.

  • Gamma parameter (v1.0.5): added \(\gamma\) to control the time-frequency resolution trade-off.

  • Kazemi window (v1.0.5): added the alternative window of Kazemi et al. [2014] alongside the original Gaussian.

  • FFTW plan caching: static FFTW plans are reused across calls, with proper cleanup at process exit to prevent segfaults.

  • Hilbert phase fix: corrected the Hilbert phase by no longer shifting the spectrum in frequency when multiplying by the window, preserving the correct instantaneous phase in the Stockwell spectrum.

  • Input validation: reject empty arrays and non-integer lo/hi parameters.

References

[1] (1,2,3)

R. G. Stockwell, L. Mansinha, and R. P. Lowe. Localization of the complex spectrum: the S transform. IEEE Transactions on Signal Processing, 44(4):998–1001, 1996. doi:10.1109/78.492555.

[2]

R. A. Brown, M. L. Lauzon, and R. Frayne. A general description of linear time-frequency transforms and formulation of a fast, invertible transform that samples the continuous S-transform spectrum nonredundantly. IEEE Transactions on Signal Processing, 58(1):281–290, 2010. doi:10.1109/TSP.2009.2028972.

[3]

R. G. Stockwell. S-transform analysis of gravity wave activity from a small scale network of airglow imagers. PhD thesis, University of Western Ontario, 1999.

[4] (1,2)

K. Kazemi, M. Amirian, and M. J. Dehghani. The S-transform using a new window to improve frequency and time resolution. Signal, Image and Video Processing, 8(3):533–541, 2014. doi:10.1007/s11760-013-0551-1.

[5]

K. S. Riedel and A. Sidorenko. Sine taper spectral estimation. IEEE Transactions on Acoustics, Speech, and Signal Processing, 27(5):534–535, 1979.