diff --git a/fooof/core/funcs.py b/fooof/core/funcs.py
index e4751c467..e65350db3 100644
--- a/fooof/core/funcs.py
+++ b/fooof/core/funcs.py
@@ -7,6 +7,8 @@
     - They are left available for easy swapping back in, if desired.
 """
 
+from inspect import isfunction
+
 import numpy as np
 
 from fooof.core.errors import InconsistentDataError
@@ -14,15 +16,21 @@
 ###################################################################################################
 ###################################################################################################
 
-def gaussian_function(xs, *params):
+def gaussian_function(xs, cf, pw, bw, *params):
     """Gaussian fitting function.
 
     Parameters
     ----------
     xs : 1d array
         Input x-axis values.
+    cf : float
+        The center of the gaussian.
+    pw : float
+        The height of the gaussian.
+    bw : float
+        The width of the gaussian.
     *params : float
-        Parameters that define gaussian function.
+        Additional centers, heights, and widths.
 
     Returns
     -------
@@ -32,16 +40,17 @@ def gaussian_function(xs, *params):
 
     ys = np.zeros_like(xs)
 
+    params = [cf, pw, bw, *params]
+
     for ii in range(0, len(params), 3):
 
         ctr, hgt, wid = params[ii:ii+3]
-
         ys = ys + hgt * np.exp(-(xs-ctr)**2 / (2*wid**2))
 
     return ys
 
 
-def expo_function(xs, *params):
+def expo_function(xs, offset, knee, exp):
     """Exponential fitting function, for fitting aperiodic component with a 'knee'.
 
     NOTE: this function requires linear frequency (not log).
@@ -50,26 +59,32 @@ 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))
+    offset : float
+        The y-intercept of the fit.
+    knee : float
+        The bend in the fit.
+    exp : float
+        The exponential slope of the fit.
 
     Returns
     -------
     ys : 1d array
         Output values for exponential function.
+
+    Notes
+    -----
+    Parameters (offset, knee, exp) that define Lorentzian function:
+    y = 10^offset * (1/(knee + x^exp))
     """
 
     ys = np.zeros_like(xs)
 
-    offset, knee, exp = params
-
     ys = ys + offset - np.log10(knee + xs**exp)
 
     return ys
 
 
-def expo_nk_function(xs, *params):
+def expo_nk_function(xs, offset, exp):
     """Exponential fitting function, for fitting aperiodic component without a 'knee'.
 
     NOTE: this function requires linear frequency (not log).
@@ -78,34 +93,40 @@ def expo_nk_function(xs, *params):
     ----------
     xs : 1d array
         Input x-axis values.
-    *params : float
-        Parameters (offset, exp) that define Lorentzian function:
-        y = 10^off * (1/(x^exp))
+    offset : float
+        The y-intercept of the fit.
+    exp : float
+        The exponential slope of the fit.
 
     Returns
     -------
     ys : 1d array
         Output values for exponential function, without a knee.
+
+    Notes
+    -----
+    Parameters (offset, exp) that define Lorentzian function:
+    y = 10^off * (1/(x^exp))
     """
 
     ys = np.zeros_like(xs)
 
-    offset, exp = params
-
     ys = ys + offset - np.log10(xs**exp)
 
     return ys
 
 
-def linear_function(xs, *params):
+def linear_function(xs, offset, slope):
     """Linear fitting function.
 
     Parameters
     ----------
     xs : 1d array
         Input x-axis values.
-    *params : float
-        Parameters that define linear function.
+    offset : float
+        The y-intercept of the fit.
+    slope : float
+        The slope of the fit.
 
     Returns
     -------
@@ -115,22 +136,24 @@ def linear_function(xs, *params):
 
     ys = np.zeros_like(xs)
 
-    offset, slope = params
-
     ys = ys + offset + (xs*slope)
 
     return ys
 
 
-def quadratic_function(xs, *params):
+def quadratic_function(xs, offset, slope, curve):
     """Quadratic fitting function.
 
     Parameters
     ----------
     xs : 1d array
         Input x-axis values.
-    *params : float
-        Parameters that define quadratic function.
+    offset : float
+        The y-intercept of the fit.
+    slope : float
+        The slope of the fit.
+    curve : float
+        The curve of the fit.
 
     Returns
     -------
@@ -140,8 +163,6 @@ def quadratic_function(xs, *params):
 
     ys = np.zeros_like(xs)
 
-    offset, slope, curve = params
-
     ys = ys + offset + (xs*slope) + ((xs**2)*curve)
 
     return ys
