Complex Demodulation & Hilbert Transform

In physical oceanography we can use 'Complex Demodulation' to depict the time series of the tidal amplitude of a time series. Despite the fancy name, what someone may want to use to impress someoneelse (likely the paper's reviwer?), it is a relatively simple method.

In essence, it separates the information from a signal carrier (as in radio signals!). The information we want is the tidal amplitude along the time, carried mainly by the semi-diurnal tide (M2).

Here you find a more detailled explanation about what is complex demodulation (Hawaii stuff, from now on)

https://currents.soest.hawaii.edu/ocn_data_analysis/_static/complex_demod.html

By the way, the Hawaii stuff was the 2.nd result from Google for 'python complex demodulation tides'!

Our goal is to apply this method to a water level record to depict the sprint-neap variations on the amplitude.

We'll use a long time series recorded off Bragança, Pará, which has a nice meso/macro tidal regime!

And what about Hilbert? Well... as I was working in this notebook, I got poor results from the 'Complex Demodulation', and found out that it can be also obtained with such 'Hilbert Transform'... which is a sort of demodulation too, but not 'complex', in the sense it work only with real numbers.

import pickle
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import numpy as np

from scipy import signal

from scipy.signal import hilbert, chirp

Bellow are two filters (Blackman and Butterworth)... it is necessary a filter for the demodulation!

# from Hawwaii stuff!

def bl_filt(y, half_width):
    """
    Simple Blackman filter.
    
    The end effects are handled by calculating the weighted
    average of however many points are available, rather than
    by zero-padding.
    """
    nf = half_width * 2 + 1
    x = np.linspace(-1, 1, nf, endpoint=True)
    x = x[1:-1]   # chop off the useless endpoints with zero weight
    w = 0.42 + 0.5 * np.cos(x * np.pi) + 0.08 * np.cos(x * 2 * np.pi)
    ytop = np.convolve(y, w, mode='same')
    ybot = np.convolve(np.ones_like(y), w, mode='same')
    
    return ytop / ybot

#####################################################################
# low pass Butterworth filter for frequency sampling fs = 1

def lowpassfilter(x):
    fs = 1
    order = 5
    cutoff = 1/30
    
    b, a = signal.butter(order, cutoff, fs=fs, btype='lowpass')
    xf = signal.filtfilt(b, a, x)

    return xf

Uploading the tidal record!

This is a hourly dataset 1.4 year long. Notice the nice spring/neap variation of the sea level!

path = r'D:\\GUTO\\1_Trabs\\2_Plataforma\\Para\\ADCP_fundeio20152017\\'

with open(path + 'Braganca_sealevel.pkl', 'rb') as io:
    sl = pickle.load(io)

time = sl[0]
sealevel = sl[1]

plt.figure(figsize=(22,4))
plt.plot(time, sealevel)
plt.show()

Context:

I start studying the 'Hawaii' stuff... getting the code and making the reverse engineering to understad each step. I do not fully understand all steps, but I must be sure I doing it right! As it is said, I don't care if the duck is male, I want the egg!

+

Below is the 'de-functionized' function from the Hawaii stuff, which I tweak here ant there... however I didn't get the desired result, as it will be showed.

t = mdates.date2num(time)*24
x = sealevel

# this is very new to me! It is a way to test the content of the array! It produces a boolean (True or False). 
# https://numpy.org/doc/stable/reference/generated/numpy.dtype.kind.html
rotary = x.dtype.kind == 'c'

# not sure, but I think this is the carrier frequency period! Changing it changes enormously the result!
central_period = 12 + 25/60
# central_period = 12

# for the filter design! 
hwidth = 2

# carrier signal. Below there is some exploration about what this is!
c = np.exp(-1j * 2 * np.pi * t / central_period)

# the product between the signal and signal carrier (plot below)
product = x * c

# the sampling interval (1 hour)
dt = t[1] - t[0]

# for the filter...
hwpts = int(round(hwidth * abs(central_period) / dt))

# filtering with Blackman and Butterworth
demod = bl_filt(product, hwpts)
demod2 = lowpassfilter(product)

# this sintax is also new to me! It is the same as 'if rotary == False:'
# in the source, there is this comment:
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
    demod2 *= 2

reconstructed = demod * np.conj(c)
reconstructed2 = demod * np.conj(c)

# this reconstruct the original signal... 
if not rotary:    
    reconstructed = reconstructed.real
    reconstructed2 = reconstructed2.real

# do this in the case we are working with complex numbers
if np.sign(central_period) < 0:
    demod = np.conj(demod)


fig, axs = plt.subplots(2, figsize=(22,4))
axs[0].plot(time, sealevel)
axs[1].plot(time, demod)
axs[1].plot(time, demod2-0.5) # both filters do the same, the the Blackman dont screw up at the extremities!

[<matplotlib.lines.Line2D at 0x2ba9260f1d0>]

This is not the expected result... the expected result for the present case would be a time series of the tidal amplitude, always positive! I changed the central_period, which defines the carrier signal, but didn't get nothing better! I think there is something wrong with this code (unlikely) or the method doesn't work with the present data (?!) or I doing something dumb!

Of course... keep trying. Further down on the Google results I found other solution!

The cells below are the 'code' exploring what are the things in the Hawaii solution ...

Further below is another way to 'demodulate' I found! Which, as I learned, is another method to do the same thing! The Hilbert stuff!

I tried also the ChatGPT, however I got the same result of Hawaii Stuff, which supports the hypothesis the code is right. So, what's the problem?

plt.figure(figsize=(22,4))
plt.plot(time, reconstructed)

[<matplotlib.lines.Line2D at 0x2ba9269ce90>]

Understanding the 'dtype.kind'

x = np.array([1.1, 2.2])
z = x.dtype.kind == 'c'
z

False

x = np.array([1.1, 2.2])
z = x.dtype.kind == 'f'
z

True

x.dtype.kind

'f'

Exploring the 'if not'

xx = False

if xx:
    print('hi')
else:
    print('by')

if not xx:
    print('HI')
else:
    print('BY')

by
HI

trying to understand the 'Signal carrier' from the Hawaii stuff... I didn't! kkkk

The signal carrier is a complex time series, which the real part is 90 degrees lagged from the imaginary... but what is curious is that neither is exacly in phase whith the signal!

# separating the real and imaginary parts of c
c_real = [v.real for v in c]
c_imag = [v.imag for v in c]

fig, (ax1, ax2) = plt.subplots(2, figsize=(22,6))

ax1.plot(sealevel)
ax1.set_xlim(0, 100)
ax1.set_title('The signal')

ax2.plot(c_real, label='Real')
ax2.plot(c_imag, label='Imaginary')
ax2.set_xlim(0, 100)
ax2.set_title('The carrier signal')
ax2.legend()

n = [4, 16, 28]
for ax in [ax1, ax2]:
    for i in n:
        ax.axvline(i)

plt.show()

fig, axs = plt.subplots(1,3, figsize=(20,4))
axs[0].plot(c_real[:100], c_imag[:100])
axs[1].plot(c_real[:100], sealevel[:100])
axs[2].plot(c_imag[:100], sealevel[:100])

for ax in axs:
    ax.axis('equal')
plt.show()

plt.figure(figsize=(22,4))
plt.plot(product)

[<matplotlib.lines.Line2D at 0x2ba998cb810>]

2nd try! and solution...

Below is the solution presented in the Scipy site! I copy from, tweak, and voilá! Got what I want!

And, for those who want to impress the reviewer, you can talk about 'Hilbert transform', 'Fourier transform' and 'Hermitian symmetry'!!

I initially thought that I was doing 'Complex Demodulation'... especially to justify such fancy name in the article! However, I learned that this a different method to achieve the same result!

This will be used too!

# https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.hilbert.html

t = mdates.date2num(time)*24

x = sealevel

fs = 1

# kkkk... the follow 2 lines are all necessary to do the job!
analytic_signal = hilbert(x)
amplitude_envelope = np.abs(analytic_signal)

# unwrap? angle? Don't know... However, the next 2 steps are unnecessary!
instantaneous_phase = np.unwrap(np.angle(analytic_signal))
instantaneous_frequency = (np.diff(instantaneous_phase) / (2.0*np.pi* fs))

plt.figure(figsize=(22,5))
plt.plot(t, x)
plt.plot(t, lowpassfilter(amplitude_envelope), lw=3)

[<matplotlib.lines.Line2D at 0x2ba8c828990>]

plt.figure(figsize=(10,5))
plt.plot(t, x)
plt.axhline(0, color='k')
plt.plot(t, lowpassfilter(amplitude_envelope), lw=3)
plt.xlim(406000, 408000)

(406000.0, 408000.0)

analytic_signal[0]

(-1.587017385147893+1.6518101469127788j)

amplitude_envelope[0]

2.290655133843956

analytic_signal[0].real

-1.587017385147893

plt.figure(figsize=(22,4))
plt.plot(amplitude_envelope)

[<matplotlib.lines.Line2D at 0x2ba9807a690>]

plt.figure(figsize=(22,4))
plt.plot(instantaneous_phase)

[<matplotlib.lines.Line2D at 0x2ba980df890>]

plt.figure(figsize=(22,4))
plt.plot(instantaneous_frequency)
plt.ylim(0.04, 0.12)

(0.04, 0.12)

Function from Hawaii struff

This is the function for demodulation from the Hawaii stuff... if you want to use it, you need to install some package (from Hawaii) that have the 'Bunch' method! It seems a wrapper... I didn't use it.

# https://currents.soest.hawaii.edu/ocn_data_analysis/_static/complex_demod.html

def complex_demod(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 Bunch; look at the code to see what it contains.
    """     
    
    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))
    
    demod = bl_filt(product, hwpts)
    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.
    
    return Bunch(t=t,
                 signal=x,  
                 hwpts=hwpts,
                 demod=demod,
                 reconstructed=reconstructed)

Adapt from ChatGPT (below!)

import numpy as np
import matplotlib.pyplot as plt

# Create a modulated complex signal
fs = 1000  # Sampling frequency
t = mdates.date2num(time)*24  # Time vector
carrier_frequency = 1/(12 + 25/60)  # Frequency of the carrier signal
modulating_signal = sealevel  # Modulating signal (baseband)
carrier_signal = np.exp(2j * np.pi * carrier_frequency * t)  # Carrier signal

modulated_signal = modulating_signal * carrier_signal

# Perform complex demodulation
demodulated_signal = modulated_signal * np.conj(carrier_signal)

# Low-pass filter the demodulated signal
cutoff_frequency = 20  # Cutoff frequency of the low-pass filter
nyquist_rate = fs / 2
num_taps = 101  # Number of taps for the filter
h = np.sinc(2 * cutoff_frequency / nyquist_rate * (np.arange(num_taps) - (num_taps - 1) / 2))
demodulated_signal_filtered = np.convolve(demodulated_signal, h, mode='same')

# Plot the original signal, modulated signal, and demodulated signal
plt.figure(figsize=(12, 6))
plt.subplot(3, 1, 1)
plt.plot(t, modulating_signal, label='Modulating Signal (Baseband)')
plt.legend()

plt.subplot(3, 1, 2)
plt.plot(t, np.real(modulated_signal), label='Real Part of Modulated Signal')
plt.plot(t, np.imag(modulated_signal), label='Imaginary Part of Modulated Signal')
plt.legend()

plt.subplot(3, 1, 3)
plt.plot(t, np.real(demodulated_signal_filtered), label='Demodulated Signal (Filtered)')
plt.plot(t, np.imag(demodulated_signal_filtered), label='Demodulated Signal Imaginary Part (Filtered)')
plt.legend()

plt.tight_layout()
plt.show()

ChatGPT

import numpy as np
import matplotlib.pyplot as plt

# Create a modulated complex signal
fs = 1000  # Sampling frequency
t = np.linspace(0, 1, fs, endpoint=False)  # Time vector
carrier_frequency = 10  # Frequency of the carrier signal
modulating_signal = np.sin(2 * np.pi * 5 * t)  # Modulating signal (baseband)
carrier_signal = np.exp(2j * np.pi * carrier_frequency * t)  # Carrier signal

modulated_signal = modulating_signal * carrier_signal

# Perform complex demodulation
demodulated_signal = modulated_signal * np.conj(carrier_signal)

# Low-pass filter the demodulated signal
cutoff_frequency = 20  # Cutoff frequency of the low-pass filter
nyquist_rate = fs / 2
num_taps = 101  # Number of taps for the filter
h = np.sinc(2 * cutoff_frequency / nyquist_rate * (np.arange(num_taps) - (num_taps - 1) / 2))
demodulated_signal_filtered = np.convolve(demodulated_signal, h, mode='same')

# Plot the original signal, modulated signal, and demodulated signal
plt.figure(figsize=(12, 6))
plt.subplot(3, 1, 1)
plt.plot(t, modulating_signal, label='Modulating Signal (Baseband)')
plt.legend()

plt.subplot(3, 1, 2)
plt.plot(t, np.real(modulated_signal), label='Real Part of Modulated Signal')
plt.plot(t, np.imag(modulated_signal), label='Imaginary Part of Modulated Signal')
plt.legend()

plt.subplot(3, 1, 3)
plt.plot(t, np.real(demodulated_signal_filtered), label='Demodulated Signal (Filtered)')
plt.plot(t, np.imag(demodulated_signal_filtered), label='Demodulated Signal Imaginary Part (Filtered)')
plt.legend()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Create a synthetic tidal signal with varying amplitude and phase
t = np.linspace(0, 30, 1000)  # Time vector
frequency = 0.04  # Tidal frequency in cycles per time unit
amplitude = 1 + 0.5 * np.sin(2 * np.pi * 0.03 * t)  # Varying amplitude
phase = 0.2 * np.sin(2 * np.pi * 0.1 * t)  # Varying phase

# Create the modulated tidal signal
tidal_signal = amplitude * np.sin(2 * np.pi * frequency * t + phase)

# Perform complex demodulation
analytic_signal = tidal_signal - np.mean(tidal_signal)  # Remove the mean to get an analytic signal
hilbert_transform = np.imag(np.fft.ifft(np.fft.fft(analytic_signal) * 2j * np.pi * np.fft.fftfreq(len(t))))

# Calculate amplitude and phase
amplitude_envelope = np.abs(analytic_signal)
instantaneous_phase = np.angle(analytic_signal)

# Plot the original signal, amplitude envelope, and instantaneous phase
plt.figure(figsize=(12, 6))

plt.subplot(3, 1, 1)
plt.plot(t, tidal_signal, label='Tidal Signal')
plt.title('Original Tidal Signal')
plt.legend()

plt.subplot(3, 1, 2)
plt.plot(t, amplitude_envelope, label='Amplitude Envelope')
plt.title('Amplitude Envelope')
plt.legend()

plt.subplot(3, 1, 3)
plt.plot(t, instantaneous_phase, label='Instantaneous Phase')
plt.title('Instantaneous Phase')
plt.legend()

plt.tight_layout()
plt.show()