Complex Demodulation

NIST/SEMATECH Section 1.3.3.8-9 Complex Demodulation

What It Is

Complex demodulation (Granger, 1964) is a time-series analysis technique that extracts the time-varying amplitude and phase of a sinusoidal component at a specified frequency. It produces two companion plots: an amplitude plot showing how the strength of the cyclic component changes over time, and a phase plot used to improve the estimate of the frequency in the sinusoidal model.

The underlying sinusoidal model is $Y_i = C + \alpha\sin(2\pi\omega t_i + \phi) + E_i$ , where $\alpha$ is the amplitude, $\phi$ is the phase shift, and $\omega$ is the dominant frequency. In this model $\alpha$ and $\phi$ are assumed constant. The amplitude plot displays estimated amplitude vs. time, and the phase plot displays phase angle vs. time. If the amplitude is not constant, the model is typically replaced with $Y_i = C + (B_0 + B_1 t_i)\sin(2\pi\omega t_i + \phi) + E_i$ , where the amplitude varies linearly with time. Quadratic amplitude models are sometimes used; higher order models are relatively rare.

Questions This Plot Answers

Does the amplitude change over time?
Are there any outliers that need to be investigated?
Is the amplitude different at the beginning of the series (start-up effect)?
Is the specified demodulation frequency correct?

Why It Matters

In the sinusoidal model $Y_i = C + \alpha\sin(2\pi\omega t_i + \phi) + E_i$ , the amplitude $\alpha$ is assumed constant. The amplitude plot checks whether this assumption is reasonable; if not, $\alpha$ should be replaced with a time-varying model. The non-linear fitting for the sinusoidal model is usually quite sensitive to the choice of good starting values. The initial estimate of $\omega$ is obtained from a spectral plot, and the phase plot assesses whether this estimate is adequate. Using the phase plot together with the spectral plot can significantly improve the quality of the non-linear fits obtained.

When to Use a Complex Demodulation

The amplitude plot is used to determine if the assumption of constant amplitude in the sinusoidal model is justifiable. If the amplitude changes over time, the constant $\alpha$ should be replaced with a time-varying model such as a linear fit ( $B_0 + B_1 t$ ) or quadratic fit. The phase plot is used to improve the estimate of the frequency $\omega$ : if the phase plot shows lines sloping from left to right, the frequency estimate should be increased; if lines slope right to left, it should be decreased; if the slope is essentially zero, the frequency estimate does not need to be modified.

How to Interpret a Complex Demodulation

In the amplitude plot, the vertical axis shows the estimated amplitude as a function of time. A flat amplitude trace indicates a stationary cyclic component with constant $\alpha$ , while trends or abrupt changes indicate that $\alpha$ varies with time and the model should be updated accordingly. In the phase plot, the vertical axis shows the phase angle as a function of time. Lines sloping from left to right indicate that the demodulation frequency should be increased. Lines sloping from right to left indicate the frequency should be decreased. An essentially flat phase trace confirms that the specified demodulation frequency is correct.

Assumptions and Limitations

Complex demodulation requires specifying the target frequency in advance, which typically comes from a preliminary spectral plot. It assumes the signal can be approximated by a sinusoidal model at the chosen frequency, and it uses a low-pass filter whose bandwidth must be selected to balance smoothing and time resolution. The mathematical computations for determining the amplitude and phase are detailed in Granger (1964).

See It In Action

This technique is demonstrated in the following case studies:

Beam Deflections Case Study

Reference: NIST/SEMATECH e-Handbook of Statistical Methods, Sections 1.3.3.8-9

Formulas

Sinusoidal Model (Constant Amplitude)

Y_i = C + \alpha\sin(2\pi\omega\, t_i + \phi) + E_i

The standard sinusoidal time series model where alpha is the amplitude, phi is the phase shift, omega is the dominant frequency, C is a constant, and E_i is the error term. Both alpha and phi are assumed constant over time.

Time-Varying Amplitude Model

Y_i = C + (B_0 + B_1\, t_i)\sin(2\pi\omega\, t_i + \phi) + E_i

When the amplitude plot shows that alpha is not constant, it is replaced by a linear function of time. This is the most common replacement; quadratic models are sometimes used, but higher order models are relatively rare.

Python Example

import numpy as np
import matplotlib.pyplot as plt

# LEW.DAT — beam deflection data (NIST Section 1.4.2.5)
# True dominant frequency ≈ 0.3026 cycles/observation.
lew = np.array([
    -213, -564,  -35,  -15,  141,  115, -420, -360,  203, -338,
    -431,  194, -220, -513,  154, -125, -559,   92,  -21, -579,
     -52,   99, -543, -175,  162, -457, -346,  204, -300, -474,
     164, -107, -572,   -8,   83, -541, -224,  180, -420, -374,
     201, -236, -531,   83,   27, -564, -112,  131, -507, -254,
     199, -311, -495,  143,  -46, -579,  -90,  136, -472, -338,
     202, -287, -477,  169, -124, -568,   17,   48, -568, -135,
     162, -430, -422,  172,  -74, -577,  -13,   92, -534, -243,
     194, -355, -465,  156,  -81, -578,  -64,  139, -449, -384,
     193, -198, -538,  110,  -44, -577,   -6,   66, -552, -164,
], dtype=float)
n = len(lew)
filter_hw = 30  # moving-average half-width

def demodulate(data, freq, hw):
    """Complex demodulation at a given frequency."""
    t = np.arange(len(data))
    kernel = np.ones(2 * hw + 1) / (2 * hw + 1)
    re = np.convolve(data * np.cos(2 * np.pi * freq * t), kernel, mode="same")
    im = np.convolve(data * np.sin(2 * np.pi * freq * t), kernel, mode="same")
    return re, im

# Amplitude: use correct frequency (0.3025) for stable envelope
re_a, im_a = demodulate(lew, 0.3025, filter_hw)
amplitude = 2 * np.sqrt(re_a**2 + im_a**2)

# Phase: use wrong frequency (0.28) to show sawtooth drift
re_p, im_p = demodulate(lew, 0.28, filter_hw)
phase = np.arctan2(im_p, re_p)  # wrapped radians

# Trim edge points where the moving average has partial windows
trim = slice(filter_hw, n - filter_hw)
t = np.arange(n)

fig, axes = plt.subplots(2, 1, figsize=(10, 6))
axes[0].plot(t, amplitude)
axes[0].set_ylabel("Estimated Amplitude")
axes[0].set_title("Complex Demodulation Amplitude Plot (freq = 0.3025)")

axes[1].scatter(t[trim], phase[trim], s=12, marker="x", color="k")
axes[1].set_xlabel("Time")
axes[1].set_ylabel("Estimated Phase")
axes[1].set_title("Complex Demodulation Phase Plot (freq = 0.28)")
plt.tight_layout()
plt.show()