Open In App

Complex Demodulation Phase and Amplitude Plot

Improve
Improve
Like Article
Like
Save
Share
Report

Complex Demodulation Phase Plot

In the frequency analysis of time series models, a common model is a sinusoidal wave:

Y_{i} = C + \alpha\sin{(2\pi\omega t_{i} + \phi)} + E_{i}

where, ∝ is the amplitude, phi is the phase shift and omega is the dominant frequency. The goal of the complex demodulation plot is to improve the frequency estimate.


 The complex demodulation plots formed by two components:

  • Vertical axis: Phase
  • Horizontal axis: Time

Importance of a Good Initial Estimate for the Frequency:

The non-linear fitting for the sinusoidal waves

Y_{i} = C + \alpha\sin{(2\pi\omega t_{i} + \phi)} + E_{i}  

The above equation is sensitive to good initial values. The initial value of frequency omega can be obtained by the spectral plot. The complex demodulation phase plot is used to assess whether this estimate is adequate. If the estimate is not adequate then whether it should be increased and decreased.

Complex Demodulation Amplitude Plot:

In the frequency analysis of time series models, a common model is a sinusoidal wave:

Y_{i} = C + \alpha\sin{(2\pi\omega t_{i} + \phi)} + E_{i}  

where, ∝ is the amplitude, phi is the phase shift and omega is the dominant frequency.

If the slope of the complex demodulation amplitude plot is not zero, then the above equation is finally replaced by the model.

Y_{i} = C + \alpha_{i}\sin{(2\pi\omega t_{i} + \phi)} + E_{i}

where, ai is some type of linear model fit with standard least squares. The most common case is linear fit, that is the model becomes as follows:

Y_{i} = C + (B_0 + B_1*t_{i})\sin{(2\pi\omega t_{i} + \phi)} + E_{i}

The complex demodulation amplitude plot is formed by

  • Vertical axis: Amplitude
  • Horizontal axis: Time

 The complex amplitude demodulation plot answers the following questions:

  • Does amplitude change over time?
  • Is there any start-up effect that led to the amplitude being unstable at the start?
  • Is there any outliers in amplitude?

Implementation:
Code: Python code for complex demodulation analysis 

python3

# code
import numpy as np
import matplotlib.pyplot as plt
 
def gen_test_data(periods, noise=0, rotary=False, npts=1000, dt=1.0/24):
    """
    Generate a simple time series for testing complex demodulation.
     
    *periods* is a sequence with the periods of one or more
        harmonics that will be added to make the test signal.
        They can be positive or negative.
    *noise* is the amplitude of independent Gaussian noise.   
    *rotary* is Boolean; if True, the test signal is complex.
    *npts* is the length of the series.
    *dt* is the time interval (default is 1.0/24)
     
    Returns t, x: ndarrays with the test times and test data values..
    """    
         
    t = np.arange(npts, dtype=float) * dt
     
    if rotary:
        x = noise * (np.random.randn(npts) + 1j * np.random.randn(npts))
    else:
        x = noise * np.random.randn(npts)
 
    for p in periods:
        if rotary:
            x += np.exp(2j * np.pi * t / p)
        else:
            x += np.cos(2 * np.pi * t / p)
     
    return t, x
   
def complex_demodulation(t, x, central_period, hwidth = 2):
    """
    Complex demodulation of a real or complex series, *x*
    of samples at times *t*, assumed to be uniformly spaced.
     
    *central_period* is the period of the central frequency
        for the demodulation.  It should be positive for real
        signals. For complex signals, a positive value will
        return the CCW rotary component, and a negative value
        will return the CW component (negative frequency).
        Period is in the same time units as are used for *t*.
 
    *hwidth* is the Blackman filter half-width in units of the
        *central_period*.  For example, the default value of 2
        makes the Blackman half-width equal to twice the
        central period.
     
    Returns a dictionary.
    """    
     
    rotary = x.dtype.kind == 'c'  # complex input
     
    # Make the complex exponential for demodulation:
    c = np.exp(-1j * 2 * np.pi * t / central_period)
     
    product = x * c
     
    # filter half-width number of points
    dt = t[1] - t[0]
    hwpts = int(round(hwidth * abs(central_period) / dt))
    nf = hwpts * 2 + 1
    x1 = np.linspace(-1, 1, nf, endpoint=True)
    x1 = x1[1:-1]   # chop off the useless endpoints with zero weight
    w1 = 0.42 + 0.5 * np.cos(x1 * np.pi) + 0.08 * np.cos(x1 * 2 * np.pi)
    ytop = np.convolve(product, w1, mode='same')
    ybot = np.convolve(np.ones_like(product), w1, mode='same')
    demod = ytop/ybot
    if not rotary:   
        # The factor of 2 below comes from fact that the
        # mean value of a squared unit sinusoid is 0.5.
        demod *= 2
     
    reconstructed = (demod * np.conj(c))
    if not rotary:
        reconstructed = reconstructed.real
         
    if np.sign(central_period) < 0:
        demod = np.conj(demod)
        # This is to make the phase increase in time
        # for both positive and negative demod frequency
        # when the frequency of the signal exceeds the
        # frequency of the demodulation.
     
    d = {'t':t,'signal' : x,'hwpts' : hwpts,'demod': demod,'reconstructed' : reconstructed}
    return d
   
def plot_demodulation(dm):
    fig, axs = plt.subplots(3, sharex=True)
    resid = dm.get('signal') - dm.get('reconstructed')
    if dm.get('signal').dtype.kind == 'c':
        axs[0].plot(dm.get('t'), dm.get('signal').real, label='signal.real')
        axs[0].plot(dm.get('t'), dm.get('signal').imag, label='signal.imag')
        axs[0].plot(dm.get('t'), resid.real, label='difference real')
        axs[0].plot(dm.get('t'), resid.imag, label='difference imag')
    else:   
        axs[0].plot(dm.get('t'), dm.get('signal'), label='signal')
        axs[0].plot(dm.get('t'), dm.get('reconstructed'), label='reconstructed')
        axs[0].plot(dm.get('t'), dm.get('signal') - dm.get('reconstructed'), label='difference')
     
    axs[0].legend(loc='upper right', fontsize='small')
     
    axs[1].plot(dm.get('t'), np.abs(dm.get('demod')), label='amplitude', color='C3')
    axs[1].legend(loc='upper right', fontsize='small')
     
    axs[2].plot(dm.get('t'), np.angle(dm.get('demod'), deg=True), '.', label='phase',
                color='C4')
    axs[2].set_ylim(-180, 180)
    axs[2].legend(loc='upper right', fontsize='small')
     
    for ax in axs:
        ax.locator_params(axis='y', nbins=5)
    return fig, axs
 
def test_demodulation(periods, central_period,
                 noise=0,
                 rotary=False,
                 hwidth = 1,
                 npts=1000,
                 dt=1.0/24):
     
    t, x = gen_test_data(periods, noise=noise, rotary=rotary,
                     npts=npts, dt=dt)
    dm = complex_demodulation(t, x, central_period, hwidth=hwidth)
    fig, axs = plot_demodulation(dm)
    return fig, axs, dm
# Example 1
test_demodulation([12.0/24], 12.0/24);
# Example 2
test_demodulation([11.0/24], 12.0/24)

                    

Example 1: Signal, amplitude demodulation, and phase demodulation

Example 2: Signal, amplitude demodulation, and phase demodulation



Last Updated : 14 Aug, 2021
Like Article
Save Article
Previous
Next
Share your thoughts in the comments
Similar Reads