[ENH] - Add gamma periodic mode definition

specparam/modes/definitions.py

@@ -6,7 +6,8 @@
 from specparam.modes.mode import Mode
 from specparam.modes.params import ParamDefinition
 from specparam.modes.funcs import (expo_function, expo_nk_function, double_expo_function,
-                                   gaussian_function, skewed_gaussian_function, cauchy_function)
+                                   gaussian_function, skewed_gaussian_function,
+                                   cauchy_function, gamma_function)
 from specparam.modes.jacobians import jacobian_gauss
 from specparam.utils.checks import check_selection
 
@@ -150,11 +151,34 @@
 )
 
 
+## PE - Gamma Mode
+
+params_gamma = ParamDefinition(OrderedDict({
+    'cf' : 'Center frequency of the peak.',
+    'pw' : 'Power of the peak, over and above the aperiodic component.',
+    'shape' : 'Shape parameter of the gamma function.',
+    'scale' : 'Scale parameter of the gamma function.',
+}))
+
+pe_gamma = Mode(
+    name='gamma',
+    component='periodic',
+    description='Fit a gamma peak function.',
+    func=gamma_function,
+    jacobian=None,
+    params=params_gamma,
+    ndim=2,
+    freq_space='linear',
+    powers_space='log10',
+)
+
+
 # Collect available periodic modes
 PE_MODES = {
     'gaussian' : pe_gaussian,
     'skewed_gaussian' : pe_skewed_gaussian,
     'cauchy' : pe_cauchy,
+    'gamma' : pe_gamma,
 }
 
 ###################################################################################################

specparam/modes/funcs.py

@@ -1,8 +1,11 @@
 """Functions that can be used for model fitting."""
 
+from math import gamma
+
 import numpy as np
 from scipy.special import erf
 
+from specparam.utils.array import normalize
 from specparam.utils.array import normalize
 
 ###################################################################################################
@@ -90,6 +93,41 @@ def cauchy_function(xs, *params):
     return ys
 
 
+def gamma_function(xs, *params):
+    """Gamma fitting function.
+
+    Parameters
+    ----------
+    xs : 1d array
+        Input x-axis values.
+    *params : float
+        Parameters that define a gamma function.
+
+    Returns
+    -------
+    ys : 1d array
+        Output values for gamma function.
+
+    Notes
+    -----
+    Parameters should come in ordered sets of 4, each including the following:
+
+    - cf: center frequency parameter
+    - pw: power (height) parameter
+    - shp: shape parameter
+    - scl: scale parameter
+    """
+
+    ys = np.zeros_like(xs)
+
+    for ctr, hgt, shp, scale in zip(*[iter(params)] * 4):
+        cxs = xs-ctr
+        cxs = cxs.clip(min=0)
+        ys = ys + hgt * normalize((1 / (gamma(shp) * scale**shp) * cxs**(shp-1) * np.exp(-cxs/scale)))
+
+    return ys
+
+
 ## APERIODIC FUNCTIONS
 
 def expo_function(xs, *params):

specparam/tests/modes/test_funcs.py

@@ -52,6 +52,15 @@ def test_cauchy_function():
 
     assert np.all(ys)
 
+def test_gamma_function():
+
+    ctr, hgt, shp, scl = 50, 5, 2, 2
+
+    xs = np.arange(1, 100, 1.)
+    ys = gamma_function(xs, ctr, hgt, shp, scl)
+
+    assert np.any(ys)
+
 ## Aperiodic functions
 
 def test_expo_function():