@@ -167,7 +188,9 @@ def get_pe_func(periodic_mode):
 
     """
 
-    if periodic_mode == 'gaussian':
+    if isfunction(periodic_mode):
+        pe_func = periodic_mode
+    elif periodic_mode == 'gaussian':
         pe_func = gaussian_function
     else:
         raise ValueError("Requested periodic mode not understood.")
@@ -194,7 +217,9 @@ def get_ap_func(aperiodic_mode):
         If the specified aperiodic mode label is not understood.
     """
 
-    if aperiodic_mode == 'fixed':
+    if isfunction(aperiodic_mode):
+        ap_func = aperiodic_mode
+    elif aperiodic_mode == 'fixed':
         ap_func = expo_nk_function
     elif aperiodic_mode == 'knee':
         ap_func = expo_function
diff --git a/fooof/data/__init__.py b/fooof/data/__init__.py
index 87fe8c368..17c0a9367 100644
--- a/fooof/data/__init__.py
+++ b/fooof/data/__init__.py
@@ -1,3 +1,3 @@
 """Data sub-module for FOOOF."""
 
-from .data import FOOOFSettings, FOOOFMetaData, FOOOFResults, SimParams
+from .data import FOOOFSettings, FOOOFMetaData, FOOOFResults, SimParams, FitParams
diff --git a/fooof/data/data.py b/fooof/data/data.py
index aef669f1a..e7d44ef78 100644
--- a/fooof/data/data.py
+++ b/fooof/data/data.py
@@ -8,7 +8,13 @@
     - this means no additional attributes can be defined (which is more memory efficient)
 """
 
-from collections import namedtuple
+from collections import namedtuple, OrderedDict
+
+import numpy as np
+
+from fooof.core.modutils import safe_import
+
+pd = safe_import('pandas')
 
 ###################################################################################################
 ###################################################################################################
@@ -78,8 +84,30 @@ class FOOOFResults(namedtuple('FOOOFResults', ['aperiodic_params', 'peak_params'
     -----
     This object is a data object, based on a NamedTuple, with immutable data attributes.
     """
+
     __slots__ = ()
 
+    def to_dict(self):
+
+        # Combined peak, aperiodic, and goodness of fit params
+        results_dict = OrderedDict()
+        results_dict.update(self.peak_params.to_dict())
+        results_dict.update(self.aperiodic_params.to_dict())
+        results_dict.update(OrderedDict(r_squared=self.r_squared, error=self.error))
+
+        return results_dict
+
+    def to_pandas(self):
+
+        if pd is None:
+            raise ValueError("Pandas is not installed.")
+
+        results_dict = self.to_dict()
+
+        results_df = pd.DataFrame(results_dict)
+
+        return results_df
+
 
 class SimParams(namedtuple('SimParams', ['aperiodic_params', 'periodic_params', 'nlv'])):
     """Parameters that define a simulated power spectrum.
@@ -98,3 +126,33 @@ class SimParams(namedtuple('SimParams', ['aperiodic_params', 'periodic_params',
     This object is a data object, based on a NamedTuple, with immutable data attributes.
     """
     __slots__ = ()
+
+
+class FitParams(np.ndarray):
+
+    def __new__(cls, params, labels):
+
+        return np.asarray(params).view(cls)
+
+    def __init__(self, params, labels):
+
+        self.params = params
+        self.labels = labels
+
+    def to_dict(self):
+
+        params = OrderedDict((k, v) for k, v in \
+            zip(self.labels, self.params.transpose()))
+
+        return params
+
+    def to_pandas(self):
+
+        if pd is None:
+            raise ValueError("Pandas is not installed.")
+
+        params = self.to_dict()
+
+        params = pd.DataFrame(params)
+
+        return params
diff --git a/fooof/objs/fit.py b/fooof/objs/fit.py
index fba745d27..78418262b 100644
--- a/fooof/objs/fit.py
+++ b/fooof/objs/fit.py
@@ -56,6 +56,7 @@
 
 import warnings
 from copy import deepcopy
+from inspect import signature
 
 import numpy as np
 from numpy.linalg import LinAlgError
@@ -67,7 +68,7 @@
 from fooof.core.reports import save_report_fm
 from fooof.core.modutils import copy_doc_func_to_method
 from fooof.core.utils import group_three, check_array_dim
