[ENH] Aperiodic Knee Fit

fooof/core/funcs.py

@@ -51,8 +51,8 @@ def expo_function(xs, *params):
     xs : 1d array
         Input x-axis values.
     *params : float
-        Parameters (offset, knee, exp) that define Lorentzian function:
-        y = 10^offset * (1/(knee + x^exp))
+        Parameters (offset, knee_freq, exp) that define Lorentzian function:
+        y = 10^offset * (1/(knee_freq^exp + x^exp))
 
     Returns
     -------
@@ -62,9 +62,56 @@ def expo_function(xs, *params):
 
     ys = np.zeros_like(xs)
 
-    offset, knee, exp = params
+    offset, knee_freq, exp = params
 
-    ys = ys + offset - np.log10(knee + xs**exp)
+    ys = ys + offset - np.log10(knee_freq**exp + xs**exp)
+
+    return ys
+
+
+def expo_const_function(xs, *params):
+    """Exponential fitting function, for fitting aperiodic component with a 'knee' and constant.
+
+    NOTE: this function requires linear frequency (not log).
+
+    Parameters
+    ----------
+    xs : 1d array
+        Input x-axis values.
+    *params : float
+        Parameters (offset, knee_freq, exp, const) that define Lorentzian function:
+        y = 10^offset * (1/(knee_freq^exp + x^exp)) + const
+
+    Returns
+    -------
+    ys : 1d array
+        Output values for exponential function.
+    """
+
+    ys = np.zeros_like(xs)
+
+    offset, knee_freq, exp, const = params
+
+    f_e = xs**exp
+    fk_e = knee_freq**exp
+
+    ys = ys + np.log10((10**offset + (const * fk_e) + (const * f_e)) / (fk_e + f_e))
+
+    return ys
+
+
+def expo_double_const_function(xs, *params):
+    """Double exponential fitting function, for fitting aperiodic component."""
+
+    offset0, fk0, exp0, const0, offset1, fk1, exp1, const1 = params
+
+    num_a = 10**offset0 + (const0 * fk0**exp0) + (const0 * xs**exp0)
+    num_b = 10**offset1 + (const1 * fk1**exp1) + (const1 * xs**exp1)
+
+    denom_a = fk0**exp0 + xs**exp0
+    denom_b =  fk1**exp1 + xs**exp1
+
+    ys =  np.log10((num_a/denom_a) + (num_b/denom_b))
 
     return ys
 
@@ -180,7 +227,7 @@ def get_ap_func(aperiodic_mode):
 
     Parameters
     ----------
-    aperiodic_mode : {'fixed', 'knee'}
+    aperiodic_mode : {'fixed', 'knee', 'knee_constant'}
         Which aperiodic fitting function to return.
 
     Returns
@@ -198,6 +245,8 @@ def get_ap_func(aperiodic_mode):
         ap_func = expo_nk_function
     elif aperiodic_mode == 'knee':
         ap_func = expo_function
+    elif aperiodic_mode == 'knee_constant':
+        ap_func = expo_const_function
     else:
         raise ValueError("Requested aperiodic mode not understood.")
 
@@ -214,7 +263,7 @@ def infer_ap_func(aperiodic_params):
 
     Returns
     -------
-    aperiodic_mode : {'fixed', 'knee'}
+    aperiodic_mode : {'fixed', 'knee', 'knee_constant'}
         Which kind of aperiodic fitting function the given parameters are consistent with.
 
     Raises
@@ -227,6 +276,8 @@ def infer_ap_func(aperiodic_params):
         aperiodic_mode = 'fixed'
     elif len(aperiodic_params) == 3:
         aperiodic_mode = 'knee'
+    elif len(aperiodic_params) == 4:
+        aperiodic_mode = 'knee_constant'
     else:
         raise InconsistentDataError("The given aperiodic parameters are "
                                     "inconsistent with available options.")

fooof/core/info.py

@@ -64,7 +64,7 @@ def get_ap_indices(aperiodic_mode):
 
     Parameters
     ----------
