From 6848319fa08888bc480f82ce2d150cc90532aec1 Mon Sep 17 00:00:00 2001 From: rdgao Date: Tue, 23 Oct 2018 16:13:22 -0700 Subject: [PATCH 1/4] added function to perform rotation about axis frequency. --- fooof/synth.py | 44 +++++++++++++++++++++++++++++++++++++++ fooof/tests/test_synth.py | 10 +++++++++ 2 files changed, 54 insertions(+) diff --git a/fooof/synth.py b/fooof/synth.py index 15a87b681..d51fd42c4 100644 --- a/fooof/synth.py +++ b/fooof/synth.py @@ -375,6 +375,50 @@ def gen_power_vals(xs, background_params, gauss_params, nlv): return ys +def rotate_powerlaw(psd, freqs, delta_f, f_rotation=30): + """Change the power law exponent of a PSD about an axis frequency. + + Parameters + ---------- + psd : 1d array + Power spectrum to be rotated. + freqs : 1d array, Hz + Frequency axis of input PSD. Must be same length as psd. + delta_f : float + Change in power law exponent to be applied. Positive is counterclockwise + rotation (flatten), negative is clockwise rotation (steepen). + f_rotation : float, Hz, optional + Axis of rotation frequency, such that power at that frequency is unchanged + by the rotation. Only matters if not further normalizing signal variance. + + Returns + ------- + x : 1d array + Rotated psd. + + """ + if np.all(f_rotationnp.abs(freqs).max()): + raise ValueError('Rotation frequency not within frequency range.') + + # make the 1/f rotation mask + f_mask = np.zeros_like(freqs) + if freqs[0] == 0.: + # If starting freq is 0Hz, default power at 0Hz to old value because log + # will return inf. Use case occurs in simulating/manipulating time series. + f_mask[0] = 1. + f_mask[1:] = 10**(np.log10(np.abs(freqs[1:])) * (delta_f)) + else: + # Otherwise, apply rotation to all frequencies. + f_mask = 10**(np.log10(np.abs(freqs)) * (delta_f)) + + # normalize power at rotation frequency + f_mask = f_mask / f_mask[np.where(freqs >= f_rotation)[0][0]] + + # apply mask + return psd * f_mask + + + ################################################################################################### ################################################################################################### diff --git a/fooof/tests/test_synth.py b/fooof/tests/test_synth.py index 00553f040..5506ecfe9 100644 --- a/fooof/tests/test_synth.py +++ b/fooof/tests/test_synth.py @@ -136,6 +136,16 @@ def test_gen_power_values(): assert np.all(ys) +def test_rotate_powerlaw(): + # Not the best test right now, just checking that a change in slope of 0 + # does nothing to a flat spectrum. + freqs = np.arange(500.) + + psd_sim = np.ones_like(freqs) + psd_rot = rotate_powerlaw(psd_sim, freqs, delta_f=0.) + + assert np.all(psd_rot==psd_sim) + def test_check_iter(): # Note: generator case not tested From cfcbe89053765f95dbf18530366b7185417c6020 Mon Sep 17 00:00:00 2001 From: rdgao Date: Wed, 24 Oct 2018 11:06:28 -0700 Subject: [PATCH 2/4] removed default at 30Hz for rotate_powerlaw --- fooof/synth.py | 12 +++++++----- fooof/tests/test_synth.py | 4 ++++ 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/fooof/synth.py b/fooof/synth.py index d51fd42c4..b17655e3e 100644 --- a/fooof/synth.py +++ b/fooof/synth.py @@ -375,7 +375,7 @@ def gen_power_vals(xs, background_params, gauss_params, nlv): return ys -def rotate_powerlaw(psd, freqs, delta_f, f_rotation=30): +def rotate_powerlaw(psd, freqs, delta_f, f_rotation=None): """Change the power law exponent of a PSD about an axis frequency. Parameters @@ -390,6 +390,7 @@ def rotate_powerlaw(psd, freqs, delta_f, f_rotation=30): f_rotation : float, Hz, optional Axis of rotation frequency, such that power at that frequency is unchanged by the rotation. Only matters if not further normalizing signal variance. + If None, the transform normalizes to power at 1Hz by defaults. Returns ------- @@ -397,9 +398,6 @@ def rotate_powerlaw(psd, freqs, delta_f, f_rotation=30): Rotated psd. """ - if np.all(f_rotationnp.abs(freqs).max()): - raise ValueError('Rotation frequency not within frequency range.') - # make the 1/f rotation mask f_mask = np.zeros_like(freqs) if freqs[0] == 0.: @@ -412,7 +410,11 @@ def rotate_powerlaw(psd, freqs, delta_f, f_rotation=30): f_mask = 10**(np.log10(np.abs(freqs)) * (delta_f)) # normalize power at rotation frequency - f_mask = f_mask / f_mask[np.where(freqs >= f_rotation)[0][0]] + if f_rotation is not None: + if np.all(f_rotationnp.abs(freqs).max()): + raise ValueError('Rotation frequency not within frequency range.') + + f_mask = f_mask / f_mask[np.where(freqs >= f_rotation)[0][0]] # apply mask return psd * f_mask diff --git a/fooof/tests/test_synth.py b/fooof/tests/test_synth.py index 5506ecfe9..dbd4c04d3 100644 --- a/fooof/tests/test_synth.py +++ b/fooof/tests/test_synth.py @@ -146,6 +146,10 @@ def test_rotate_powerlaw(): assert np.all(psd_rot==psd_sim) + psd_rot = rotate_powerlaw(psd_sim, freqs, delta_f=0., f_rotation=30.) + + assert np.all(psd_rot==psd_sim) + def test_check_iter(): # Note: generator case not tested From a920117bb7d11fe738099f10d2d9596c24335c88 Mon Sep 17 00:00:00 2001 From: TomDonoghue Date: Wed, 24 Oct 2018 17:49:49 -0700 Subject: [PATCH 3/4] Updates to code org & strip back scope, update docs --- fooof/synth.py | 47 +++++++++++++++++++++++------------------------ 1 file changed, 23 insertions(+), 24 deletions(-) diff --git a/fooof/synth.py b/fooof/synth.py index b17655e3e..ea5b38b11 100644 --- a/fooof/synth.py +++ b/fooof/synth.py @@ -375,50 +375,49 @@ def gen_power_vals(xs, background_params, gauss_params, nlv): return ys -def rotate_powerlaw(psd, freqs, delta_f, f_rotation=None): - """Change the power law exponent of a PSD about an axis frequency. +def rotate_spectrum(freqs, power_spectrum, delta_f, f_rotation): + """Rotate a power spectrum about a frequency point, changing the power law exponent. Parameters ---------- - psd : 1d array - Power spectrum to be rotated. - freqs : 1d array, Hz - Frequency axis of input PSD. Must be same length as psd. + freqs : 1d array + Frequency axis of input power spectrum, in Hz. + power_spectrum : 1d array + Power values of the spectrum that is to be rotated. delta_f : float - Change in power law exponent to be applied. Positive is counterclockwise - rotation (flatten), negative is clockwise rotation (steepen). - f_rotation : float, Hz, optional - Axis of rotation frequency, such that power at that frequency is unchanged - by the rotation. Only matters if not further normalizing signal variance. - If None, the transform normalizes to power at 1Hz by defaults. + Change in power law exponent to be applied. + Positive is counterclockwise rotation (flatten) + Negative is clockwise rotation (steepen). + f_rotation : float + Frequency, in Hz, at which to do the rotation, such that power at that frequency is unchanged. Returns ------- - x : 1d array + 1d array Rotated psd. """ - # make the 1/f rotation mask + + # Check that the requested frequency rotation value is within the given range + if f_rotation < freqs.min() or f_rotation > freqs.max(): + raise ValueError('Rotation frequency not within frequency range.') + f_mask = np.zeros_like(freqs) + + + # NOTE: Update this to use data checking from fooof.fit if freqs[0] == 0.: # If starting freq is 0Hz, default power at 0Hz to old value because log - # will return inf. Use case occurs in simulating/manipulating time series. + # will return inf. Use case occurs in simulating/manipulating time series. f_mask[0] = 1. f_mask[1:] = 10**(np.log10(np.abs(freqs[1:])) * (delta_f)) else: # Otherwise, apply rotation to all frequencies. f_mask = 10**(np.log10(np.abs(freqs)) * (delta_f)) - # normalize power at rotation frequency - if f_rotation is not None: - if np.all(f_rotationnp.abs(freqs).max()): - raise ValueError('Rotation frequency not within frequency range.') - - f_mask = f_mask / f_mask[np.where(freqs >= f_rotation)[0][0]] - - # apply mask - return psd * f_mask + f_mask = f_mask / f_mask[np.where(freqs >= f_rotation)[0][0]] + return f_mask * power_spectrum ################################################################################################### From 5fd0a4905fcebbaadf6a635bcb9e05ce4e8125a6 Mon Sep 17 00:00:00 2001 From: TomDonoghue Date: Wed, 24 Oct 2018 19:05:55 -0700 Subject: [PATCH 4/4] Update to fixed smoke tests --- fooof/tests/test_synth.py | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/fooof/tests/test_synth.py b/fooof/tests/test_synth.py index dbd4c04d3..03ef03a93 100644 --- a/fooof/tests/test_synth.py +++ b/fooof/tests/test_synth.py @@ -55,6 +55,8 @@ def test_stepper(): assert Stepper(8,12,.1) + # TODO: add more tests of Stepper + def test_param_sampler(): @@ -136,19 +138,19 @@ def test_gen_power_values(): assert np.all(ys) -def test_rotate_powerlaw(): - # Not the best test right now, just checking that a change in slope of 0 - # does nothing to a flat spectrum. - freqs = np.arange(500.) +def test_rotate_spectrum(): - psd_sim = np.ones_like(freqs) - psd_rot = rotate_powerlaw(psd_sim, freqs, delta_f=0.) + # Create a spectrum to use for test rotations + freqs, spectrum = gen_power_spectrum([1, 100], [1, 1], []) - assert np.all(psd_rot==psd_sim) + # Check that rotation transforms the power spectrum + rotated_spectrum = rotate_spectrum(freqs, spectrum, delta_f=0.5, f_rotation=25.) + assert not np.all(rotated_spectrum == spectrum) - psd_rot = rotate_powerlaw(psd_sim, freqs, delta_f=0., f_rotation=30.) + # Check that 0 rotation returns the same spectrum + rotated_spectrum = rotate_spectrum(freqs, spectrum, delta_f=0., f_rotation=25.) + assert np.all(rotated_spectrum == spectrum) - assert np.all(psd_rot==psd_sim) def test_check_iter():