-from fooof.core.funcs import gaussian_function, get_ap_func, infer_ap_func
+from fooof.core.funcs import get_pe_func, get_ap_func, infer_ap_func
 from fooof.core.errors import (FitError, NoModelError, DataError,
                                NoDataError, InconsistentDataError)
 from fooof.core.strings import (gen_settings_str, gen_results_fm_str,
@@ -77,7 +78,7 @@
 from fooof.plts.style import style_spectrum_plot
 from fooof.utils.data import trim_spectrum
 from fooof.utils.params import compute_gauss_std
-from fooof.data import FOOOFResults, FOOOFSettings, FOOOFMetaData
+from fooof.data import FOOOFResults, FOOOFSettings, FOOOFMetaData, FitParams
 from fooof.sim.gen import gen_freqs, gen_aperiodic, gen_periodic, gen_model
 
 ###################################################################################################
@@ -154,8 +155,9 @@ class FOOOF():
     """
     # 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):
+    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',
+                 periodic_mode='gaussian', verbose=True):
         """Initialize object with desired settings."""
 
         # Set input settings
@@ -164,6 +166,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
@@ -439,6 +442,9 @@ def fit(self, freqs=None, power_spectrum=None, freq_range=None):
         if self.verbose:
             self._check_width_limits()
 
+        # Determine the aperiodic and periodic fit funcs
+        self._set_fit_funcs()
+
         # In rare cases, the model fails to fit, and so uses try / except
         try:
 
@@ -479,6 +485,14 @@ def fit(self, freqs=None, power_spectrum=None, freq_range=None):
             # Convert gaussian definitions to peak parameters
             self.peak_params_ = self._create_peak_params(self.gaussian_params_)
 
+            # Create a flexible datatype object for ap/pe parameters
+            pe_labels = list(signature(self._pe_func).parameters.keys())[1:]
+            ap_labels = list(signature(self._ap_func).parameters.keys())[1:]
+
+            self.peak_params_ = FitParams(self.peak_params_, pe_labels)
+            self.gaussian_params_ = FitParams(self.gaussian_params_, pe_labels)
+            self.aperiodic_params_ = FitParams(self.aperiodic_params_, ap_labels)
+
             # Calculate R^2 and error of the model fit
             self._calc_r_squared()
             self._calc_error()
@@ -631,7 +645,6 @@ def get_results(self):
         return FOOOFResults(**{key.strip('_') : getattr(self, key) \
             for key in OBJ_DESC['results']})
 
-
     @copy_doc_func_to_method(plot_fm)
     def plot(self, plot_peaks=None, plot_aperiodic=True, plt_log=False,
              add_legend=True, save_fig=False, file_name=None, file_path=None,
@@ -715,6 +728,11 @@ def set_check_data_mode(self, check_data):
 
         self._check_data = check_data
 
+    def _set_fit_funcs(self):
+        """Set the requested aperiodic and periodic fit functions."""
+
+        self._pe_func = get_pe_func(self.periodic_mode)
+        self._ap_func = get_ap_func(self.aperiodic_mode)
 
     def _check_width_limits(self):
         """Check and warn about peak width limits / frequency resolution interaction."""
@@ -762,8 +780,7 @@ def _simple_ap_fit(self, freqs, power_spectrum):
         try:
             with warnings.catch_warnings():
                 warnings.simplefilter("ignore")
-                aperiodic_params, _ = curve_fit(get_ap_func(self.aperiodic_mode),
-                                                freqs, power_spectrum, p0=guess,
+                aperiodic_params, _ = curve_fit(self._ap_func, freqs, power_spectrum, p0=guess,
                                                 maxfev=self._maxfev, bounds=ap_bounds)
         except RuntimeError:
             raise FitError("Model fitting failed due to not finding parameters in "
@@ -818,9 +835,8 @@ def _robust_ap_fit(self, freqs, power_spectrum):
         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)
+                aperiodic_params, _ = curve_fit(self._ap_func, freqs_ignore, spectrum_ignore,
+                                                p0=popt, maxfev=self._maxfev, bounds=ap_bounds)
         except RuntimeError:
             raise FitError("Model fitting failed due to not finding "
                            "parameters in the robust aperiodic fit.")
@@ -904,7 +920,7 @@ def _fit_peaks(self, flat_iter):
 
             # Collect guess parameters and subtract this guess gaussian from the data
             guess = np.vstack((guess, (guess_freq, guess_height, guess_std)))
-            peak_gauss = gaussian_function(self.freqs, guess_freq, guess_height, guess_std)
+            peak_gauss = self._pe_func(self.freqs, guess_freq, guess_height, guess_std)
             flat_iter = flat_iter - peak_gauss
 
         # Check peaks based on edges, and on overlap, dropping any that violate requirements
@@ -963,7 +979,7 @@ def _fit_peak_guess(self, guess):
 
         # Fit the peaks
         try:
-            gaussian_params, _ = curve_fit(gaussian_function, self.freqs, self._spectrum_flat,
+            gaussian_params, _ = curve_fit(self._pe_func, self.freqs, self._spectrum_flat,
                                            p0=guess, maxfev=self._maxfev, bounds=gaus_param_bounds)
         except RuntimeError:
             raise FitError("Model fitting failed due to not finding "