-    aperiodic_mode : {'fixed', 'knee'}
+    aperiodic_mode : {'fixed', 'knee', 'knee_constant'}
         Which mode was used for the aperiodic component.
 
     Returns
@@ -77,6 +77,8 @@ def get_ap_indices(aperiodic_mode):
         labels = ('offset', 'exponent')
     elif aperiodic_mode == 'knee':
         labels = ('offset', 'knee', 'exponent')
+    elif aperiodic_mode == 'knee_constant':
+        labels = ('offset', 'knee', 'exponent', 'constant')
     else:
         raise ValueError("Aperiodic mode not understood.")
 
@@ -90,7 +92,7 @@ def get_indices(aperiodic_mode):
 
     Parameters
     ----------
-    aperiodic_mode : {'fixed', 'knee'}
+    aperiodic_mode : {'fixed', 'knee', 'knee_constant'}
         Which mode was used for the aperiodic component.
 
     Returns

fooof/core/strings.py

@@ -282,6 +282,13 @@ def gen_results_fm_str(fm, concise=False):
         return _no_model_str(concise)
 
     # Create the formatted strings for printing
+    ap_params = ['offset', 'exponent']
+
+    if fm.aperiodic_mode != 'fixed':
+        ap_params.insert(1, 'knee')
+    if fm.aperiodic_mode == 'knee_constant':
+        ap_params.insert(3, 'constant')
+
     str_lst = [
 
         # Header
@@ -297,8 +304,7 @@ def gen_results_fm_str(fm, concise=False):
         '',
 
         # Aperiodic parameters
-        ('Aperiodic Parameters (offset, ' + ('knee, ' if fm.aperiodic_mode == 'knee' else '') + \
-         'exponent): '),
+        ('Aperiodic Parameters (' + ', '.join(ap_params) + '): '),
         ', '.join(['{:2.4f}'] * len(fm.aperiodic_params_)).format(*fm.aperiodic_params_),
         '',
 
@@ -355,6 +361,8 @@ def gen_results_fg_str(fg, concise=False):
     exps = fg.get_params('aperiodic_params', 'exponent')
     kns = fg.get_params('aperiodic_params', 'knee') \
         if fg.aperiodic_mode == 'knee' else np.array([0])
+    consts = fg.get_params('aperiodic_params', 'constant') \
+        if fg.aperiodic_mode == 'knee_constant' else np.array([0])
 
     # Check if there are any power spectra that failed to fit
     n_failed = sum(np.isnan(exps))
@@ -379,15 +387,19 @@ def gen_results_fg_str(fg, concise=False):
         '',
 
         # Aperiodic parameters - knee fit status, and quick exponent description
-        'Power spectra were fit {} a knee.'.format(\
-            'with' if fg.aperiodic_mode == 'knee' else 'without'),
+        'Power spectra were fit {w} a knee{e}'.format(
+            w='with' if fg.aperiodic_mode != 'fixed' else 'without',
+            e=' and constant.' if fg.aperiodic_mode == 'knee_constant' else '.'),
         '',
         'Aperiodic Fit Values:',
-        *[el for el in ['    Knees - Min: {:6.2f}, Max: {:6.2f}, Mean: {:5.2f}'
+        *[el for el in ['    Knees - Min: {:8.2f}, Max: {:8.2f}, Mean: {:8.2f}'
                         .format(np.nanmin(kns), np.nanmax(kns), np.nanmean(kns)),
-                       ] if fg.aperiodic_mode == 'knee'],
-        'Exponents - Min: {:6.3f}, Max: {:6.3f}, Mean: {:5.3f}'
+                       ] if fg.aperiodic_mode != 'fixed'],
+        'Exponents - Min: {:8.3f}, Max: {:8.3f}, Mean: {:8.3f}'
         .format(np.nanmin(exps), np.nanmax(exps), np.nanmean(exps)),
+        *[el for el in [' Constants - Min: {:8.1e}, Max: {:8.1e}, Mean: {:8.1e}'
+                        .format(np.nanmin(consts), np.nanmax(consts), np.nanmean(consts)),
+                       ] if fg.aperiodic_mode == 'knee_constant'],
         '',
 
         # Peak Parameters
@@ -397,9 +409,9 @@ def gen_results_fg_str(fg, concise=False):
 
         # Goodness if fit
         'Goodness of fit metrics:',
-        '   R2s -  Min: {:6.3f}, Max: {:6.3f}, Mean: {:5.3f}'
+        '   R2s -  Min: {:8.3f}, Max: {:8.3f}, Mean: {:8.3f}'
         .format(np.nanmin(r2s), np.nanmax(r2s), np.nanmean(r2s)),
-        'Errors -  Min: {:6.3f}, Max: {:6.3f}, Mean: {:5.3f}'
+        'Errors -  Min: {:8.3f}, Max: {:8.3f}, Mean: {:8.3f}'
         .format(np.nanmin(errors), np.nanmax(errors), np.nanmean(errors)),
         '',
 

fooof/data/data.py

@@ -28,7 +28,7 @@ class FOOOFSettings(namedtuple('FOOOFSettings', ['peak_width_limits', 'max_n_pea
         Absolute threshold for detecting peaks, in units of the input data.
     peak_threshold : float
         Relative threshold for detecting peaks, in units of standard deviation of the input data.
-    aperiodic_mode : {'fixed', 'knee'}
+    aperiodic_mode : {'fixed', 'knee', 'knee_constant'}
         Which approach to take for fitting the aperiodic component.
 
     Notes
@@ -62,8 +62,9 @@ class FOOOFResults(namedtuple('FOOOFResults', ['aperiodic_params', 'peak_params'
     Parameters
     ----------
     aperiodic_params : 1d array
-        Parameters that define the aperiodic fit. As [Offset, (Knee), Exponent].
+        Parameters that define the aperiodic fit. As [Offset, (Knee), Exponent, Constant].
         The knee parameter is only included if aperiodic is fit with knee.
+        The constant paramter is only included is aperiodic is fit with knee_constant.
     peak_params : 2d array
         Fitted parameter values for the peaks. Each row is a peak, as [CF, PW, BW].
     r_squared : float

fooof/objs/fit.py

@@ -100,7 +100,7 @@ class FOOOF():
         Absolute threshold for detecting peaks, in units of the input data.
     peak_threshold : float, optional, default: 2.0
         Relative threshold for detecting peaks, in units of standard deviation of the input data.
-    aperiodic_mode : {'fixed', 'knee'}
+    aperiodic_mode : {'fixed', 'knee', 'knee_constant'}
         Which approach to take for fitting the aperiodic component.
     verbose : bool, optional, default: True
         Verbosity mode. If True, prints out warnings and general status updates.
@@ -172,12 +172,18 @@ def __init__(self, peak_width_limits=(0.5, 12.0), max_n_peaks=np.inf, min_peak_h
         # Guess parameters for aperiodic fitting, [offset, knee, exponent]
         #   If offset guess is None, the first value of the power spectrum is used as offset guess
         #   If exponent guess is None, the abs(log-log slope) of first & last points is used
-        self._ap_guess = (None, 0, None)
-        # Bounds for aperiodic fitting, as: ((offset_low_bound, knee_low_bound, exp_low_bound),
-        #                                    (offset_high_bound, knee_high_bound, exp_high_bound))
+        self._ap_guess = (None, 1, None, 0)
+        # Bounds for aperiodic fitting, as:
+        #   ((offset_low_bound, knee_low_bound, exp_low_bound, const_low_bound),
+        #    (offset_high_bound, knee_high_bound, exp_high_bound, const_high_bounds))
         # By default, aperiodic fitting is unbound, but can be restricted here, if desired
-        #   Even if fitting without knee, leave bounds for knee (they are dropped later)
-        self._ap_bounds = ((-np.inf, -np.inf, -np.inf), (np.inf, np.inf, np.inf))
+        self._ap_bounds = ((-np.inf, -np.inf, -np.inf, 0), (np.inf, np.inf, np.inf, np.inf))
+
+        if self.aperiodic_mode == 'knee':
+            self._ap_bounds = tuple(bound[:-1] for bound in self._ap_bounds)
+        elif self.aperiodic_mode == 'fixed':
+            self._ap_bounds = tuple(bound[0::2] for bound in self._ap_bounds)
+
         # Threshold for how far a peak has to be from edge to keep.
         #   This is defined in units of gaussian standard deviation
         self._bw_std_edge = 1.0
@@ -276,8 +282,7 @@ def _reset_data_results(self, clear_freqs=False, clear_spectrum=False, clear_res
 
         if clear_results:
 
-            self.aperiodic_params_ = np.array([np.nan] * \
-                (2 if self.aperiodic_mode == 'fixed' else 3))
+            self.aperiodic_params_ = np.array([np.nan] * len(self._ap_bounds[0]))
             self.gaussian_params_ = np.empty([0, 3])
             self.peak_params_ = np.empty([0, 3])
             self.r_squared_ = np.nan
@@ -741,17 +746,14 @@ def _simple_ap_fit(self, freqs, power_spectrum):
         # Get the guess parameters and/or calculate from the data, as needed
         #   Note that these are collected as lists, to concatenate with or without knee later
         off_guess = [power_spectrum[0] if not self._ap_guess[0] else self._ap_guess[0]]
-        kne_guess = [self._ap_guess[1]] if self.aperiodic_mode == 'knee' else []
+        kne_guess = [self._ap_guess[1]] if self.aperiodic_mode != 'fixed' else []
         exp_guess = [np.abs((self.power_spectrum[-1] - self.power_spectrum[0]) /
                             (np.log10(self.freqs[-1]) - np.log10(self.freqs[0])))
                      if not self._ap_guess[2] else self._ap_guess[2]]
-
-        # Get bounds for aperiodic fitting, dropping knee bound if not set to fit knee
-        ap_bounds = self._ap_bounds if self.aperiodic_mode == 'knee' \
-            else tuple(bound[0::2] for bound in self._ap_bounds)
+        const_guess = [self._ap_guess[-1]] if self.aperiodic_mode == 'knee_constant' else []
 
         # Collect together guess parameters
-        guess = np.array(off_guess + kne_guess + exp_guess)
+        guess = np.array(off_guess + kne_guess + exp_guess + const_guess)
 
         # Ignore warnings that are raised in curve_fit
         #   A runtime warning can occur while exploring parameters in curve fitting
@@ -762,7 +764,7 @@ def _simple_ap_fit(self, freqs, power_spectrum):
                 warnings.simplefilter("ignore")
                 aperiodic_params, _ = curve_fit(get_ap_func(self.aperiodic_mode),
                                                 freqs, power_spectrum, p0=guess,
-                                                maxfev=self._maxfev, bounds=ap_bounds)
+                                                maxfev=self._maxfev, bounds=self._ap_bounds)
         except RuntimeError:
             raise FitError("Model fitting failed due to not finding parameters in "
                            "the simple aperiodic component fit.")
@@ -807,18 +809,14 @@ def _robust_ap_fit(self, freqs, power_spectrum):
         freqs_ignore = freqs[perc_mask]
         spectrum_ignore = power_spectrum[perc_mask]
 
-        # Get bounds for aperiodic fitting, dropping knee bound if not set to fit knee
-        ap_bounds = self._ap_bounds if self.aperiodic_mode == 'knee' \
-            else tuple(bound[0::2] for bound in self._ap_bounds)
-
         # Second aperiodic fit - using results of first fit as guess parameters
         #  See note in _simple_ap_fit about warnings
         try:
             with warnings.catch_warnings():
                 warnings.simplefilter("ignore")
                 aperiodic_params, _ = curve_fit(get_ap_func(self.aperiodic_mode),
                                                 freqs_ignore, spectrum_ignore, p0=popt,
-                                                maxfev=self._maxfev, bounds=ap_bounds)
+                                                maxfev=self._maxfev, bounds=self._ap_bounds)
         except RuntimeError:
             raise FitError("Model fitting failed due to not finding "
                            "parameters in the robust aperiodic fit.")

fooof/objs/group.py

@@ -43,7 +43,7 @@ class FOOOFGroup(FOOOF):
         Absolute threshold for detecting peaks, in units of the input data.
     peak_threshold : float, optional, default: 2.0
         Relative threshold for detecting peaks, in units of standard deviation of the input data.
-    aperiodic_mode : {'fixed', 'knee'}
+    aperiodic_mode : {'fixed', 'knee', 'knee_constant'}
         Which approach to take for fitting the aperiodic component.
     verbose : bool, optional, default: True
         Verbosity mode. If True, prints out warnings and general status updates.
@@ -336,7 +336,7 @@ def get_params(self, name, col=None):
         ----------
         name : {'aperiodic_params', 'peak_params', 'gaussian_params', 'error', 'r_squared'}
             Name of the data field to extract across the group.
-        col : {'CF', 'PW', 'BW', 'offset', 'knee', 'exponent'} or int, optional
+        col : {'CF', 'PW', 'BW', 'offset', 'knee', 'exponent', 'constant'} or int, optional
             Column name / index to extract from selected data, if requested.
             Only used for name of {'aperiodic_params', 'peak_params', 'gaussian_params'}.
 
@@ -369,7 +369,7 @@ def get_params(self, name, col=None):
         if isinstance(col, str):
             col = get_indices(self.aperiodic_mode)[col]
         elif isinstance(col, int):
-            if col not in [0, 1, 2]:
+            if col not in [0, 1, 2, 3]:
                 raise ValueError("Input value for `col` not valid.")
 
         # Pull out the requested data field from the group data

fooof/sim/gen.py

@@ -51,7 +51,7 @@ def gen_power_spectrum(freq_range, aperiodic_params, periodic_params, nlv=0.005,
         Frequency range to simulate power spectrum across, as [f_low, f_high], inclusive.
     aperiodic_params : list of float
         Parameters to create the aperiodic component of a power spectrum.
-        Length should be 2 or 3 (see note).
+        Length should be 2, 3, or 4 (see note).
     periodic_params : list of float or list of list of float
         Parameters to create the periodic component of a power spectrum.
         Total length of n_peaks * 3 (see note).
@@ -80,7 +80,8 @@ def gen_power_spectrum(freq_range, aperiodic_params, periodic_params, nlv=0.005,
     Aperiodic Parameters:
 
     - The function for the aperiodic process to use is inferred from the provided parameters.
-    - If length of 2, the 'fixed' aperiodic mode is used, if length of 3, 'knee' is used.
+    - If length of 2, the 'fixed' aperiodic mode is used, if length of 3, 'knee' is used, and
+      if length of 4, 'knee_constant' is used.
 
     Periodic Parameters:
 
@@ -140,7 +141,7 @@ def gen_power_spectrum(freq_range, aperiodic_params, periodic_params, nlv=0.005,
 
         # The rotation changes the offset, so recalculate it's value & update params
         new_offset = compute_rotation_offset(aperiodic_params[1], f_rotation)
-        aperiodic_params = [new_offset, aperiodic_params[1]]
+        aperiodic_params[0] = new_offset
 
     else:
 
@@ -299,7 +300,7 @@ def gen_aperiodic(freqs, aperiodic_params, aperiodic_mode=None):
         Frequency vector to create aperiodic component for.
     aperiodic_params : list of float
         Parameters that define the aperiodic component.
-    aperiodic_mode : {'fixed', 'knee'}, optional
+    aperiodic_mode : {'fixed', 'knee', 'knee_constant'}, optional
         Which kind of aperiodic component to generate.
         If not provided, is inferred from the parameters.
 

fooof/sim/params.py

@@ -44,7 +44,7 @@ def update_sim_ap_params(sim_params, delta, field=None):
         Object storing the current parameter definition.
     delta : float or list of float
         Value(s) by which to update the parameters.
-    field : {'offset', 'knee', 'exponent'} or list of string
+    field : {'offset', 'knee', 'exponent', 'constant'} or list of string
         Field of the aperiodic parameter(s) to update.
 
     Returns