Repository: LinkedEarth/Pyleoclim_util (language: python)

Issue:
Clean up nomenclature for coherence/cross-wavelet transform
**Is your feature request related to a problem? Please describe.**
Nomenclature in developers API is confusing and do not track to the user APIS. 

**Describe the solution you'd like**
Clean up and allow access to both coherence and xwt. 



Discussion / hints:
Suggestion: mirror Grinsted's nomenclature and call `xwt` what is currently `cross_wt`  and `wtc` what is currently `wavelet_coherence`. Also, `xwt` needs to be exposed at the user API level (currently only `wwz`, and `xwc` are. )

Fix patch:
diff --git a/README.md b/README.md
index 0b95a09b..770d09ad 100644
--- a/README.md
+++ b/README.md
@@ -13,7 +13,7 @@
 Paleoclimate data, whether from observations or model simulations, offer unique challenges to the analyst.
 Pyleoclim is a Python package primarily geared towards the analysis and visualization of paleoclimate data. Such data usually come in the form of timeseries with missing values and age uncertainties, so the package includes several low-level methods to deal with these issues to simplify the user's life, with intuitive, high-level analysis and plotting methods that support publication-quality scientific workflows.
 
-There are many entry points to Pyleoclim, thanks to its underlying [data structures](https://pyleoclim-util.readthedocs.io/en/stable/core/ui.html). The package makes use of the Linked Paleo Data ([LiPD](http://www.clim-past.net/12/1093/2016/)) standard container and its associated [utilities](http://nickmckay.github.io/LiPD-utilities/). The package is aware of age ensembles stored via LiPD and uses them for time-uncertain analyses, very much like [GeoChronR](https://doi.org/10.5194/gchron-2020-25).
+There are many entry points to Pyleoclim, thanks to its underlying [data structures](https://pyleoclim-util.readthedocs.io/en/master/core/ui.html). The package makes use of the Linked Paleo Data ([LiPD](http://www.clim-past.net/12/1093/2016/)) standard container and its associated [utilities](http://nickmckay.github.io/LiPD-utilities/). The package is aware of age ensembles stored via LiPD and uses them for time-uncertain analyses, very much like [GeoChronR](https://doi.org/10.5194/gchron-2020-25).
 
 LiPD, however, is not an obligatory entry point to Pyleoclim. Low-level modules are primarily based on [NumPy](http://www.numpy.org) arrays or [Pandas](https://pandas.pydata.org) dataframes, so Pyleoclim contains a lot of timeseries analysis code (e.g. spectral analysis, singular spectrum analysis, wavelet analysis, correlation analysis) that can apply to these more common types as well, including those generated by numerical models (via [xarray](http://xarray.pydata.org)). This makes the package suitable for rigorous model-data comparisons, like [this one](https://www.pnas.org/content/116/18/8728.short).
 
diff --git a/doc_build/installation.rst b/doc_build/installation.rst
index cd1fbf55..4b2d2c21 100644
--- a/doc_build/installation.rst
+++ b/doc_build/installation.rst
@@ -103,3 +103,13 @@ To do that, a Fortran compiler (e.g. :code:`gfortran` or :code:`ifort`) is requi
   print(pyleo)
 
 Again, unless you are planning to make heavy use of the WWZ functionality, we recommend using the default installation.
+
+Docker Container
+""""""""""""""""
+Docker containers with various versions of Pyleoclim are available at: `https://quay.io/repository/2i2c/paleohack-2021?tab=tags <https://quay.io/repository/2i2c/paleohack-2021?tab=tags>`_.
+
+To pull an image:
+
+.. code-block:: bash
+
+  docker pull quay.io/2i2c/paleohack-2021:latest
diff --git a/doc_build/savefig/coh.png b/doc_build/savefig/coh.png
index 47a25f13..1126cb31 100644
Binary files a/doc_build/savefig/coh.png and b/doc_build/savefig/coh.png differ
diff --git a/doc_build/savefig/coh_ntau.png b/doc_build/savefig/coh_ntau.png
index 81e3d601..df6e0a87 100644
Binary files a/doc_build/savefig/coh_ntau.png and b/doc_build/savefig/coh_ntau.png differ
diff --git a/doc_build/savefig/coh_tau.png b/doc_build/savefig/coh_tau.png
index 51e490f2..e1540f93 100644
Binary files a/doc_build/savefig/coh_tau.png and b/doc_build/savefig/coh_tau.png differ
diff --git a/doc_build/savefig/gen_ar1_t0.png b/doc_build/savefig/gen_ar1_t0.png
index 2cf462cd..290e56cf 100644
Binary files a/doc_build/savefig/gen_ar1_t0.png and b/doc_build/savefig/gen_ar1_t0.png differ
diff --git a/doc_build/savefig/gen_ar1_t1.png b/doc_build/savefig/gen_ar1_t1.png
index 2f2a3d81..21d7f9a7 100644
Binary files a/doc_build/savefig/gen_ar1_t1.png and b/doc_build/savefig/gen_ar1_t1.png differ
diff --git a/doc_build/savefig/gen_ar1_t2.png b/doc_build/savefig/gen_ar1_t2.png
index 9791f931..684a19ad 100644
Binary files a/doc_build/savefig/gen_ar1_t2.png and b/doc_build/savefig/gen_ar1_t2.png differ
diff --git a/doc_build/savefig/gen_cn2_t0.png b/doc_build/savefig/gen_cn2_t0.png
index 58317b03..84779b77 100644
Binary files a/doc_build/savefig/gen_cn2_t0.png and b/doc_build/savefig/gen_cn2_t0.png differ
diff --git a/doc_build/savefig/gen_cn2_t1.png b/doc_build/savefig/gen_cn2_t1.png
index f3853cb8..561f8c8b 100644
Binary files a/doc_build/savefig/gen_cn2_t1.png and b/doc_build/savefig/gen_cn2_t1.png differ
diff --git a/doc_build/savefig/gen_cn_t0.png b/doc_build/savefig/gen_cn_t0.png
index 8e155c8d..9ac29f78 100644
Binary files a/doc_build/savefig/gen_cn_t0.png and b/doc_build/savefig/gen_cn_t0.png differ
diff --git a/doc_build/savefig/gen_cn_t1.png b/doc_build/savefig/gen_cn_t1.png
index eec29cca..3705d605 100644
Binary files a/doc_build/savefig/gen_cn_t1.png and b/doc_build/savefig/gen_cn_t1.png differ
diff --git a/doc_build/savefig/mapallarchive.png b/doc_build/savefig/mapallarchive.png
index 104d5497..36713155 100644
Binary files a/doc_build/savefig/mapallarchive.png and b/doc_build/savefig/mapallarchive.png differ
diff --git a/doc_build/savefig/mapallarchive_marker.png b/doc_build/savefig/mapallarchive_marker.png
index c1ef0ced..524ccf03 100644
Binary files a/doc_build/savefig/mapallarchive_marker.png and b/doc_build/savefig/mapallarchive_marker.png differ
diff --git a/doc_build/savefig/mapone.png b/doc_build/savefig/mapone.png
index 4c6ab456..0bbac54a 100644
Binary files a/doc_build/savefig/mapone.png and b/doc_build/savefig/mapone.png differ
diff --git a/doc_build/savefig/ms_ct.png b/doc_build/savefig/ms_ct.png
index 14e27ebb..40cb2989 100644
Binary files a/doc_build/savefig/ms_ct.png and b/doc_build/savefig/ms_ct.png differ
diff --git a/doc_build/savefig/mts_stackplot.png b/doc_build/savefig/mts_stackplot.png
index c8bf8844..bd266832 100644
Binary files a/doc_build/savefig/mts_stackplot.png and b/doc_build/savefig/mts_stackplot.png differ
diff --git a/doc_build/savefig/random_series.png b/doc_build/savefig/random_series.png
index 34518727..0976bb6c 100644
Binary files a/doc_build/savefig/random_series.png and b/doc_build/savefig/random_series.png differ
diff --git a/doc_build/savefig/scree_nmc.png b/doc_build/savefig/scree_nmc.png
index 928f27a1..2364b1a7 100644
Binary files a/doc_build/savefig/scree_nmc.png and b/doc_build/savefig/scree_nmc.png differ
diff --git a/doc_build/savefig/spec_ls.png b/doc_build/savefig/spec_ls.png
index 27291b6b..d1c6c0d0 100644
Binary files a/doc_build/savefig/spec_ls.png and b/doc_build/savefig/spec_ls.png differ
diff --git a/doc_build/savefig/spec_ls_freq.png b/doc_build/savefig/spec_ls_freq.png
index 077e03ea..07ea8117 100644
Binary files a/doc_build/savefig/spec_ls_freq.png and b/doc_build/savefig/spec_ls_freq.png differ
diff --git a/doc_build/savefig/spec_ls_n50.png b/doc_build/savefig/spec_ls_n50.png
index 2a2619e5..61bc57ed 100644
Binary files a/doc_build/savefig/spec_ls_n50.png and b/doc_build/savefig/spec_ls_n50.png differ
diff --git a/doc_build/savefig/spec_mtm.png b/doc_build/savefig/spec_mtm.png
index 87a21596..14043ec5 100644
Binary files a/doc_build/savefig/spec_mtm.png and b/doc_build/savefig/spec_mtm.png differ
diff --git a/doc_build/savefig/spec_perio.png b/doc_build/savefig/spec_perio.png
index e9116bc5..6e16af6b 100644
Binary files a/doc_build/savefig/spec_perio.png and b/doc_build/savefig/spec_perio.png differ
diff --git a/doc_build/savefig/spec_welch.png b/doc_build/savefig/spec_welch.png
index c7dbb416..a5f2030d 100644
Binary files a/doc_build/savefig/spec_welch.png and b/doc_build/savefig/spec_welch.png differ
diff --git a/doc_build/savefig/spec_wwz.png b/doc_build/savefig/spec_wwz.png
index dff510e2..638e132a 100644
Binary files a/doc_build/savefig/spec_wwz.png and b/doc_build/savefig/spec_wwz.png differ
diff --git a/doc_build/savefig/ssa_recon.png b/doc_build/savefig/ssa_recon.png
index 1147a5db..4534ff7b 100644
Binary files a/doc_build/savefig/ssa_recon.png and b/doc_build/savefig/ssa_recon.png differ
diff --git a/doc_build/savefig/ts_air.png b/doc_build/savefig/ts_air.png
index f45f62e8..9059f1c0 100644
Binary files a/doc_build/savefig/ts_air.png and b/doc_build/savefig/ts_air.png differ
diff --git a/doc_build/savefig/ts_dashboard.png b/doc_build/savefig/ts_dashboard.png
index c06c6eb5..595c66fd 100644
Binary files a/doc_build/savefig/ts_dashboard.png and b/doc_build/savefig/ts_dashboard.png differ
diff --git a/doc_build/savefig/ts_dist.png b/doc_build/savefig/ts_dist.png
index 38577fc3..b9a9d903 100644
Binary files a/doc_build/savefig/ts_dist.png and b/doc_build/savefig/ts_dist.png differ
diff --git a/doc_build/savefig/ts_filter1.png b/doc_build/savefig/ts_filter1.png
index 57abf2eb..d82b069b 100644
Binary files a/doc_build/savefig/ts_filter1.png and b/doc_build/savefig/ts_filter1.png differ
diff --git a/doc_build/savefig/ts_filter2.png b/doc_build/savefig/ts_filter2.png
index 246b0a72..622c131d 100644
Binary files a/doc_build/savefig/ts_filter2.png and b/doc_build/savefig/ts_filter2.png differ
diff --git a/doc_build/savefig/ts_filter3.png b/doc_build/savefig/ts_filter3.png
index 0d017b09..d4d88205 100644
Binary files a/doc_build/savefig/ts_filter3.png and b/doc_build/savefig/ts_filter3.png differ
diff --git a/doc_build/savefig/ts_filter4.png b/doc_build/savefig/ts_filter4.png
index 1ec1d0aa..3dc4f801 100644
Binary files a/doc_build/savefig/ts_filter4.png and b/doc_build/savefig/ts_filter4.png differ
diff --git a/doc_build/savefig/ts_filter5.png b/doc_build/savefig/ts_filter5.png
index fb1e070f..343dc6d4 100644
Binary files a/doc_build/savefig/ts_filter5.png and b/doc_build/savefig/ts_filter5.png differ
diff --git a/doc_build/savefig/ts_nino.png b/doc_build/savefig/ts_nino.png
index db748c5e..c080aea4 100644
Binary files a/doc_build/savefig/ts_nino.png and b/doc_build/savefig/ts_nino.png differ
diff --git a/doc_build/savefig/ts_plot.png b/doc_build/savefig/ts_plot.png
index 968d0578..e657dda1 100644
Binary files a/doc_build/savefig/ts_plot.png and b/doc_build/savefig/ts_plot.png differ
diff --git a/doc_build/savefig/ts_plot2.png b/doc_build/savefig/ts_plot2.png
index e1c2ee28..335ecfb2 100644
Binary files a/doc_build/savefig/ts_plot2.png and b/doc_build/savefig/ts_plot2.png differ
diff --git a/doc_build/savefig/ts_plot4.png b/doc_build/savefig/ts_plot4.png
index 968d0578..e657dda1 100644
Binary files a/doc_build/savefig/ts_plot4.png and b/doc_build/savefig/ts_plot4.png differ
diff --git a/doc_build/savefig/ts_plot5.png b/doc_build/savefig/ts_plot5.png
index 968d0578..e657dda1 100644
Binary files a/doc_build/savefig/ts_plot5.png and b/doc_build/savefig/ts_plot5.png differ
diff --git a/doc_build/savefig/ts_sg.png b/doc_build/savefig/ts_sg.png
index 6b6ee37c..4e917b27 100644
Binary files a/doc_build/savefig/ts_sg.png and b/doc_build/savefig/ts_sg.png differ
diff --git a/doc_build/savefig/wwa_wwz.png b/doc_build/savefig/wwa_wwz.png
index fdc61647..3f0e4a85 100644
Binary files a/doc_build/savefig/wwa_wwz.png and b/doc_build/savefig/wwa_wwz.png differ
diff --git a/doc_build/ts_plot3.png b/doc_build/ts_plot3.png
index 774b0f58..6d063d20 100644
Binary files a/doc_build/ts_plot3.png and b/doc_build/ts_plot3.png differ
diff --git a/doc_build/utils/introduction.rst b/doc_build/utils/introduction.rst
index 2e420bb1..d21f2cf2 100644
--- a/doc_build/utils/introduction.rst
+++ b/doc_build/utils/introduction.rst
@@ -144,8 +144,14 @@ Functions for wavelet analysis. Includes some pre-processing and post-processing
 
    * - :ref:`wwz <utils_wwz>`
      - Weighted wavelet amplitude (WWA) for unevenly-spaced data
-   * - :ref:`xwc <utils_xwc>`
-     - Cross-wavelet analysis for unevenly-spaced data.
+   * - :ref:`wavelet_coherence <utils_wavelet_coherence>`
+     - Cross-wavelet coherence for unevenly-spaced data.
+   * - :ref:`xwt <utils_xwt`
+     - Cross-Wavelet transform for unevenly-spaced data.
+   * - :ref:`wtc <utils_wtc`
+     - Cross-Wavelet coherence for unevenly-spaced data.
+   * - :ref:`cwt <utils_cwt`
+     - Continuous wavelet transform after Torrence and Compo
 
 
 Tsutils
@@ -162,19 +168,17 @@ This modules contain pre-processing functions for time series analysis.
    * - :ref:`bin <utils_bin>`
      - Bin the values into evenly-spaced bins
    * - :ref:`gkernel <utils_gkernel>`
-     - Coarsen time resolution using a Gaussian Kernel
-   * - :ref:`grid_properties <utils_grid_properties>`
-     - Establishes the grid properties of a numerical array
+     - Coarsens time resolution using a Gaussian Kernel
+   * - :ref:`increments <utils_increments>`
+     - Assesses the increments of a numerical array
    * - :ref:`interp <utils_interp>`
      - Interpolation function based on `scipy.interpolate.interp1d <https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html>_`
-   * - :ref:`on_common_axis <utils_on_common_axis>`
-     - Places two timeseries on a common time axis
    * - :ref:`standardize <utils_standardize>`
      - Standardizes a timeseries
    * - :ref:`ts2segments <utils_ts2segments>`
      - Chop a timeseries into several segments based on gap detection
    * - :ref:`annualize <utils_annualize>`
-     - Annualize a time series whose time resolution is finer than 1 year
+     - Annualizes a time series whose time resolution is finer than 1 year
    * - :ref:`gaussianize <utils_gaussianize>`
      - Maps a (proxy) timeseries to a Gaussian distribution
    * - :ref:`gaussianize_single <utils_gaussianize_single>`
@@ -182,9 +186,9 @@ This modules contain pre-processing functions for time series analysis.
    * - :ref:`detrend <utils_detrend>`
      - Applies linear, constant, low-pass filter, or decomposition-based detrending
    * - :ref:`detect_outliers <utils_detect_outliers>`
-     - Detect outliers in a timeseries
+     - Detects outliers in a timeseries
    * - :ref:`remove_outliers <utils_remove_outliers>`
-     - Remove outliers in a timeseries
+     - Removes outliers in a timeseries
 
 
 Tsbase
@@ -233,7 +237,7 @@ The most relevant functions concern querying the LinkedEarth wiki. The first 5 f
 jsonutils
 """""""""
 
-This module converts Pyleoclim objects to and from JSON files. Useful for obtaining a human-readable output and keeping the results of an analysis.
+This module converts Pyleoclim objects to and from JSON files. Useful for obtaining a human-readable output and keeping the results of an analysis. The JSON file can also be used to swap analysis results between programming language. Please note that this utilities are maintained on a as-needed basis and that not all objects are currently available.
 
 .. list-table::
    :widths: 30 70
@@ -243,7 +247,13 @@ This module converts Pyleoclim objects to and from JSON files. Useful for obtain
      - Saves a Pyleoclim Object (e.g., Series, PSD, Scalogram) to a json file
    * - :ref:`json_to_Series <utils_json_to_Series>`
      - Load a pyleoclim Series object from a JSON file
+   * - :ref:`json_to_LipdSeries <utils_json_to_LipdSeries>`
+     - Load a pyleoclim LipdSeries object from a JSON file
    * - :ref:`json_to_PSD <utils_json_to_PSD>`
      - Load a pyleoclim PSD object from a JSON file
    * - :ref:`json_to_Scalogram <utils_json_to_Scalogram>`
      - Load a pyleoclim Scalogram object from a JSON file
+   * - :ref:`json_to_Coherence <utils_json_to_Coherence>`
+     - Load a pyleoclim Coherence object from a JSON file
+   * - :ref:`json_to_MultiplePSD <utils_json_to_MultiplePSD>`
+     - Load a pyleoclim Multiple object from a JSON file
diff --git a/doc_build/utils/jsonutils/json_to_Coherence.rst b/doc_build/utils/jsonutils/json_to_Coherence.rst
new file mode 100644
index 00000000..5d6d520c
--- /dev/null
+++ b/doc_build/utils/jsonutils/json_to_Coherence.rst
@@ -0,0 +1,6 @@
+.. _utils_json_to_Coherence:
+
+json_to_Coherence (pyleoclim.utils.jsonutils.json_to_Coherence)
+===============================================================
+
+.. autofunction:: pyleoclim.utils.jsonutils.json_to_Coherence
diff --git a/doc_build/utils/jsonutils/json_to_LipdSeries.rst b/doc_build/utils/jsonutils/json_to_LipdSeries.rst
new file mode 100644
index 00000000..7e132d29
--- /dev/null
+++ b/doc_build/utils/jsonutils/json_to_LipdSeries.rst
@@ -0,0 +1,6 @@
+.. _utils_json_to_LipdSeries:
+
+json_to_LipdSeries (pyleoclim.utils.jsonutils.json_to_LipdSeries)
+=================================================================
+
+.. autofunction:: pyleoclim.utils.jsonutils.json_to_LipdSeries
diff --git a/doc_build/utils/jsonutils/json_to_MultiplePSD.rst b/doc_build/utils/jsonutils/json_to_MultiplePSD.rst
new file mode 100644
index 00000000..c02eff3b
--- /dev/null
+++ b/doc_build/utils/jsonutils/json_to_MultiplePSD.rst
@@ -0,0 +1,6 @@
+.. _utils_json_to_MultiplePSD:
+
+json_to_MultiplePSD (pyleoclim.utils.jsonutils.json_to_MultiplePSD)
+==================================================================
+
+.. autofunction:: pyleoclim.utils.jsonutils.json_to_MultiplePSD
diff --git a/doc_build/utils/tsutils/grid_properties.rst b/doc_build/utils/tsutils/grid_properties.rst
deleted file mode 100644
index d5032b99..00000000
--- a/doc_build/utils/tsutils/grid_properties.rst
+++ /dev/null
@@ -1,6 +0,0 @@
-.. _utils_grid_properties:
-
-grid_properties (pyleoclim.utils.tsutils.grid_properties)
-========================================================
-
-.. autofunction:: pyleoclim.utils.tsutils.grid_properties
diff --git a/doc_build/utils/tsutils/increments.rst b/doc_build/utils/tsutils/increments.rst
new file mode 100644
index 00000000..537897eb
--- /dev/null
+++ b/doc_build/utils/tsutils/increments.rst
@@ -0,0 +1,6 @@
+.. _utils_increments:
+
+increments (pyleoclim.utils.tsutils.increments)
+========================================================
+
+.. autofunction:: pyleoclim.utils.tsutils.increments
diff --git a/doc_build/utils/tsutils/on_common_axis.rst b/doc_build/utils/tsutils/on_common_axis.rst
deleted file mode 100644
index 77dd20ce..00000000
--- a/doc_build/utils/tsutils/on_common_axis.rst
+++ /dev/null
@@ -1,6 +0,0 @@
-.. _utils_on_common_axis:
-
-on_common_axis (pyleoclim.utils.tsutils.on_common_axis)
-===============================================
-
-.. autofunction:: pyleoclim.utils.tsutils.on_common_axis
diff --git a/doc_build/utils/wavelet/chisquare_inv.rst b/doc_build/utils/wavelet/chisquare_inv.rst
new file mode 100644
index 00000000..c9904ea0
--- /dev/null
+++ b/doc_build/utils/wavelet/chisquare_inv.rst
@@ -0,0 +1,6 @@
+.. _utils_chisquare_inv:
+
+chisquare_inv (pyleoclim.utils.wavelet.chisquare_inv)
+=====================================================
+
+.. autofunction:: pyleoclim.utils.wavelet.chisquare_inv
diff --git a/doc_build/utils/wavelet/chisquare_solve.rst b/doc_build/utils/wavelet/chisquare_solve.rst
new file mode 100644
index 00000000..ba59e72a
--- /dev/null
+++ b/doc_build/utils/wavelet/chisquare_solve.rst
@@ -0,0 +1,6 @@
+.. _utils_chisquare_solve:
+
+chisquare_solve (pyleoclim.utils.wavelet.chisquare_solve)
+========================================================
+
+.. autofunction:: pyleoclim.utils.wavelet.chisquare_solve
diff --git a/doc_build/utils/wavelet/clean_ts.rst b/doc_build/utils/wavelet/clean_ts.rst
new file mode 100644
index 00000000..c18b7253
--- /dev/null
+++ b/doc_build/utils/wavelet/clean_ts.rst
@@ -0,0 +1,6 @@
+.. _wavelet_clean_ts:
+
+clean_ts (pyleoclim.utils.wavelet.clean_ts)
+===========================================
+
+.. autofunction:: pyleoclim.utils.wavelet.clean_ts
diff --git a/doc_build/utils/wavelet/cross_wt.rst b/doc_build/utils/wavelet/cross_wt.rst
deleted file mode 100644
index 053b5e22..00000000
--- a/doc_build/utils/wavelet/cross_wt.rst
+++ /dev/null
@@ -1,6 +0,0 @@
-.. _wavelet_cross_wt:
-
-cross_wt (pyleoclim.utils.wavelet.cross_wt)
-===========================================
-
-.. autofunction:: pyleoclim.utils.wavelet.cross_wt
diff --git a/doc_build/utils/wavelet/cwt.rst b/doc_build/utils/wavelet/cwt.rst
new file mode 100644
index 00000000..d2a6ad4d
--- /dev/null
+++ b/doc_build/utils/wavelet/cwt.rst
@@ -0,0 +1,6 @@
+.. _utils_cwt:
+
+cwt (pyleoclim.utils.wavelet.cwt)
+=================================
+
+.. autofunction:: pyleoclim.utils.wavelet.cwt
diff --git a/doc_build/utils/wavelet/is_evenly_spaced.rst b/doc_build/utils/wavelet/is_evenly_spaced.rst
new file mode 100644
index 00000000..630a0d8a
--- /dev/null
+++ b/doc_build/utils/wavelet/is_evenly_spaced.rst
@@ -0,0 +1,6 @@
+.. _wavelet_is_evenly_spaced:
+
+is_evenly_spaced (pyleoclim.utils.wavelet.is_evenly_spaced)
+==========================================================
+
+.. autofunction:: pyleoclim.utils.wavelet.is_evenly_spaced
diff --git a/doc_build/utils/wavelet/tc_wave_bases.rst b/doc_build/utils/wavelet/tc_wave_bases.rst
new file mode 100644
index 00000000..62e30ec0
--- /dev/null
+++ b/doc_build/utils/wavelet/tc_wave_bases.rst
@@ -0,0 +1,6 @@
+.. _utils_tc_wave_bases:
+
+tc_wave_bases (pyleoclim.utils.wavelet.tc_wave_bases)
+=====================================================
+
+.. autofunction:: pyleoclim.utils.wavelet.tc_wave_bases
diff --git a/doc_build/utils/wavelet/tc_wave_signif.rst b/doc_build/utils/wavelet/tc_wave_signif.rst
new file mode 100644
index 00000000..513d8259
--- /dev/null
+++ b/doc_build/utils/wavelet/tc_wave_signif.rst
@@ -0,0 +1,6 @@
+.. _utils_tc_wave_signif:
+
+tc_wave_signif (pyleoclim.utils.wavelet.tc_wave_signif)
+=====================================================
+
+.. autofunction:: pyleoclim.utils.wavelet.tc_wave_signif
diff --git a/doc_build/utils/wavelet/tc_wavelet.rst b/doc_build/utils/wavelet/tc_wavelet.rst
new file mode 100644
index 00000000..2be29303
--- /dev/null
+++ b/doc_build/utils/wavelet/tc_wavelet.rst
@@ -0,0 +1,6 @@
+.. _utils_tc_wavelet:
+
+tc_wavelet (pyleoclim.utils.wavelet.tc_wavelet)
+===============================================
+
+.. autofunction:: pyleoclim.utils.wavelet.tc_wavelet
diff --git a/doc_build/utils/wavelet/wtc.rst b/doc_build/utils/wavelet/wtc.rst
new file mode 100644
index 00000000..48246d1e
--- /dev/null
+++ b/doc_build/utils/wavelet/wtc.rst
@@ -0,0 +1,6 @@
+.. _utils_wtc:
+
+wtc (pyleoclim.utils.wavelet.wtc)
+==================================
+
+.. autofunction:: pyleoclim.utils.wavelet.wtc
diff --git a/doc_build/utils/wavelet/xwt.rst b/doc_build/utils/wavelet/xwt.rst
new file mode 100644
index 00000000..d86ed81c
--- /dev/null
+++ b/doc_build/utils/wavelet/xwt.rst
@@ -0,0 +1,6 @@
+.. _wavelet_xwt:
+
+xwt (pyleoclim.utils.wavelet.xwt)
+===========================================
+
+.. autofunction:: pyleoclim.utils.wavelet.xwt
diff --git a/pyleoclim/core/ui.py b/pyleoclim/core/ui.py
index 21580e51..0b968b2b 100644
--- a/pyleoclim/core/ui.py
+++ b/pyleoclim/core/ui.py
@@ -32,7 +32,6 @@
 import cartopy.crs as ccrs
 import cartopy.feature as cfeature
 
-
 from tqdm import tqdm
 from scipy.stats.mstats import mquantiles
 from scipy import stats
@@ -42,6 +41,8 @@
 
 import lipd as lpd
 
+import collections
+
 def pval_format(p, threshold=0.01, style='exp'):
     ''' Print p-value with proper format when p is close to 0
     '''
@@ -1878,7 +1879,7 @@ def spectral(self, method='lomb_scargle', freq_method='log', freq_kwargs=None, s
         ----------
 
         method : str
-            {'wwz', 'mtm', 'lomb_scargle', 'welch', 'periodogram'}
+            {'wwz', 'mtm', 'lomb_scargle', 'welch', 'periodogram', 'cwt'}
 
         freq_method : str
             {'log','scale', 'nfft', 'lomb_scargle', 'welch'}
@@ -1893,7 +1894,7 @@ def spectral(self, method='lomb_scargle', freq_method='log', freq_kwargs=None, s
             Label for the PSD object
 
         scalogram : pyleoclim.core.ui.Series.Scalogram
-            The return of the wavelet analysis; effective only when the method is 'wwz'
+            The return of the wavelet analysis; effective only when the method is 'wwz' or 'cwt'
 
         verbose : bool
             If True, will print warning messages if there is any
@@ -1916,6 +1917,8 @@ def spectral(self, method='lomb_scargle', freq_method='log', freq_kwargs=None, s
 
         pyleoclim.utils.spectral.wwz_psd : Spectral analysis using the Wavelet Weighted Z transform
 
+        pyleoclim.utils.spectral.cwt_psd : Spectral analysis using the continuous Wavelet Transform as implemented by Torrence and Compo
+
         pyleoclim.utils.wavelet.make_freq_vector : Functions to create the frequency vector
 
         pyleoclim.utils.tsutils.detrend : Detrending function
@@ -2020,7 +2023,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) #in practice, need more AR1 simulations
+            psd_perio_signif = psd_perio.signif_test(number=20, method='ar1sim') #in practice, need more AR1 simulations
             @savefig spec_perio.png
             fig, ax = psd_perio_signif.plot(title='PSD using Periodogram method')
             pyleo.closefig(fig)
@@ -2033,7 +2036,7 @@ def spectral(self, method='lomb_scargle', freq_method='log', freq_kwargs=None, s
 
             ts_interp = ts_std.interp()
             psd_welch = ts_interp.spectral(method='welch')
-            psd_welch_signif = psd_welch.signif_test(number=20) #in practice, need more AR1 simulations
+            psd_welch_signif = psd_welch.signif_test(number=20, method='ar1sim') #in practice, need more AR1 simulations
             @savefig spec_welch.png
             fig, ax = psd_welch_signif.plot(title='PSD using Welch method')
             pyleo.closefig(fig)
@@ -2046,11 +2049,24 @@ def spectral(self, method='lomb_scargle', freq_method='log', freq_kwargs=None, s
 
             ts_interp = ts_std.interp()
             psd_mtm = ts_interp.spectral(method='mtm')
-            psd_mtm_signif = psd_mtm.signif_test(number=20) #in practice, need more AR1 simulations
+            psd_mtm_signif = psd_mtm.signif_test(number=20, method='ar1sim') #in practice, need more AR1 simulations
             @savefig spec_mtm.png
             fig, ax = psd_mtm_signif.plot(title='PSD using MTM method')
             pyleo.closefig(fig)
 
+        - Continuous Wavelet Transform
+
+        .. ipython:: python
+            :okwarning:
+            :okexcept:
+
+            ts_interp = ts_std.interp()
+            psd_cwt = ts_interp.spectral(method='cwt')
+            psd_cwt_signif = psd_cwt.signif_test()
+            @savefig spec_cwt.png
+            fig, ax = psd_cwt_signif.plot(title='PSD using CWT method')
+            pyleo.closefig(fig)
+
         '''
         if not verbose:
             warnings.simplefilter('ignore')
@@ -2061,13 +2077,15 @@ def spectral(self, method='lomb_scargle', freq_method='log', freq_kwargs=None, s
             'mtm': specutils.mtm,
             'lomb_scargle': specutils.lomb_scargle,
             'welch': specutils.welch,
-            'periodogram': specutils.periodogram
+            'periodogram': specutils.periodogram,
+            'cwt':specutils.cwt_psd
         }
         args = {}
         freq_kwargs = {} if freq_kwargs is None else freq_kwargs.copy()
         freq = waveutils.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'] = {}
@@ -2083,6 +2101,16 @@ def spectral(self, method='lomb_scargle', freq_method='log', freq_kwargs=None, s
                 }
             )
 
+        if method == 'cwt' and scalogram is not None:
+            Results = collections.namedtuple('Results', ['amplitude', 'coi', 'freq', 'time', 'scale', 'mother','param'])
+            res = Results(amplitude=scalogram.amplitude, coi=scalogram.coi,
+                          freq=scalogram.frequency, time=scalogram.time,
+                          scale=scalogram.wave_args['scale'],
+                          mother=scalogram.wave_args['mother'],
+                          param=scalogram.wave_args['param'])
+            args['cwt'].update({'cwt_res':res})
+
+
         spec_res = spec_func[method](self.value, self.time, **args[method])
         if type(spec_res) is dict:
             spec_res = dict2namedtuple(spec_res)
@@ -2095,6 +2123,12 @@ def spectral(self, method='lomb_scargle', freq_method='log', freq_kwargs=None, s
             args['wwz'].pop('wwz_Neffs')
             args['wwz'].pop('wwz_freq')
 
+        if method == 'cwt':
+            args['cwt'].update({'scale':spec_res.scale,'mother':spec_res.mother,'param':spec_res.param})
+            if scalogram is not None:
+                args['cwt'].pop('cwt_res')
+
+
         psd = PSD(
             frequency=spec_res.freq,
             amplitude=spec_res.psd,
@@ -2109,13 +2143,11 @@ def spectral(self, method='lomb_scargle', freq_method='log', freq_kwargs=None, s
     def wavelet(self, method='wwz', settings=None, freq_method='log', ntau=None, freq_kwargs=None, verbose=False):
         ''' Perform wavelet analysis on the timeseries
 
-        cwt wavelets documented on https://pywavelets.readthedocs.io/en/latest/ref/cwt.html
-
         Parameters
         ----------
 
         method : {wwz, cwt}
-            Whether to use the wwz method for unevenly spaced timeseries or traditional cwt (from pywavelets)
+            Whether to use the wwz method for unevenly spaced timeseries or traditional cwt (from Torrence and Compo)
 
         freq_method : str
             {'log', 'scale', 'nfft', 'lomb_scargle', 'welch'}
@@ -2124,11 +2156,12 @@ def wavelet(self, method='wwz', settings=None, freq_method='log', ntau=None, fre
             Arguments for frequency vector
 
         ntau : int
-            The length of the time shift points that determins the temporal resolution of the result.
+            The length of the time shift points that determines the temporal resolution of the result.
             If None, it will be either the length of the input time axis, or at most 50.
+            Not relevant for cwt
 
         settings : dict
-            Arguments for the specific spectral method
+            Arguments for the specific wavelet method
 
         verbose : bool
             If True, will print warning messages if there is any
@@ -2153,6 +2186,12 @@ def wavelet(self, method='wwz', settings=None, freq_method='log', ntau=None, fre
 
         pyleoclim.core.ui.MultipleScalogram : Multiple Scalogram object
 
+        References
+        ----------
+
+        Torrence, C. and G. P. Compo, 1998: A Practical Guide to Wavelet Analysis. Bull. Amer. Meteor. Soc., 79, 61-78.
+        Python routines available at http://paos.colorado.edu/research/wavelets/
+
         Examples
         --------
 
@@ -2175,21 +2214,18 @@ def wavelet(self, method='wwz', settings=None, freq_method='log', ntau=None, fre
             fig, ax = scal_signif.plot()
             pyleo.closefig(fig)
 
+        See the Series.spectral() functionality as an example to change the wavelet method.
+
         '''
         if not verbose:
             warnings.simplefilter('ignore')
 
         settings = {} if settings is None else settings.copy()
         wave_func = {
-            'wwz': waveutils.wwz
-            # 'cwt': waveutils.cwt,
+            'wwz': waveutils.wwz,
+            'cwt': waveutils.cwt
         }
 
-        if method == 'cwt' and 'freq' in settings.keys():
-            scales=1/np.array(settings['freq'])
-            settings.update({'scales':scales})
-            del settings['freq']
-
         freq_kwargs = {} if freq_kwargs is None else freq_kwargs.copy()
         freq = waveutils.make_freq_vector(self.time, method=freq_method, **freq_kwargs)
 
@@ -2201,15 +2237,15 @@ def wavelet(self, method='wwz', settings=None, freq_method='log', ntau=None, fre
         tau = np.linspace(np.min(self.time), np.max(self.time), ntau)
 
         args['wwz'] = {'tau': tau, 'freq': freq}
-        args['cwt'] = {'wavelet' : 'morl', 'scales':1/freq}
-
+        args['cwt'] = {'freq':freq}
 
         args[method].update(settings)
         wave_res = wave_func[method](self.value, self.time, **args[method])
         if method == 'wwz':
             wwz_Neffs = wave_res.Neffs
-        else:
+        elif method=='cwt':
             wwz_Neffs = None
+            args[method].update({'scale':wave_res.scale,'mother':wave_res.mother,'param':wave_res.param})
 
         scal = Scalogram(
             frequency=wave_res.freq,
@@ -2245,11 +2281,11 @@ def wavelet_coherence(self, target_series, method='wwz', settings=None, freq_met
             Arguments for frequency vector
 
         tau : array
-            The time shift points that determins the temporal resolution of the result.
+            The time shift points that determines the temporal resolution of the result.
             If None, it will be calculated using ntau.
 
         ntau : int
-            The length of the time shift points that determins the temporal resolution of the result.
+            The length of the time shift points that determines the temporal resolution of the result.
             If None, it will be either the length of the input time axis, or at most 50.
 
         settings : dict
@@ -2329,7 +2365,7 @@ def wavelet_coherence(self, target_series, method='wwz', settings=None, freq_met
 
         settings = {} if settings is None else settings.copy()
         xwc_func = {
-            'wwz': waveutils.xwc,
+            'wwz': waveutils.wavelet_coherence,
         }
 
         freq_kwargs = {} if freq_kwargs is None else freq_kwargs.copy()
@@ -2578,13 +2614,13 @@ def causality(self, target_series, method='liang', settings=None):
         causal_res = spec_func[method](sorted_self.value, sorted_target.value, **args[method])
         return causal_res
 
-    def surrogates(self, method='ar1', number=1, length=None, seed=None, settings=None):
+    def surrogates(self, method='ar1sim', number=1, length=None, seed=None, settings=None):
         ''' Generate surrogates with increasing time axis
 
         Parameters
         ----------
 
-        method : {ar1}
+        method : {ar1sim}
             Uses an AR1 model to generate surrogates of the timeseries
 
         number : int
@@ -2610,10 +2646,10 @@ def surrogates(self, method='ar1', number=1, length=None, seed=None, settings=No
         '''
         settings = {} if settings is None else settings.copy()
         surrogate_func = {
-            'ar1': tsmodel.ar1_sim,
+            'ar1sim': tsmodel.ar1_sim,
         }
         args = {}
-        args['ar1'] = {'t': self.time}
+        args['ar1sim'] = {'t': self.time}
         args[method].update(settings)
 
         if seed is not None:
@@ -2636,7 +2672,7 @@ def outliers(self, auto=True, remove=True, fig_outliers=True,fig_knee=True,
                  plot_outliers_kwargs=None,plot_knee_kwargs=None,figsize=[10,4],
                  saveknee_settings=None,saveoutliers_settings=None, mute=False):
         '''
-        Detects outliers in a timeseries and removes if specified
+        Detects outliers in a timeseries and removes if specified. The method uses clustering to locate outliers.
 
         Parameters
         ----------
@@ -2678,11 +2714,49 @@ def outliers(self, auto=True, remove=True, fig_outliers=True,fig_knee=True,
 
         pyleoclim.utils.plotting.plot_scatter_xy : Scatter plot on top of a line plot
 
+        Examples
+        --------
+
+        Let's create a "perfect" sinusoidal signal and add outliers
+
+        .. ipython:: python
+            :okwarning:
+            :okexcept:
+
+            import pyleoclim as pyleo
+            import numpy as np
+
+            # create the signal
+            freqs=[1/20,1/80]
+            time=np.arange(2001)
+            signals=[]
+            for freq in freqs:
+                signals.append(np.cos(2*np.pi*freq*time))
+            signal=sum(signals)
+
+            #add outliers
+            outliers_start = np.mean(signal)+5*np.std(signal)
+            outliers_end = np.mean(signal)+7*np.std(signal)
+            outlier_values = np.arange(outliers_start,outliers_end,0.1)
+            index = np.random.randint(0,len(signal),6)
+            signal_out = signal
+            for i,ind in enumerate(index):
+                signal_out[ind] = outlier_values[i]
+
+            #Make a Series object
+            ts = pyleo.Series(time=time,value=signal_out)
+            @savefig outliers.png
+            fig, ax = ts.plot()
+            pyleo.closefig(fig)
+
+            #Detect and remove outliers
+            ts_new=ts.outliers()
+            @savefig outliers_remove.png
+            fig, ax = ts_new.plot()
+            pyleo.closefig(fig)
+
         '''
         new = self.copy()
-
-        #outlier_indices,fig1,ax1,fig2,ax2 = tsutils.detect_outliers(self.time, self.value, auto=auto, plot_knee=fig_knee,plot_outliers=fig_outliers,\
-        #                                                   figsize=figsize,save_knee=save_knee,save_outliers=save_outliers,plot_outliers_kwargs=plot_outliers_kwargs,plot_knee_kwargs=plot_knee_kwargs)
         outlier_indices = tsutils.detect_outliers(
             self.time, self.value, auto=auto, plot_knee=fig_knee,plot_outliers=fig_outliers,
             figsize=figsize,saveknee_settings=saveknee_settings,saveoutliers_settings=saveoutliers_settings,
@@ -2826,7 +2900,7 @@ def __str__(self):
         msg = print(tabulate(table, headers='keys'))
         return f'Length: {np.size(self.frequency)}'
 
-    def signif_test(self, number=None, method='ar1', seed=None, qs=[0.95],
+    def signif_test(self, method='ar1sim', number=None, seed=None, qs=[0.95],
                     settings=None, scalogram = None):
         '''
 
@@ -2835,14 +2909,14 @@ def signif_test(self, number=None, method='ar1', seed=None, qs=[0.95],
         ----------
         number : int, optional
             Number of surrogate series to generate for significance testing. The default is None.
-        method : {ar1}, optional
-            Method to generate surrogates. The default is 'ar1'.
+        method : {ar1asym,'ar1sim'}
+            Method to generate surrogates. AR1sim uses simulated timeseries with similar persistence. AR1asymp represents the closed form solution. The default is AR1sim
         seed : int, optional
             Option to set the seed for reproducibility. The default is None.
         qs : list, optional
             Singificance levels to return. The default is [0.95].
         settings : dict, optional
-            Parameters. The default is None.
+            Parameters for the specific significance test. The default is None. Note that the default value for the asymptotic solution is `time-average`
         scalogram : Pyleoclim Scalogram object, optional
             Scalogram containing signif_scals exported during significance testing of scalogram.
             If number is None and signif_scals are present, will use length of scalogram list as number of significance tests
@@ -2879,36 +2953,94 @@ def signif_test(self, number=None, method='ar1', seed=None, qs=[0.95],
 
             pyleo.closefig(fig)
 
+        See also
+        --------
+
+        pyleoclim.utils.wavelet.tc_wave_signif : asymptotic significance calculation
+
         '''
 
-        signif_scals = None
-        if scalogram:
-            try:
-                signif_scals = scalogram.signif_scals
-            except:
-                return ValueError('Could not find signif_scals in passed object, make sure this is a scalogram with signif_scals that were saved during significance testing')
+        if self.spec_method == 'wwz' 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 number is None and signif_scals:
-            number = len(signif_scals.scalogram_list)
-        elif number is None and signif_scals is None:
-            number = 200
-        elif number == 0:
-            return self
+        if method not in ['ar1sim', 'ar1asym']:
+                raise ValueError("The available methods are ar1sim'and 'ar1asym'")
+
+        if method == 'ar1sim':
+            signif_scals = None
+            if scalogram:
+                try:
+                    signif_scals = scalogram.signif_scals
+                except:
+                    return ValueError('Could not find signif_scals in passed object, make sure this is a scalogram with signif_scals that were saved during significance testing')
 
-        new = self.copy()
-        surr = self.timeseries.surrogates(
-            number=number, seed=seed, method=method, settings=settings
-        )
 
-        if signif_scals:
-            surr_psd = surr.spectral(
-                method=self.spec_method, settings=self.spec_args, scalogram_list=signif_scals
+            if number is None and signif_scals:
+                number = len(signif_scals.scalogram_list)
+            elif number is None and signif_scals is None:
+                number = 200
+            elif number == 0:
+                return self
+
+            new = self.copy()
+            surr = self.timeseries.surrogates(
+                number=number, seed=seed, method=method, settings=settings
             )
-        else:
-            surr_psd = surr.spectral(method=self.spec_method, settings=self.spec_args)
-        new.signif_qs = surr_psd.quantiles(qs=qs)
-        new.signif_method = method
+
+            if signif_scals:
+                surr_psd = surr.spectral(
+                    method=self.spec_method, settings=self.spec_args, scalogram_list=signif_scals
+                )
+            else:
+                surr_psd = surr.spectral(method=self.spec_method, settings=self.spec_args)
+            new.signif_qs = surr_psd.quantiles(qs=qs)
+            new.signif_method = method
+
+        elif method == 'ar1asym':
+            new=self.copy()
+
+            if type(qs) is not list:
+                raise TypeError('qs should be a list')
+
+            settings = {'sigtest':'time-average'} if settings is None else settings.copy()
+
+            if self.spec_method=='cwt':
+                if 'dof' not in settings.keys():
+                    dof = len(self.timeseries.value) - self.spec_args['scale']
+                    settings.update({'dof':dof})
+                signif_levels=waveutils.tc_wave_signif(self.timeseries.value,
+                                                       self.timeseries.time,
+                                                       self.spec_args['scale'],
+                                                       self.spec_args['mother'],
+                                                       self.spec_args['param'],
+                                                       qs=qs, **settings)
+            else:
+                # hard code Mortlet values to obtain the spectrum
+                param = 6
+                fourier_factor = 4 * np.pi / (param + np.sqrt(2 + param**2))
+                scale = 1/(fourier_factor*self.frequency)
+                if 'dof' not in settings.keys():
+                    dof = len(self.timeseries.value) - scale
+                    settings.update({'dof':dof})
+                signif_levels=waveutils.tc_wave_signif(self.timeseries.value,
+                                                       self.timeseries.time,
+                                                       scale,
+                                                       'MORLET',
+                                                       param,
+                                                       qs=qs, **settings)
+
+            # get it back into the object
+            new.signif_method = method
+
+            ms_base = []
+            for idx, item in enumerate(signif_levels):
+                label = str(int(qs[idx]*100))+'%'
+                s = PSD(frequency=self.frequency, amplitude = item, label=label)
+                ms_base.append(s)
+                new.signif_qs = MultiplePSD(ms_base)
 
         return new
 
@@ -3148,7 +3280,8 @@ def plot(self, in_loglog=True, in_period=True, label=None, xlabel=None, ylabel='
         # plot significance levels
         if self.signif_qs is not None:
             signif_method_label = {
-                'ar1': 'AR(1)',
+                'ar1sim': 'AR(1) simulations',
+                'ar1asym': 'AR(1) asymptotic solution'
             }
             nqs = np.size(self.signif_qs.psd_list)
 
@@ -3317,7 +3450,7 @@ def __str__(self):
         msg = print(tabulate(table, headers='keys'))
         return f'Dimension: {np.size(self.frequency)} x {np.size(self.time)}'
 
-    def plot(self, in_period=True, xlabel=None, ylabel=None, title=None,
+    def plot(self, variable = 'amplitude', in_period=True, xlabel=None, ylabel=None, title=None,
              ylim=None, xlim=None, yticks=None, figsize=[10, 8], mute=False,
              signif_clr='white', signif_linestyles='-', signif_linewidths=1,
              contourf_style={}, cbar_style={}, savefig_settings={}, ax=None):
@@ -3325,6 +3458,8 @@ def plot(self, in_period=True, xlabel=None, ylabel=None, title=None,
 
         Parameters
         ----------
+        variable : {'amplitude','power'}
+            Whether to plot the amplitude or power. Default is amplitude
         in_period : bool, optional
             Plot the in period instead of frequency space. The default is True.
         xlabel : str, optional
@@ -3391,7 +3526,12 @@ def plot(self, in_period=True, xlabel=None, ylabel=None, title=None,
             if ylabel is None:
                 ylabel = f'Frequency [1/{self.period_unit}]' if self.period_unit is not None else 'Frequency'
 
-        cont = ax.contourf(self.time, y_axis, self.amplitude.T, **contourf_args)
+        if variable == 'amplitude':
+            cont = ax.contourf(self.time, y_axis, self.amplitude.T, **contourf_args)
+        elif variable=='power':
+            cont = ax.contourf(self.time, y_axis, self.amplitude.T**2, **contourf_args)
+        else:
+            raise ValueError('Variable should be either "amplitude" or "power"')
         ax.set_yscale('log')
 
         # plot colorbar
@@ -3455,16 +3595,15 @@ def plot(self, in_period=True, xlabel=None, ylabel=None, title=None,
         else:
             return ax
 
-    def signif_test(self, number=None, method='ar1', seed=None, qs=[0.95],
+    def signif_test(self, method='ar1sim', number=None, seed=None, qs=[0.95],
                     settings=None, export_scal = False):
         '''Significance test for wavelet analysis
 
         Parameters
         ----------
-        number : int, optional
-            Number of surrogates to generate for significance analysis. The default is 200.
-        method : {'ar1'}, optional
-            Method to use to generate the surrogates. The default is 'ar1'.
+        method : {'ar1asym', 'ar1sim'}
+            Method to use to generate the surrogates.  AR1sim uses simulated timeseries with similar persistence. AR1asymp represents the closed form solution.The default is AR1sim
+            Number of surrogates to generate for significance analysis based on simulations. The default is 200.
         seed : int, optional
             Set the seed for the random number generator. Useful for reproducibility The default is None.
         qs : list, optional
@@ -3472,7 +3611,7 @@ def signif_test(self, number=None, method='ar1', seed=None, qs=[0.95],
         settings : dict, optional
             Parameters for the model. The default is None.
         export_scal : bool
-            Whether or not to export the scalograms used in the noise realizations. Note: The scalograms used for wavelet analysis are slightly different
+            Whether or not to export the scalograms used in the noise realizations. Note: For the wwz method, the scalograms used for wavelet analysis are slightly different
             than those used for spectral analysis (different decay constant). As such, this functionality should be used only to expedite exploratory analysis.
 
         Raises
@@ -3507,53 +3646,91 @@ def signif_test(self, number=None, method='ar1', seed=None, qs=[0.95],
 
         pyleoclim.core.ui.Series.wavelet : wavelet analysis
 
+        pyleoclim.utils.wavelet.tc_wave_signif : asymptotic significance calculation
+
         '''
+        if self.wave_method == 'wwz' and method == 'ar1asym':
+            raise ValueError('Asymptotic solution is not supported for the wwz method')
 
-        if hasattr(self,'signif_scals'):
-            signif_scals = self.signif_scals
+        if method not in ['ar1sim', 'ar1asym']:
+                raise ValueError("The available methods are ar1sim'and 'ar1asym'")
 
-        #Allow for a few different configurations of passed number of signif tests, default behavior is to set number = 200
-        if number is None and signif_scals is not None:
-            number = len(signif_scals.scalogram_list)
-        elif number is None and signif_scals is None:
-            number = 200
-        elif number == 0:
-            return self
+        if method == 'ar1sim':
 
-        new = self.copy()
+            if hasattr(self,'signif_scals'):
+                signif_scals = self.signif_scals
+
+            #Allow for a few different configurations of passed number of signif tests, default behavior is to set number = 200
+            if number is None and signif_scals is not None:
+                number = len(signif_scals.scalogram_list)
+            elif number is None and signif_scals is None:
+                number = 200
+            elif number == 0:
+                return self
+
+            new = self.copy()
 
-        if signif_scals:
-            scalogram_list = signif_scals.scalogram_list
-            #If signif_scals already in scalogram object are more than those requested for significance testing, use as many of them as required
-            if len(scalogram_list) > number:
-                surr_scal = MultipleScalogram(scalogram_list=scalogram_list[:number])
-            #If number is the same as the length of signif_scals, just use signif_scals
-            elif len(scalogram_list) == number:
-                surr_scal = signif_scals
-            #If the number is more than the length of signif_scals, reuse what is available and calculate the rest
-            elif len(scalogram_list) < number:
-                number -= len(scalogram_list)
-                surr_scal_tmp = []
-                surr_scal_tmp.extend(scalogram_list)
+            if signif_scals:
+                scalogram_list = signif_scals.scalogram_list
+                #If signif_scals already in scalogram object are more than those requested for significance testing, use as many of them as required
+                if len(scalogram_list) > number:
+                    surr_scal = MultipleScalogram(scalogram_list=scalogram_list[:number])
+                #If number is the same as the length of signif_scals, just use signif_scals
+                elif len(scalogram_list) == number:
+                    surr_scal = signif_scals
+                #If the number is more than the length of signif_scals, reuse what is available and calculate the rest
+                elif len(scalogram_list) < number:
+                    number -= len(scalogram_list)
+                    surr_scal_tmp = []
+                    surr_scal_tmp.extend(scalogram_list)
+                    surr = self.timeseries.surrogates(
+                number=number, seed=seed, method=method, settings=settings
+            )
+                    surr_scal_tmp.extend(surr.wavelet(method=self.wave_method, settings=self.wave_args,).scalogram_list)
+                    surr_scal = MultipleScalogram(scalogram_list=surr_scal_tmp)
+            else:
                 surr = self.timeseries.surrogates(
-            number=number, seed=seed, method=method, settings=settings
-        )
-                surr_scal_tmp.extend(surr.wavelet(method=self.wave_method, settings=self.wave_args,).scalogram_list)
-                surr_scal = MultipleScalogram(scalogram_list=surr_scal_tmp)
-        else:
-            surr = self.timeseries.surrogates(
-            number=number, seed=seed, method=method, settings=settings
-        )
-            surr_scal = surr.wavelet(method=self.wave_method, settings=self.wave_args,)
+                number=number, seed=seed, method=method, settings=settings
+            )
+                surr_scal = surr.wavelet(method=self.wave_method, settings=self.wave_args,)
 
-        if len(qs) > 1:
-            raise ValueError('qs should be a list with size 1!')
+            if type(qs) is not list:
+                raise TypeError('qs should be a list')
 
-        new.signif_qs = surr_scal.quantiles(qs=qs)
-        new.signif_method = method
+            new.signif_qs = surr_scal.quantiles(qs=qs)
+            new.signif_method = method
 
-        if export_scal == True:
-            new.signif_scals = surr_scal
+            if export_scal == True:
+                new.signif_scals = surr_scal
+
+        elif method == 'ar1asym':
+
+            new = self.copy()
+
+            if type(qs) is not list:
+                raise TypeError('qs should be a list')
+
+            settings = {} if settings is None else settings.copy()
+
+            signif_levels=waveutils.tc_wave_signif(self.timeseries.value,
+                                                   self.timeseries.time,
+                                                   self.wave_args['scale'],
+                                                   self.wave_args['mother'],
+                                                   self.wave_args['param'],
+                                                   qs=qs, **settings)
+
+            #Create a scalogram for each of the significance levels
+            ms_base =[]
+            for idx, item in enumerate(signif_levels):
+                label = str(int(qs[idx]*100))+'%'
+                # expand
+                signif = item[:, np.newaxis].dot(np.ones(len(self.timeseries.value))[np.newaxis, :])
+                s = Scalogram(frequency=self.frequency, time =self.time,
+                              amplitude = signif.T, label=label)
+                ms_base.append(s)
+
+            new.signif_qs = MultipleScalogram(ms_base)
+            new.signif_method = method
 
         return new
 
@@ -4081,7 +4258,7 @@ def standardize(self):
             ms.series_list[idx]=s
         return ms
 
-    def grid_properties(self, step_style='median'):
+    def increments(self, step_style='median'):
         '''
         Extract grid properties (start, stop, step) of all the Series objects in
         a collection.
@@ -4099,18 +4276,21 @@ def grid_properties(self, step_style='median'):
         Returns
         -------
 
-        grid_properties : numpy array
-            n x 3 array, where n is the number of series
+        increments : numpy array
+            n x 3 array, where n is the number of series,
+            index 0 is the earliest time
+            index 1 is the latest time
+            index 2 is the step chosen according to step_style
 
         '''
         gp = np.empty((len(self.series_list),3)) # obtain grid parameters
         for idx,item in enumerate(self.series_list):
             item      = item.clean(verbose=idx==0)
-            gp[idx,:] = tsutils.grid_properties(item.time, step_style=step_style)
+            gp[idx,:] = tsutils.increments(item.time, step_style=step_style)
 
         return gp
 
-    def common_time(self, method='interp', common_step = 'max', start=None, stop = None, step=None, step_style = None, **kwargs):
+    def common_time(self, method='interp', step = None, start = None, stop = None, step_style = None, **kwargs):
         ''' Aligns the time axes of a MultipleSeries object, via binning
         interpolation., or Gaussian kernel. Alignmentis critical for workflows
         that need to assume a common time axis for the group of series under consideration.
@@ -4122,32 +4302,30 @@ def common_time(self, method='interp', common_step = 'max', start=None, stop = N
         stop  : the earliest stop date of the bunch (minimum of the maxima)
         step  : The representative spacing between consecutive values
 
-        Optional arguments for binning or interpolation are those of the underling functions.
+        Optional arguments for binning, Gaussian kernel (gkernel) interpolation are those of the underling functions.
 
-        If the time axis are retrograde, this step makes them prograde.
+        If any of the time axes are retrograde, this step makes them prograde.
 
         Parameters
         ----------
         method :  string
             either 'bin', 'interp' [default] or 'gkernel'
-        common_step : string
-            Method to obtain a representative step among all Series
-            Valid entries: 'median' [default], 'mean', 'mode' or 'max'
+        step : float
+            common step for all time axes  [default = None]
         start : float
             starting point of the common time axis [default = None]
         stop : float
             end point of the common time axis [default = None]
-        step : float
-            increment the common time axis  [default = None]
-
-        if not provided, `pyleoclim` will use `grid_properties()` to determine these parameters
-
-        step_style : 'string'
-            step style to be applied from `grid_properties` [default = None]
-
+        
+        step_style : string
+            Method to obtain a representative step among all Series
+            (using tsutils.increments)
+            Valid entries: 'median', 'mean', 'mode' or 'max'
+            Default value is None, so that it will be chosen according
+            to the method: 'max' for bin and gkernel, 'mean' for interp. 
 
         kwargs: dict
-            keyword arguments (dictionary) of the various methods
+            keyword arguments (dictionary) of the (bin, gkernel or interp) methods
 
         Returns
         -------
@@ -4161,7 +4339,7 @@ def common_time(self, method='interp', common_step = 'max', start=None, stop = N
         pyleoclim.utils.tsutils.bin : put timeseries values into bins of equal size (possibly leaving NaNs in).
         pyleoclim.utils.tsutils.gkernel : coarse-graining using a Gaussian kernel
         pyleoclim.utils.tsutils.interp : interpolation onto a regular grid (default = linear interpolation)
-        pyleoclim.utils.tsutils.grid_properties : infer grid properties
+        pyleoclim.utils.tsutils.increments : infer grid properties
 
         Examples
         --------
@@ -4192,7 +4370,7 @@ def common_time(self, method='interp', common_step = 'max', start=None, stop = N
             # create MS object from the list
             ms = pyleo.MultipleSeries(serieslist)
 
-            fig, ax = plt.subplots(2,2)
+            fig, ax = plt.subplots(2,2,sharex=True,sharey=True)
             ax = ax.flatten()
             # apply common_time with default parameters
             msc = ms.common_time()
@@ -4217,37 +4395,46 @@ def common_time(self, method='interp', common_step = 'max', start=None, stop = N
             pyleo.closefig(fig)
 
         '''
-
-        if step_style == None:
+        # specify stepping style
+        if step_style == None: # if step style isn't specified, pick a robust choice according to method
             if method == 'bin' or method == 'gkernel':
                step_style = 'max'
             elif  method == 'interp':
                step_style = 'mean'
-
-        gp = self.grid_properties(step_style=step_style)
-
-        # define parameters for common time axis
-        start = gp[:,0].max()
-        stop  = gp[:,1].min()
+               
+        # obtain grid properties with given step_style
+        gp = self.increments(step_style=step_style)  
+        
+        # define grid step     
+        if step is not None and step > 0:
+            common_step = step 
+        else:
+            if step_style == 'mean':
+                common_step = gp[:,2].mean()
+            elif step_style == 'max':
+                common_step = gp[:,2].max()
+            elif step_style == 'mode':
+                common_step = stats.mode(gp[:,2])[0][0]
+            else:
+                common_step = np.median(gp[:,2])
+   
+        # define start and stop
+        if start is None: 
+            start = gp[:,0].max() # pick the latest of the start times
+        
+        if stop is None:
+            stop  = gp[:,1].min() # pick the earliest of the stop times
+        
         if start > stop:
             raise ValueError('At least one series has no common time interval with others. Please check the time axis of the series.')
-
-        if step is None:
-            if common_step == 'mean':
-                step = gp[:,2].mean()
-            elif common_step == 'max':
-                step = gp[:,2].max()
-            elif common_step == 'mode':
-                step = stats.mode(gp[:,2])[0][0]
-            else:
-                step = np.median(gp[:,2])
-
+    
         ms = self.copy()
 
+        # apply each method
         if method == 'bin':
             for idx,item in enumerate(self.series_list):
                 ts = item.copy()
-                d = tsutils.bin(ts.time, ts.value, bin_size=step, start=start, stop=stop, evenly_spaced = False, **kwargs)
+                d = tsutils.bin(ts.time, ts.value, bin_size=common_step, start=start, stop=stop, evenly_spaced = False, **kwargs)
                 ts.time  = d['bins']
                 ts.value = d['binned_values']
                 ms.series_list[idx] = ts
@@ -4255,7 +4442,7 @@ def common_time(self, method='interp', common_step = 'max', start=None, stop = N
         elif method == 'interp':
             for idx,item in enumerate(self.series_list):
                 ts = item.copy()
-                ti, vi = tsutils.interp(ts.time, ts.value, step=step, start=start, stop=stop,**kwargs)
+                ti, vi = tsutils.interp(ts.time, ts.value, step=common_step, start=start, stop=stop,**kwargs)
                 ts.time  = ti
                 ts.value = vi
                 ms.series_list[idx] = ts
@@ -4263,7 +4450,7 @@ def common_time(self, method='interp', common_step = 'max', start=None, stop = N
         elif method == 'gkernel':
             for idx,item in enumerate(self.series_list):
                 ts = item.copy()
-                ti, vi = tsutils.gkernel(ts.time,ts.value,step=step, start=start, stop=stop, **kwargs)
+                ti, vi = tsutils.gkernel(ts.time,ts.value,step=common_step, start=start, stop=stop, **kwargs)
                 ts.time  = ti
                 ts.value = vi
                 ms.series_list[idx] = ts
@@ -4785,7 +4972,7 @@ def spectral(self, method='lomb_scargle', settings=None, mute_pbar=False, freq_m
         ----------
 
         method : str
-            {'wwz', 'mtm', 'lomb_scargle', 'welch', 'periodogram'}
+            {'wwz', 'mtm', 'lomb_scargle', 'welch', 'periodogram', 'cwt'}
 
         freq_method : str
             {'log','scale', 'nfft', 'lomb_scargle', 'welch'}
@@ -4806,7 +4993,7 @@ def spectral(self, method='lomb_scargle', settings=None, mute_pbar=False, freq_m
             Mute the progress bar. Default is False.
 
         scalogram_list : pyleoclim.MultipleScalogram object, optional
-            Multiple scalogram object containing pre-computed scalograms to use when calculating spectra, only works with wwz
+            Multiple scalogram object containing pre-computed scalograms to use when calculating spectra, only works with wwz or cwt
 
         Returns
         -------
@@ -4826,6 +5013,8 @@ def spectral(self, method='lomb_scargle', settings=None, mute_pbar=False, freq_m
 
         pyleoclim.utils.spectral.wwz_psd : Spectral analysis using the Wavelet Weighted Z transform
 
+        pyleoclim.utils.spectral.cwt_psd : Spectral analysis using the continuous Wavelet Transform as implemented by Torrence and Compo
+
         pyleoclim.utils.wavelet.make_freq_vector : Functions to create the frequency vector
 
         pyleoclim.utils.tsutils.detrend : Detrending function
@@ -4839,7 +5028,7 @@ def spectral(self, method='lomb_scargle', settings=None, mute_pbar=False, freq_m
         settings = {} if settings is None else settings.copy()
 
         psd_list = []
-        if method == 'wwz' and scalogram_list:
+        if method in ['wwz','cwt'] and scalogram_list:
             scalogram_list_len = len(scalogram_list.scalogram_list)
             series_len = len(self.series_list)
 
@@ -4873,7 +5062,7 @@ def wavelet(self, method='wwz', settings={}, freq_method='log', ntau=None, freq_
         Parameters
         ----------
         method : {wwz, cwt}
-            Whether to use the wwz method for unevenly spaced timeseries or traditional cwt (from pywavelets)
+            Whether to use the wwz method for unevenly spaced timeseries or traditional cwt (from Torrence and Compo)
 
         settings : dict, optional
             Settings for the particular method. The default is {}.
@@ -4906,6 +5095,8 @@ def wavelet(self, method='wwz', settings={}, freq_method='log', ntau=None, freq_
         --------
         pyleoclim.utils.wavelet.wwz : wwz function
 
+        pyleoclim.utils.wavelet.cwt : cwt function
+
         pyleoclim.utils.wavelet.make_freq_vector : Functions to create the frequency vector
 
         pyleoclim.utils.tsutils.detrend : Detrending function
@@ -4914,6 +5105,12 @@ def wavelet(self, method='wwz', settings={}, freq_method='log', ntau=None, freq_
 
         pyleoclim.core.ui.MultipleScalogram : Multiple Scalogram object
 
+        References
+        ----------
+
+        Torrence, C. and G. P. Compo, 1998: A Practical Guide to Wavelet Analysis. Bull. Amer. Meteor. Soc., 79, 61-78.
+        Python routines available at http://paos.colorado.edu/research/wavelets/
+
         '''
         settings = {} if settings is None else settings.copy()
 
diff --git a/pyleoclim/utils/jsonutils.py b/pyleoclim/utils/jsonutils.py
index 73e139c5..c5fa1fa4 100644
--- a/pyleoclim/utils/jsonutils.py
+++ b/pyleoclim/utils/jsonutils.py
@@ -8,11 +8,14 @@
 Utilities to import/export Pyleoclim objects from JSON files
 """
 
-__all__ =['PyleoObj_to_json', 'json_to_Series','json_to_PSD','json_to_Scalogram']
+__all__ =['PyleoObj_to_json', 'json_to_Series', 'json_to_LipdSeries',
+          'json_to_PSD','json_to_Scalogram', 'json_to_Coherence', 'json_to_MultiplePSD']
 
 import numpy as np
 import json
 import pyleoclim as pyleo
+from urllib.request import urlopen
+import re
 
 def transform(obj_dict):
     '''
@@ -32,7 +35,7 @@ def transform(obj_dict):
     for k in obj_dict.keys():
         if isinstance(obj_dict[k],(np.ndarray)):
             obj_dict[k] = obj_dict[k].tolist()
-        elif isinstance(obj_dict[k],(pyleo.core.ui.Series,pyleo.core.ui.Scalogram)):
+        elif isinstance(obj_dict[k],(pyleo.core.ui.Series,pyleo.core.ui.Scalogram,pyleo.core.ui.PSD)):
             obj_dict[k]=PyleoObj_to_json(obj_dict[k],dict_return=True)
         elif isinstance(obj_dict[k],pyleo.core.ui.MultiplePSD):
             obj_dict[k]=PyleoObj_to_json(obj_dict[k],dict_return=True)
@@ -52,6 +55,11 @@ def transform(obj_dict):
             obj_dict[k]['scalogram_list']=c            
         elif isinstance(obj_dict[k],(dict)):
             obj_dict[k]=transform(obj_dict[k])
+        elif isinstance(obj_dict[k],(list)):
+             new_list = []
+             for item in obj_dict[k]:
+                 new_list.append(PyleoObj_to_json(item,dict_return=True))
+             obj_dict[k] = new_list
     return obj_dict
 
 def list_to_array(obj_dict):
@@ -78,6 +86,8 @@ def list_to_array(obj_dict):
             obj_dict[k]=obj_dict[k]
     return obj_dict
 
+ 
+
 def PyleoObj_to_json(PyleoObj,filename=None,dict_return=False):
     '''
     
@@ -116,16 +126,29 @@ def open_json(filename):
     Parameters
     ----------
     filename : str
-        Path to the json file
+        Path to the json file or URL
 
     Returns
     -------
     t : dict
-        A Python diciotanary from the json
+        A Python dictionary from the json
 
     '''
-    with open(filename,'r') as f:
-        t=json.load(f)
+    
+    regex = re.compile(
+        r'^(?:http|ftp)s?://' # http:// or https://
+        r'(?:(?:[A-Z0-9](?:[A-Z0-9-]{0,61}[A-Z0-9])?\.)+(?:[A-Z]{2,6}\.?|[A-Z0-9-]{2,}\.?)|' #domain...
+        r'localhost|' #localhost...
+        r'\d{1,3}\.\d{1,3}\.\d{1,3}\.\d{1,3})' # ...or ip
+        r'(?::\d+)?' # optional port
+        r'(?:/?|[/?]\S+)$', re.IGNORECASE)
+    
+    if (re.match(regex, filename) is not None)==True:
+        response = urlopen(filename)
+        t = json.loads(response.read())
+    else:    
+        with open(filename,'r') as f:
+            t=json.load(f)
     return t
 
 def json_to_Series(filename):
@@ -135,7 +158,7 @@ def json_to_Series(filename):
     Parameters
     ----------
     filename : str
-        The name of the JSON file containing the Series information. 
+        The name of the JSON file/URL containing the Series information. 
 
     Returns
     -------
@@ -163,7 +186,7 @@ def json_to_LipdSeries(filename):
     Parameters
     ----------
     filename : str
-        The name of the JSON file containing the Series information
+        The name of the JSON file/URL containing the Series information
 
     Returns
     -------
@@ -171,6 +194,7 @@ def json_to_LipdSeries(filename):
         A pyleoclim.Series object
 
     '''
+        
     t = open_json(filename)
     ts = pyleo.LipdSeries(tso=t['lipd_ts'],clean_ts=t['clean_ts'],verbose=t['verbose']) 
     return ts
@@ -229,7 +253,7 @@ def json_to_PSD(filename):
     Parameters
     ----------
     filename : str
-        Name of the json file containing the necessary information
+        Name of the json file/URL containing the necessary information
 
     Returns
     -------
@@ -237,8 +261,7 @@ def json_to_PSD(filename):
         a pyleoclim PSD object
 
     '''
-    with  open(filename,'r') as f:
-        t = json.load(f)
+    t = open_json(filename)
     t = list_to_array(t)
     
     #Deal with significance testing 
@@ -312,7 +335,7 @@ def Scalogram_to_MultipleScalogram(series_list):
                                freq_kwargs=t['freq_kwargs'],
                                period_unit=t['period_unit'],
                                time_label=t['time_label'],
-                               wwz_Neffs=t['wwz_Neffs'],
+                               wwz_Neffs=np.array(t['wwz_Neffs']),
                                signif_scals=t['signif_scals'])
 
         d.append(scalogram_obj)
@@ -320,9 +343,21 @@ def Scalogram_to_MultipleScalogram(series_list):
     return mscalogram
 
 def json_to_Scalogram(filename):
-    with open(filename,'r') as f:
-        t = json.load(f)
-    t = list_to_array(t)
+    '''
+    Transform a Scalogram object stored in JSON format back to a Pyleoclim object
+
+    Parameters
+    ----------
+    filename : str
+        Path of the JSON file/URL
+
+    Returns
+    -------
+    scalogram : pyleoclim.Scalogram
+        A Pyleoclim Scalogram object.
+
+    '''
+    t = open_json(filename)
     temp = t['timeseries']
     ts = pyleo.Series(time=np.array(temp['time']),
                      value=np.array(temp['value']),
@@ -356,7 +391,7 @@ def json_to_Scalogram(filename):
                                freq_kwargs=t['freq_kwargs'],
                                period_unit=t['period_unit'],
                                time_label=t['time_label'],
-                               wwz_Neffs=t['wwz_Neffs'],
+                               wwz_Neffs=np.array(t['wwz_Neffs']),
                                signif_scals=d)
     
     return scalogram
@@ -369,7 +404,7 @@ def json_to_Coherence(filename):
     Parameters
     ----------
     filename : str
-        json file to unpack.
+        Path/URL to json file to unpack.
 
     Returns
     -------
@@ -377,8 +412,7 @@ def json_to_Coherence(filename):
         A coherence object
 
     '''
-    with open(filename,'r') as f:
-        t = json.load(f)
+    t = open_json(filename)
     t = list_to_array(t)
     temp1 = t['timeseries1']
     ts1 = pyleo.Series(time=np.array(temp1['time']),
@@ -418,4 +452,58 @@ def json_to_Coherence(filename):
                                time_label=t['time_label'],
                                signif_qs = c,
                                signif_method=t['signif_method'])
-    return coherence
\ No newline at end of file
+    return coherence
+
+def json_to_MultiplePSD(filename):
+    '''
+    Transform a MultiplePSD object stored in JSON format back to a Pyleoclim object
+
+    Parameters
+    ----------
+    filename : str
+        Filename or URL for the JSON file
+
+    Returns
+    -------
+    mpsd : pyleoclim.MultiplePSD
+        The MultiplePSD object
+
+    '''
+    
+    t = open_json(filename)
+    psd_list = []
+    for item in t['psd_list']:
+        item=list_to_array(item)
+        if type(item['signif_qs']) is dict:
+            c = PSD_to_MultiplePSD(item['signif_qs']['psd_list'])
+        else:
+            c = item['signif_qs']
+        psd = pyleo.PSD(frequency=np.array(item['frequency']),
+                        amplitude=np.array(item['amplitude']),
+                        label=item['label'],
+                        timeseries = pyleo.Series(time=np.array(item['timeseries']['time']),
+                                                 value=np.array(item['timeseries']['value']),
+                                                 time_name=item['timeseries']['time_name'],
+                                                 time_unit=item['timeseries']['time_unit'],
+                                                 value_name=item['timeseries']['value_name'],
+                                                 value_unit=item['timeseries']['value_unit'],
+                                                 label=item['timeseries']['label'],
+                                                 clean_ts=item['timeseries']['clean_ts'], 
+                                                 verbose=item['timeseries']['verbose']),
+                        spec_method = item['spec_method'],
+                        spec_args =  item['spec_args'],
+                        signif_qs = c,
+                        signif_method = item['signif_method'],
+                        plot_kwargs=item['plot_kwargs'],
+                        period_unit=item['period_unit'],
+                        beta_est_res=item['beta_est_res'])
+        psd_list.append(psd)
+    
+    if t['beta_est_res']==None:
+        pass
+    else:
+        t['beta_est_res']=np.array(t['beta_est_res'])
+    
+    mpsd=pyleo.MultiplePSD(psd_list=psd_list,beta_est_res=t['beta_est_res'])
+    
+    return mpsd
\ No newline at end of file
diff --git a/pyleoclim/utils/spectral.py b/pyleoclim/utils/spectral.py
index 9da09da6..2ff39d16 100644
--- a/pyleoclim/utils/spectral.py
+++ b/pyleoclim/utils/spectral.py
@@ -34,6 +34,7 @@
     prepare_wwz,
     wwz,
     wwa2psd,
+    cwt,
 )
 #from .tsutils import clean_ts, interp, bin
 
@@ -741,38 +742,95 @@ def wwz_psd(ys, ts, freq=None, freq_method='log', freq_kwargs=None,
         wwz_freq = res_wwz.freq
 
     psd = wwa2psd(wwa, ts_cut, wwz_Neffs, freq=wwz_freq, Neff=Neff, anti_alias=anti_alias, avgs=avgs)
-    #  psd[1/freqs > np.max(coi)] = np.nan  # cut off the unreliable part out of the coi
-    #  psd = psd[1/freqs <= np.max(coi)] # cut off the unreliable part out of the coi
-    #  freqs = freqs[1/freqs <= np.max(coi)]
-
-    # Monte-Carlo simulations of AR1 process
-    #nf = np.size(freq)
-
-    #  psd_ar1 = np.ndarray(shape=(nMC, nf))
-
-    #  if nMC >= 1:
-        #  #  tauest = wa.tau_estimation(ys_cut, ts_cut, detrend=detrend)
-
-        #  for i in tqdm(range(nMC), desc='Monte-Carlo simulations'):
-            #  #  r = wa.ar1_model(ts_cut, tauest)
-            #  r = ar1_sim(ys_cut, np.size(ts_cut), 1, ts=ts_cut)
-            #  res_red = wwz(r, ts_cut, freq=freq, tau=tau, c=c, nproc=nproc, nMC=0,
-                                                                     #  detrend=detrend, params=params,
-                                                                     #  gaussianize=gaussianize, standardize=standardize,
-                                                                     #  method=method)
-            #  psd_ar1[i, :] = wa.wwa2psd(res_red.wwa, ts_cut, res_red.Neffs,
-                                       #  freq=res_red.freq, Neff=Neff, anti_alias=anti_alias, avgs=avgs)
-            #  #  psd_ar1[i, 1/freqs_red > np.max(coi_red)] = np.nan  # cut off the unreliable part out of the coi
-            #  #  psd_ar1 = psd_ar1[1/freqs_red <= np.max(coi_red)] # cut off the unreliable part out of the coi
-
-        #  psd_ar1_q95 = mquantiles(psd_ar1, 0.95, axis=0)[0]
-
-    #  else:
-        #  psd_ar1_q95 = None
-
-    # Results = collections.namedtuple('Results', ['psd', 'freq', 'psd_ar1_q95', 'psd_ar1'])
-    # res = Results(psd=psd, freq=freq, psd_ar1_q95=psd_ar1_q95, psd_ar1=psd_ar1)
     Results = collections.namedtuple('Results', ['psd', 'freq'])
     res = Results(psd=psd, freq=freq)
 
     return res
+
+def cwt_psd(ys, ts, freq=None, freq_method='log', freq_kwargs=None,scale = None, 
+            detrend=False,sg_kwargs={}, gaussianize=False, standardize =False, pad=False, 
+            mother='MORLET',param=None, cwt_res=None):
+    '''
+    Wrapper function to implement Torrence and Compo continuous wavelet transform
+
+    Parameters
+    ----------
+    ys : numpy.array
+        the time series.
+    ts : numpy.array
+        the time axis.
+    freq : numpy.array, optional
+        The frequency vector. The default is None, which will prompt the use of one the underlying functions
+    freq_method : string, optional
+        The method by which to obtain the frequency vector. The default is 'log'.
+        Options are 'log' (default), 'nfft', 'lomb_scargle', 'welch', and 'scale'
+    freq_kwargs : dict, optional
+        Optional parameters for the choice of the frequency vector. See make_freq_vector and additional methods for details. The default is {}.
+    scale : numpy.array
+        Optional scale vector in place of a frequency vector. Default is None. If scale is not None, frequency method and attached arguments will be ignored. 
+    detrend : bool, string, {'linear', 'constant', 'savitzy-golay', 'emd'}
+        Whether to detrend and with which option. The default is False.
+    sg_kwargs : dict, optional
+        Additional parameters for the savitzy-golay method. The default is {}.
+    gaussianize : bool, optional
+        Whether to gaussianize. The default is False.
+    standardize : bool, optional
+        Whether to standardize. The default is False.     
+    pad : bool, optional
+        Whether or not to pad the timeseries. with zeroes to get N up to the next higher power of 2. 
+        This prevents wraparound from the end of the time series to the beginning, and also speeds up the FFT's used to do the wavelet transform.
+        This will not eliminate all edge effects. The default is False.
+    mother : string, optional
+        the mother wavelet function. The default is 'MORLET'. Options are: 'MORLET', 'PAUL', or 'DOG'
+    param : flaot, optional
+        the mother wavelet parameter. The default is None since it varies for each mother
+            - For 'MORLET' this is k0 (wavenumber), default is 6.
+            - For 'PAUL' this is m (order), default is 4.
+            - For 'DOG' this is m (m-th derivative), default is 2.
+    cwt_res : dict
+        Results from pyleoclim.utils.wavelet.cwt
+
+    Returns
+    -------
+    res : dict
+        Dictionary containing:
+            - psd: the power density function
+            - freq: frequency vector
+            - scale: the scale vector
+            - mother: the mother wavelet
+            - param : the wavelet parameter
+            
+    See also
+    --------
+    
+    pyleoclim.utils.wavelet.make_freq_vector : make the frequency vector with various methods
+    
+    pyleoclim.utils.wavelet.cwt: Torrence and Compo implementation of the continuous wavelet transform 
+    
+    pyleoclim.utils.tsutils.detrend : detrending functionalities in Pyleoclim
+    
+    References
+    ----------
+    
+    Torrence, C. and G. P. Compo, 1998: A Practical Guide to Wavelet Analysis. Bull. Amer. Meteor. Soc., 79, 61-78.
+    Python routines available at http://paos.colorado.edu/research/wavelets/
+    
+    '''
+    
+    
+        #get the wavelet:
+    if cwt_res is None:
+        cwt_res = cwt(ys,ts,freq=freq, freq_method=freq_method, freq_kwargs=freq_kwargs,
+              scale = scale, detrend=detrend,sg_kwargs=sg_kwargs, gaussianize=gaussianize, 
+              standardize = standardize, pad=pad, mother=mother, param=param) 
+        n= len(ts)
+    else:
+        n=len(cwt_res.time)
+    
+    psd = np.sum(cwt_res.amplitude.T**2,axis=1)/n
+    
+    
+    Results = collections.namedtuple('Results', ['psd', 'freq','scale','mother','param'])
+    res = Results(psd=psd, freq=cwt_res.freq, scale=cwt_res.scale, mother=cwt_res.mother,param=cwt_res.param)
+
+    return res
\ No newline at end of file
diff --git a/pyleoclim/utils/tsmodel.py b/pyleoclim/utils/tsmodel.py
index e4baccfc..12d5b9b0 100644
--- a/pyleoclim/utils/tsmodel.py
+++ b/pyleoclim/utils/tsmodel.py
@@ -14,6 +14,7 @@
 
 __all__ = [
     'ar1_sim',
+    'ar1_fit',
     'colored_noise',
     'colored_noise_2regimes',
     'gen_ar1_evenly',
diff --git a/pyleoclim/utils/tsutils.py b/pyleoclim/utils/tsutils.py
index aff83513..39dea22b 100644
--- a/pyleoclim/utils/tsutils.py
+++ b/pyleoclim/utils/tsutils.py
@@ -13,7 +13,7 @@
     'bin',
     'interp',
     'gkernel',
-    'grid_properties',
+    'increments',
     'standardize',
     'ts2segments',
     'annualize',
@@ -193,7 +193,7 @@ def gkernel(t,y, h = 3.0, step=None,start=None,stop=None, step_style = 'max'):
         The interpolation step. Default is max spacing between consecutive points.
    
     step_style : 'string'
-            step style to be applied from `grid_properties` [default = 'max']
+            step style to be applied from `increments` [default = 'max']
     
     start : float
         where/when to start the interpolation. Default is min(t).
@@ -224,7 +224,7 @@ def gkernel(t,y, h = 3.0, step=None,start=None,stop=None, step_style = 'max'):
         
     # get the interpolation step if not provided
     if step is None:
-        _, _, step = grid_properties(np.asarray(t), step_style = step_style)
+        _, _, step = increments(np.asarray(t), step_style = step_style)
         # Get the start and end point if not given
     if start is None:
         start = np.nanmin(np.asarray(t))
@@ -254,8 +254,8 @@ def gkernel(t,y, h = 3.0, step=None,start=None,stop=None, step_style = 'max'):
     return tc, yc
 
 
-def grid_properties(x,step_style='median'):
-    ''' Establishes the grid properties of a numerical array:
+def increments(x,step_style='median'):
+    ''' Establishes the increments of a numerical array:
         start, stop, and representative step.
 
     Parameters
@@ -331,7 +331,7 @@ def interp(x,y, interp_type='linear', step=None,start=None,stop=None, step_style
 
     # get the interpolation step if not available
     if step is None:
-        _, _, step = grid_properties(np.asarray(x), step_style = step_style)
+        _, _, step = increments(np.asarray(x), step_style = step_style)
 
         # Get the start and end point if not given
     if start is None:
@@ -352,86 +352,86 @@ def interp(x,y, interp_type='linear', step=None,start=None,stop=None, step_style
     return xi, yi
 
 
-def on_common_axis(x1, y1, x2, y2, method = 'interpolation', step=None, start=None, stop=None):
-    """Places two timeseries on a common axis
-
-    Note this function assumes that the time representation and units are the same (e.g., BP vs CE)
-
-    Parameters
-    ----------
-    x1 : array
-        x-axis values of the first timeseries
-    y1 : array
-        y-axis values of the first timeseries
-    x2 : array
-        x-axis values of the second timeseries
-    y2 : array
-        y-axis values of the second timeseries
-    method : str
-        Which method to use to get the timeseries on the same x axis.
-        Valid entries: 'interpolation' (default, linear interpolation),
-        'bin', 'None'. 'None' only cuts the timeseries to the common
-        period but does not attempt to generate a common time axis
-    step : float
-        The interpolation step. Default is mean resolution
-        of lowest resolution series
-    start : float
-        where/when to start. Default is the maximum of the minima of
-        the two timeseries
-    stop : float
-        Where/when to stop. Default is the minimum of the maxima of
-        the two timeseries
-
-    Returns
-    -------
-
-    xi1, xi2 : array
-        The interpolated x-axis
-    interp_values1, interp_values2 : array
-        the interpolated y-values
-    """
-    # make sure that x1, y1, x2, y2 are numpy arrays
-    x1 = np.array(x1, dtype='float64')
-    y1 = np.array(y1, dtype='float64')
-    x2 = np.array(x2, dtype='float64')
-    y2 = np.array(y2, dtype='float64')
-
-    # Find the mean/max x-axis is not provided
-    if start is None:
-        start = np.nanmax([np.nanmin(x1), np.nanmin(x2)])
-    if stop is None:
-        stop = np.nanmin([np.nanmax(x1), np.nanmax(x2)])
-
-    # Get the interp_step
-    if step is None:
-        step = np.nanmin([np.nanmean(np.diff(x1)), np.nanmean(np.diff(x2))])
-
-    if method == 'interpolation':
-    # perform the interpolation
-        xi1, interp_values1 = interp(x1, y1, step=step, start=start,
-                                stop=stop)
-        xi2, interp_values2 = interp(x2, y2, step=step, start=start,
-                                stop=stop)
-    elif method == 'bin':
-        xi1, interp_values1, _ , _ = bin(x1, y1, bin_size=step, start=start,
-                                stop=stop)
-        xi2, interp_values2, _ , _ = bin(x2, y2, bin_size=step, start=start,
-                                stop=stop)
-    elif method == None:
-        min_idx1 = np.where(x1>=start)[0][0]
-        min_idx2 = np.where(x2>=start)[0][0]
-        max_idx1 = np.where(x1<=stop)[0][-1]
-        max_idx2 = np.where(x2<=stop)[0][-1]
-
-        xi1 = x1[min_idx1:max_idx1+1]
-        xi2 = x2[min_idx2:max_idx2+1]
-        interp_values1 = y1[min_idx1:max_idx1+1]
-        interp_values2 = y2[min_idx2:max_idx2+1]
-
-    else:
-        raise KeyError('Not a valid interpolation method')
-
-    return xi1, xi2, interp_values1, interp_values2
+# def on_common_axis(x1, y1, x2, y2, method = 'interpolation', step=None, start=None, stop=None):
+#     """Places two timeseries on a common axis
+
+#     Note this function assumes that the time representation and units are the same (e.g., BP vs CE)
+
+#     Parameters
+#     ----------
+#     x1 : array
+#         x-axis values of the first timeseries
+#     y1 : array
+#         y-axis values of the first timeseries
+#     x2 : array
+#         x-axis values of the second timeseries
+#     y2 : array
+#         y-axis values of the second timeseries
+#     method : str
+#         Which method to use to get the timeseries on the same x axis.
+#         Valid entries: 'interpolation' (default, linear interpolation),
+#         'bin', 'None'. 'None' only cuts the timeseries to the common
+#         period but does not attempt to generate a common time axis
+#     step : float
+#         The interpolation step. Default is mean resolution
+#         of lowest resolution series
+#     start : float
+#         where/when to start. Default is the maximum of the minima of
+#         the two timeseries
+#     stop : float
+#         Where/when to stop. Default is the minimum of the maxima of
+#         the two timeseries
+
+#     Returns
+#     -------
+
+#     xi1, xi2 : array
+#         The interpolated x-axis
+#     interp_values1, interp_values2 : array
+#         the interpolated y-values
+#     """
+#     # make sure that x1, y1, x2, y2 are numpy arrays
+#     x1 = np.array(x1, dtype='float64')
+#     y1 = np.array(y1, dtype='float64')
+#     x2 = np.array(x2, dtype='float64')
+#     y2 = np.array(y2, dtype='float64')
+
+#     # Find the mean/max x-axis is not provided
+#     if start is None:
+#         start = np.nanmax([np.nanmin(x1), np.nanmin(x2)])
+#     if stop is None:
+#         stop = np.nanmin([np.nanmax(x1), np.nanmax(x2)])
+
+#     # Get the interp_step
+#     if step is None:
+#         step = np.nanmin([np.nanmean(np.diff(x1)), np.nanmean(np.diff(x2))])
+
+#     if method == 'interpolation':
+#     # perform the interpolation
+#         xi1, interp_values1 = interp(x1, y1, step=step, start=start,
+#                                 stop=stop)
+#         xi2, interp_values2 = interp(x2, y2, step=step, start=start,
+#                                 stop=stop)
+#     elif method == 'bin':
+#         xi1, interp_values1, _ , _ = bin(x1, y1, bin_size=step, start=start,
+#                                 stop=stop)
+#         xi2, interp_values2, _ , _ = bin(x2, y2, bin_size=step, start=start,
+#                                 stop=stop)
+#     elif method == None:
+#         min_idx1 = np.where(x1>=start)[0][0]
+#         min_idx2 = np.where(x2>=start)[0][0]
+#         max_idx1 = np.where(x1<=stop)[0][-1]
+#         max_idx2 = np.where(x2<=stop)[0][-1]
+
+#         xi1 = x1[min_idx1:max_idx1+1]
+#         xi2 = x2[min_idx2:max_idx2+1]
+#         interp_values1 = y1[min_idx1:max_idx1+1]
+#         interp_values2 = y2[min_idx2:max_idx2+1]
+
+#     else:
+#         raise KeyError('Not a valid interpolation method')
+
+#     return xi1, xi2, interp_values1, interp_values2
 
 
 def standardize(x, scale=1, axis=0, ddof=0, eps=1e-3):
diff --git a/pyleoclim/utils/wavelet.py b/pyleoclim/utils/wavelet.py
index 29a997d1..f30286cf 100644
--- a/pyleoclim/utils/wavelet.py
+++ b/pyleoclim/utils/wavelet.py
@@ -10,9 +10,11 @@
 
 
 __all__ = [
-    #'cwt',
+    'cwt',
     'wwz',
-    'xwc',
+    'xwt',
+    'wtc',
+    'wavelet_coherence',
 ]
 
 import numpy as np
@@ -32,10 +34,8 @@
 from .tsutils import preprocess
 from .tsbase import (
     clean_ts,
-    is_evenly_spaced,
-)
-
-from .filter import ts_pad
+    is_evenly_spaced)
+from .tsmodel import(ar1_fit)
 
 warnings.filterwarnings("ignore", category=NumbaPerformanceWarning)
 
@@ -1201,14 +1201,14 @@ def wwz(ys, ts, tau=None, ntau=None, freq=None, freq_method='log', freq_kwargs={
     nproc : int
         the number of processes for multiprocessing
 
-    detrend : string, {None, 'linear', 'constant', 'savitzy-golay'}
+    detrend : string, {None, 'linear', 'constant', 'savitzy-golay', 'emd'}
         available methods for detrending, including
 
         - None: the original time series is assumed to have no trend;
         - 'linear': a linear least-squares fit to `ys` is subtracted;
         - 'constant': the mean of `ys` is subtracted
         - 'savitzy-golay': ys is filtered using the Savitzky-Golay filters and the resulting filtered series is subtracted from y.
-        Empirical mode decomposition. The last mode is assumed to be the trend and removed from the series
+        - 'emd': Empirical mode decomposition. The last mode is assumed to be the trend and removed from the series
 
     sg_kwargs : dict
         The parameters for the Savitzky-Golay filters. See :func:`pyleoclim.utils.filter.savitzky_golay()` for details.
@@ -1271,7 +1271,7 @@ def wwz(ys, ts, tau=None, ntau=None, freq=None, freq_method='log', freq_kwargs={
 
     pyleoclim.utils.wavelet.kirchner_f2py : Returns the weighted wavelet amplitude (WWA) modified by Kirchner. Uses Fortran. Fastest method but requires a compiler.
 
-    pyleoclim.utils.filter.savitzky_golay : Smooth (and optionally differentiate) data with a Savitzky-Golay filter.
+    pyleoclim.utils.tsutils.detrend : detrending functionalities in Pyleoclim
 
     pyleoclim.utils.wavelet.make_freq_vector : Make frequency vector
 
@@ -1361,7 +1361,7 @@ def wwz(ys, ts, tau=None, ntau=None, freq=None, freq_method='log', freq_kwargs={
 
     return res
 
-def xwc(ys1, ts1, ys2, ts2, smooth_factor=0.25,
+def wavelet_coherence(ys1, ts1, ys2, ts2, smooth_factor=0.25,
         tau=None, freq=None, freq_method='log', freq_kwargs=None,
         c=1/(8*np.pi**2), Neff=3, nproc=8, detrend=False, sg_kwargs=None,
         nMC=200,
@@ -1491,8 +1491,8 @@ def xwc(ys1, ts1, ys2, ts2, smooth_factor=0.25,
     wt_coeff1 = res_wwz1.coeff[1] - res_wwz1.coeff[2]*1j
     wt_coeff2 = res_wwz2.coeff[1] - res_wwz2.coeff[2]*1j
 
-    xw_coherence, xw_phase = wavelet_coherence(wt_coeff1, wt_coeff2, freq, tau, smooth_factor=smooth_factor)
-    xwt, xw_amplitude, _ = cross_wt(wt_coeff1, wt_coeff2)
+    xw_coherence, xw_phase = wtc(wt_coeff1, wt_coeff2, freq, tau, smooth_factor=smooth_factor)
+    xw_t, xw_amplitude, _ = xwt(wt_coeff1, wt_coeff2)
 
     # Monte-Carlo simulations of AR1 process
     nt = np.size(tau)
@@ -1527,8 +1527,8 @@ def xwc(ys1, ts1, ys2, ts2, smooth_factor=0.25,
         #  AR1_q = None
 
     coi = make_coi(tau, Neff=Neff)
-    Results = collections.namedtuple('Results', ['xw_coherence', 'xw_amplitude', 'xw_phase', 'xwt', 'freq', 'time', 'AR1_q', 'coi'])
-    res = Results(xw_coherence=xw_coherence, xw_amplitude=xw_amplitude, xw_phase=xw_phase, xwt=xwt,
+    Results = collections.namedtuple('Results', ['xw_coherence', 'xw_amplitude', 'xw_phase', 'xw_t', 'freq', 'time', 'AR1_q', 'coi'])
+    res = Results(xw_coherence=xw_coherence, xw_amplitude=xw_amplitude, xw_phase=xw_phase, xw_t=xw_t,
                   freq=freq, time=tau, AR1_q=AR1_q, coi=coi)
 
     return res
@@ -1679,17 +1679,29 @@ def freq_vector_nfft(ts):
 
     return freq
 
-def freq_vector_scale(ts, nv=12, fourier_factor=1):
-    ''' Return the frequency vector based on scales for wavelet analysis
+def freq_vector_scale(ts, dj=0.25, s0=None,j1=None, mother='MORLET',param=None):
+    ''' Return the frequency vector based on scales for wavelet analysis. 
+    This function is adapted from Torrence and Compo
 
     Parameters
     ----------
 
-    ts : array
-        time axis of the time series
-
-    nv : int
-        the parameter that controls the number of freq points
+    ts: numpy.array
+        The time axis for the timeseries    
+    dj : float, optional
+        The spacing between discrete scales. The default is 0.25. A smaller number will give better scale resolution, but be slower to plot.
+    s0 : float, optional
+        the smallest scale of the wavelet. The default is None, representing 2*dT.
+    j1 : float, optional
+        the number of scales minus one. Scales range from S0 up to S0*2**(J1*DJ),
+        to give a total of (J1+1) scales. The default is None, which represents (LOG2(N DT/S0))/DJ.
+    mother : string, optional
+        the mother wavelet function. The default is 'MORLET'. Options are: 'MORLET', 'PAUL', or 'DOG'
+    param : flaot, optional
+        the mother wavelet parameter. The default is None since it varies for each mother
+            - For 'MORLET' this is k0 (wavenumber), default is 6.
+            - For 'PAUL' this is m (order), default is 4.
+            - For 'DOG' this is m (m-th derivative), default is 2.
 
     Returns
     -------
@@ -1710,14 +1722,44 @@ def freq_vector_scale(ts, nv=12, fourier_factor=1):
     pyleoclim.utils.wavelet.freq_vector_log : Return the frequency vector based on logspace
 
     pyleoclim.utils.wavelet.make_freq_vector : Make frequency vector
-
+    
+    References
+    ----------
+    
+    Torrence, C. and G. P. Compo, 1998: A Practical Guide to Wavelet Analysis. Bull. Amer. Meteor. Soc., 79, 61-78.
+    Python routines available at http://paos.colorado.edu/research/wavelets/
     '''
+    
+    if mother.upper() not in ['MORLET','DOG','PAUL']:
+        raise ValueError('The mother wavelet should be either "MORLET","PAUL", or "DOG"')
+    
+    dt = np.diff(ts).mean()
+    n1=len(ts)
+    
+    if s0 is None:
+        s0 = 2 * dt
+    if j1 is None:
+        j1 = np.fix((np.log(n1 * dt / s0) / np.log(2)) / dj)
+    
+    # construct SCALE array & empty PERIOD & WAVE arrays
+    if mother.upper() == 'MORLET':
+        if param == None:
+            param = 6.
+        fourier_factor = 4 * np.pi / (param + np.sqrt(2 + param**2))
+    elif mother.upper() == 'PAUL':
+        if param == None:
+            param = 4.
+        fourier_factor = 4 * np.pi / (2 * param + 1)
+    elif mother.upper() == 'DOG':
+        if param == None:
+            param = 2.
+        fourier_factor = 2 * np.pi * np.sqrt(2. / (2 * param + 1))
+    else:
+        fourier_factor = np.nan
 
-    s0 = 2*np.median(np.diff(ts))
-    a0 = 2**(1/nv)
-    noct = np.floor(np.log2(np.size(ts)))-1  # number of octave
-    scale = s0*a0**(np.arange(noct*nv+1))
-    freq = 1/(scale[::-1]*fourier_factor)
+    j = np.arange(0, j1 + 1)
+    scale = s0 * 2. ** (j * dj)
+    freq = 1. / (fourier_factor * scale)       
 
     return freq
 
@@ -2249,7 +2291,7 @@ def prepare_wwz(ys, ts, freq=None, freq_method='log', freq_kwargs=None, tau=None
 
     return ys_cut, ts_cut, freq, tau
 
-def cross_wt(coeff1, coeff2):
+def xwt(coeff1, coeff2):
     ''' Return the cross wavelet transform.
 
     Parameters
@@ -2267,6 +2309,8 @@ def cross_wt(coeff1, coeff2):
     Returns
     -------
 
+    xw_t : array
+        the cross wavelet transform complex number
     xw_amplitude : array
         the cross wavelet amplitude
     xw_phase : array
@@ -2279,13 +2323,13 @@ def cross_wt(coeff1, coeff2):
     wavelet coherence to geophysical time series. Nonlin. Processes Geophys. 11, 561–566 (2004).
 
     '''
-    xwt = coeff1 * np.conj(coeff2)
-    xw_amplitude = np.sqrt(xwt.real**2 + xwt.imag**2)
-    xw_phase = np.arctan2(xwt.imag, xwt.real)
+    xw_t = coeff1 * np.conj(coeff2)
+    xw_amplitude = np.sqrt(xw_t.real**2 + xw_t.imag**2)
+    xw_phase = np.arctan2(xw_t.imag, xw_t.real)
 
-    return xwt, xw_amplitude, xw_phase
+    return xw_t, xw_amplitude, xw_phase
 
-def wavelet_coherence(coeff1, coeff2, freq, tau, smooth_factor=0.25):
+def wtc(coeff1, coeff2, freq, tau, smooth_factor=0.25):
     ''' Return the cross wavelet coherence.
 
     Parameters
@@ -2459,214 +2503,553 @@ def reconstruct_ts(coeff, freq, tau, t, len_bd=0):
 
     return rec_ts, t
 
-# ## Methods for Torrence and compo
-
-# # This is the main function, which has been rewritten to work with functionalities in Pyleoclim
+############ Methods for Torrence and Compo#############
 
-# def cwt(ys,ts,mother='morlet',param=None,freq=None,freq_method='scale',
-#         freq_kwargs={},detrend=False, sg_kwargs={}, gaussianize=False,
-#         standardize=False,pad=False,pad_kwargs={}):
-
-#     ys=np.array(ys)
-#     ts=np.array(ts)
-
-#     ys, ts = clean_ts(ys, ts) #clean up time
+def cwt(ys,ts,freq=None,freq_method='log',freq_kwargs={}, scale = None, detrend=False,sg_kwargs={},
+        gaussianize=False, standardize=False, pad=False, mother='MORLET',param=None):
+    '''
+    Wrapper function to implement Torrence and Compo continuous wavelet transform
 
-#     #make sure that the time series is evenly-spaced
-#     if is_evenly_spaced(ts) == True:
-#         dt = np.mean(np.diff(ts))
-#     else:
-#         raise ValueError('Time series must be evenly spaced in time')
-
-#     # prepare the time series
-#     pd_ys = preprocess(ys, ts, detrend=detrend, sg_kwargs=sg_kwargs,
-#                        gaussianize=gaussianize, standardize=standardize)
-
-#     # Get the fourier factor
-#     if mother.lower() == 'morlet':
-#         if param is None:
-#             param = 6.
-#         fourier_factor = 4 * np.pi / (param + np.sqrt(2 + param**2))
-#     elif mother.lower() == 'paul':
-#         if param is None:
-#             param = 4.
-#         fourier_factor = 4 * np.pi / (2 * param + 1)
-#     elif mother.lower() == 'dog':
-#         if param is None:
-#             param = 2.
-#         fourier_factor = 2 * np.pi * np.sqrt(2. / (2 * param + 1))
-#     else:
-#         fourier_factor = 1
-
-#     #get the frequency/scale information
-#     if freq is None:
-#         freq_kwargs = {} if freq_kwargs is None else freq_kwargs.copy()
-#         if freq_method == 'scale':
-#             freq_kwargs.update({'fourier_factor':fourier_factor})
-#         freq = make_freq_vector(ts, method=freq_method, **freq_kwargs)
-#     # Use scales
-#     scale = np.sort(1/(freq*fourier_factor))
-
-#     #Normalize
-#     #n_ys = pd_ys-np.mean(pd_ys)
-
-#     #pad if wanted
-#     if pad == True:
-#         pad_kwargs = {} if pad_kwargs is None else pad_kwargs.copy()
-#         yp,tp = ts_pad(pd_ys,ts,**pad_kwargs)
-#     else:
-#         yp=pd_ys
-#         tp=ts
+    Parameters
+    ----------
+    ys : numpy.array
+        the time series.
+    ts : numpy.array
+        the time axis.
+    freq : numpy.array, optional
+        The frequency vector. The default is None, which will prompt the use of one the underlying functions
+    freq_method : string, optional
+        The method by which to obtain the frequency vector. The default is 'log'.
+        Options are 'log' (default), 'nfft', 'lomb_scargle', 'welch', and 'scale'
+    freq_kwargs : dict, optional
+        Optional parameters for the choice of the frequency vector. See make_freq_vector and additional methods for details. The default is {}.
+    scale : numpy.array
+        Optional scale vector in place of a frequency vector. Default is None. If scale is not None, frequency method and attached arguments will be ignored. 
+    detrend : bool, string, {'linear', 'constant', 'savitzy-golay', 'emd'}
+        Whether to detrend and with which option. The default is False.
+    sg_kwargs : dict, optional
+        Additional parameters for the savitzy-golay method. The default is {}.
+    gaussianize : bool, optional
+        Whether to gaussianize. The default is False.
+    standardize : bool, optional
+        Whether to standardize. The default is False.
+    pad : bool, optional
+        Whether or not to pad the timeseries. with zeroes to get N up to the next higher power of 2. 
+        This prevents wraparound from the end of the time series to the beginning, and also speeds up the FFT's used to do the wavelet transform.
+        This will not eliminate all edge effects. The default is False.
+    mother : string, optional
+        the mother wavelet function. The default is 'MORLET'. Options are: 'MORLET', 'PAUL', or 'DOG'
+    param : flaot, optional
+        the mother wavelet parameter. The default is None since it varies for each mother
+            - For 'MORLET' this is k0 (wavenumber), default is 6.
+            - For 'PAUL' this is m (order), default is 4.
+            - For 'DOG' this is m (m-th derivative), default is 2.
 
-#     # Wave calculation
-#     n = len(yp)
+    Returns
+    -------
+    res : dict
+        Dictionary containing:
+            - amplitude: the wavelet amplitude
+            - coi: cone of influence
+            - freq: frequency vector
+            - coeff: the wavelet coefficients
+            - scale: the scale vector
+            - time: the time vector
+            - mother: the mother wavelet
+            - param : the wavelet parameter
+            
+    See also
+    --------
+    
+    pyleoclim.utils.wavelet.make_freq_vector : make the frequency vector with various methods
+    
+    pyleoclim.utils.wavelet.tc_wavelet: the underlying wavelet function by Torrence and Compo
+    
+    pyleoclim.utils.tsutils.detrend : detrending functionalities in Pyleoclim
+    
+    References
+    ----------
+    
+    Torrence, C. and G. P. Compo, 1998: A Practical Guide to Wavelet Analysis. Bull. Amer. Meteor. Soc., 79, 61-78.
+    Python routines available at http://paos.colorado.edu/research/wavelets/
+
+    '''    
+    
+    ts = np.array(ts)
+    ys = np.array(ys)
+    
+    ys, ts = clean_ts(ys,ts)
+
+    if len(ts) != len(ys):
+        raise ValueError('Time and value axis should be the same length')
+    
+    if is_evenly_spaced(ts) == False:
+        raise ValueError('Time vector should be evenly spaced for this method. Interpolate or use WWZ.')
+    
+    if mother.upper() not in ['MORLET','DOG','PAUL']:
+        raise ValueError('The mother wavelet should be either "MORLET","PAUL", or "DOG"')
+        
+    #preprocessing
+    # remove NaNs
+    
+    dt = np.diff(ts).mean()
+
+    ys = preprocess(ys, ts, detrend=detrend, sg_kwargs=sg_kwargs,
+               gaussianize=gaussianize, standardize=standardize) #TC seems to require standardization
+    
+    # fourier factor determination
+    if mother.upper() == 'MORLET':
+        if param == None:
+            param = 6.
+        fourier_factor = 4 * np.pi / (param + np.sqrt(2 + param**2))
+    elif mother.upper() == 'PAUL':
+        if param == None:
+            param = 4.
+        fourier_factor = 4 * np.pi / (2 * param + 1)
+    elif mother.upper() == 'DOG':
+        if param == None:
+            param = 2.
+        fourier_factor = 2 * np.pi * np.sqrt(2. / (2 * param + 1))
+    else:
+        fourier_factor = np.nan
+    
+    #get the scale
+    if scale is None:
+        if freq is None:
+            if freq_method == 'scale':
+                freq_kwargs = {} if freq_kwargs is None else freq_kwargs.copy()
+                freq_kwargs.update({'mother':mother,'param':param})
+            freq = make_freq_vector(ts,method=freq_method,**freq_kwargs)
+        scale = 1. / (fourier_factor * freq)
+        
+    #calculate wavelet
+    wave, coi = tc_wavelet(ys, dt, scale, mother, param, pad)
+    amplitude=np.abs(wave)
+    
+    Results = collections.namedtuple('Results', ['amplitude', 'coi', 'freq', 'time', 'scale', 'coeff', 'mother','param'])
+    res = Results(amplitude=amplitude.T, coi=coi, freq=freq, time=ts, scale=scale, coeff=wave, mother=mother,param=param)
 
-#     # construct wavenumber array used in transform [Eqn(5)]
-#     kplus = np.arange(1, int(n / 2) + 1)
-#     kplus = (kplus * 2 * np.pi / (n * dt))
-#     kminus = np.arange(1, int((n - 1) / 2) + 1)
-#     kminus = np.sort((-kminus * 2 * np.pi / (n * dt)))
-#     k = np.concatenate(([0.], kplus, kminus))
+    return res
+    
 
-#     # compute FFT of the (padded) time series
-#     f = np.fft.fft(yp)
+def tc_wavelet(Y, dt, scale, mother, param, pad=False):
+    '''
+    WAVELET  1D Wavelet transform. Adapted from Torrence and Compo to fit existing Pyleoclim functionalities
+    
+    Computes the wavelet transform of the vector Y (length N),
+    with sampling rate DT.
 
-#     # define the wavelet array
-#     wave = np.zeros(shape=(len(scale), n), dtype=complex)
+    By default, the Morlet wavelet (k0=6) is used.
+    The wavelet basis is normalized to have total energy=1 at all scales.
 
-#     # loop through all scales and compute transform
-#     for a1 in range(0, len(scale)):
-#         daughter, fourier_factor, coi, _ = \
-#             wave_bases(mother, k, scale[a1], param)
-#         wave[a1, :] = np.fft.ifft(f * daughter)  # wavelet transform[Eqn(4)]
+    Parameters
+    ----------
+    Y : numpy.array
+        the time series of length N.
+    dt : float
+        the sampling time
+    
+    mother : string, optional
+        the mother wavelet function. The default is 'MORLET'. Options are: 'MORLET', 'PAUL', or 'DOG'
+    param : flaot, optional
+        the mother wavelet parameter. The default is None since it varies for each mother
+            - For 'MORLET' this is k0 (wavenumber), default is 6.
+            - For 'PAUL' this is m (order), default is 4.
+            - For 'DOG' this is m (m-th derivative), default is 2.
+    pad : {True,False}, optional
+        Whether or not to pad the timeseries. with zeroes to get N up to the next higher power of 2. 
+        This prevents wraparound from the end of the time series to the beginning, and also speeds up the FFT's used to do the wavelet transform.
+        This will not eliminate all edge effects. The default is False.
 
-#     #COI
-#     coi = coi * dt * np.concatenate((
-#         np.insert(np.arange(int((len(ys) + 1) / 2) - 1), [0], [1E-5]),
-#         np.insert(np.flipud(np.arange(0, int(len(ys) / 2) - 1)), [-1], [1E-5])))
+    Returns
+    -------
+    wave : numpy.array
+        The wavelet coefficients
+    coi : numpy.array
+        The cone of influence. Periods greater than this are subject to edge effects.
+        
+    See also
+    --------
+    
+    pyleoclim.utils.wavelet.tc_wave_bases: 1D wavelet functions Morlet, Paul or Dog
+    
+    References
+    ----------
+    
+    Torrence, C. and G. P. Compo, 1998: A Practical Guide to Wavelet Analysis. Bull. Amer. Meteor. Soc., 79, 61-78.
+    Python routines available at http://paos.colorado.edu/research/wavelets/
+    
+    '''
+    
+    n1 = len(Y)
+
+    # construct time series to analyze, pad if necessary
+    x = Y - np.mean(Y)
+    if pad == True:
+        # power of 2 nearest to N
+        base2 = np.fix(np.log(n1) / np.log(2) + 0.4999)
+        nzeroes = (2 ** (base2 + 1) - n1).astype(np.int64)
+        x = np.concatenate((x, np.zeros(nzeroes)))
+
+    n = len(x)
+
+    # construct wavenumber array used in transform [Eqn(5)]
+    kplus = np.arange(1, int(n / 2) + 1)
+    kplus = (kplus * 2 * np.pi / (n * dt))
+    kminus = np.arange(1, int((n - 1) / 2) + 1)
+    kminus = np.sort((-kminus * 2 * np.pi / (n * dt)))
+    k = np.concatenate(([0.], kplus, kminus))
+
+    # compute FFT of the (padded) time series
+    f = np.fft.fft(x)  # [Eqn(3)]
+    
+    # define the wavelet array
+    wave = np.zeros(shape=(len(scale), n), dtype=complex)
+
+    # loop through all scales and compute transform
+    for a1 in range(0, len(scale)):
+        daughter, fourier_factor, coi, _ = \
+            tc_wave_bases(mother, k, scale[a1], param)
+        wave[a1, :] = np.fft.ifft(f * daughter)  # wavelet transform[Eqn(4)]
+
+    # COI [Sec.3g]
+    coi = coi * dt * np.concatenate((
+        np.insert(np.arange(int((n1 + 1) / 2) - 1), [0], [1E-5]),
+        np.insert(np.flipud(np.arange(0, int(n1 / 2) - 1)), [-1], [1E-5])))
+    wave = wave[:, :n1]  # get rid of padding before returning
+
+    return wave, coi
+
+def tc_wave_bases(mother, k, scale, param):
+    '''
+    WAVE_BASES  1D Wavelet functions Morlet, Paul, or DOG
 
-#     #Remove the padding
-#     if pad == True:
-#         idx = np.in1d(tp,ts)
-#         wave = wave[:,idx]
+    Parameters
+    ----------
+    mother : string 
+        equal to 'MORLET' or 'PAUL' or 'DOG'
+    k : numpy.array
+        the Fourier frequencies at which to calculate the wavelet
+    scale : float
+        The wavelet scale
+    param : float
+        the nondimensional parameter for the wavelet function
 
-#     res = {}
+    Returns
+    -------
+    daughter : numpy.array
+        a vector, the wavelet function
+    fourier_factor : float
+        the ratio of Fourier period to scale
+    coi : float
+        the cone-of-influence size at the scale
+    dofmin : float
+        degrees of freedom for each point in the wavelet power
+             (either 2 for Morlet and Paul, or 1 for the DOG)
+    
+    References
+    ----------
+    
+    Torrence, C. and G. P. Compo, 1998: A Practical Guide to Wavelet Analysis. Bull. Amer. Meteor. Soc., 79, 61-78.
+    Python routines available at http://paos.colorado.edu/research/wavelets/
 
-#     return res
+    '''
+    
+    
+    n = len(k)
+    kplus = np.array(k > 0., dtype=float)
+
+    if mother == 'MORLET':  # -----------------------------------  Morlet
+
+        if param == -1:
+            param = 6.
+
+        k0 = np.copy(param)
+        # calc psi_0(s omega) from Table 1
+        expnt = -(scale * k - k0) ** 2 / 2. * kplus
+        norm = np.sqrt(scale * k[1]) * (np.pi ** (-0.25)) * np.sqrt(n)
+        daughter = norm * np.exp(expnt)
+        daughter = daughter * kplus  # Heaviside step function
+        # Scale-->Fourier [Sec.3h]
+        fourier_factor = (4 * np.pi) / (k0 + np.sqrt(2 + k0 ** 2))
+        coi = fourier_factor / np.sqrt(2)  # Cone-of-influence [Sec.3g]
+        dofmin = 2  # Degrees of freedom
+        
+    elif mother == 'PAUL':  # --------------------------------  Paul
+        if param == -1:
+            param = 4.
+        m = param
+        # calc psi_0(s omega) from Table 1
+        expnt = -scale * k * kplus
+        norm_bottom = np.sqrt(m * np.prod(np.arange(1, (2 * m))))
+        norm = np.sqrt(scale * k[1]) * (2 ** m / norm_bottom) * np.sqrt(n)
+        daughter = norm * ((scale * k) ** m) * np.exp(expnt) * kplus
+        fourier_factor = 4 * np.pi / (2 * m + 1)
+        coi = fourier_factor * np.sqrt(2)
+        dofmin = 2
+        
+    elif mother == 'DOG':  # --------------------------------  DOG
+        if param == -1:
+            param = 2.
+        m = param
+        # calc psi_0(s omega) from Table 1
+        expnt = -(scale * k) ** 2 / 2.0
+        norm = np.sqrt(scale * k[1] / gamma(m + 0.5)) * np.sqrt(n)
+        daughter = -norm * (1j ** m) * ((scale * k) ** m) * np.exp(expnt)
+        fourier_factor = 2 * np.pi * np.sqrt(2. / (2 * m + 1))
+        coi = fourier_factor / np.sqrt(2)
+        dofmin = 1
+    else:
+        print('Mother must be one of MORLET, PAUL, DOG')
 
+    return daughter, fourier_factor, coi, dofmin
 
-# def wave_bases(mother, k, scale, param):
-#     '''
+def tc_wave_signif(ys, ts, scale, mother, param, sigtest='chi-square', qs=[0.95],
+                dof=None, gws=None):
+    '''
+    Asymptotic singificance testing.
 
+    Parameters
+    ----------
+    ys : numpy.array
+        Values for the timeseries
+    ts : numpy.array
+        time vector.
+    scale : numpy.array
+        vector of scale
+    mother : str
+        Type of mother wavelet
+    param : int
+        mother wavelet parameter
+    sigtest : {'chi-square','time-average','scale-average'}, optional
+        Type of significance test to perform . The default is 'chi-square'.
+        - chi-square: a regular chi-square test Eq(18) from Torrence and Compo
+        - time-average: DOF should be set to NA, the number of local wavelet spectra that were averaged together.
+             For the Global Wavelet Spectrum, this would be NA=N, where N is the number of points in your time series. Eq23 in Torrence and Compo
+        -scale-average: In this case, DOF should be set to a two-element vector [S1,S2], which gives the scale range that was averaged together.
+             e.g. if one scale-averaged scales between 2 and 8, then DOF=[2,8].
+    qs : list, optional
+        Significance level. The default is [0.95].
+    dof : None, optional
+        Degrees of freedon for signif test. The default is None, which will automatically assign:
+            - chi-square: DOF = 2 (or 1 for MOTHER='DOG')
+            - time-average: DOF = NA, the number of times averaged together.
+            -scale-average: DOF = [S1,S2], the range of scales averaged.
+    gws : np.array, optional
+        Global wavelet spectrum. a vector of the same length as scale. If input then this is used as the theoretical background spectrum, rather than white or red noise. The default is None.
 
-#     Parameters
-#     ----------
-#     mother : string, {}
-#         DESCRIPTION.
-#     k : TYPE
-#         DESCRIPTION.
-#     scale : TYPE
-#         DESCRIPTION.
-#     param : TYPE
-#         DESCRIPTION.
-
-#     Raises
-#     ------
-#     KeyError
-#         DESCRIPTION.
+    Returns
+    -------
+    signif_level : numpy.array
+        Array of values for significance level 
+    
+    References
+    ----------
+    
+    Torrence, C. and G. P. Compo, 1998: A Practical Guide to Wavelet Analysis. Bull. Amer. Meteor. Soc., 79, 61-78.
+    Python routines available at http://paos.colorado.edu/research/wavelets/
+    
+    See also
+    --------
+    
+    pyleoclim.utils.wavelet.chisquare_inv : inverse of chi-square CDF
+    
+    pyleoclim.utils.wavelet.chisquare_solve : return the difference between calculated percentile and true P
 
-#     Returns
-#     -------
-#     daughter : TYPE
-#         DESCRIPTION.
-#     fourier_factor : TYPE
-#         DESCRIPTION.
-#     coi : TYPE
-#         DESCRIPTION.
-#     dofmin : TYPE
-#         DESCRIPTION.
 
-#     '''
+    '''
+    
+    if mother.upper() not in ['MORLET','DOG','PAUL']:
+        raise ValueError('The mother wavelet should be either "MORLET","PAUL", or "DOG"')
+    
+    if sigtest not in ['chi-square','time-average','scale-average']:
+        raise ValueError("The type of significance test should be either 'chi-square','time-average','scale-average'")
+    
+    
+    J1 = len(scale) - 1
+    dj = np.log2(scale[1] / scale[0])
+
+    variance = np.std(ys,ddof=1) ** 2
+
+    # get the appropriate parameters [see Table(2)]
+    if mother.upper() == 'MORLET':  # ----------------------------------  Morlet
+        empir = ([2., -1, -1, -1])
+        if param == 6:
+            empir[1:] = ([0.776, 2.32, 0.60])
+        k0 = param
+        # Scale-->Fourier [Sec.3h]
+        fourier_factor = (4 * np.pi) / (k0 + np.sqrt(2 + k0 ** 2))
+        
+    elif mother.upper() == 'PAUL':
+        empir = ([2, -1, -1, -1])
+        if param == 4:
+            empir[1:] = ([1.132, 1.17, 1.5])
+        m = param
+        fourier_factor = (4 * np.pi) / (2 * m + 1)
+        
+    elif mother.upper() == 'DOG':  # -------------------------------------Paul
+        empir = ([1., -1, -1, -1])
+        if param ==2:
+            empir[1:] = ([3.541, 1.43, 1.4])
+        elif param == 6:  # --------------------------------------DOG
+            empir[1:] = ([1.966, 1.37, 0.97])
+        m = param
+        fourier_factor = 2 * np.pi * np.sqrt(2. / (2 * m + 1))
+
+
+    period = scale * fourier_factor
+    dofmin = empir[0]  # Degrees of freedom with no smoothing
+    Cdelta = empir[1]  # reconstruction factor
+    gamma_fac = empir[2]  # time-decorrelation factor
+    dj0 = empir[3]  # scale-decorrelation factor
+
+    dt = float(np.mean(np.diff(ts)))
+    freq =  dt / period  # normalized frequency
+
+    if gws is not None:   # use global-wavelet as background spectrum
+        fft_theor = gws
+    else:
+        # [Eqn(16)]
+        lag1 = ar1_fit(ys,ts)
+        fft_theor = (1 - lag1 ** 2) / \
+            (1 - 2 * lag1 * np.cos(freq * 2 * np.pi) + lag1 ** 2)
+        fft_theor = variance * fft_theor  # include time-series variance
+
+    signif = fft_theor
+    
+    if dof is None:
+        dof = dofmin
+
+    signif_level = []
+    
+    for siglvl in qs:
+            
+        if sigtest == 'chi-square':  # no smoothing, DOF=dofmin [Sec.4]
+            dof = dofmin
+            chisquare = chisquare_inv(siglvl, dof) / dof
+            signif = fft_theor * chisquare  # [Eqn(18)]
+            
+        #expand
+            #signif = signif[:, np.newaxis].dot(np.ones(len(ys))[np.newaxis, :])
+            
+        
+        elif sigtest == 'time-average':  # time-averaged significance
+            if len(np.atleast_1d(dof)) == 1:
+                dof = np.zeros(J1) + dof
+            dof[dof < 1] = 1
+            # [Eqn(23)]
+            dof = dofmin * np.sqrt(1 + (dof * dt / gamma_fac / scale) ** 2)
+            dof[dof < dofmin] = dofmin   # minimum DOF is dofmin
+            for a1 in range(0, J1 + 1):
+                chisquare = chisquare_inv(siglvl, dof[a1]) / dof[a1]
+                signif[a1] = fft_theor[a1] * chisquare
+        
+        elif sigtest == 'scale-average':  # time-averaged significance
+            if len(dof) != 2:
+                raise ValueError('DOF must be set to [S1,S2],'
+                    ' the range of scale-averages')
+            if Cdelta == -1:
+                raise ValueError('Cdelta & dj0 not defined'
+                      ' for ' + mother + ' with param = ' + str(param))
+    
+            s1 = dof[0]
+            s2 = dof[1]
+            avg = np.logical_and(scale >= 2, scale < 8)  # scales between S1 & S2
+            navg = np.sum(np.array(np.logical_and(scale >= 2, scale < 8),
+                dtype=int))
+            if navg == 0:
+                raise ValueError('No valid scales between ' + s1 + ' and ' + s2)
+            Savg = 1. / np.sum(1. / scale[avg])  # [Eqn(25)]
+            Smid = np.exp((np.log(s1) + np.log(s2)) / 2.)  # power-of-two midpoint
+            dof = (dofmin * navg * Savg / Smid) * \
+                np.sqrt(1 + (navg * dj / dj0) ** 2)  # [Eqn(28)]
+            fft_theor = Savg * np.sum(fft_theor[avg] / scale[avg])  # [Eqn(27)]
+            chisquare = chisquare_inv(siglvl, dof) / dof
+            signif = (dj * dt / Cdelta / Savg) * fft_theor * chisquare  # [Eqn(26)]
+        
+        signif_level.append(np.sqrt(signif.T))
+        
+    return signif_level
+
+def chisquare_inv(P, V):
+    '''
+    Returns the inverse of chi-square CDF with V degrees of freedom at fraction P
+    
 
-#     n = len(k)
-#     kplus = np.array(k > 0., dtype=float)
-
-#     if mother == 'morlet':  # -----------------------------------  Morlet
-
-#         if param == -1:
-#             param = 6.
-
-#         k0 = np.copy(param)
-#         # calc psi_0(s omega) from Table 1
-#         expnt = -(scale * k - k0) ** 2 / 2. * kplus
-#         norm = np.sqrt(scale * k[1]) * (np.pi ** (-0.25)) * np.sqrt(n)
-#         daughter = norm * np.exp(expnt)
-#         daughter = daughter * kplus  # Heaviside step function
-#         # Scale-->Fourier [Sec.3h]
-#         fourier_factor = (4 * np.pi) / (k0 + np.sqrt(2 + k0 ** 2))
-#         coi = fourier_factor / np.sqrt(2)  # Cone-of-influence [Sec.3g]
-#         dofmin = 2  # Degrees of freedom
-#     elif mother == 'paul':  # --------------------------------  Paul
-#         if param == -1:
-#             param = 4.
-#         m = param
-#         # calc psi_0(s omega) from Table 1
-#         expnt = -scale * k * kplus
-#         norm_bottom = np.sqrt(m * np.prod(np.arange(1, (2 * m))))
-#         norm = np.sqrt(scale * k[1]) * (2 ** m / norm_bottom) * np.sqrt(n)
-#         daughter = norm * ((scale * k) ** m) * np.exp(expnt) * kplus
-#         fourier_factor = 4 * np.pi / (2 * m + 1)
-#         coi = fourier_factor * np.sqrt(2)
-#         dofmin = 2
-#     elif mother == 'dog':  # --------------------------------  DOG
-#         if param == -1:
-#             param = 2.
-#         m = param
-#         # calc psi_0(s omega) from Table 1
-#         expnt = -(scale * k) ** 2 / 2.0
-#         norm = np.sqrt(scale * k[1] / gamma(m + 0.5)) * np.sqrt(n)
-#         daughter = -norm * (1j ** m) * ((scale * k) ** m) * np.exp(expnt)
-#         fourier_factor = 2 * np.pi * np.sqrt(2. / (2 * m + 1))
-#         coi = fourier_factor / np.sqrt(2)
-#         dofmin = 1
-#     else:
-#         raise KeyError('Mother must be one of "morlet", "paul", "dog"')
+    Parameters
+    ----------
+    P : float
+        fraction
+    V : float
+        degress of freedom
 
-#     return daughter, fourier_factor, coi, dofmin
+    Returns
+    -------
+    X : float
+        Inverse chi-square
+    
+    References
+    ----------
+    
+    Torrence, C. and G. P. Compo, 1998: A Practical Guide to Wavelet Analysis. Bull. Amer. Meteor. Soc., 79, 61-78.
+    Python routines available at http://paos.colorado.edu/research/wavelets/
 
+    '''
+    
+    if (1 - P) < 1E-4:
+        raise ValueError('P must be < 0.9999')
 
-# def chisquare_inv(P, V):
+    if P == 0.95 and V == 2:  # this is a no-brainer
+        X = 5.9915
+        return X
 
-#     if (1 - P) < 1E-4:
-#         print('P must be < 0.9999')
+    MINN = 0.01  # hopefully this is small enough
+    MAXX = 1  # actually starts at 10 (see while loop below)
+    X = 1
+    TOLERANCE = 1E-4  # this should be accurate enough
 
-#     if P == 0.95 and V == 2:  # this is a no-brainer
-#         X = 5.9915
-#         return X
+    while (X + TOLERANCE) >= MAXX:  # should only need to loop thru once
+        MAXX = MAXX * 10.
+    # this calculates value for X, NORMALIZED by V
+        X = fminbound(chisquare_solve, MINN, MAXX, args=(P, V), xtol=TOLERANCE)
+        MINN = MAXX
 
-#     MINN = 0.01  # hopefully this is small enough
-#     MAXX = 1  # actually starts at 10 (see while loop below)
-#     X = 1
-#     TOLERANCE = 1E-4  # this should be accurate enough
+    X = X * V  # put back in the goofy V factor
 
-#     while (X + TOLERANCE) >= MAXX:  # should only need to loop thru once
-#         MAXX = MAXX * 10.
-#     # this calculates value for X, NORMALIZED by V
-#         X = fminbound(chisquare_solve, MINN, MAXX, args=(P, V), xtol=TOLERANCE)
-#         MINN = MAXX
+    return X 
 
-#     X = X * V  # put back in the goofy V factor
+def chisquare_solve(XGUESS, P, V):
+    '''
+    Given XGUESS, a percentile P, and degrees-of-freedom V, return the difference between calculated percentile and P.
 
-#     return X
+    Parameters
+    ----------
+    XGUESS : float
+        sig level
+    P : float
+        percentile
+    V : float
+        degrees of freedom
 
-# def chisquare_solve(XGUESS, P, V):
+    Returns
+    -------
+    PDIFF : float
+        difference between calculated percentile and P
+        
+    References
+    ----------
+    
+    Torrence, C. and G. P. Compo, 1998: A Practical Guide to Wavelet Analysis. Bull. Amer. Meteor. Soc., 79, 61-78.
+    Python routines available at http://paos.colorado.edu/research/wavelets/
 
-#     PGUESS = gammainc(V / 2, V * XGUESS / 2)  # incomplete Gamma function
+    '''
+    
+    PGUESS = gammainc(V / 2, V * XGUESS / 2)  # incomplete Gamma function
 
-#     PDIFF = np.abs(PGUESS - P)            # error in calculated P
+    PDIFF = np.abs(PGUESS - P)            # error in calculated P
 
-#     TOL = 1E-4
-#     if PGUESS >= 1 - TOL:  # if P is very close to 1 (i.e. a bad guess)
-#         PDIFF = XGUESS   # then just assign some big number like XGUESS
+    TOL = 1E-4
+    if PGUESS >= 1 - TOL:  # if P is very close to 1 (i.e. a bad guess)
+        PDIFF = XGUESS   # then just assign some big number like XGUESS
 
-#     return PDIFF
+    return PDIFF


Test patch:
diff --git a/pyleoclim/tests/test_ui_MultipleSeries.py b/pyleoclim/tests/test_ui_MultipleSeries.py
index 789eb663..f26791f9 100644
--- a/pyleoclim/tests/test_ui_MultipleSeries.py
+++ b/pyleoclim/tests/test_ui_MultipleSeries.py
@@ -316,13 +316,13 @@ def test_pca_t3(self):
         assert abs(res.pctvar.sum() - 100)<0.001 
         
  
-class TestMultipleSeriesGridProperties:
-    '''Test for MultipleSeries.grid_properties()
+class TestMultipleSeriesIncrements:
+    '''Test for MultipleSeries.increments()
     
     '''
     @pytest.mark.parametrize('step_style', ['min', 'max', 'mean', 'median'])
-    def test_grid_properties(self, step_style):
-        p = 10; n = 100
+    def test_increments(self, step_style):
+        p = 2; n = 100
         signal = pyleo.gen_ts(model='colored_noise',nt=n,alpha=1.0).standardize() 
         X = signal.value[:,None] + np.random.randn(n,p)
         t = np.arange(n)
@@ -332,7 +332,7 @@ def test_grid_properties(self, step_style):
             mslist.append(pyleo.Series(time = t, value = X[:,i]))
         ms = pyleo.MultipleSeries(mslist)
         
-        gp = ms.grid_properties(step_style=step_style)
+        gp = ms.increments(step_style=step_style)
         
         assert (gp[0,:] == np.array((t.min(), t.max(), 1.))).all()
                
@@ -420,17 +420,20 @@ def test_StackPlot_t1(self, plot_kwargs):
 class TestMultipleSeriesSpectral():
     ''' Test for MultipleSeries.spectral
     '''
-    
-    def test_spectral_t0(self):
+    @pytest.mark.parametrize('spec_method',['wwz','cwt'])
+    def test_spectral_t0(self,spec_method):
         '''Test the spectral function with pre-generated scalogram objects
         '''
         
         d=load_data()
         sst = d.to_LipdSeries(number=5)
         d18Osw = d.to_LipdSeries(number=3)
+        if spec_method == 'cwt':
+            sst=sst.interp()
+            d18Osw=d18Osw.interp()
         ms = pyleo.MultipleSeries([sst,d18Osw])
-        scals = ms.wavelet()
-        psds = ms.spectral(method='wwz',scalogram_list=scals)
+        scals = ms.wavelet(method=spec_method)
+        psds = ms.spectral(method=spec_method,scalogram_list=scals)
         
         
         
diff --git a/pyleoclim/tests/test_ui_Scalogram.py b/pyleoclim/tests/test_ui_Scalogram.py
index f6900fd5..cccda4db 100644
--- a/pyleoclim/tests/test_ui_Scalogram.py
+++ b/pyleoclim/tests/test_ui_Scalogram.py
@@ -47,12 +47,23 @@ def gen_colored_noise(alpha=1, nt=100, f0=None, m=None, seed=None):
 class TestUiScalogramSignifTest:
     ''' Tests for Scalogram.signif_test()
     '''
-
-    def test_signif_test_t0(self):
+    
+    @pytest.mark.parametrize('wave_method',['wwz','cwt'])
+    def test_signif_test_t0(self, wave_method):
         ''' Test scalogram.signif_test() with default parameters
         '''
         alpha = 1
         t, v = gen_colored_noise(nt=100, alpha=alpha)
         ts = pyleo.Series(time=t, value=v)
-        scal = ts.wavelet()
+        scal = ts.wavelet(method=wave_method)
         scal_signif = scal.signif_test(number=1)
+    
+    @pytest.mark.parametrize('ar1_method',['ar1asym', 'ar1sim'])
+    def test_signif_test_t1(self,ar1_method):
+        ''' Test scalogram.signif_test() with default parameters
+        '''
+        alpha = 1
+        t, v = gen_colored_noise(nt=100, alpha=alpha)
+        ts = pyleo.Series(time=t, value=v)
+        scal = ts.wavelet(method='cwt')
+        scal_signif = scal.signif_test(method=ar1_method,number=1)
diff --git a/pyleoclim/tests/test_ui_Series.py b/pyleoclim/tests/test_ui_Series.py
index 1ed43a18..802d5237 100644
--- a/pyleoclim/tests/test_ui_Series.py
+++ b/pyleoclim/tests/test_ui_Series.py
@@ -54,14 +54,11 @@ def gen_colored_noise(alpha=1, nt=100, f0=None, m=None, seed=None):
     return t, v
     
 def load_data():
-    # note: JEG swapped the two choices on 03/07/2022 to get through PaleoHack 3. 
-    # Awaits fixing of jsonutils (issue 209)
     try:
-        d = pyleo.utils.jsonutils.json_to_Scalogram('../../example_data/scal_signif_benthic.json')
-    except:
         url = 'https://raw.githubusercontent.com/LinkedEarth/Pyleoclim_util/Development/example_data/scal_signif_benthic.json'
-        response = urlopen(url)
-        d = json.loads(response.read())
+        d = pyleo.utils.jsonutils.json_to_Scalogram(url)
+    except:
+        d = pyleo.utils.jsonutils.json_to_Scalogram('../../example_data/scal_signif_benthic.json')
     return d
 
 # Tests below
@@ -136,7 +133,7 @@ class TestUiSeriesSpectral:
     [parametrizing tests](https://docs.pytest.org/en/stable/example/parametrize.html).
     '''
 
-    @pytest.mark.parametrize('spec_method', ['wwz', 'mtm', 'lomb_scargle', 'welch', 'periodogram'])
+    @pytest.mark.parametrize('spec_method', ['wwz', 'mtm', 'lomb_scargle', 'welch', 'periodogram','cwt'])
     def test_spectral_t0(self, spec_method, eps=0.5):
         ''' Test Series.spectral() with available methods using default arguments
 
@@ -176,8 +173,8 @@ def test_spectral_t2(self, nfreq, eps=0.3):
         beta = psd.beta_est().beta_est_res['beta']
         assert np.abs(beta-alpha) < eps
 
-    @pytest.mark.parametrize('nv', [10, 20, 30])
-    def test_spectral_t3(self, nv, eps=0.3):
+    @pytest.mark.parametrize('dj', [0.25, 0.5, 1])
+    def test_spectral_t3(self, dj, eps=0.3):
         ''' Test Series.spectral() with MTM using `freq_method='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.
@@ -185,7 +182,7 @@ def test_spectral_t3(self, nv, eps=0.3):
         alpha = 1
         t, v = gen_colored_noise(nt=500, alpha=alpha)
         ts = pyleo.Series(time=t, value=v)
-        psd = ts.spectral(method='mtm', freq_method='scale', freq_kwargs={'nv': nv})
+        psd = ts.spectral(method='mtm', freq_method='scale', freq_kwargs={'dj': dj})
         beta = psd.beta_est().beta_est_res['beta']
         assert np.abs(beta-alpha) < eps
 
@@ -236,15 +233,16 @@ def test_spectral_t6(self, spec_method, eps=0.3):
         psd = ts.spectral(method=spec_method)
         beta = psd.beta_est().beta_est_res['beta']
         assert np.abs(beta-alpha) < eps
-        
-    def test_spectral_t7(self):
+    
+    @pytest.mark.parametrize('spec_method', ['wwz','cwt'])
+    def test_spectral_t7(self, spec_method,):
         '''Test the spectral significance testing with pre-generated scalogram objects
         '''
         
         ts = pyleo.gen_ts(model='colored_noise')
-        scal = ts.wavelet()
+        scal = ts.wavelet(method=spec_method)
         signif = scal.signif_test(number=2,export_scal = True)
-        sig_psd = ts.spectral(method='wwz',scalogram=scal)
+        sig_psd = ts.spectral(method=spec_method,scalogram=scal)
         sig_psd.signif_test(number=2,scalogram=signif).plot()
 
 class TestUiSeriesBin:
@@ -487,7 +485,7 @@ def test_summary_plot_t1(self):
         Passing just a pre generated psd.
         '''
         scal = load_data()
-        ts = scal.__dict__['timeseries']
+        ts = scal.timeseries
         fig, ax = ts.summary_plot()
     
         plt.close(fig)  
@@ -499,7 +497,7 @@ def test_summary_plot_t2(self):
         Passing just a pre generated psd.
         '''
         scal = load_data()
-        ts = scal.__dict__['timeseries']
+        ts = scal.timeseries
         fig, ax = ts.summary_plot(
             scalogram = scal
         )
@@ -513,7 +511,7 @@ def test_summary_plot_t3(self):
         Passing just a pre generated psd.
         '''
         scal = load_data()
-        ts = scal.__dict__['timeseries']
+        ts = scal.timeseries
         psd = ts.spectral(scalogram=scal)
         fig, ax = ts.summary_plot(
             psd = psd
@@ -528,7 +526,7 @@ def test_summary_plot_t4(self):
         Passing just a pre generated psd.
         '''
         scal = load_data()
-        ts = scal.__dict__['timeseries']
+        ts = scal.timeseries
         fig, ax = ts.summary_plot(
             scalogram = scal, psd_method='lomb_scargle'
         )
@@ -788,8 +786,7 @@ class TestUISeriesWavelet():
     ''' Test the wavelet functionalities
     '''
 
-    #@pytest.mark.parametrize('wave_method',['wwz','cwt'])
-    @pytest.mark.parametrize('wave_method',['wwz'])
+    @pytest.mark.parametrize('wave_method',['wwz','cwt'])
     def test_wave_t0(self, wave_method):
         ''' Test Series.wavelet() with available methods using default arguments
         '''
@@ -798,8 +795,7 @@ def test_wave_t0(self, wave_method):
         ts = pyleo.Series(time=t, value=v)
         scal = ts.wavelet(method=wave_method)
 
-    #@pytest.mark.parametrize('wave_method',['wwz','cwt'])
-    @pytest.mark.parametrize('wave_method',['wwz'])
+    @pytest.mark.parametrize('wave_method',['wwz','cwt'])
     def test_wave_t1(self,wave_method):
         '''Test Series.spectral() with WWZ/cwt with specified frequency vector passed via `settings`
         '''