diff --git a/pyleoclim/core/geoseries.py b/pyleoclim/core/geoseries.py
index be8dd444..fd62fc2d 100644
--- a/pyleoclim/core/geoseries.py
+++ b/pyleoclim/core/geoseries.py
@@ -713,11 +713,11 @@ def dashboard(self, figsize=[11, 8], gs=None, plt_kwargs=None, histplt_kwargs=No
         if 'method' in spectral_kwargs.keys():
             pass
         else:
-            spectral_kwargs.update({'method': 'lomb_scargle'})
-        if 'freq_method' in spectral_kwargs.keys():
+            spectral_kwargs.update({'method': 'lomb_scargle'}) # unneeded as it is already the default 
+        if 'freq' in spectral_kwargs.keys():
             pass
         else:
-            spectral_kwargs.update({'freq_method': 'lomb_scargle'})
+            spectral_kwargs.update({'freq': 'lomb_scargle'})
 
         ax['spec'] = fig.add_subplot(gs[-1, -2:])
         spectralfig_kwargs = {} if spectralfig_kwargs is None else spectralfig_kwargs.copy()
diff --git a/pyleoclim/core/lipdseries.py b/pyleoclim/core/lipdseries.py
index fa4e8b2c..21694858 100644
--- a/pyleoclim/core/lipdseries.py
+++ b/pyleoclim/core/lipdseries.py
@@ -871,11 +871,11 @@ def dashboard(self, figsize=[11, 8], plt_kwargs=None, histplt_kwargs=None, spect
             pass
         else:
             spectral_kwargs.update({'method': 'lomb_scargle'})
-        if 'freq_method' in spectral_kwargs.keys():
+        if 'freq' in spectral_kwargs.keys():
             pass
         else:
             if ensemble == False:
-                spectral_kwargs.update({'freq_method': 'lomb_scargle'})
+                spectral_kwargs.update({'freq': 'lomb_scargle'})
             elif ensemble == True:
                 pass
 
diff --git a/pyleoclim/core/multipleseries.py b/pyleoclim/core/multipleseries.py
index f9c73bbf..e3c1aabb 100644
--- a/pyleoclim/core/multipleseries.py
+++ b/pyleoclim/core/multipleseries.py
@@ -1303,7 +1303,7 @@ def detrend(self,method='emd',**kwargs):
             ms.series_list[idx]=s
         return ms
 
-    def spectral(self, method='lomb_scargle', settings=None, mute_pbar=False, freq_method='log', 
+    def spectral(self, method='lomb_scargle', freq=None, settings=None, mute_pbar=False, 
                 freq_kwargs=None, label=None, verbose=False, scalogram_list=None):
         ''' Perform spectral analysis on the timeseries
 
@@ -1312,7 +1312,14 @@ def spectral(self, method='lomb_scargle', settings=None, mute_pbar=False, freq_m
 
         method : str; {'wwz', 'mtm', 'lomb_scargle', 'welch', 'periodogram', 'cwt'}
 
-        freq_method : str; {'log','scale', 'nfft', 'lomb_scargle', 'welch'}
+        freq : str or array, optional
+           Information to produce the frequency vector. 
+           This can be 'log','scale', 'nfft', 'lomb_scargle' or 'welch' or a NumPy array.
+           If a string, will use make_freq_vector with the specified frequency-generating method.
+           If an array, this will be passed directly to the spectral method.
+           If None (default), will use 'log' for WWZ and 'lomb_scargle' for Lomb-Scargle. 
+           This parameter is highly consequential for the WWZ and Lomb-Scargle methods, 
+           but is otherwise ignored, as other spectral methods generate their frequency vector internally.
 
         freq_kwargs : dict
         
@@ -1375,14 +1382,12 @@ def spectral(self, method='lomb_scargle', settings=None, mute_pbar=False, freq_m
 
         .. jupyter-execute::
 
-            import pyleoclim as pyleo
-            url = 'http://wiki.linked.earth/wiki/index.php/Special:WTLiPD?op=export&lipdid=MD982176.Stott.2004'
-            data = pyleo.Lipd(usr_path = url)
-            tslist = data.to_LipdSeriesList()
-            tslist = tslist[2:] # drop the first two series which only concerns age and depth
-            ms = pyleo.MultipleSeries(tslist)
-            ms_psd = ms.spectral()
-
+            soi = pyleo.utils.load_dataset('SOI')
+            nino = pyleo.utils.load_dataset('NINO3')
+            ms = soi & nino
+            ms_psd = ms.spectral(method='mtm')
+            ms_psd.plot()
+                    
         '''
         settings = {} if settings is None else settings.copy()
 
@@ -1395,20 +1400,20 @@ def spectral(self, method='lomb_scargle', settings=None, mute_pbar=False, freq_m
             #OR if the scalogram list is longer than the series list we use as many scalograms from the scalogram list as we need
             if scalogram_list_len >= series_len:
                 for idx, s in enumerate(tqdm(self.series_list, desc='Performing spectral analysis on individual series', position=0, leave=True, disable=mute_pbar)):
-                    psd_tmp = s.spectral(method=method, settings=settings, freq_method=freq_method, freq_kwargs=freq_kwargs, label=label, verbose=verbose,scalogram = scalogram_list.scalogram_list[idx])
+                    psd_tmp = s.spectral(method=method, settings=settings, freq=freq, freq_kwargs=freq_kwargs, label=label, verbose=verbose,scalogram = scalogram_list.scalogram_list[idx])
                     psd_list.append(psd_tmp)
             #If the scalogram list isn't as long as the series list, we re-use all the scalograms we can and then calculate the rest
             elif scalogram_list_len < series_len:
                 for idx, s in enumerate(tqdm(self.series_list, desc='Performing spectral analysis on individual series', position=0, leave=True, disable=mute_pbar)):
                     if idx < scalogram_list_len:
-                        psd_tmp = s.spectral(method=method, settings=settings, freq_method=freq_method, freq_kwargs=freq_kwargs, label=label, verbose=verbose,scalogram = scalogram_list.scalogram_list[idx])
+                        psd_tmp = s.spectral(method=method, settings=settings, freq=freq, freq_kwargs=freq_kwargs, label=label, verbose=verbose,scalogram = scalogram_list.scalogram_list[idx])
                         psd_list.append(psd_tmp)
                     else:
-                        psd_tmp = s.spectral(method=method, settings=settings, freq_method=freq_method, freq_kwargs=freq_kwargs, label=label, verbose=verbose)
+                        psd_tmp = s.spectral(method=method, settings=settings, freq=freq, freq_kwargs=freq_kwargs, label=label, verbose=verbose)
                         psd_list.append(psd_tmp)
         else:
             for s in tqdm(self.series_list, desc='Performing spectral analysis on individual series', position=0, leave=True, disable=mute_pbar):
-                psd_tmp = s.spectral(method=method, settings=settings, freq_method=freq_method, freq_kwargs=freq_kwargs, label=label, verbose=verbose)
+                psd_tmp = s.spectral(method=method, settings=settings, freq=freq, freq_kwargs=freq_kwargs, label=label, verbose=verbose)
                 psd_list.append(psd_tmp)
 
         psds = MultiplePSD(psd_list=psd_list)
diff --git a/pyleoclim/core/psds.py b/pyleoclim/core/psds.py
index 2dcf217b..50b98200 100644
--- a/pyleoclim/core/psds.py
+++ b/pyleoclim/core/psds.py
@@ -249,14 +249,11 @@ def signif_test(self, method='ar1sim', number=None, seed=None, qs=[0.95],
         '''
         from ..core.surrogateseries import SurrogateSeries
 
-        if self.spec_method == 'wwz' and method == 'ar1asym':
+        if self.spec_method in ['wwz','lomb_scargle'] and method == 'ar1asym':
             raise ValueError('Asymptotic solution is not supported for the wwz method')
 
-        if self.spec_method == 'lomb_scargle' and method == 'ar1asym':
-            raise ValueError('Asymptotic solution is not supported for the Lomb-Scargle method')
-
         if method not in ['ar1sim', 'uar1','ar1asym']:
-                raise ValueError("The available methods are 'ar1sim' and 'ar1asym'")
+                raise ValueError("The available methods are 'ar1sim', 'uar1' and 'ar1asym'")
 
         if method in ['ar1sim', 'uar1']:
             signif_scals = None
@@ -310,7 +307,7 @@ def signif_test(self, method='ar1sim', number=None, seed=None, qs=[0.95],
                                                        self.spec_args['param'],
                                                        qs=qs, **settings)
             else:
-                # hard code Mortlet values to obtain the spectrum
+                # hard code Morlet values to obtain the spectrum
                 param = 6
                 fourier_factor = 4 * np.pi / (param + np.sqrt(2 + param**2))
                 scale = 1/(fourier_factor*self.frequency)
diff --git a/pyleoclim/core/series.py b/pyleoclim/core/series.py
index b5cb9d5f..98d18309 100644
--- a/pyleoclim/core/series.py
+++ b/pyleoclim/core/series.py
@@ -2757,7 +2757,7 @@ def detrend(self, method='emd', keep_log=False, preserve_mean = False, **kwargs)
             new.log += ({len(new.log): 'detrend','method': method, 'args': kwargs, 'previous_trend': trend},)
         return new
 
-    def spectral(self, method='lomb_scargle', freq_method='log', freq_kwargs=None, settings=None, label=None, scalogram=None, verbose=False):
+    def spectral(self, method='lomb_scargle', freq=None, freq_kwargs=None, settings=None, label=None, scalogram=None, verbose=False):
         ''' Perform spectral analysis on the timeseries
 
         Parameters
@@ -2765,12 +2765,20 @@ def spectral(self, method='lomb_scargle', freq_method='log', freq_kwargs=None, s
 
         method : str;
             {'wwz', 'mtm', 'lomb_scargle', 'welch', 'periodogram', 'cwt'}
-
-        freq_method : str
-            {'log','scale', 'nfft', 'lomb_scargle', 'welch'}
-
+            Default is Lomb-Scargle, because it can handle unevenly spaced series, and is fast. 
+            However, for evenly spaced data mtm would almost surely be a better choice. 
+            
+        freq : str or array, optional
+           Information to produce the frequency vector. 
+           This can be 'log','scale', 'nfft', 'lomb_scargle' or 'welch' or a NumPy array.
+           If a string, will use make_freq_vector with the specified frequency-generating method.
+           If an array, this will be passed directly to the spectral method.
+           If None (default), will use 'log' for WWZ and 'lomb_scargle' for Lomb-Scargle. 
+           This parameter is highly consequential for the WWZ and Lomb-Scargle methods, 
+           but is otherwise ignored, as other spectral methods generate their frequency vector internally.
+            
         freq_kwargs : dict
-            Arguments for frequency vector
+            Arguments for frequency vector 
 
         settings : dict
             Arguments for the specific spectral method
@@ -2829,28 +2837,33 @@ def spectral(self, method='lomb_scargle', freq_method='log', freq_kwargs=None, s
         .. jupyter-execute::
 
             psd_ls = ts_std.spectral(method='lomb_scargle')
-            psd_ls_signif = psd_ls.signif_test(number=20) #in practice, need more AR1 simulations
+            psd_ls_signif = psd_ls.signif_test(number=20) #in practice, need more AR(1) simulations
             fig, ax = psd_ls_signif.plot(title='PSD using Lomb-Scargle method')
 
         We may pass in method-specific arguments via "settings", which is a dictionary.
         For instance, to adjust the number of overlapping segment for Lomb-Scargle, we may specify the method-specific argument "n50";
-        to adjust the frequency vector, we may modify the "freq_method" or modify the method-specific argument "freq".
+        to adjust the frequency vector, we may modify the "freq" or modify the method-specific argument "freq".
 
         .. jupyter-execute::
 
             import numpy as np
             psd_LS_n50 = ts_std.spectral(method='lomb_scargle', settings={'n50': 4})  # c=1e-2 yields lower frequency resolution
-            psd_LS_freq = ts_std.spectral(method='lomb_scargle', settings={'freq': np.linspace(1/20, 1/0.2, 51)})
-            psd_LS_LS = ts_std.spectral(method='lomb_scargle', freq_method='lomb_scargle')  # with frequency vector generated using REDFIT method
+            psd_LS_freq = ts_std.spectral(method='lomb_scargle',freq=np.linspace(1/20, 1/0.2, 51))
+            psd_LS_log = ts_std.spectral(method='lomb_scargle', freq='log')  # with frequency vector generated using REDFIT method
             fig, ax = psd_LS_n50.plot(
-                title='PSD using Lomb-Scargle method with 4 overlapping segments',
+                title='Lomb-Scargle PSD with 4 overlapping segments',
                 label='settings={"n50": 4}')
-            psd_ls.plot(ax=ax, label='settings={"n50": 3}', marker='o')
+            psd_ls.plot(ax=ax, label='settings={"n50": 3}', marker='o', alpha=.2)
 
             fig, ax = psd_LS_freq.plot(
-                title='PSD using Lomb-Scargle method with different frequency vectors',
+                title='Lomb-Scargle PSD with specified frequency vector',
                 label='freq=np.linspace(1/20, 1/0.2, 51)', marker='o')
-            psd_ls.plot(ax=ax, label='freq_method="log"', marker='o')
+            psd_ls.plot(ax=ax, label='Lomb-Scargle frequency vector', marker='o', alpha=.2)
+            
+            fig, ax = psd_LS_log.plot(
+                title='Lomb-Scargle PSD with different frequency vectors',
+                label='log frequency method', marker='o')
+            psd_ls.plot(ax=ax, label='Lomb-Scargle frequency vector', marker='o', alpha=.2)
 
         You may notice the differences in the PSD curves regarding smoothness and the locations of the analyzed period points.
 
@@ -2879,7 +2892,7 @@ def spectral(self, method='lomb_scargle', freq_method='log', freq_kwargs=None, s
 
             ts_interp = ts_std.interp()
             psd_perio = ts_interp.spectral(method='periodogram')
-            psd_perio_signif = psd_perio.signif_test(number=20, method='ar1sim') #in practice, need more AR1 simulations
+            psd_perio_signif = psd_perio.signif_test(number=20, method='ar1sim') #in practice, need more AR(1) simulations
             fig, ax = psd_perio_signif.plot(title='PSD using Periodogram method')
 
 
@@ -2888,7 +2901,7 @@ def spectral(self, method='lomb_scargle', freq_method='log', freq_kwargs=None, s
         .. jupyter-execute::
 
             psd_welch = ts_interp.spectral(method='welch')
-            psd_welch_signif = psd_welch.signif_test(number=20, method='ar1sim') #in practice, need more AR1 simulations
+            psd_welch_signif = psd_welch.signif_test(number=20, method='ar1sim') #in practice, need more AR(1) simulations
             fig, ax = psd_welch_signif.plot(title='PSD using Welch method')
 
 
@@ -2897,7 +2910,7 @@ def spectral(self, method='lomb_scargle', freq_method='log', freq_kwargs=None, s
         .. jupyter-execute::
 
             psd_mtm = ts_interp.spectral(method='mtm', label='MTM, NW=4')
-            psd_mtm_signif = psd_mtm.signif_test(number=20, method='ar1sim') #in practice, need more AR1 simulations
+            psd_mtm_signif = psd_mtm.signif_test(number=20, method='ar1sim') #in practice, need more AR(1) simulations
             fig, ax = psd_mtm_signif.plot(title='PSD using the multitaper method')
 
 
@@ -2909,7 +2922,6 @@ def spectral(self, method='lomb_scargle', freq_method='log', freq_kwargs=None, s
             psd_mtm2 = ts_interp.spectral(method='mtm', settings={'NW':2}, label='MTM, NW=2')
             fig, ax = psd_mtm2.plot(title='MTM with NW=2')
 
-
         - Continuous Wavelet Transform
 
         .. jupyter-execute::
@@ -2919,7 +2931,6 @@ def spectral(self, method='lomb_scargle', freq_method='log', freq_kwargs=None, s
             psd_cwt_signif = psd_cwt.signif_test(number=20)
             fig, ax = psd_cwt_signif.plot(title='PSD using the CWT method')
 
-
         '''
         if not verbose:
             warnings.simplefilter('ignore')
@@ -2935,15 +2946,35 @@ def spectral(self, method='lomb_scargle', freq_method='log', freq_kwargs=None, s
         }
         args = {}
         freq_kwargs = {} if freq_kwargs is None else freq_kwargs.copy()
-        freq = specutils.make_freq_vector(self.time, method=freq_method, **freq_kwargs)
-
-        args['wwz'] = {'freq': freq}
-        args['cwt'] = {'freq': freq}
-        args['mtm'] = {}
-        args['lomb_scargle'] = {'freq': freq}
-        args['welch'] = {}
-        args['periodogram'] = {}
-        args[method].update(settings)
+        
+        if 'freq' in settings.keys():
+            freq_vec = settings['freq']
+        else:
+            if freq is None: # assign the frequency method automatically based on context
+                if method in ['wwz','cwt']:
+                    freq_vec = specutils.make_freq_vector(self.time, method='log', **freq_kwargs) 
+                elif method=='lomb_scargle':
+                    freq_vec = specutils.make_freq_vector(self.time, method='lomb_scargle', **freq_kwargs) 
+                else:
+                    warnings.warn(f'freq argument ignored; it is determined automatically by {method}')
+            elif isinstance(freq, str):   # apply the specified method     
+                freq_vec = specutils.make_freq_vector(self.time, method=freq, **freq_kwargs) 
+            elif isinstance(freq,np.ndarray): # use the specified vector if dimensions check out
+                freq_vec = np.squeeze(freq)
+                if freq.ndim != 1:
+                    raise ValueError("freq should be a 1-dimensional array")
+        
+        # pass frequency
+        if method in ['wwz','cwt','lomb_scargle']:
+            args[method] = {'freq': freq_vec}
+        else:
+            args[method] = {}
+            
+        # args['mtm'] = {'freq': freq_vec}
+        # args['welch'] = {}
+        # args['periodogram'] = {}
+        
+        args[method].update(settings) # this overrides the frequency vector
 
         if method == 'wwz' and scalogram is not None:
             args['wwz'].update(
@@ -3006,9 +3037,9 @@ def wavelet(self, method='cwt', settings=None, freq_method='log', freq_kwargs=No
                 is appropriate for unevenly-spaced series.
             Default is cwt, returning an error if the Series is unevenly-spaced.
 
-        freq_method : str
-            {'log', 'scale', 'nfft', 'lomb_scargle', 'welch'}
-
+        freq_method : str, optional
+            Can be one of 'log', 'scale', 'nfft', 'lomb_scargle', 'welch'. 
+            
         freq_kwargs : dict
             Arguments for the frequency vector
 
diff --git a/pyleoclim/tests/test_core_Series.py b/pyleoclim/tests/test_core_Series.py
index 32a32ff5..d7a3c9a3 100644
--- a/pyleoclim/tests/test_core_Series.py
+++ b/pyleoclim/tests/test_core_Series.py
@@ -213,52 +213,52 @@ def test_spectral_t0(self, pinkseries, spec_method, eps=0.5):
         beta = psd.beta_est().beta_est_res['beta']
         assert np.abs(beta-1.0) < eps
 
-    @pytest.mark.parametrize('freq_method', ['log', 'scale', 'nfft', 'lomb_scargle', 'welch'])
-    def test_spectral_t1(self, pinkseries, freq_method, eps=0.3):
-        ''' Test Series.spectral() with MTM using available `freq_method` options with other arguments being default
+    @pytest.mark.parametrize('freq', ['log', 'scale', 'nfft', 'lomb_scargle', 'welch'])
+    def test_spectral_t1(self, pinkseries, freq, eps=0.3):
+        ''' Test Series.spectral() with MTM using available `freq` options with other arguments being default
 
         We will estimate the scaling slope of an ideal colored noise to make sure the result is reasonable.
         '''
         ts = pinkseries
-        psd = ts.spectral(method='mtm', freq_method=freq_method)
+        psd = ts.spectral(method='mtm', freq=freq)
         beta = psd.beta_est().beta_est_res['beta']
         assert np.abs(beta-1.0) < eps
 
     @pytest.mark.parametrize('nf', [10, 20, 30])
     def test_spectral_t2(self, pinkseries, nf, eps=0.3):
-        ''' Test Series.spectral() with MTM using `freq_method='log'` with different values for its keyword argument `nfreq`
+        ''' Test Series.spectral() with MTM using `freq='log'` with different values for its keyword argument `nfreq`
 
         We will estimate the scaling slope of an ideal colored noise to make sure the result is reasonable.
         '''
         ts = pinkseries
-        psd = ts.spectral(method='mtm', freq_method='log', freq_kwargs={'nf': nf})
+        psd = ts.spectral(method='mtm', freq='log', freq_kwargs={'nf': nf})
         beta = psd.beta_est().beta_est_res['beta']
         assert np.abs(beta-1.0) < eps
 
     @pytest.mark.parametrize('dj', [0.25, 0.5, 1])
     def test_spectral_t3(self, pinkseries, dj, eps=0.3):
-        ''' Test Series.spectral() with MTM using `freq_method='scale'` with different values for its keyword argument `nv`
+        ''' Test Series.spectral() with MTM using `freq='scale'` with different values for its keyword argument `nv`
 
         We will estimate the scaling slope of an ideal colored noise to make sure the result is reasonable.
         '''
         ts = pinkseries
-        psd = ts.spectral(method='mtm', freq_method='scale', freq_kwargs={'dj': dj})
+        psd = ts.spectral(method='mtm', freq='scale', freq_kwargs={'dj': dj})
         beta = psd.beta_est().beta_est_res['beta']
         assert np.abs(beta-1.0) < eps
 
     @pytest.mark.parametrize('dt, nf, ofac, hifac', [(None, 20, 1, 1), (None, None, 2, 0.5)])
     def test_spectral_t4(self, pinkseries, dt, nf, ofac, hifac, eps=0.5):
-        ''' Test Series.spectral() with Lomb_Scargle using `freq_method=lomb_scargle` with different values for its keyword arguments
+        ''' Test Series.spectral() with Lomb_Scargle using `freq=lomb_scargle` with different values for its keyword arguments
 
         We will estimate the scaling slope of an ideal colored noise to make sure the result is reasonable.
         '''
         ts = pinkseries
-        psd = ts.spectral(method='mtm', freq_method='lomb_scargle', freq_kwargs={'dt': dt, 'nf': nf, 'ofac': ofac, 'hifac': hifac})
+        psd = ts.spectral(method='mtm', freq='lomb_scargle', freq_kwargs={'dt': dt, 'nf': nf, 'ofac': ofac, 'hifac': hifac})
         beta = psd.beta_est().beta_est_res['beta']
         assert np.abs(beta-1.0) < eps
 
     def test_spectral_t5(self, pinkseries, eps=0.6):
-        ''' Test Series.spectral() with WWZ with specified frequency vector passed via `settings`
+        ''' Test Series.spectral() with LS with specified frequency vector passed via `settings`
 
         We will estimate the scaling slope of an ideal colored noise to make sure the result is reasonable.
         Note `asser_array_equal(psd.frequency, freq)` is used to make sure the specified frequency vector is really working.
@@ -266,7 +266,7 @@ def test_spectral_t5(self, pinkseries, eps=0.6):
         '''
         ts = pinkseries
         freq = np.linspace(1/500, 1/2, 100)
-        psd = ts.spectral(method='wwz', settings={'freq': freq}, label='WWZ')
+        psd = ts.spectral(method='lomb_scargle', freq=freq)
         beta = psd.beta_est(fmin=1/200, fmax=1/10).beta_est_res['beta']
         assert_array_equal(psd.frequency, freq)
         assert np.abs(beta-1.0) < eps
@@ -284,7 +284,7 @@ def test_spectral_t6(self, pinkseries, spec_method, eps=0.3):
         t_unevenly =  np.delete(ts.time, deleted_idx)
         v_unevenly =  np.delete(ts.value, deleted_idx)
 
-        ts2 = pyleo.Series(time=t_unevenly, value=v_unevenly)
+        ts2 = pyleo.Series(time=t_unevenly, value=v_unevenly,auto_time_params=True)
         psd = ts2.spectral(method=spec_method)
         beta = psd.beta_est().beta_est_res['beta']
         assert np.abs(beta-1.0) < eps
diff --git a/pyleoclim/utils/spectral.py b/pyleoclim/utils/spectral.py
index 9df79105..efed16ae 100644
--- a/pyleoclim/utils/spectral.py
+++ b/pyleoclim/utils/spectral.py
@@ -346,10 +346,10 @@ def lomb_scargle(ys, ts, freq=None, freq_method='lomb_scargle',
                  gaussianize=False,
                  standardize=True,
                  average='mean'):
-    """ Lomb-scargle priodogram
+    """ Lomb-scargle periodogram
 
     Appropriate for unevenly-spaced arrays.
-    Uses the lomb-scargle implementation from scipy.signal: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lombscargle.html
+    Uses the lomb-scargle implementation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lombscargle.html
 
     Parameters
     ----------
@@ -472,8 +472,6 @@ def lomb_scargle(ys, ts, freq=None, freq_method='lomb_scargle',
 
     Scargle, J. D. (1982). Studies in astronomical time series analysis. II. Statistical aspects of spectral analyis of unvenly spaced data. The Astrophysical Journal, 263(2), 835-853.
 
-    Scargle, J. D. (1982). Studies in astronomical time series analysis. II. Statistical aspects of spectral analyis of unvenly spaced data. The Astrophysical Journal, 263(2), 835-853.
-
     """
     
     if standardize == True:
@@ -1281,6 +1279,6 @@ def psd_fBM(freq, ts, H):
 
     for k in range(nf):
         tmp = 2 * omega[k] * T
-        psd[k] = (1 - 2**(1 - 2*H)*np.sin(tmp)/tmp) / np.abs(omega[k])**(1 + 2*H)
+        psd[k] = (1 - 2**(1 - 2*H)*np.sin(tmp)/tmp) / np.abs(omega[k])**(1 + 2*H) # Eq 6 from Flandrin, 1989
 
     return psd
diff --git a/pyleoclim/utils/wavelet.py b/pyleoclim/utils/wavelet.py
index 39ecf190..cefff8d7 100644
--- a/pyleoclim/utils/wavelet.py
+++ b/pyleoclim/utils/wavelet.py
@@ -1735,7 +1735,7 @@ def wwz_coherence(y1, t1, y2, t2, smooth_factor=0.25,
 
     return res
 
-def freq_vector_lomb_scargle(ts, dt= None, nf=None, ofac=4, hifac=1):
+def freq_vector_lomb_scargle(ts, dt= None, nf=None, ofac=4, hifac=0.95):
     ''' Return the frequency vector based on the REDFIT recommendation.
 
     Parameters
@@ -1774,8 +1774,9 @@ def freq_vector_lomb_scargle(ts, dt= None, nf=None, ofac=4, hifac=1):
     References
     ----------
 
-    Trauth, M. H. MATLAB® Recipes for Earth Sciences. (Springer, 2015). pp 181.
-
+    - Trauth, M. H. MATLAB® Recipes for Earth Sciences. (Springer, 2015). pp 181.
+    
+    - https://www.mathworks.com/matlabcentral/fileexchange/22215-lomb-normalized-periodogram?focused=5108122&tab=function
 
     See also
     --------
@@ -1800,6 +1801,7 @@ def freq_vector_lomb_scargle(ts, dt= None, nf=None, ofac=4, hifac=1):
     if nf is None:
         df = flo
         nf = int((fhi - flo) / df + 1)
+        
 
     freq = np.linspace(flo, fhi, nf)
 
@@ -2032,17 +2034,15 @@ def freq_vector_log(ts, fmin=None, fmax= None, nf=None):
     if nf is None:
         nf = nt//10 + 1
     if fmin is None:    
-        fmin = 2/(np.max(ts)-np.min(ts))
+        fmin = 1/np.ptp(ts)  # used to be 2/np.ptp(ts)
     if fmax is None:
         fmax = fs/2
         
-    start = np.log2(fmin)
-    stop = np.log2(fmax)
-    freq = np.logspace(start, stop, nf, base=2)
+    freq = np.logspace(np.log2(fmin), np.log2(fmax), nf, base=2)
 
     return freq
 
-def make_freq_vector(ts, method='log', **kwargs):
+def make_freq_vector(ts, method='log', rayleigh_bound=True, **kwargs):
     ''' Make frequency vector
 
     This function selects among five methods to obtain the frequency
@@ -2059,23 +2059,13 @@ def make_freq_vector(ts, method='log', **kwargs):
 
         The method to use. Options are 'log' (default), 'nfft', 'lomb_scargle', 'welch', and 'scale'
         
+    rayleigh_bound : boolean
+    
+        If True (default), sets the minimum frequency to the Rayleigh frequency 1/L, where L is the length of the series
+        
     kwargs : dict, optional
-            -fmin : float
-                minimum frequency. If None is provided (default), inferred by the method.
-                
-            - fmax : float
-                maximum frequency. If None is provided (default), inferred by the method. 
-                
-            - nf (int): number of frequency points
     
-            For Lomb_Scargle, additional parameters may be passed:
-
-            - ofac (float): Oversampling rate that influences the resolution of the frequency axis,
-                 when equals to 1, it means no oversamling (should be >= 1).
-                 The default value 4 is usaually a good value.
-            - hifac (float): fhi/fnyq (should be >= 1), where fhi is the highest frequency that
-                  can be analyzed by the Lomb-Scargle algorithm and fnyq is the Nyquist frequency.
-
+            See underlying methods for optional parameters
     Returns
     -------
 
@@ -2099,19 +2089,24 @@ def make_freq_vector(ts, method='log', **kwargs):
     '''
 
     if method == 'lomb_scargle':
-        freq = freq_vector_lomb_scargle(ts,**kwargs)
+        fr = freq_vector_lomb_scargle(ts,**kwargs)
     elif method == 'welch':
-        freq = freq_vector_welch(ts)
+        fr = freq_vector_welch(ts)
     elif method == 'nfft':
-        freq = freq_vector_nfft(ts)
+        fr = freq_vector_nfft(ts)
     elif method == 'scale':
-        freq = freq_vector_scale(ts, **kwargs)
+        fr = freq_vector_scale(ts, **kwargs)
     elif method == 'log':
-        freq = freq_vector_log(ts, **kwargs)
+        fr = freq_vector_log(ts, **kwargs)
     else:
         raise ValueError('This method is not supported')
-    #  freq = freq[1:]  # discard the first element 0
-
+        
+    if rayleigh_bound:
+        fR = 1/np.ptp(ts) # Rayleigh frequency
+        freq = fr[fr>=fR]
+    else:
+        freq = fr
+        
     return freq