Source code for neurodsp.timefrequency.hilbert

"""Tools for estimating time frequency properties, using Hilbert methods."""

import numpy as np
from scipy.signal import hilbert

from neurodsp.filt import filter_signal
from neurodsp.utils.decorators import multidim
from neurodsp.utils.outliers import remove_nans, restore_nans
from neurodsp.filt.utils import infer_passtype, remove_filter_edges

###################################################################################################
###################################################################################################

[docs]@multidim() def robust_hilbert(sig): """Compute the Hilbert transform, ignoring any boundaries that are NaN. Parameters ---------- sig : 1d array Time series. Returns ------- sig_hilb : 1d array The analytic signal, of which the imaginary part is the Hilbert transform of the input. Examples -------- Compute the analytic signal, using zero padding: >>> from neurodsp.sim import sim_combined >>> sig = sim_combined(n_seconds=10, fs=500, ... components={'sim_powerlaw': {}, 'sim_oscillation': {'freq': 10}}) >>> sig_hilb = robust_hilbert(sig) """ # Extract the signal that is not nan sig_nonan, sig_nans = remove_nans(sig) # Compute Hilbert transform of signal without nans sig_hilb_nonan = hilbert(sig_nonan) # Fill in output hilbert with nans on edges sig_hilb = restore_nans(sig_hilb_nonan, sig_nans, dtype=complex) return sig_hilb
[docs]@multidim() def phase_by_time(sig, fs, f_range=None, remove_edges=True, **filter_kwargs): """Compute the instantaneous phase of a time series. Parameters ---------- sig : 1d array Time series. fs : float Sampling rate, in Hz. f_range : tuple of float or None, optional default: None Filter range, in Hz, as (low, high). If None, no filtering is applied. remove_edges : bool, optional, default: True If True, replace samples that are within half of the filter's length to the edge with nan. This removes edge artifacts from the filtered signal. Only used if `f_range` is defined. **filter_kwargs Keyword parameters to pass to `filter_signal`. Returns ------- pha : 1d array Instantaneous phase time series. Examples -------- Compute the instantaneous phase, for the alpha range: >>> from neurodsp.sim import sim_combined >>> sig = sim_combined(n_seconds=10, fs=500, ... components={'sim_powerlaw': {}, 'sim_oscillation': {'freq': 10}}) >>> pha = phase_by_time(sig, fs=500, f_range=(8, 12)) """ if f_range: sig, filter_kernel = filter_signal(sig, fs, infer_passtype(f_range), f_range=f_range, remove_edges=False, return_filter=True, **filter_kwargs) pha = np.angle(robust_hilbert(sig)) if f_range and remove_edges: pha = remove_filter_edges(pha, len(filter_kernel)) return pha
[docs]@multidim() def amp_by_time(sig, fs, f_range=None, remove_edges=True, **filter_kwargs): """Compute the instantaneous amplitude of a time series. Parameters ---------- sig : 1d array Time series. fs : float Sampling rate, in Hz. f_range : tuple of float or None, optional default: None Filter range, in Hz, as (low, high). If None, no filtering is applied. remove_edges : bool, optional, default: True If True, replace samples that are within half of the filter's length to the edge with nan. This removes edge artifacts from the filtered signal. Only used if `f_range` is defined. **filter_kwargs Keyword parameters to pass to `filter_signal`. Returns ------- amp : 1d array Instantaneous amplitude time series. Examples -------- Compute the instantaneous amplitude, for the alpha range: >>> from neurodsp.sim import sim_combined >>> sig = sim_combined(n_seconds=10, fs=500, ... components={'sim_powerlaw': {}, 'sim_oscillation' : {'freq': 10}}) >>> amp = amp_by_time(sig, fs=500, f_range=(8, 12)) """ if f_range: sig, filter_kernel = filter_signal(sig, fs, infer_passtype(f_range), f_range=f_range, remove_edges=False, return_filter=True, **filter_kwargs) amp = np.abs(robust_hilbert(sig)) if f_range and remove_edges: amp = remove_filter_edges(amp, len(filter_kernel)) return amp
[docs]@multidim() def freq_by_time(sig, fs, f_range=None, remove_edges=True, **filter_kwargs): """Compute the instantaneous frequency of a time series. Parameters ---------- sig : 1d array Time series. fs : float Sampling rate, in Hz. f_range : tuple of float or None, optional default: None Filter range, in Hz, as (low, high). If None, no filtering is applied. remove_edges : bool, optional, default: True If True, replace samples that are within half of the filter's length to the edge with nan. This removes edge artifacts from the filtered signal. Only used if `f_range` is defined. **filter_kwargs Keyword parameters to pass to `filter_signal`. Returns ------- i_f : 1d array Instantaneous frequency time series. Notes ----- This function assumes monotonic phase, so phase slips will be processed as high frequencies. Examples -------- Compute the instantaneous frequency, for the alpha range: >>> from neurodsp.sim import sim_combined >>> sig = sim_combined(n_seconds=10, fs=500, ... components={'sim_powerlaw': {}, 'sim_oscillation' : {'freq': 10}}) >>> instant_freq = freq_by_time(sig, fs=500, f_range=(8, 12)) """ # Calculate instantaneous phase, from which we can compute instantaneous frequency pha = phase_by_time(sig, fs, f_range, remove_edges, **filter_kwargs) # If filter edges were removed, temporarily drop nans for subsequent operations pha, nans = remove_nans(pha) # Differentiate the unwrapped phase, to estimate frequency as the rate of change of phase pha_unwrapped = np.unwrap(pha) pha_diff = np.diff(pha_unwrapped) # Convert differentiated phase to frequency i_f = fs * pha_diff / (2 * np.pi) # Prepend nan value to re-align with original signal (necessary due to differentiation) i_f = np.insert(i_f, 0, np.nan) # Restore nans, to re-align signal if filter edges were removed i_f = restore_nans(i_f, nans) return i_f