diff --git a/specparam/modes/definitions.py b/specparam/modes/definitions.py
index 2384d3e2..d53d2b58 100644
--- a/specparam/modes/definitions.py
+++ b/specparam/modes/definitions.py
@@ -7,7 +7,7 @@
 from specparam.modes.params import ParamDefinition
 from specparam.modes.funcs import (powerlaw_function, lorentzian_function, double_expo_function,
                                    gaussian_function, skewed_gaussian_function,
-                                   cauchy_function, triangle_function)
+                                   cauchy_function, gamma_function, triangle_function)
 from specparam.modes.jacobians import jacobian_gauss
 from specparam.utils.checks import check_selection
 
@@ -156,6 +156,30 @@
     powers_space='log10',
 )
 
+
+## 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.',
+    formula=r'P(F)_n = a * \frac{1}{\Gamma (k)\theta^{k}}F_c^{k-1}e^{-\frac{F_c}{\theta}}',
+    func=gamma_function,
+    jacobian=None,
+    params=params_gamma,
+    ndim=2,
+    freq_space='linear',
+    powers_space='log10',
+)
+
+
 ## PE - Triangle Mode
 
 params_triangle = ParamDefinition(OrderedDict({
@@ -183,6 +207,7 @@
     'gaussian' : pe_gaussian,
     'skewed_gaussian' : pe_skewed_gaussian,
     'cauchy' : pe_cauchy,
+    'gamma' : pe_gamma,
     'triangle' : pe_triangle,
 }
 
diff --git a/specparam/modes/funcs.py b/specparam/modes/funcs.py
index 14752879..ec5ac6ac 100644
--- a/specparam/modes/funcs.py
+++ b/specparam/modes/funcs.py
@@ -13,9 +13,12 @@
 - `P` : relating to the periodic component.
 """
 
+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
 
 ###################################################################################################
@@ -127,6 +130,39 @@ def cauchy_function(xs, *params):
     return ys
 
 
+def gamma_function(xs, *params):
+    r"""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
+    -----
+    Defines a gamma fit function as:
+
+        P(F)_n = a * \frac{1}{\Gamma (k)\theta^{k}}F_c^{k-1}e^{-\frac{F_c}{\theta}}
+
+    where k is the shape parameter, theta is the scale parameter, and F_c is the
+    frequency values scaled to center the peak at the desired center frequency.
+    """
+
+    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)))
+
+
 def triangle_function(xs, *params):
     r"""Triangle fitting function.
 
diff --git a/specparam/tests/modes/test_funcs.py b/specparam/tests/modes/test_funcs.py
index 4af73131..e8dec36d 100644
--- a/specparam/tests/modes/test_funcs.py
+++ b/specparam/tests/modes/test_funcs.py
@@ -49,6 +49,13 @@ 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)
+
 def test_triangle_function():
 
     ctr, hgt, wid = 50, 5, 10