Add 'compute_rotation_frequency'

doc/api.rst

@@ -181,8 +181,9 @@ Functions for transforming power spectra.
     rotate_spectrum
     rotate_sim_spectrum
     compute_rotation_offset
+    compute_rotation_frequency
 
-Utitilies
+Utilities
 ~~~~~~~~~
 
 Utilities for simulating power spectra.

fooof/sim/transform.py

@@ -182,3 +182,47 @@ def compute_rotation_offset(delta_exponent, f_rotation):
     """
 
     return -np.log10(f_rotation) * -delta_exponent
+
+
+def compute_rotation_frequency(delta_exponent_b, f_rotation_b, delta_exponent_c, f_rotation_c):
+    """Calculate the rotation frequency between two rotated power spectra.
+
+    Parameters
+    ----------
+    delta_exponent_b : float
+        The applied change in exponent value for power spectrum 'B'.
+    f_rotation_b : float
+        The rotation frequency applied to power spectrum 'B'.
+    delta_exponent_c : float
+        The applied change in exponent value for power spectrum 'C'.
+    f_rotation_c : float
+        The rotation frequency applied to power spectrum 'C'.
+
+    Returns
+    -------
+    float
+        The frequency rotation point between spectra 'B' & 'C'.
+
+    Notes
+    -----
+    **Code Notes**
+
+    This computes the rotation frequency for two power spectra 'B' & 'C',
+    under the assumption that they are both rotated versions of a the
+    same original power spectrum 'A'.
+
+    **Derivation**
+
+    Given an original power spectrum A, then:
+
+    - B = A*(f_rotation_b/freqs)^delta_exponent_b
+    - C = A*(f_rotation_c/freqs)^delta_exponent_c
+
+    Therefore, what you want is f_rotation_bc, which is the frequency where B==C.
+
+    To find this, we can plug everything back into the equation, to find where
+    B[freqs] == C[freqs], which is how we arrive at the solution below.
+    """
+
+    return (((f_rotation_c**delta_exponent_c) / (f_rotation_b**delta_exponent_b))) ** \
+         (1/(delta_exponent_c-delta_exponent_b))

fooof/tests/sim/test_transform.py

@@ -58,3 +58,13 @@ def test_translate_sim_spectrum():
 def test_compute_rotation_offset():
 
     assert compute_rotation_offset(20, 0.5)
+
+def test_compute_rotation_frequency():
+
+    delta_exp_b, delta_exp_c = 0.5, 0.75
+    f_rot_b, f_rot_c = 5, 10
+
+    f_rot_bc = compute_rotation_frequency(delta_exp_b, f_rot_b, delta_exp_c, f_rot_c)
+
+    assert isinstance(f_rot_bc, float)
+    assert np.isclose(f_rot_bc, 40.)