diff --git a/fooof_mat/fooof.m b/fooof_mat/fooof.m
index 30636c5..bb9fd3e 100644
--- a/fooof_mat/fooof.m
+++ b/fooof_mat/fooof.m
@@ -1,7 +1,7 @@
 % fooof() - Fit the FOOOF model on a neural power spectrum.
 %
 % Usage:
-%   >> fooof_results = fooof(freqs, power_spectrum, f_range, settings);
+%   >> fooof_results = fooof(freqs, power_spectrum, f_range, settings, return_model);
 %
 % Inputs:
 %   freqs           = row vector of frequency values
@@ -14,7 +14,7 @@
 %       settings.peak_threshold
 %       settings.aperiodic_mode
 %       settings.verbose
-%   return_model    = boolean of whether to return actual model, optional
+%   return_model    = boolean of whether to return the FOOOF model fit, optional
 %
 % Outputs:
 %   fooof_results   = fooof model ouputs, in a struct, including:
@@ -30,9 +30,9 @@
 %            fooof_results.ap_fit
 %
 % Notes
-%   Not all settings need to be set. Any settings that are not
-%     provided as set to default values. To run with all defaults,
-%     input settings as an empty struct.
+%   Not all settings need to be defined by the user.
+%     Any settings that are not provided are set to default values.
+%     To run with all defaults, input settings as an empty struct.
 
 function fooof_results = fooof(freqs, power_spectrum, f_range, settings, return_model)
 
@@ -63,7 +63,7 @@
     %   This will default to not return model, if variable not set
     if exist('return_model', 'var') && return_model
 
-        % Get the model, and add outputs to foof_results
+        % Get the model, and add outputs to fooof_results
         model_out = fooof_get_model(fm);
         for field = fieldnames(model_out)'
             fooof_results.(field{1}) = model_out.(field{1});
diff --git a/fooof_mat/fooof_check_settings.m b/fooof_mat/fooof_check_settings.m
index df5c6e0..0563ca4 100644
--- a/fooof_mat/fooof_check_settings.m
+++ b/fooof_mat/fooof_check_settings.m
@@ -12,7 +12,7 @@
 %       settings.aperiodic_mode
 %       settings.verbose
 %
-% Returns:
+% Outputs:
 %   settings        = struct, with all settings defined:
 %       settings.peak_width_limts
 %       settings.max_n_peaks
@@ -20,7 +20,7 @@
 %       settings.peak_threshold
 %       settings.aperiodic_mode
 %       settings.verbose
-
+%
 % Notes:
 %   This is a helper function, probably not called directly by the user.
 %   Any settings not specified are set to default values
@@ -36,7 +36,7 @@
         'aperiodic_mode', 'fixed', ...
         'verbose', true);
 
-    % Overwrite any non-existent of nan settings with defaults
+    % Overwrite any non-existent or nan settings with defaults
     for field = fieldnames(defaults)'
         if ~isfield(settings, field) || all(isnan(settings.(field{1})))
             settings.(field{1}) = defaults.(field{1});
diff --git a/fooof_mat/fooof_get_model.m b/fooof_mat/fooof_get_model.m
index 209bd09..1470546 100644
--- a/fooof_mat/fooof_get_model.m
+++ b/fooof_mat/fooof_get_model.m
@@ -14,11 +14,9 @@
 %       model_fit.ap_fit
 %
 % Notes
-%   This function is mostly an internal function, but can be called
-%     directly by the user if you are interacting with FOOOF objects
-%     directly.
+%   This function is mostly an internal function.
+%     It can be called directly by the user if you are interacting with FOOOF objects directly.
 
-% Get the actual fit model from a FOOOF object
 function model_fit = fooof_get_model(fm)
 
     model_fit = struct();
diff --git a/fooof_mat/fooof_plot.m b/fooof_mat/fooof_plot.m
index 1920f3a..798c055 100644
--- a/fooof_mat/fooof_plot.m
+++ b/fooof_mat/fooof_plot.m
@@ -13,11 +13,13 @@
 function fooof_plot(fooof_results, log_freqs)
 
     %% Data Checking
+
     if ~isfield(fooof_results, 'freqs')
        error('FOOOF results struct does not contain model output.')
     end
 
     %% Set Up
+
     if exist('log_freqs', 'var') && log_freqs
         plt_freqs = log10(fooof_results.freqs);
     else
