From 0da8c6a82228322cc3e935bdedcbef90e15bcb89 Mon Sep 17 00:00:00 2001 From: Eric Lybrand Date: Thu, 25 Jun 2020 15:59:26 -0700 Subject: [PATCH 1/9] Transpose wavelet transform output. Time is now plotted on the x-axis, and frequencies along the y-axis. --- neurodsp/timefrequency/wavelets.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/neurodsp/timefrequency/wavelets.py b/neurodsp/timefrequency/wavelets.py index 63892915..934c95dc 100644 --- a/neurodsp/timefrequency/wavelets.py +++ b/neurodsp/timefrequency/wavelets.py @@ -48,9 +48,9 @@ def compute_wavelet_transform(sig, fs, freqs, n_cycles=7, scaling=0.5): freqs = create_freqs(*freqs) n_cycles = check_n_cycles(n_cycles, len(freqs)) - mwt = np.zeros([len(sig), len(freqs)], dtype=complex) + mwt = np.zeros([len(freqs), len(sig)], dtype=complex) for ind, (freq, n_cycle) in enumerate(zip(freqs, n_cycles)): - mwt[:, ind] = convolve_wavelet(sig, fs, freq, n_cycle, scaling) + mwt[ind, :] = convolve_wavelet(sig, fs, freq, n_cycle, scaling) return mwt From 7296b0bbf54e0487a5355eef2321b297e2c86db2 Mon Sep 17 00:00:00 2001 From: Eric Lybrand Date: Thu, 25 Jun 2020 16:52:49 -0700 Subject: [PATCH 2/9] use np.unwrap to wrap negative phases in [0, 2pi] --- neurodsp/timefrequency/hilbert.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/neurodsp/timefrequency/hilbert.py b/neurodsp/timefrequency/hilbert.py index bc821479..cd0e756c 100644 --- a/neurodsp/timefrequency/hilbert.py +++ b/neurodsp/timefrequency/hilbert.py @@ -193,12 +193,10 @@ def freq_by_time(sig, fs, f_range=None, hilbert_increase_n=False, ... components={'sim_powerlaw': {}, 'sim_oscillation' : {'freq': 10}}) >>> instant_freq = freq_by_time(sig, fs=500, f_range=(8, 12)) """ - - pha = phase_by_time(sig, fs, f_range, hilbert_increase_n, - remove_edges, **filter_kwargs) + pha = np.unwrap(phase_by_time(sig, fs, f_range, hilbert_increase_n, + remove_edges, **filter_kwargs)) phadiff = np.diff(pha) - phadiff[phadiff < 0] = phadiff[phadiff < 0] + 2 * np.pi i_f = fs * phadiff / (2 * np.pi) i_f = np.insert(i_f, 0, np.nan) From 61f4b26d55345308659eede9d882e57793e3b2a9 Mon Sep 17 00:00:00 2001 From: Eric Lybrand Date: Thu, 25 Jun 2020 16:59:29 -0700 Subject: [PATCH 3/9] Average wavelet transform array entries within rows, since each row corresponds to a fixed frequency. --- neurodsp/spectral/power.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neurodsp/spectral/power.py b/neurodsp/spectral/power.py index 0f54f359..a5a88662 100644 --- a/neurodsp/spectral/power.py +++ b/neurodsp/spectral/power.py @@ -106,7 +106,7 @@ def compute_spectrum_wavelet(sig, fs, freqs, avg_type='mean', **kwargs): freqs = create_freqs(*freqs) mwt = compute_wavelet_transform(sig, fs, freqs, **kwargs) - spectrum = get_avg_func(avg_type)(mwt, axis=0) + spectrum = get_avg_func(avg_type)(mwt, axis=1) return freqs, spectrum From cabc436b52a13b2b0ed4cff8a82ce915ff8b1344 Mon Sep 17 00:00:00 2001 From: Eric Lybrand Date: Thu, 25 Jun 2020 17:03:53 -0700 Subject: [PATCH 4/9] Add norm kwarg for compute_wavelet_transform --- neurodsp/timefrequency/wavelets.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/neurodsp/timefrequency/wavelets.py b/neurodsp/timefrequency/wavelets.py index 934c95dc..2d73ff7c 100644 --- a/neurodsp/timefrequency/wavelets.py +++ b/neurodsp/timefrequency/wavelets.py @@ -11,7 +11,7 @@ ################################################################################################### @multidim() -def compute_wavelet_transform(sig, fs, freqs, n_cycles=7, scaling=0.5): +def compute_wavelet_transform(sig, fs, freqs, n_cycles=7, scaling=0.5, norm='sss'): """Compute the time-frequency representation of a signal using morlet wavelets. Parameters @@ -50,7 +50,7 @@ def compute_wavelet_transform(sig, fs, freqs, n_cycles=7, scaling=0.5): mwt = np.zeros([len(freqs), len(sig)], dtype=complex) for ind, (freq, n_cycle) in enumerate(zip(freqs, n_cycles)): - mwt[ind, :] = convolve_wavelet(sig, fs, freq, n_cycle, scaling) + mwt[ind, :] = convolve_wavelet(sig, fs, freq, n_cycle, scaling, norm=norm) return mwt From cee912fd31ba4fbc7f51916b98d0e96ae479ccac Mon Sep 17 00:00:00 2001 From: Eric Lybrand Date: Sun, 28 Jun 2020 20:11:17 -0700 Subject: [PATCH 5/9] Update documentation to reflect return of analytic signal, not hilbert transform. --- neurodsp/timefrequency/hilbert.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neurodsp/timefrequency/hilbert.py b/neurodsp/timefrequency/hilbert.py index cd0e756c..b3126af6 100644 --- a/neurodsp/timefrequency/hilbert.py +++ b/neurodsp/timefrequency/hilbert.py @@ -26,7 +26,7 @@ def robust_hilbert(sig, increase_n=False): Returns ------- sig_hilb : 1d array - The Hilbert transform of the input signal. + The analytic signal, of which the imaginary part is the Hilbert transform of the input signal. Examples -------- From 2443269c3245c51b6a7f9797353ea1d22de3b7e1 Mon Sep 17 00:00:00 2001 From: Eric Lybrand Date: Thu, 2 Jul 2020 13:51:57 -0700 Subject: [PATCH 6/9] Add docstring for norm kwarg. --- neurodsp/timefrequency/wavelets.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/neurodsp/timefrequency/wavelets.py b/neurodsp/timefrequency/wavelets.py index 2d73ff7c..2f14ee00 100644 --- a/neurodsp/timefrequency/wavelets.py +++ b/neurodsp/timefrequency/wavelets.py @@ -29,6 +29,12 @@ def compute_wavelet_transform(sig, fs, freqs, n_cycles=7, scaling=0.5, norm='sss scaling : float Scaling factor. + norm : {'sss', 'amp'}, optional + Normalization method: + + * 'sss' - divide by the square root of the sum of squares + * 'amp' - divide by the sum of amplitudes + Returns ------- mwt : 2d array From b48fc10aa6f73a5218353c2fe90a877da72ec8cc Mon Sep 17 00:00:00 2001 From: Eric Lybrand Date: Thu, 2 Jul 2020 16:11:20 -0700 Subject: [PATCH 7/9] Take modulus squared of morelet wavelet transform, then average. --- neurodsp/spectral/power.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neurodsp/spectral/power.py b/neurodsp/spectral/power.py index a5a88662..b77efe46 100644 --- a/neurodsp/spectral/power.py +++ b/neurodsp/spectral/power.py @@ -106,7 +106,7 @@ def compute_spectrum_wavelet(sig, fs, freqs, avg_type='mean', **kwargs): freqs = create_freqs(*freqs) mwt = compute_wavelet_transform(sig, fs, freqs, **kwargs) - spectrum = get_avg_func(avg_type)(mwt, axis=1) + spectrum = get_avg_func(avg_type)(abs(mwt)**2, axis=1) return freqs, spectrum From 2201fd8d8fede2edf2e31c891260a1dd6da1a9d5 Mon Sep 17 00:00:00 2001 From: Eric Lybrand Date: Fri, 3 Jul 2020 12:33:35 -0700 Subject: [PATCH 8/9] Update docstring of compute_wavelet_transform to reflect that you can feed in an array of n_cycles. --- neurodsp/timefrequency/wavelets.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/neurodsp/timefrequency/wavelets.py b/neurodsp/timefrequency/wavelets.py index 2f14ee00..0268bca4 100644 --- a/neurodsp/timefrequency/wavelets.py +++ b/neurodsp/timefrequency/wavelets.py @@ -24,8 +24,9 @@ def compute_wavelet_transform(sig, fs, freqs, n_cycles=7, scaling=0.5, norm='sss If array, frequency values to estimate with morlet wavelets. If list, define the frequency range, as [freq_start, freq_stop, freq_step]. The `freq_step` is optional, and defaults to 1. Range is inclusive of `freq_stop` value. - n_cycles : float + n_cycles : float, or 1d array. Length of the filter, as the number of cycles for each frequency. + If 1d array, this defines n_cycles for each frequency. scaling : float Scaling factor. From 12337110645f009d168b8798b9f214c430391b21 Mon Sep 17 00:00:00 2001 From: TomDonoghue Date: Mon, 6 Jul 2020 12:46:41 -0700 Subject: [PATCH 9/9] updates & makes i_f handle nans --- neurodsp/timefrequency/hilbert.py | 22 +++++++++++++++++----- neurodsp/timefrequency/wavelets.py | 3 +-- 2 files changed, 18 insertions(+), 7 deletions(-) diff --git a/neurodsp/timefrequency/hilbert.py b/neurodsp/timefrequency/hilbert.py index b3126af6..13c92dc3 100644 --- a/neurodsp/timefrequency/hilbert.py +++ b/neurodsp/timefrequency/hilbert.py @@ -30,7 +30,7 @@ def robust_hilbert(sig, increase_n=False): Examples -------- - Compute a Hilbert transform of a signal, using zero padding: + Compute the analytic signal, using zero padding: >>> from neurodsp.sim import sim_combined >>> sig = sim_combined(n_seconds=10, fs=500, @@ -193,12 +193,24 @@ def freq_by_time(sig, fs, f_range=None, hilbert_increase_n=False, ... components={'sim_powerlaw': {}, 'sim_oscillation' : {'freq': 10}}) >>> instant_freq = freq_by_time(sig, fs=500, f_range=(8, 12)) """ - pha = np.unwrap(phase_by_time(sig, fs, f_range, hilbert_increase_n, - remove_edges, **filter_kwargs)) - phadiff = np.diff(pha) + # Calculate instantaneous phase, from which we can compute instantaneous frequency + pha = phase_by_time(sig, fs, f_range, hilbert_increase_n, remove_edges, **filter_kwargs) - i_f = fs * phadiff / (2 * np.pi) + # 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 diff --git a/neurodsp/timefrequency/wavelets.py b/neurodsp/timefrequency/wavelets.py index 0268bca4..3cd76749 100644 --- a/neurodsp/timefrequency/wavelets.py +++ b/neurodsp/timefrequency/wavelets.py @@ -24,12 +24,11 @@ def compute_wavelet_transform(sig, fs, freqs, n_cycles=7, scaling=0.5, norm='sss If array, frequency values to estimate with morlet wavelets. If list, define the frequency range, as [freq_start, freq_stop, freq_step]. The `freq_step` is optional, and defaults to 1. Range is inclusive of `freq_stop` value. - n_cycles : float, or 1d array. + n_cycles : float or 1d array Length of the filter, as the number of cycles for each frequency. If 1d array, this defines n_cycles for each frequency. scaling : float Scaling factor. - norm : {'sss', 'amp'}, optional Normalization method: