diff --git a/fooof/synth.py b/fooof/synth.py index 15a87b681..ea5b38b11 100644 --- a/fooof/synth.py +++ b/fooof/synth.py @@ -375,6 +375,51 @@ def gen_power_vals(xs, background_params, gauss_params, nlv): return ys +def rotate_spectrum(freqs, power_spectrum, delta_f, f_rotation): + """Rotate a power spectrum about a frequency point, changing the power law exponent. + + Parameters + ---------- + 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 + Frequency, in Hz, at which to do the rotation, such that power at that frequency is unchanged. + + Returns + ------- + 1d array + Rotated psd. + + """ + + # 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. + 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)) + + f_mask = f_mask / f_mask[np.where(freqs >= f_rotation)[0][0]] + + return f_mask * power_spectrum + + ################################################################################################### ################################################################################################### diff --git a/fooof/tests/test_synth.py b/fooof/tests/test_synth.py index 00553f040..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,6 +138,20 @@ def test_gen_power_values(): assert np.all(ys) +def test_rotate_spectrum(): + + # Create a spectrum to use for test rotations + freqs, spectrum = gen_power_spectrum([1, 100], [1, 1], []) + + # 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) + + # 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) + + def test_check_iter(): # Note: generator case not tested