Ошибка вычисления фазового угла между двумя временными рядами с использованием преобразования Гильберта ⇐ Python

[Anonymous]

Я пытаюсь вычислить фазовый угол между двумя временными рядами действительных чисел. Чтобы проверить, работает ли моя функция без ошибок, я создал две синусоидальные волны с фазой 17 градусов. Однако, когда я вычисляю фазовый угол между этими двумя синусоидами, я не получаю 17 градусов. Вот мой сценарий:

Код: Выделить всё

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()

Заранее благодарим за любую помощь.

Подробнее здесь: https://stackoverflow.com/questions/791 ... -transform