I'm trying to compute the phase angle between two time-series of real numbers. To check if my function is working without errors I have created two sine waves with a phase of 17 degrees. However, when I compute the phase angle between those two sine waves I do not get the 17 degrees. Here's my script:
import numpy as np
from scipy.signal import hilbert
import matplotlib.pyplot as plt
def coupling_angle_hilbert(x, y, datatype, center=True, pad=True):
"""
Compute the phase angle between two time series using the Hilbert transform.
Parameters:
- x: numpy array
Time series data for the first signal.
- y: numpy array
Time series data for the second signal.
- center: bool, optional
If True, center the amplitude of the data around zero. Default is True.
- pad: bool, optional
If True, perform data reflection to address issues arising with data distortion. Default is True.
- unwrap: bool, optional
If True, unwrap the phase angle to avoid phase wrapping. Default is True.
Returns:
- phase_angle: numpy array
Phase angle between the two signals.
"""
# Convert input data to radians if specified as degrees
if datatype.lower().strip() == 'degs':
x = np.radians(x)
y = np.radians(y)
# Center the signals if the 'center' option is enabled
if center:
# Adjust x to be centered around zero: subtract minimum, then offset by half the range
x = x - np.min(x) - ((np.max(x) - np.min(x))/2)
# Adjust y to be centered around zero: subtract minimum, then offset by half the range
y = y - np.min(y) - ((np.max(y) - np.min(y))/2)
# Reflect and pad the data if padding is enabled
if pad:
# Number of padding samples equal to signal length
# Ensure that the number of pads is even
npads = x.shape[0] // 2 * 2 # Ensure npads is even
# Reflect data at the beginning and end to create padding for 'x' and 'y'
x_padded = np.concatenate((x[:npads][::-1], x, x[-npads:][::-1]))
y_padded = np.concatenate((y[:npads][::-1], y, y[-npads:][::-1]))
else:
# If padding not enabled, use original signals without modification
x_padded = x
y_padded = y
# Apply the Hilbert transform to the time series data
hilbert_x = hilbert(x_padded)
hilbert_y = hilbert(y_padded)
# Calculate the phase of each signal by using arctan2 on imaginary and real parts
phase_angle_x = np.arctan2(hilbert_x.imag, x_padded)
phase_angle_y = np.arctan2(hilbert_y.imag, y_padded)
# Calculate the phase difference between y and x
phase_angle = phase_angle_y - phase_angle_x
# Trim the phase_angle to match the shape of x or y
if pad:
# Remove initial and ending padding to return only the original signal's phase angle difference
phase_angle = phase_angle[npads:npads + x.shape[0]]
return phase_angle
# input data
angles = np.radians(np.arange(0, 360, 1))
phase_offset = np.radians(17)
wav1 = np.sin(angles)
wav2 = np.sin(angles + phase_offset)
# Compute phase_angle usig Hilbert transform
ca_hilbert = coupling_angle_hilbert(wav1,
wav2,
'rads',
center=True,
pad=True)
plt.plot(np.degrees(ca_hilbert))
plt.show()
Thank you n advance for any help.
Instead of using arctan2(hilbert.imag, x), you can use np.angle(), which returns the argument (angle) of a complex number (always in the range [-π, π]). It's essentially doing arctan2(y, x) for a complex number x + iy.
Also, after phase_angle = phase_y - phase_x, we need to ensure again that it lies in [-π, π], so we do phase_angle = np.angle(np.exp(1j * phase_angle)), following the documentation.
Thus, your function becomes:
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import hilbert
def coupling_angle_hilbert(x, y, datatype="rads", center=True, pad=True):
"""
Compute the phase angle between two time series using the Hilbert transform.
Parameters:
- x: numpy array
Time series data for the first signal.
- y: numpy array
Time series data for the second signal.
- datatype: str, optional
Specify if input is in 'rads' or 'degs'. Default is 'rads'.
- center: bool, optional
If True, center the amplitude of the data around zero. Default is True.
- pad: bool, optional
If True, perform data reflection to address issues arising with data distortion. Default is True.
Returns:
- phase_angle: numpy array
Phase angle between the two signals in radians.
"""
if datatype.lower().strip() == "degs":
x = np.radians(x)
y = np.radians(y)
# Center the signals if the 'center' option is enabled
if center:
# Adjust x to be centered around zero: subtract minimum, then offset by half the range
x = x - np.min(x) - ((np.max(x) - np.min(x)) / 2)
# Adjust y to be centered around zero: subtract minimum, then offset by half the range
y = y - np.min(y) - ((np.max(y) - np.min(y)) / 2)
# Reflect and pad the data if padding is enabled
if pad:
# Number of padding samples equal to signal length
# Ensure that the number of pads is even
npads = x.shape[0] // 2 * 2 # Ensure npads is even
# Reflect data at the beginning and end to create padding for 'x' and 'y'
x_padded = np.concatenate((x[:npads][::-1], x, x[-npads:][::-1]))
y_padded = np.concatenate((y[:npads][::-1], y, y[-npads:][::-1]))
else:
# If padding not enabled, use original signals without modification
x_padded = x
y_padded = y
# Apply the Hilbert transform to the time series data
hilbert_x = hilbert(x_padded)
hilbert_y = hilbert(y_padded)
# Calculate the instantaneous phases
phase_x = np.angle(hilbert_x)
phase_y = np.angle(hilbert_y)
# Calculate the phase difference
phase_angle = phase_y - phase_x
# Ensure phase angle is in [-π, π]
phase_angle = np.angle(np.exp(1j * phase_angle))
# Trim the phase_angle to match the shape of x or y
if pad:
# Remove initial and ending padding to return only the original signal's phase angle difference
phase_angle = phase_angle[npads : npads + x.shape[0]]
return phase_angle
With that, I get 17.41 as phase difference.