@@ -27,17 +29,17 @@ function fooof_plot(fooof_results, log_freqs)
     % Plot settings
     lw = 2.5;
 
-    %% Create the Plots
+    %% Create the plots
 
     figure()
     hold on
+
     % Plot the original data
     data = plot(plt_freqs, fooof_results.power_spectrum, 'black');
 
     % Plot the full model fit
     model = plot(plt_freqs, fooof_results.fooofed_spectrum, 'red');
 
-
     % Plot the aperiodic fit
     ap_fit = plot(plt_freqs, fooof_results.ap_fit, 'b--');
 
diff --git a/fooof_mat/fooof_unpack_results.m b/fooof_mat/fooof_unpack_results.m
index d71c40b..0412c7c 100644
--- a/fooof_mat/fooof_unpack_results.m
+++ b/fooof_mat/fooof_unpack_results.m
@@ -1,10 +1,11 @@
 % fooof_unpack_results() - Extract model fit results from FOOOFResults.
 %
 % Usage:
-%   >> results_out = fooof_unpack_results(fooof_results);
+%   >> results_out = fooof_unpack_results(results_in);
 %
 % Inputs:
 %   fooof_results   = FOOOFResults object
+%
 % Outputs:
 %   results_out     = fooof model results, in a struct, including:
 %       results_out.aperiodic_params
@@ -13,10 +14,9 @@
 %       results_out.error
 %       results_out.r_squared
 %
-% Notes:
-%   This function is mostly an internal function, and doesn't need to be
-%   called directly by the user - but can be if you are interacting
-%   directly with FOOOF objects.
+% Notes
+%   This function is mostly an internal function.
+%     It can be called directly by the user if you are interacting with FOOOF objects directly.
 
 function results_out = fooof_unpack_results(results_in)
 
@@ -36,10 +36,11 @@
     results_out.error = ...
         double(py.array.array('d', py.numpy.nditer(results_in.error)));
 
-    % Note: r_squared seems to come out as float
-    %   It's not clear why this value is different, but doesn't seem
-    %   to require the type casting like other (code commented below)
     results_out.r_squared = results_in.r_squared;
+
+    % Note: r_squared seems to come out as float, and does not need type casting
+    %   It is not clear why this value is different
+    %   Just in case, the code for type casting is commented out below:
     %results_out.r_squared = ...
     %    double(py.array.array('d', py.numpy.nditer(results_in.r_squared)));
 
diff --git a/fooof_mat/fooof_version.m b/fooof_mat/fooof_version.m
index 3c1a63b..b8488ec 100644
--- a/fooof_mat/fooof_version.m
+++ b/fooof_mat/fooof_version.m
@@ -1,7 +1,11 @@
-% fooof_version() - Get FOOOF version information, of both Python & Wrapper.
+% fooof_version() - Get FOOOF version information.
 %
 % Usage:
 %   >> fooof_version()
+%
+% Notes
+%   This function returns versions of both the installed FOOOF in Python, and of this Matlab wrapper.
+%     The version number of the Matlab FOOOF wrapper is defined in this function.
 
 function fooof_version()
 
diff --git a/fooof_mat/new_fooof_group.m b/fooof_mat/new_fooof_group.m
deleted file mode 100644
index 6e908c8..0000000
--- a/fooof_mat/new_fooof_group.m
+++ /dev/null
@@ -1,26 +0,0 @@
-% NOTE: possible WIP for using FOOOFGroup directly
-%   Current problem: can't convert 2d arrays...
-
-function out = new_fooof_group(freqs, psds, f_range, settings)
-
-    % Check settings - get defaults for those not provided
-    settings = fooof_check_settings(settings);
-
-    % Initialize object to collect FOOOF results
-    fooof_results = [];
-
-    % Convert inputs
-    freqs = py.numpy.array(freqs);
-    psds = py.numpy.array(psds);
-    f_range = py.list(f_range);
-
-    % Initialize FOOOF object
-    fg = py.fooof.FOOOFGroup(settings.peak_width_limits, ...
-                             settings.max_n_peaks, ...
-                             settings.min_peak_amplitude, ...
-                             settings.peak_threshold, ...
-                             settings.aperiodic_mode, ...
-                             settings.verbose);
-    fg.fit(freqs, psds, f_range)
-
-    out = 1
\ No newline at end of file