diff --git a/specparam/core/funcs.py b/specparam/core/funcs.py
index eef4d81b..7f42c3a8 100644
--- a/specparam/core/funcs.py
+++ b/specparam/core/funcs.py
@@ -38,6 +38,29 @@ def gaussian_function(xs, *params):
 
     return ys
 
+def lorentzian_function(xs, *params):
+    """Lorentzian fitting function.
+
+    Parameters
+    ----------
+    xs : 1d array
+        Input x-axis values.
+    *params : float
+        Parameters that define a lorentzian function.
+
+    Returns
+    -------
+    ys : 1d array
+        Output values for lorentzian function.
+    """
+    
+    ys = np.zeros_like(xs)
+
+    for ctr, hgt, wid in zip(*[iter(params)] * 3):
+
+        ys = ys + hgt*wid**2/((xs-ctr)**2+wid**2)
+
+    return ys
 
 def expo_function(xs, *params):
     """Exponential fitting function, for fitting aperiodic component with a 'knee'.
@@ -138,7 +161,7 @@ def get_pe_func(periodic_mode):
 
     Parameters
     ----------
-    periodic_mode : {'gaussian'}
+    periodic_mode : {'gaussian','lorentzian'}
         Which periodic fitting function to return.
 
     Returns
@@ -155,6 +178,10 @@ def get_pe_func(periodic_mode):
 
     if periodic_mode == 'gaussian':
         pe_func = gaussian_function
+
+    elif periodic_mode == 'lorentzian':
+        pe_func = lorentzian_function
+
     else:
         raise ValueError("Requested periodic mode not understood.")
 
diff --git a/specparam/objs/fit.py b/specparam/objs/fit.py
index bb2146f2..c7c99753 100644
--- a/specparam/objs/fit.py
+++ b/specparam/objs/fit.py
@@ -160,7 +160,7 @@ class SpectralModel():
     # pylint: disable=attribute-defined-outside-init
 
     def __init__(self, peak_width_limits=(0.5, 12.0), max_n_peaks=np.inf, min_peak_height=0.0,
-                 peak_threshold=2.0, aperiodic_mode='fixed', verbose=True):
+                 peak_threshold=2.0, periodic_mode:str='gaussian', aperiodic_mode='fixed', verbose=True):
         """Initialize model object."""
 
         # Set input settings
@@ -169,6 +169,7 @@ def __init__(self, peak_width_limits=(0.5, 12.0), max_n_peaks=np.inf, min_peak_h
         self.min_peak_height = min_peak_height
         self.peak_threshold = peak_threshold
         self.aperiodic_mode = aperiodic_mode
+        self.periodic_mode = periodic_mode
         self.verbose = verbose
 
         ## PRIVATE SETTINGS
diff --git a/specparam/sim/gen.py b/specparam/sim/gen.py
index f63b2421..b81aba74 100644
--- a/specparam/sim/gen.py
+++ b/specparam/sim/gen.py
@@ -191,7 +191,7 @@ def gen_rotated_power_vals(freqs, aperiodic_params, periodic_params, nlv, f_rota
     return powers
 
 
-def gen_model(freqs, aperiodic_params, periodic_params, return_components=False):
+def gen_model(freqs, aperiodic_params, periodic_params, return_components=False, periodic_mode:str='gaussian'):
     """Generate a power spectrum model for a given parameter definition.
 
     Parameters
@@ -225,7 +225,7 @@ def gen_model(freqs, aperiodic_params, periodic_params, return_components=False)
     """
 
     ap_fit = gen_aperiodic(freqs, aperiodic_params)
-    pe_fit = gen_periodic(freqs, np.ndarray.flatten(periodic_params))
+    pe_fit = gen_periodic(freqs, np.ndarray.flatten(periodic_params), periodic_mode=periodic_mode)
     full_model = pe_fit + ap_fit
 
     if return_components: