Repository: dodge-research-group/thztools (language: python)

Issue:
41 doc unclarity in the documentation
Resolves #41.


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/draft-pdf.yml b/.github/workflows/draft-pdf.yml
new file mode 100644
index 0000000..c4f5f2f
--- /dev/null
+++ b/.github/workflows/draft-pdf.yml
@@ -0,0 +1,23 @@
+on: [push]
+
+jobs:
+  paper:
+    runs-on: ubuntu-latest
+    name: Paper Draft
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v4
+      - name: Build draft PDF
+        uses: openjournals/openjournals-draft-action@master
+        with:
+          journal: joss
+          # This should be the path to the paper within your repo.
+          paper-path: paper/paper.md
+      - name: Upload
+        uses: actions/upload-artifact@v4
+        with:
+          name: paper
+          # This is the output path where Pandoc will write the compiled
+          # PDF. Note, this should be the same directory as the input
+          # paper.md
+          path: paper/paper.pdf
\ No newline at end of file
diff --git a/.gitignore b/.gitignore
index d81ebdf..557dd9a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,7 +4,11 @@
 
 # Local .gitignore
 Untitled*.ipynb
-**/correspondence/
+**/doctest/**
+**/tmp/**
+
+# NumPy data files
+**/*.npz
 
 # Python .gitignore:
 #	https://github.com/github/gitignore/blob/main/Python.gitignore
diff --git a/.idea/thztools.iml b/.idea/thztools.iml
index d38d46c..7cf3cbb 100644
--- a/.idea/thztools.iml
+++ b/.idea/thztools.iml
@@ -7,12 +7,14 @@
       <excludeFolder url="file://$MODULE_DIR$/.pytest_cache" />
       <excludeFolder url="file://$MODULE_DIR$/dist" />
       <excludeFolder url="file://$MODULE_DIR$/docs/build" />
-      <excludeFolder url="file://$MODULE_DIR$/examples/.ipynb_checkpoints" />
       <excludeFolder url="file://$MODULE_DIR$/thztools.egg-info" />
       <excludeFolder url="file://$MODULE_DIR$/.idea/dictionaries" />
       <excludeFolder url="file://$MODULE_DIR$/.mypy_cache/3.11" />
       <excludeFolder url="file://$MODULE_DIR$/.ruff_cache/content" />
       <excludeFolder url="file://$MODULE_DIR$/htmlcov" />
+      <excludeFolder url="file://$MODULE_DIR$/docs/doctest/.doctrees" />
+      <excludeFolder url="file://$MODULE_DIR$/docs/doctest/plot_directive" />
+      <excludeFolder url="file://$MODULE_DIR$/docs/source/examples/.ipynb_checkpoints" />
     </content>
     <orderEntry type="jdk" jdkName="thztools" jdkType="Python SDK" />
     <orderEntry type="sourceFolder" forTests="false" />
diff --git a/CHANGELOG.md b/CHANGELOG.md
index f332b95..1135f92 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,8 +2,62 @@
 
 [![Common Changelog](https://common-changelog.org/badge.svg)](https://common-changelog.org)
 
-Please follow the guidelines at [Common Changelog](https://common-changelog.org) for maintaining
-this file.
+Please follow the guidelines at [Common Changelog](https://common-changelog.org)
+for maintaining this file.
+
+## [0.4.0] - 2024-05-09
+
+### Changed
+
+- Rename module constant `global_options` as `options`
+- Rename `tdnoisefit` as `noisefit`
+- Rename `NoiseModel` attributes `alpha`, `beta` and `tau` as `sigma_alpha`,
+  `sigma_beta` and `sigma_tau`, respectively
+- Rename `NoiseModel` methods `amplitude`, `noise`, and `variance` as
+  `noise_amp`, `noise_sim`, and `noise_var`, respectively
+- Revise `noisefit` and associated functions to treat each of the 3 noise
+  parameters separately instead of as a 3-component array, so that each may
+  fixed independently in the optimization process
+- Adjust internal scaling in `noisefit`
+- Reorganize `NoiseResult` class as output for `noisefit`
+- Revise `NoiseResult` docstring
+- Rename `transfer_out` as `transfer`
+- Rename `fft_sign` parameter as `numpy_sign_convention` in `transfer`
+- Make `dt` and `t0` keyword-only arguments in `wave` and set defaults
+- Change position of `dt` in `scaleshift`, `tdnoisefit`, `transfer` docstrings
+  and change defaults
+- Revise defaults for `wave`
+- Change default value for `dt` in `fit`
+- Improve type annotations
+- Update tests
+- Update dependencies in `environment-dev.yml` and `pyproject.toml`
+- Revise *Getting Started* page
+- Revise documentation, including Sphinx format and project-specific layout
+
+### Added
+
+- Add R. P. Hall, Laleh Mohtashemi and Naod Yimam to the author list
+- Add missing ORCID IDs
+- Add `get_option` and `set_option` functions to handle global options
+- Add attributes to `NoiseResult` class
+- Add `timebase` function
+- Add `seed` parameter in `NoiseModel.noise_sim`
+- Add `NoiseResult` class to documentation
+- Add doctests
+- Add warnings
+- Add `paper` directory with bib-file and stub for JOSS paper
+- Add `draft-pdf.yml` GitHub Action to autogenerate JOSS paper; update
+  deprecated `upload-artifact` Action in workflow
+- Add Contributing and Examples pages to documentation
+- Add Code of Conduct
+- Add gallery of Jupyter notebook examples to documentation using `nbsphinx`
+
+### Removed
+
+### Fixed
+
+- Fix error in units of `NoiseModel.sigma_tau` in `noisefit`
+- Fix bug related to `eta` scaling in `noisefit`
 
 ## [0.3.6] - 2023-11-12
 
@@ -22,8 +76,8 @@ this file.
 
 ### Added
 
-- `GlobalOptions` class
-- `global_options` module constant
+- Add `GlobalOptions` class
+- Add `global_options` module constant
 
 ## [0.3.4] - 2023-11-10
 
@@ -41,7 +95,8 @@ this file.
 ### Changed
 
 - Prune and update notebooks
-- Fix error in `noiseamp` docstring example (now `amplitude` method of `NoiseModel` class)
+- Fix error in `noiseamp` docstring example (now `amplitude` method of
+  `NoiseModel` class)
 - Darken text color for logo in dark mode
 - Move make-logo script into default hatch environment
 - Rename sampling time variable from `ts` to `dt`
@@ -81,14 +136,16 @@ _Changelog introduced._
 
 - Add changelog (Steve)
 - Add frequency bounds to fit function (#27) (Alireza)
-- Add Notes, References, and Examples docstring sections to `noisevar`, `noiseamp`, and `thzgen` (Steve)
+- Add Notes, References, and Examples docstring sections to `noisevar`,
+  `noiseamp`, and `thzgen` (Steve)
 - Add scaling parameters to `tdnll` (now `_tdnll_scaled`) (Steve)
 - Add logo (Steve)
 
 ### Removed
 
 - Remove `costfunlsq`
-- Remove `costfuntls`, `tdnll` from public functions (rename as `_costfuntls` and `_tdnll_scaled`)
+- Remove `costfuntls`, `tdnll` from public functions (rename as `_costfuntls`
+  and `_tdnll_scaled`)
 
 ## [0.3.0] - 2023-09-07
 
@@ -101,6 +158,7 @@ _Changelog introduced._
 
 - Remove `tdtf`
 
+[0.4.0]: https://github.com/dodge-research-group/thztools/releases/tag/v0.4.0
 [0.3.6]: https://github.com/dodge-research-group/thztools/releases/tag/v0.3.6
 [0.3.5]: https://github.com/dodge-research-group/thztools/releases/tag/v0.3.5
 [0.3.4]: https://github.com/dodge-research-group/thztools/releases/tag/v0.3.4
diff --git a/CODE_OF_CONDUCT.md b/CODE_OF_CONDUCT.md
new file mode 100644
index 0000000..a3d4108
--- /dev/null
+++ b/CODE_OF_CONDUCT.md
@@ -0,0 +1,134 @@
+
+# Contributor Covenant Code of Conduct
+
+## Our Pledge
+
+We as members, contributors, and leaders pledge to make participation in our
+community a harassment-free experience for everyone, regardless of age, body
+size, visible or invisible disability, ethnicity, sex characteristics, gender
+identity and expression, level of experience, education, socio-economic status,
+nationality, personal appearance, race, caste, color, religion, or sexual
+identity and orientation.
+
+We pledge to act and interact in ways that contribute to an open, welcoming,
+diverse, inclusive, and healthy community.
+
+## Our Standards
+
+Examples of behavior that contributes to a positive environment for our
+community include:
+
+* Demonstrating empathy and kindness toward other people
+* Being respectful of differing opinions, viewpoints, and experiences
+* Giving and gracefully accepting constructive feedback
+* Accepting responsibility and apologizing to those affected by our mistakes,
+  and learning from the experience
+* Focusing on what is best not just for us as individuals, but for the overall
+  community
+
+Examples of unacceptable behavior include:
+
+* The use of sexualized language or imagery, and sexual attention or advances of
+  any kind
+* Trolling, insulting or derogatory comments, and personal or political attacks
+* Public or private harassment
+* Publishing others' private information, such as a physical or email address,
+  without their explicit permission
+* Other conduct which could reasonably be considered inappropriate in a
+  professional setting
+
+## Enforcement Responsibilities
+
+Community leaders are responsible for clarifying and enforcing our standards of
+acceptable behavior and will take appropriate and fair corrective action in
+response to any behavior that they deem inappropriate, threatening, offensive,
+or harmful.
+
+Community leaders have the right and responsibility to remove, edit, or reject
+comments, commits, code, wiki edits, issues, and other contributions that are
+not aligned to this Code of Conduct, and will communicate reasons for moderation
+decisions when appropriate.
+
+## Scope
+
+This Code of Conduct applies within all community spaces, and also applies when
+an individual is officially representing the community in public spaces.
+Examples of representing our community include using an official email address,
+posting via an official social media account, or acting as an appointed
+representative at an online or offline event.
+
+## Enforcement
+
+Instances of abusive, harassing, or otherwise unacceptable behavior may be
+reported to Steve Dodge at jsdodge@sfu.ca, who is currently the sole community
+leader responsible for enforcement. All complaints will be reviewed and 
+investigated promptly and fairly.
+
+All community leaders are obligated to respect the privacy and security of the
+reporter of any incident.
+
+## Enforcement Guidelines
+
+Community leaders will follow these Community Impact Guidelines in determining
+the consequences for any action they deem in violation of this Code of Conduct:
+
+### 1. Correction
+
+**Community Impact**: Use of inappropriate language or other behavior deemed
+unprofessional or unwelcome in the community.
+
+**Consequence**: A private, written warning from community leaders, providing
+clarity around the nature of the violation and an explanation of why the
+behavior was inappropriate. A public apology may be requested.
+
+### 2. Warning
+
+**Community Impact**: A violation through a single incident or series of
+actions.
+
+**Consequence**: A warning with consequences for continued behavior. No
+interaction with the people involved, including unsolicited interaction with
+those enforcing the Code of Conduct, for a specified period of time. This
+includes avoiding interactions in community spaces as well as external channels
+like social media. Violating these terms may lead to a temporary or permanent
+ban.
+
+### 3. Temporary Ban
+
+**Community Impact**: A serious violation of community standards, including
+sustained inappropriate behavior.
+
+**Consequence**: A temporary ban from any sort of interaction or public
+communication with the community for a specified period of time. No public or
+private interaction with the people involved, including unsolicited interaction
+with those enforcing the Code of Conduct, is allowed during this period.
+Violating these terms may lead to a permanent ban.
+
+### 4. Permanent Ban
+
+**Community Impact**: Demonstrating a pattern of violation of community
+standards, including sustained inappropriate behavior, harassment of an
+individual, or aggression toward or disparagement of classes of individuals.
+
+**Consequence**: A permanent ban from any sort of public interaction within the
+community.
+
+## Attribution
+
+This Code of Conduct is adapted from the [Contributor Covenant][homepage],
+version 2.1, available at
+[https://www.contributor-covenant.org/version/2/1/code_of_conduct.html][v2.1].
+
+Community Impact Guidelines were inspired by
+[Mozilla's code of conduct enforcement ladder][Mozilla CoC].
+
+For answers to common questions about this code of conduct, see the FAQ at
+[https://www.contributor-covenant.org/faq][FAQ]. Translations are available at
+[https://www.contributor-covenant.org/translations][translations].
+
+[homepage]: https://www.contributor-covenant.org
+[v2.1]: https://www.contributor-covenant.org/version/2/1/code_of_conduct.html
+[Mozilla CoC]: https://github.com/mozilla/diversity
+[FAQ]: https://www.contributor-covenant.org/faq
+[translations]: https://www.contributor-covenant.org/translations
+
diff --git a/README.md b/README.md
index fa452b0..7767f92 100644
--- a/README.md
+++ b/README.md
@@ -15,11 +15,11 @@ to use statistically optimal methods for doing this analysis, as described in [L
 
 This is *alpha* software that is currently under development.
 
-| Information       | Links                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         |
-|:------------------|:----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
-| **Project**       | [![DOI](https://zenodo.org/badge/569133241.svg)](https://zenodo.org/doi/10.5281/zenodo.10100093) ![PyPI - Status](https://img.shields.io/pypi/status/thztools) ![PyPI - Python Version](https://img.shields.io/pypi/pyversions/thztools) ![GitHub](https://img.shields.io/github/license/dodge-research-group/thztools) [![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/psf/black) [![Common Changelog](https://common-changelog.org/badge.svg)](https://common-changelog.org)                                                                                                          |
-| **Build**         | ![GitHub Workflow Status](https://img.shields.io/github/actions/workflow/status/dodge-research-group/thztools/sphinx.yml?label=build%3Adocs) ![GitHub Workflow Status](https://img.shields.io/github/actions/workflow/status/dodge-research-group/thztools/pytest-with-coverage.yml?label=build%3Atests%20(conda)) ![GitHub Workflow Status](https://img.shields.io/github/actions/workflow/status/dodge-research-group/thztools/test-pip.yml?label=build%3Atests%20(pip)) [![codecov](https://codecov.io/gh/dodge-research-group/thztools/branch/dev/graph/badge.svg?token=U8PLKTQ7AH)](https://codecov.io/gh/dodge-research-group/thztools) |
-| **Documentation** | https://dodge-research-group.github.io/thztools/                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              |
-| **Cite**          | L. Mohtashemi et al., *Opt. Express* **29**, 4912 (2021). [(DOI)](https://doi.org/10.1364/OE.417724)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          |
+| Information       | Links                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
+|:------------------|:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
+| **Project**       | [![DOI](https://zenodo.org/badge/569133241.svg)](https://zenodo.org/doi/10.5281/zenodo.10100093) ![PyPI - Status](https://img.shields.io/pypi/status/thztools) ![PyPI - Python Version](https://img.shields.io/pypi/pyversions/thztools) ![GitHub](https://img.shields.io/github/license/dodge-research-group/thztools) [![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/psf/black) [![Common Changelog](https://common-changelog.org/badge.svg)](https://common-changelog.org) [![Contributor Covenant](https://img.shields.io/badge/Contributor%20Covenant-2.1-4baaaa.svg)](code_of_conduct.md) |
+| **Build**         | ![GitHub Workflow Status](https://img.shields.io/github/actions/workflow/status/dodge-research-group/thztools/sphinx.yml?label=build%3Adocs) ![GitHub Workflow Status](https://img.shields.io/github/actions/workflow/status/dodge-research-group/thztools/pytest-with-coverage.yml?label=build%3Atests%20(conda)) ![GitHub Workflow Status](https://img.shields.io/github/actions/workflow/status/dodge-research-group/thztools/test-pip.yml?label=build%3Atests%20(pip)) [![codecov](https://codecov.io/gh/dodge-research-group/thztools/branch/dev/graph/badge.svg?token=U8PLKTQ7AH)](https://codecov.io/gh/dodge-research-group/thztools)          |
+| **Documentation** | https://dodge-research-group.github.io/thztools/                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       |
+| **Cite**          | L. Mohtashemi et al., *Opt. Express* **29**, 4912 (2021). [(DOI)](https://doi.org/10.1364/OE.417724)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   |
 
 The original MATLAB code is available at [Zenodo](https://zenodo.org/record/4876388).
diff --git a/docs/source/_templates/autosummary/attribute.rst b/docs/source/_templates/autosummary/attribute.rst
index a176e99..ccba94d 100644
--- a/docs/source/_templates/autosummary/attribute.rst
+++ b/docs/source/_templates/autosummary/attribute.rst
@@ -1,8 +1,9 @@
 :orphan:
 
-{{ fullname }}
-{{ underline }}
+{{ fullname | escape | underline}}
 
 .. currentmodule:: {{ module }}
 
-.. autoattribute:: {{ objname }}
\ No newline at end of file
+attribute
+
+.. auto{{ objtype }}:: {{ fullname }}
diff --git a/docs/source/_templates/autosummary/class.rst b/docs/source/_templates/autosummary/class.rst
index 68c38cc..0ae65c9 100644
--- a/docs/source/_templates/autosummary/class.rst
+++ b/docs/source/_templates/autosummary/class.rst
@@ -1,31 +1,20 @@
-{{ fullname }}
-{{ underline }}
+{% extends "!autosummary/class.rst" %}
 
-.. currentmodule:: {{ module }}
-
-.. autoclass:: {{ objname }}
-   :no-members:
-   :no-inherited-members:
-   :no-special-members:
-
-  {% block methods %}
+{% block methods %}
+{% if methods %}
    .. HACK -- the point here is that we don't want this to appear in the output, but the autosummary should still generate the pages.
       .. autosummary::
          :toctree:
       {% for item in all_methods %}
-         {%- if not item.startswith('_') or item in ['__call__', '__mul__', '__getitem__', '__len__', '__pow__'] %}
-         {{ name }}.{{ item }}
-         {%- endif -%}
-      {%- endfor %}
-      {% for item in inherited_members %}
-         {%- if item in ['__call__', '__mul__', '__getitem__', '__len__', '__pow__'] %}
+         {%- if not item.startswith('_') or item in ['__call__'] %}
          {{ name }}.{{ item }}
          {%- endif -%}
       {%- endfor %}
-  {% endblock %}
+{% endif %}
+{% endblock %}
 
-  {% block attributes %}
-  {% if attributes %}
+{% block attributes %}
+{% if attributes %}
    .. HACK -- the point here is that we don't want this to appear in the output, but the autosummary should still generate the pages.
       .. autosummary::
          :toctree:
@@ -34,5 +23,5 @@
          {{ name }}.{{ item }}
          {%- endif -%}
       {%- endfor %}
-  {% endif %}
-  {% endblock %}
\ No newline at end of file
+{% endif %}
+{% endblock %}
\ No newline at end of file
diff --git a/docs/source/_templates/autosummary/method.rst b/docs/source/_templates/autosummary/method.rst
index 6dce222..4b591ab 100644
--- a/docs/source/_templates/autosummary/method.rst
+++ b/docs/source/_templates/autosummary/method.rst
@@ -1,8 +1,9 @@
 :orphan:
 
-{{ fullname }}
-{{ underline }}
+{{ fullname | escape | underline}}
 
 .. currentmodule:: {{ module }}
 
-.. automethod:: {{ objname }}
\ No newline at end of file
+method
+
+.. auto{{ objtype }}:: {{ fullname }}
diff --git a/docs/source/api.rst b/docs/source/api.rst
index b167942..f099e42 100644
--- a/docs/source/api.rst
+++ b/docs/source/api.rst
@@ -3,36 +3,44 @@ THzTools
 
 .. currentmodule:: thztools
 
-Noise model
-^^^^^^^^^^^
+Global options
+--------------
 
 .. autosummary::
     :toctree: generated/
+    :caption: Global options
 
-    NoiseModel
-    NoiseResult
+    get_option
+    set_option
+
+Noise model
+-----------
 
-Waveform simulation and transformation
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 .. autosummary::
     :toctree: generated/
+    :caption: Noise model
 
-    scaleshift
-    transfer_out
-    wave
+    NoiseModel
+    NoiseResult
+    noisefit
 
 Parameter estimation
-^^^^^^^^^^^^^^^^^^^^
+--------------------
+
 .. autosummary::
     :toctree: generated/
+    :caption: Parameter estimation
 
     fit
-    tdnoisefit
 
-Global options
-^^^^^^^^^^^^^^
+Waveform utilities
+------------------
 
 .. autosummary::
     :toctree: generated/
+    :caption: Waveform utilities
 
-    GlobalOptions
+    scaleshift
+    timebase
+    transfer
+    wave
diff --git a/docs/source/conf.py b/docs/source/conf.py
index ac97cf9..3337e82 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -17,6 +17,7 @@
 # The runpy package is used to import the version number from
 # __about__.py using OPTION 1 below. Not needed for OPTION 2.
 import math
+import os
 import runpy
 import sys
 from pathlib import Path
@@ -46,36 +47,100 @@
 # version = importlib.metadata.version(project)
 # release = version
 
-# -- General configuration ---------------------------------------------------
+# -----------------------------------------------------------------------------
+# General configuration
+# -----------------------------------------------------------------------------
 
 # Add any Sphinx extension module names here, as strings. They can be
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
 # ones.
 extensions = [
-    "matplotlib.sphinxext.plot_directive",
     "sphinx.ext.autodoc",
-    "sphinx.ext.autosummary",
-    "sphinx.ext.coverage",
-    "sphinx.ext.extlinks",
+    "nbsphinx",
+    "numpydoc",
+    "sphinx_design",
+    "sphinx_gallery.load_style",
     "sphinx.ext.intersphinx",
+    "sphinx.ext.coverage",
+    "sphinx.ext.doctest",
+    "sphinx.ext.autosummary",
+    "matplotlib.sphinxext.plot_directive",
     "sphinx.ext.mathjax",
     "sphinx.ext.viewcode",
-    "numpydoc",
+    "sphinx.ext.extlinks",
 ]
 
-# -----------------------------------------------------------------------------
-# Matplotlib plot_directive options (adapted from SciPy docs)
-# -----------------------------------------------------------------------------
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ["_templates"]
 
-plot_pre_code = """
+# The suffix(es) of source filenames.
+# You can specify multiple suffix as a list of string:
+#
+# source_suffix = ['.rst', '.md']
+source_suffix = ".rst"
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+# This pattern also affects html_static_path and html_extra_path.
+exclude_patterns = [
+    "_build",
+    "Thumbs.db",
+    ".DS_Store",
+    "**.ipynb_checkpoints",
+]
+
+# Show successful tests
+doctest_show_successes = True
+
+# Setup code for doctests
+doctest_global_setup = """
 import numpy as np
-np.random.seed(123)
+import thztools as thz
+rng = np.random.default_rng(0)
 """
 
+# -----------------------------------------------------------------------------
+# nbsphinx configuration
+# -----------------------------------------------------------------------------
+
+nbsphinx_allow_errors = True
+
+# Set PYTHONPATH to local src directory (otherwise installs from PyPI)
+src = os.path.abspath('../../src')
+os.environ["PYTHONPATH"] = str(src)
+
+# -----------------------------------------------------------------------------
+# NumPy extensions
+# -----------------------------------------------------------------------------
+
+# Make numpydoc to generate plots for example sections
+numpydoc_use_plots = True
+
+numpydoc_class_members_toctree = True
+numpydoc_attributes_as_param_list = False
+
+# Cross-reference parameter types when True
+numpydoc_xref_param_type = False
+
+# Report warnings for all validation checks
+numpydoc_validation_checks = {"all", "SA01", "ES01", "RT02", "EX01", "PR01"}
+
+# Ensure all our internal links work
+nitpicky = True
+
+# -----------------------------------------------------------------------------
+# Matplotlib plot_directive options (adapted from SciPy docs)
+# -----------------------------------------------------------------------------
+
+# Display source code for Matplotlib plot directives
 plot_include_source = True
-plot_formats = [("png", 96)]
+plot_formats = [("png", 100), "pdf"]
 plot_html_show_formats = False
 plot_html_show_source_link = False
+plot_pre_code = """
+import numpy as np
+"""
+plot_apply_rcparams = True
 
 phi = (math.sqrt(5) + 1) / 2
 
@@ -97,79 +162,38 @@
     "text.usetex": False,
 }
 
-autodoc_type_aliases = {
-    "ArrayLike": "ArrayLike",
-    "Boolean": "bool",
-    "BooleanOrArrayLike": "BooleanOrArrayLike",
-    "BooleanOrNDArray": "BooleanOrNDArray",
-    "DType": "DType",
-    "DTypeBoolean": "DTypeBoolean",
-    "DTypeComplex": "DTypeComplex",
-    "DTypeFloating": "DTypeFloating",
-    "DTypeInteger": "DTypeInteger",
-    "DTypeNumber": "DTypeNumber",
-    "Floating": "float",
-    "FloatingOrArrayLike": "FloatingOrArrayLike",
-    "FloatingOrNDArray": "FloatingOrNDArray",
-    "Integer": "int",
-    "IntegerOrArrayLike": "IntegerOrArrayLike",
-    "IntegerOrNDArray": "IntegerOrNDArray",
-    "NestedSequence": "NestedSequence",
-    "Number": "Number",
-    "NumberOrArrayLike": "NumberOrArrayLike",
-    "NumberOrNDArray": "NumberOrNDArray",
-    "StrOrArrayLike": "StrOrArrayLike",
-    "StrOrNDArray": "StrOrNDArray",
-}
+# -----------------------------------------------------------------------------
+# Autodoc
+# -----------------------------------------------------------------------------
 
-autosummary_generate = True
 autodoc_typehints = "none"
 
-# Report warnings for all validation checks
-numpydoc_validation_checks = {"all", "SA01", "ES01", "RT02", "EX01"}
-
-# Add any paths that contain templates here, relative to this directory.
-templates_path = ["_templates"]
-
-# The suffix(es) of source filenames.
-# You can specify multiple suffix as a list of string:
-#
-# source_suffix = ['.rst', '.md']
-source_suffix = ".rst"
-
-# The master toctree document.
-master_doc = "index"
+# If true, '()' will be appended to :func: etc. cross-reference text.
+add_function_parentheses = False
 
-# Ensure all our internal links work
-nitpicky = True
 
-# List of patterns, relative to source directory, that match files and
-# directories to ignore when looking for source files.
-# This pattern also affects html_static_path and html_extra_path.
-exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"]
+# -----------------------------------------------------------------------------
+# Autosummary
+# -----------------------------------------------------------------------------
 
-# The name of the Pygments (syntax highlighting) style to use.
-pygments_style = "default"
+autosummary_generate = True
 
-mathjax_path = "https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"
+# -----------------------------------------------------------------------------
+# Intersphinx configuration
+# -----------------------------------------------------------------------------
 
 intersphinx_mapping = {
-    "python": (
-        f"https://docs.python.org/{sys.version_info.major}",
-        None,
-    ),
-    "numpy": ("https://numpy.org/doc/stable", None),
-    "scipy": ("https://docs.scipy.org/doc/scipy/", None),
-    "matplotlib": ("https://matplotlib.org/", None),
+    "numpy": ("https://numpy.org/devdocs", None),
+    "python": ("https://docs.python.org/3", None),
+    "scipy": ("https://docs.scipy.org/doc/scipy", None),
+    "matplotlib": ("https://matplotlib.org/stable", None),
 }
 
-
 # -- Options for HTML output -------------------------------------------------
 
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.
 #
-# html_theme = 'sphinx_rtd_theme'
 html_theme = "pydata_sphinx_theme"
 
 # Add any paths that contain custom static files (such as style sheets) here,
@@ -197,6 +221,11 @@
     ],
 }
 
+mathjax_path = "https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"
+
+# The name of the Pygments (syntax highlighting) style to use.
+pygments_style = "default"
+
 # -- Options for HTMLHelp output ---------------------------------------------
 
 # Output file base name for HTML help builder.
@@ -204,6 +233,9 @@
 
 # -- Options for LaTeX output ------------------------------------------------
 
+# The master toctree document.
+master_doc = "index"
+
 latex_elements = {
     # The paper size ('letterpaper' or 'a4paper').
     #
diff --git a/docs/source/contributing.rst b/docs/source/contributing.rst
new file mode 100644
index 0000000..f129599
--- /dev/null
+++ b/docs/source/contributing.rst
@@ -0,0 +1,113 @@
+How to Contribute
+=================
+
+Contributions to THzTools are welcome!
+Here are some ways you can contribute:
+
+* Create a
+  `discussion topic <https://github.com/dodge-research-group/thztools/discussions>`_
+  if you have an idea for a new feature or a general topic for
+  discussion.
+* Create an `issue <https://github.com/dodge-research-group/thztools/issues>`_
+  if you find a bug with THzTools, an issue with its documentation, or have
+  a suggestion for improving the package. New feature requests are also welcome!
+* If you would like to make improvements to the source code or
+  documentation then you may do so directly by opening a
+  `pull request <https://github.com/dodge-research-group/thztools/pulls>`_.
+  Please review the following section for details about the development process.
+
+Please review our `code of conduct <https://github.com/dodge-research-group/thztools/CODE_OF_CONDUCT>`_
+and follow the guidelines described there to help us maintain clear and
+respectful communication.
+
+Development Process
+-------------------
+
+Please review the `Contributing to a project
+<https://docs.github.com/en/get-started/exploring-projects-on-github
+/contributing-to-a-project>`_ instructions on GitHub if you would like to
+contribute to the THzTools code base and documentation.
+
+Development Tools
+-----------------
+
+THzTools development relies on the following tools:
+
+- Unit tests: `pytest <https://pytest.org>`_
+- Code coverage: `coverage <https://coverage.readthedocs.io>`_
+- Documentation builder: `Sphinx <https://www.sphinx-doc.org>`_,
+  using the `PyData Sphinx Theme
+  <https://pydata-sphinx-theme.readthedocs.io>`_ and the
+  `numpydoc <https://numpydoc.readthedocs.io>`_ extension
+- Documentation tests:
+  `doctest <https://docs.python.org/3/library/doctest.html>`_
+- Linter: `ruff <https://docs.astral.sh/ruff/linter/>`_
+- Code format: `black <https://black.readthedocs.io>`_
+- Type checker: `mypy <https://mypy.readthedocs.io>`_
+
+Configuration metadata for each of these tools is stored in ``pyproject.toml``.
+
+Development Environment
+-----------------------
+
+THzTools has been developed with the `Hatch <https://hatch.pypa.io>`_ and
+`Conda <https://docs.conda.io/en/latest/>`_. Brief instructions for their
+installation and use are below.
+
+Hatch
+^^^^^
+
+Follow the instructions at this `link <https://hatch.pypa.io/latest/install/>`__
+to install Hatch.
+
+Enter the environment by typing the following command from the top
+directory of the project (ie, ``~/thztools``)::
+
+  hatch shell
+
+This will update any changes to the dependencies and enable several script
+shortcuts that can be used for testing, linting, and building documentation.
+
+To execute tests with ``pytest``::
+
+  hatch run test
+
+To lint the code base with ``ruff`` and apply ``black`` format rules::
+
+  hatch run lint:all
+
+To build the documentation with ``sphinx`` and test examples with ``doctest``::
+
+  hatch run docs:build
+
+To exit the ``hatch`` environment::
+
+  exit
+
+Conda
+^^^^^
+
+Follow the instructions at this
+`link <https://docs.anaconda.com/free/miniconda/>`__ to install Miniconda, a
+lightweight version of the Conda package and environment manager.
+
+Create the environment by typing the following command from the top directory
+of the project (ie, ``~/thztools``). The first line of the
+``environment-dev.yml`` file sets the environment name to ``thztools``, and the
+last line uses ``pip`` to install the ``thztools`` package in
+`editable <https://pip.pypa.io/en/stable/cli/pip_install/#install-editable/>`_
+mode for development. ::
+
+  conda env create -f envs/environment-dev.yml
+
+To activate the ``thztools`` environment::
+
+  conda activate thztools
+
+To deactivate the environment when you are finished using ``thztools``::
+
+  conda deactivate
+
+More detailed instructions on using Conda to manage environments are available
+`here <https://docs.conda.io/projects/conda/en/stable/user-guide/tasks
+/manage-environments.html>`_.
diff --git a/docs/source/examples.rst b/docs/source/examples.rst
new file mode 100644
index 0000000..1c2c0f7
--- /dev/null
+++ b/docs/source/examples.rst
@@ -0,0 +1,11 @@
+Examples
+========
+
+.. toctree::
+   :maxdepth: 1
+
+.. nbgallery::
+    :name: example-gallery
+    :glob:
+
+    examples/*
\ No newline at end of file
diff --git a/docs/source/examples/simulate-ideal-waveform.ipynb b/docs/source/examples/simulate-ideal-waveform.ipynb
new file mode 100644
index 0000000..63f68f1
--- /dev/null
+++ b/docs/source/examples/simulate-ideal-waveform.ipynb
@@ -0,0 +1,97 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "source": [
+    "# Simulate waveform\n",
+    "\n",
+    "Objects from `thztools`: `thzgen`"
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "code",
+   "source": [
+    "from matplotlib import pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "import thztools as thz"
+   ],
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "markdown",
+   "source": "The first cell below sets the parameters for the simulation, and the subsequent cell uses the function `thzgen` to generate a terahertz time-domain waveform `mu` with `n = 256` samples.",
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "code",
+   "metadata": {},
+   "source": [
+    "n = 256  # Number of samples\n",
+    "dt = 0.05  # Sampling time [ps]"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "code",
+   "metadata": {},
+   "source": [
+    "# Compute the timebase and waveform\n",
+    "t = thz.timebase(n, dt=dt)\n",
+    "mu = thz.wave(n, dt=dt)"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "_, ax = plt.subplots()\n",
+    "\n",
+    "ax.plot(t, mu)\n",
+    "\n",
+    "ax.set_xlabel('Time (ps)')\n",
+    "ax.set_ylabel(r'Amplitude (units of $\\mu_{p})$')\n",
+    "\n",
+    "ax.set_xticks(np.arange(0, 11, 5))\n",
+    "ax.set_xlim(0, 10)\n",
+    "\n",
+    "plt.show()"
+   ],
+   "outputs": [],
+   "execution_count": null
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.16"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 1
+}
diff --git a/docs/source/examples/simulate-noise.ipynb b/docs/source/examples/simulate-noise.ipynb
new file mode 100644
index 0000000..91e064b
--- /dev/null
+++ b/docs/source/examples/simulate-noise.ipynb
@@ -0,0 +1,106 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "source": [
+    "# Simulate noise\n",
+    "\n",
+    "Objects from `thztools`: `thzgen`, `NoiseModel`, `NoiseModel.noise_sim`, `NoiseModel.noise_amp`"
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "code",
+   "source": [
+    "from matplotlib import pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "import thztools as thz"
+   ],
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "markdown",
+   "source": "The first cell below sets the parameters for the simulation, and the subsequent cell uses the function `thzgen` to generate a terahertz time-domain waveform `mu` with `n = 256` samples.",
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "code",
+   "metadata": {},
+   "source": [
+    "n = 256  # Number of samples\n",
+    "dt = 0.05  # Sampling time [ps]\n",
+    "\n",
+    "sigma_alpha = 1e-4  # Additive noise amplitude [signal units]\n",
+    "sigma_beta = 1e-2  # Multiplicative noise amplitude [dimensionless]\n",
+    "sigma_tau = 1e-3  # Time base noise amplitude [ps]\n",
+    "noise_model = thz.NoiseModel(sigma_alpha, sigma_beta, sigma_tau, dt=dt)"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "code",
+   "metadata": {},
+   "source": [
+    "# Compute the timebase and ideal waveform\n",
+    "t = thz.timebase(n, dt=dt)\n",
+    "mu = thz.wave(n, dt=dt)\n",
+    "delta = noise_model.noise_sim(mu, seed=0)\n",
+    "sigma_mu = noise_model.noise_amp(mu)"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "_, ax = plt.subplots()\n",
+    "\n",
+    "ax.plot(t, delta)\n",
+    "ax.plot(t, sigma_mu, 'k:')\n",
+    "ax.plot(t, -sigma_mu, 'k:')\n",
+    "\n",
+    "ax.set_xlabel('Time (ps)')\n",
+    "ax.set_ylabel(r'Amplitude (units of $\\mu_{p})$')\n",
+    "\n",
+    "ax.set_xticks(np.arange(0, 11, 5))\n",
+    "ax.set_xlim(0, 10)\n",
+    "\n",
+    "plt.show()"
+   ],
+   "outputs": [],
+   "execution_count": null
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.16"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 1
+}
diff --git a/docs/source/generated/thztools.GlobalOptions.rst b/docs/source/generated/thztools.GlobalOptions.rst
deleted file mode 100644
index cd39dff..0000000
--- a/docs/source/generated/thztools.GlobalOptions.rst
+++ /dev/null
@@ -1,27 +0,0 @@
-﻿thztools.GlobalOptions
-======================
-
-.. currentmodule:: thztools
-
-.. autoclass:: GlobalOptions
-   :no-members:
-   :no-inherited-members:
-   :no-special-members:
-
-  
-   .. HACK -- the point here is that we don't want this to appear in the output, but the autosummary should still generate the pages.
-      .. autosummary::
-         :toctree:
-      
-      
-  
-
-  
-  
-   .. HACK -- the point here is that we don't want this to appear in the output, but the autosummary should still generate the pages.
-      .. autosummary::
-         :toctree:
-      
-         GlobalOptions.sampling_time
-  
-  
\ No newline at end of file
diff --git a/docs/source/generated/thztools.GlobalOptions.sampling_time.rst b/docs/source/generated/thztools.GlobalOptions.sampling_time.rst
deleted file mode 100644
index e140abb..0000000
--- a/docs/source/generated/thztools.GlobalOptions.sampling_time.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-:orphan:
-
-thztools.GlobalOptions.sampling_time
-====================================
-
-.. currentmodule:: thztools
-
-.. autoattribute:: GlobalOptions.sampling_time
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseModel.alpha.rst b/docs/source/generated/thztools.NoiseModel.alpha.rst
deleted file mode 100644
index 47d40ae..0000000
--- a/docs/source/generated/thztools.NoiseModel.alpha.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-:orphan:
-
-thztools.NoiseModel.alpha
-=========================
-
-.. currentmodule:: thztools
-
-.. autoattribute:: NoiseModel.alpha
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseModel.amplitude.rst b/docs/source/generated/thztools.NoiseModel.amplitude.rst
deleted file mode 100644
index aaab7a5..0000000
--- a/docs/source/generated/thztools.NoiseModel.amplitude.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-:orphan:
-
-thztools.NoiseModel.amplitude
-=============================
-
-.. currentmodule:: thztools
-
-.. automethod:: NoiseModel.amplitude
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseModel.beta.rst b/docs/source/generated/thztools.NoiseModel.beta.rst
deleted file mode 100644
index 74b1014..0000000
--- a/docs/source/generated/thztools.NoiseModel.beta.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-:orphan:
-
-thztools.NoiseModel.beta
-========================
-
-.. currentmodule:: thztools
-
-.. autoattribute:: NoiseModel.beta
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseModel.dt.rst b/docs/source/generated/thztools.NoiseModel.dt.rst
index e68cc19..6634d7f 100644
--- a/docs/source/generated/thztools.NoiseModel.dt.rst
+++ b/docs/source/generated/thztools.NoiseModel.dt.rst
@@ -5,4 +5,6 @@ thztools.NoiseModel.dt
 
 .. currentmodule:: thztools
 
-.. autoattribute:: NoiseModel.dt
\ No newline at end of file
+attribute
+
+.. autoattribute:: thztools.NoiseModel.dt
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseModel.noise.rst b/docs/source/generated/thztools.NoiseModel.noise.rst
deleted file mode 100644
index 854bafa..0000000
--- a/docs/source/generated/thztools.NoiseModel.noise.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-:orphan:
-
-thztools.NoiseModel.noise
-=========================
-
-.. currentmodule:: thztools
-
-.. automethod:: NoiseModel.noise
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseModel.noise_amp.rst b/docs/source/generated/thztools.NoiseModel.noise_amp.rst
new file mode 100644
index 0000000..88a8764
--- /dev/null
+++ b/docs/source/generated/thztools.NoiseModel.noise_amp.rst
@@ -0,0 +1,10 @@
+:orphan:
+
+thztools.NoiseModel.noise\_amp
+==============================
+
+.. currentmodule:: thztools
+
+method
+
+.. automethod:: thztools.NoiseModel.noise_amp
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseModel.noise_sim.rst b/docs/source/generated/thztools.NoiseModel.noise_sim.rst
new file mode 100644
index 0000000..232baae
--- /dev/null
+++ b/docs/source/generated/thztools.NoiseModel.noise_sim.rst
@@ -0,0 +1,10 @@
+:orphan:
+
+thztools.NoiseModel.noise\_sim
+==============================
+
+.. currentmodule:: thztools
+
+method
+
+.. automethod:: thztools.NoiseModel.noise_sim
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseModel.noise_var.rst b/docs/source/generated/thztools.NoiseModel.noise_var.rst
new file mode 100644
index 0000000..04612ea
--- /dev/null
+++ b/docs/source/generated/thztools.NoiseModel.noise_var.rst
@@ -0,0 +1,10 @@
+:orphan:
+
+thztools.NoiseModel.noise\_var
+==============================
+
+.. currentmodule:: thztools
+
+method
+
+.. automethod:: thztools.NoiseModel.noise_var
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseModel.rst b/docs/source/generated/thztools.NoiseModel.rst
index 8a4aaa1..49aaf02 100644
--- a/docs/source/generated/thztools.NoiseModel.rst
+++ b/docs/source/generated/thztools.NoiseModel.rst
@@ -4,30 +4,27 @@
 .. currentmodule:: thztools
 
 .. autoclass:: NoiseModel
-   :no-members:
-   :no-inherited-members:
-   :no-special-members:
 
-  
+   
+
    .. HACK -- the point here is that we don't want this to appear in the output, but the autosummary should still generate the pages.
       .. autosummary::
          :toctree:
       
-         NoiseModel.amplitude
-         NoiseModel.noise
-         NoiseModel.variance
-      
-  
+         NoiseModel.noise_amp
+         NoiseModel.noise_sim
+         NoiseModel.noise_var
+
+
+
+   
 
-  
-  
    .. HACK -- the point here is that we don't want this to appear in the output, but the autosummary should still generate the pages.
       .. autosummary::
          :toctree:
       
          NoiseModel.dt
-         NoiseModel.alpha
-         NoiseModel.beta
-         NoiseModel.tau
-  
-  
\ No newline at end of file
+         NoiseModel.sigma_alpha
+         NoiseModel.sigma_beta
+         NoiseModel.sigma_tau
+
diff --git a/docs/source/generated/thztools.NoiseModel.sigma_alpha.rst b/docs/source/generated/thztools.NoiseModel.sigma_alpha.rst
new file mode 100644
index 0000000..9bd4601
--- /dev/null
+++ b/docs/source/generated/thztools.NoiseModel.sigma_alpha.rst
@@ -0,0 +1,10 @@
+:orphan:
+
+thztools.NoiseModel.sigma\_alpha
+================================
+
+.. currentmodule:: thztools
+
+attribute
+
+.. autoattribute:: thztools.NoiseModel.sigma_alpha
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseModel.sigma_beta.rst b/docs/source/generated/thztools.NoiseModel.sigma_beta.rst
new file mode 100644
index 0000000..9e2d092
--- /dev/null
+++ b/docs/source/generated/thztools.NoiseModel.sigma_beta.rst
@@ -0,0 +1,10 @@
+:orphan:
+
+thztools.NoiseModel.sigma\_beta
+===============================
+
+.. currentmodule:: thztools
+
+attribute
+
+.. autoattribute:: thztools.NoiseModel.sigma_beta
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseModel.sigma_tau.rst b/docs/source/generated/thztools.NoiseModel.sigma_tau.rst
new file mode 100644
index 0000000..7436dab
--- /dev/null
+++ b/docs/source/generated/thztools.NoiseModel.sigma_tau.rst
@@ -0,0 +1,10 @@
+:orphan:
+
+thztools.NoiseModel.sigma\_tau
+==============================
+
+.. currentmodule:: thztools
+
+attribute
+
+.. autoattribute:: thztools.NoiseModel.sigma_tau
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseModel.tau.rst b/docs/source/generated/thztools.NoiseModel.tau.rst
deleted file mode 100644
index 29c2b9c..0000000
--- a/docs/source/generated/thztools.NoiseModel.tau.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-:orphan:
-
-thztools.NoiseModel.tau
-=======================
-
-.. currentmodule:: thztools
-
-.. autoattribute:: NoiseModel.tau
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseModel.variance.rst b/docs/source/generated/thztools.NoiseModel.variance.rst
deleted file mode 100644
index 10a268e..0000000
--- a/docs/source/generated/thztools.NoiseModel.variance.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-:orphan:
-
-thztools.NoiseModel.variance
-============================
-
-.. currentmodule:: thztools
-
-.. automethod:: NoiseModel.variance
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseResult.a.rst b/docs/source/generated/thztools.NoiseResult.a.rst
index b6f77e0..c2c7fc0 100644
--- a/docs/source/generated/thztools.NoiseResult.a.rst
+++ b/docs/source/generated/thztools.NoiseResult.a.rst
@@ -5,4 +5,6 @@ thztools.NoiseResult.a
 
 .. currentmodule:: thztools
 
-.. autoattribute:: NoiseResult.a
\ No newline at end of file
+attribute
+
+.. autoattribute:: thztools.NoiseResult.a
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseResult.diagnostic.rst b/docs/source/generated/thztools.NoiseResult.diagnostic.rst
index caa4343..153faf7 100644
--- a/docs/source/generated/thztools.NoiseResult.diagnostic.rst
+++ b/docs/source/generated/thztools.NoiseResult.diagnostic.rst
@@ -5,4 +5,6 @@ thztools.NoiseResult.diagnostic
 
 .. currentmodule:: thztools
 
-.. autoattribute:: NoiseResult.diagnostic
\ No newline at end of file
+attribute
+
+.. autoattribute:: thztools.NoiseResult.diagnostic
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseResult.err_a.rst b/docs/source/generated/thztools.NoiseResult.err_a.rst
new file mode 100644
index 0000000..d1233dc
--- /dev/null
+++ b/docs/source/generated/thztools.NoiseResult.err_a.rst
@@ -0,0 +1,10 @@
+:orphan:
+
+thztools.NoiseResult.err\_a
+===========================
+
+.. currentmodule:: thztools
+
+attribute
+
+.. autoattribute:: thztools.NoiseResult.err_a
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseResult.err_eta.rst b/docs/source/generated/thztools.NoiseResult.err_eta.rst
new file mode 100644
index 0000000..fa32774
--- /dev/null
+++ b/docs/source/generated/thztools.NoiseResult.err_eta.rst
@@ -0,0 +1,10 @@
+:orphan:
+
+thztools.NoiseResult.err\_eta
+=============================
+
+.. currentmodule:: thztools
+
+attribute
+
+.. autoattribute:: thztools.NoiseResult.err_eta
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseResult.err_mu.rst b/docs/source/generated/thztools.NoiseResult.err_mu.rst
new file mode 100644
index 0000000..dad9def
--- /dev/null
+++ b/docs/source/generated/thztools.NoiseResult.err_mu.rst
@@ -0,0 +1,10 @@
+:orphan:
+
+thztools.NoiseResult.err\_mu
+============================
+
+.. currentmodule:: thztools
+
+attribute
+
+.. autoattribute:: thztools.NoiseResult.err_mu
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseResult.err_sigma_alpha.rst b/docs/source/generated/thztools.NoiseResult.err_sigma_alpha.rst
new file mode 100644
index 0000000..5bdb2d9
--- /dev/null
+++ b/docs/source/generated/thztools.NoiseResult.err_sigma_alpha.rst
@@ -0,0 +1,10 @@
+:orphan:
+
+thztools.NoiseResult.err\_sigma\_alpha
+======================================
+
+.. currentmodule:: thztools
+
+attribute
+
+.. autoattribute:: thztools.NoiseResult.err_sigma_alpha
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseResult.err_sigma_beta.rst b/docs/source/generated/thztools.NoiseResult.err_sigma_beta.rst
new file mode 100644
index 0000000..2d5bde8
--- /dev/null
+++ b/docs/source/generated/thztools.NoiseResult.err_sigma_beta.rst
@@ -0,0 +1,10 @@
+:orphan:
+
+thztools.NoiseResult.err\_sigma\_beta
+=====================================
+
+.. currentmodule:: thztools
+
+attribute
+
+.. autoattribute:: thztools.NoiseResult.err_sigma_beta
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseResult.err_sigma_tau.rst b/docs/source/generated/thztools.NoiseResult.err_sigma_tau.rst
new file mode 100644
index 0000000..9c0980c
--- /dev/null
+++ b/docs/source/generated/thztools.NoiseResult.err_sigma_tau.rst
@@ -0,0 +1,10 @@
+:orphan:
+
+thztools.NoiseResult.err\_sigma\_tau
+====================================
+
+.. currentmodule:: thztools
+
+attribute
+
+.. autoattribute:: thztools.NoiseResult.err_sigma_tau
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseResult.eta.rst b/docs/source/generated/thztools.NoiseResult.eta.rst
index 5674ccf..40c6d2e 100644
--- a/docs/source/generated/thztools.NoiseResult.eta.rst
+++ b/docs/source/generated/thztools.NoiseResult.eta.rst
@@ -5,4 +5,6 @@ thztools.NoiseResult.eta
 
 .. currentmodule:: thztools
 
-.. autoattribute:: NoiseResult.eta
\ No newline at end of file
+attribute
+
+.. autoattribute:: thztools.NoiseResult.eta
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseResult.fval.rst b/docs/source/generated/thztools.NoiseResult.fval.rst
index 75d15d4..4d4bcb8 100644
--- a/docs/source/generated/thztools.NoiseResult.fval.rst
+++ b/docs/source/generated/thztools.NoiseResult.fval.rst
@@ -5,4 +5,6 @@ thztools.NoiseResult.fval
 
 .. currentmodule:: thztools
 
-.. autoattribute:: NoiseResult.fval
\ No newline at end of file
+attribute
+
+.. autoattribute:: thztools.NoiseResult.fval
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseResult.hess_inv.rst b/docs/source/generated/thztools.NoiseResult.hess_inv.rst
new file mode 100644
index 0000000..b58a607
--- /dev/null
+++ b/docs/source/generated/thztools.NoiseResult.hess_inv.rst
@@ -0,0 +1,10 @@
+:orphan:
+
+thztools.NoiseResult.hess\_inv
+==============================
+
+.. currentmodule:: thztools
+
+attribute
+
+.. autoattribute:: thztools.NoiseResult.hess_inv
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseResult.mu.rst b/docs/source/generated/thztools.NoiseResult.mu.rst
index 762859b..1ae7f94 100644
--- a/docs/source/generated/thztools.NoiseResult.mu.rst
+++ b/docs/source/generated/thztools.NoiseResult.mu.rst
@@ -5,4 +5,6 @@ thztools.NoiseResult.mu
 
 .. currentmodule:: thztools
 
-.. autoattribute:: NoiseResult.mu
\ No newline at end of file
+attribute
+
+.. autoattribute:: thztools.NoiseResult.mu
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseResult.noise_model.rst b/docs/source/generated/thztools.NoiseResult.noise_model.rst
index 1bf3c8a..30e1dc4 100644
--- a/docs/source/generated/thztools.NoiseResult.noise_model.rst
+++ b/docs/source/generated/thztools.NoiseResult.noise_model.rst
@@ -1,8 +1,10 @@
 :orphan:
 
-thztools.NoiseResult.noise_model
-================================
+thztools.NoiseResult.noise\_model
+=================================
 
 .. currentmodule:: thztools
 
-.. autoattribute:: NoiseResult.noise_model
\ No newline at end of file
+attribute
+
+.. autoattribute:: thztools.NoiseResult.noise_model
\ No newline at end of file
diff --git a/docs/source/generated/thztools.NoiseResult.rst b/docs/source/generated/thztools.NoiseResult.rst
index 77bc99d..a592499 100644
--- a/docs/source/generated/thztools.NoiseResult.rst
+++ b/docs/source/generated/thztools.NoiseResult.rst
@@ -4,20 +4,18 @@
 .. currentmodule:: thztools
 
 .. autoclass:: NoiseResult
-   :no-members:
-   :no-inherited-members:
-   :no-special-members:
 
-  
+   
+
    .. HACK -- the point here is that we don't want this to appear in the output, but the autosummary should still generate the pages.
       .. autosummary::
          :toctree:
       
-      
-  
 
-  
-  
+
+
+   
+
    .. HACK -- the point here is that we don't want this to appear in the output, but the autosummary should still generate the pages.
       .. autosummary::
          :toctree:
@@ -27,6 +25,12 @@
          NoiseResult.a
          NoiseResult.eta
          NoiseResult.fval
+         NoiseResult.hess_inv
+         NoiseResult.err_sigma_alpha
+         NoiseResult.err_sigma_beta
+         NoiseResult.err_sigma_tau
+         NoiseResult.err_mu
+         NoiseResult.err_a
+         NoiseResult.err_eta
          NoiseResult.diagnostic
-  
-  
\ No newline at end of file
+
diff --git a/docs/source/generated/thztools.get_option.rst b/docs/source/generated/thztools.get_option.rst
new file mode 100644
index 0000000..0998225
--- /dev/null
+++ b/docs/source/generated/thztools.get_option.rst
@@ -0,0 +1,6 @@
+﻿thztools.get\_option
+====================
+
+.. currentmodule:: thztools
+
+.. autofunction:: get_option
\ No newline at end of file
diff --git a/docs/source/generated/thztools.noisefit.rst b/docs/source/generated/thztools.noisefit.rst
new file mode 100644
index 0000000..0b4fe82
--- /dev/null
+++ b/docs/source/generated/thztools.noisefit.rst
@@ -0,0 +1,6 @@
+﻿thztools.noisefit
+=================
+
+.. currentmodule:: thztools
+
+.. autofunction:: noisefit
\ No newline at end of file
diff --git a/docs/source/generated/thztools.set_option.rst b/docs/source/generated/thztools.set_option.rst
new file mode 100644
index 0000000..9120ce5
--- /dev/null
+++ b/docs/source/generated/thztools.set_option.rst
@@ -0,0 +1,6 @@
+﻿thztools.set\_option
+====================
+
+.. currentmodule:: thztools
+
+.. autofunction:: set_option
\ No newline at end of file
diff --git a/docs/source/generated/thztools.tdnoisefit.rst b/docs/source/generated/thztools.tdnoisefit.rst
deleted file mode 100644
index d6400c7..0000000
--- a/docs/source/generated/thztools.tdnoisefit.rst
+++ /dev/null
@@ -1,6 +0,0 @@
-﻿thztools.tdnoisefit
-===================
-
-.. currentmodule:: thztools
-
-.. autofunction:: tdnoisefit
\ No newline at end of file
diff --git a/docs/source/generated/thztools.timebase.rst b/docs/source/generated/thztools.timebase.rst
new file mode 100644
index 0000000..77c8a63
--- /dev/null
+++ b/docs/source/generated/thztools.timebase.rst
@@ -0,0 +1,6 @@
+﻿thztools.timebase
+=================
+
+.. currentmodule:: thztools
+
+.. autofunction:: timebase
\ No newline at end of file
diff --git a/docs/source/generated/thztools.transfer.rst b/docs/source/generated/thztools.transfer.rst
new file mode 100644
index 0000000..98a0b6f
--- /dev/null
+++ b/docs/source/generated/thztools.transfer.rst
@@ -0,0 +1,6 @@
+﻿thztools.transfer
+=================
+
+.. currentmodule:: thztools
+
+.. autofunction:: transfer
\ No newline at end of file
diff --git a/docs/source/generated/thztools.transfer_out.rst b/docs/source/generated/thztools.transfer_out.rst
deleted file mode 100644
index 8ea4d67..0000000
--- a/docs/source/generated/thztools.transfer_out.rst
+++ /dev/null
@@ -1,6 +0,0 @@
-﻿thztools.transfer\_out
-======================
-
-.. currentmodule:: thztools
-
-.. autofunction:: transfer_out
\ No newline at end of file
diff --git a/docs/source/getting_started.rst b/docs/source/getting_started.rst
new file mode 100644
index 0000000..1726287
--- /dev/null
+++ b/docs/source/getting_started.rst
@@ -0,0 +1,28 @@
+Getting Started
+===================
+
+.. toctree::
+   :maxdepth: 2
+
+Installation
+------------
+
+Using pip: ::
+
+    pip install thztools
+
+The package name is ``thztools``, and the import convention is: ::
+
+    import thztools as thz
+
+Recommended imports
+-------------------
+
+The following packages are not required to use ``thztools`` but make it more
+convenient to use.
+
+- `matplotlib <https://matplotlib.org/>`_
+
+- `numpy <https://numpy.org/>`_
+
+- `scipy <https://docs.scipy.org/>`_
diff --git a/docs/source/gettingstarted.rst b/docs/source/gettingstarted.rst
deleted file mode 100644
index 820e51c..0000000
--- a/docs/source/gettingstarted.rst
+++ /dev/null
@@ -1,29 +0,0 @@
-Getting Started
-===================
-
-.. toctree::
-   :maxdepth: 2
-
-Installation
-+++++++++++++++++++
-
-The instructions below are for developers using the `Conda <https://docs.conda.io/en/latest/>`_ package and environment manager.
-
-1. Create the environment by typing the following command from the top directory. The first line of the **environment-dev.yml** file sets the environment name to **thztools**, and the last line uses **pip** to install the **thztools** package in `editable <https://pip.pypa.io/en/stable/cli/pip_install/#install-editable/>`_ mode for development. ::
-
-    conda env create -f envs/environment-dev.yml
-
-2. Activate the **thztools** environment. ::
-
-    conda activate thztools
-
-3. Use the following command to run a Jupyter notebook server from within this environment. ::
-
-    jupyter notebook
-
-4. Deactivate the environment when you are finished using **thztools**. ::
-
-    conda deactivate
-
-More detailed instructions on using Conda to manage environments are available `here <https://docs.conda.io/projects/conda/en/latest/user-guide/tasks/manage-environments.html/>`_.
-
diff --git a/docs/source/index.rst b/docs/source/index.rst
index a4e4489..d299e30 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -1,3 +1,5 @@
+:html_theme.sidebar_secondary.remove:
+
 THzTools Documentation
 ====================================
 
@@ -10,10 +12,47 @@ represent as a function of *frequency* for analysis. THzTools makes it easier fo
 to use statistically optimal methods for doing this analysis, as described in Laleh Mohtashemi
 et al., Opt. Express **29**, 4912 (2021), `<https://doi.org/10.1364/OE.417724>`_.
 
-Table of Contents
--------------------
 .. toctree::
-   :maxdepth: 3
+   :hidden:
 
-   Getting started <gettingstarted>
+   Home <self>
+   Getting started <getting_started>
+   Examples <examples>
    API Reference <api>
+   Contributing <contributing>
+
+.. grid::
+
+   .. grid-item-card::
+       :link: getting_started
+       :link-type: doc
+       :text-align: center
+       :class-body: sd-fs-3 sd-text-primary
+
+       Getting started
+
+   .. grid-item-card::
+       :link: examples
+       :link-type: doc
+       :text-align: center
+       :class-body: sd-fs-3 sd-text-primary
+
+       Examples
+
+.. grid::
+
+   .. grid-item-card::
+       :link: api
+       :link-type: doc
+       :text-align: center
+       :class-body: sd-fs-3 sd-text-primary
+
+       API Reference
+
+   .. grid-item-card::
+       :link: contributing
+       :link-type: doc
+       :text-align: center
+       :class-body: sd-fs-3 sd-text-primary
+
+       Contributing
diff --git a/envs/environment-dev.yml b/envs/environment-dev.yml
index 4868f2f..8bba99e 100644
--- a/envs/environment-dev.yml
+++ b/envs/environment-dev.yml
@@ -4,19 +4,19 @@ channels:
   - nodefaults
 dependencies:
   - black
+  - coverage
+  - jupyter
+  - matplotlib
   - mypy
-  - scipy
   - numdifftools
+  - numpydoc
   - pandas
-  - jupyter
-  - matplotlib
-  - sphinx
+  - pip
   - pydata-sphinx-theme
-  - numpydoc
-  - tomli
-  - ruff
   - pytest-cov
-  - coverage
-  - pip
+  - ruff
+  - scipy
+  - sphinx
+  - tomli
   - pip:
     - -e ../
diff --git a/examples/figgen/calcnoise.py b/examples/figgen/calcnoise.py
index 281ef6c..ceb42b4 100644
--- a/examples/figgen/calcnoise.py
+++ b/examples/figgen/calcnoise.py
@@ -2,7 +2,7 @@
 
 import numpy as np
 
-from src.thztools import airscancorrect, shiftmtx, tdnoisefit, tdtf
+from src.thztools import airscancorrect, noisefit, shiftmtx, tdtf
 
 
 def calcnoise(t, x):
@@ -35,9 +35,7 @@ def calcnoise(t, x):
 
     mu0 = np.mean(x, axis=1)
     param["delay"] = {"v0": v0, "mu0": mu0, "ts": ts}
-    eta0 = tdnoisefit(x, param["delay"], fix["delay"], ignore["delay"])[0][
-        "eta"
-    ]
+    eta0 = noisefit(x, param["delay"], fix["delay"], ignore["delay"])[0]["eta"]
 
     print("Elapsed time is", time.time() - start_time)
 
@@ -49,7 +47,7 @@ def calcnoise(t, x):
     fix["amp"] = {"logv": True, "mu": True, "a": False, "eta": True}
     ignore["amp"] = {"a": False, "eta": True}
     param["amp"] = {"v0": v0, "mu0": mu0, "eta0": eta0, "ts": ts}
-    a0 = tdnoisefit(x, param["amp"], fix["amp"], ignore["amp"])[0]["a"]
+    a0 = noisefit(x, param["amp"], fix["amp"], ignore["amp"])[0]["a"]
 
     print("Elapsed time is", time.time() - start_time)
 
@@ -69,7 +67,7 @@ def calcnoise(t, x):
     fix["var"] = {"logv": False, "mu": True, "a": True, "eta": True}
     ignore["var"] = {"a": False, "eta": False}
     param["var"] = {"v0": v0, "mu0": mu0, "a0": a0, "eta": eta0, "ts": ts}
-    v0 = tdnoisefit(x, param["var"], fix["var"], ignore["var"])[0]["var"]
+    v0 = noisefit(x, param["var"], fix["var"], ignore["var"])[0]["var"]
 
     print("Elapsed time is", time.time() - start_time)
 
@@ -80,7 +78,7 @@ def calcnoise(t, x):
     fix["all"] = {"logv": False, "mu": False, "a": False, "eta": False}
     ignore["all"] = {"a": False, "eta": False}
     param["all"] = {"v0": v0, "mu0": mu0, "a": a0, "eta": eta0, "ts": ts}
-    all_fit = tdnoisefit(x, param["all"], fix["all"], ignore["all"])
+    all_fit = noisefit(x, param["all"], fix["all"], ignore["all"])
 
     # Return results through Output Structure
     output["x"] = x
diff --git a/examples/figgen/datamc.py b/examples/figgen/datamc.py
index 4e91ad5..d3c5a2b 100644
--- a/examples/figgen/datamc.py
+++ b/examples/figgen/datamc.py
@@ -65,7 +65,7 @@ def datamc(**kwargs):
         "Vi": vi,
         "V": v,
         "Nf": nf,
-        "sigma_t": sigma_t,
+        "noise_amp": sigma_t,
     }
 
     return data
diff --git a/examples/figgen/plotnoiseest.py b/examples/figgen/plotnoiseest.py
index cfd1004..4826083 100644
--- a/examples/figgen/plotnoiseest.py
+++ b/examples/figgen/plotnoiseest.py
@@ -12,7 +12,7 @@
 # using DataPulse class get the amplitude
 air_scans = []
 for file in os.listdir(path):
-    obj = DataPulse(path + "/" + file).amplitude
+    obj = DataPulse(path + "/" + file).noise_amp
     air_scans.append(obj)
 
 # get the time array
diff --git a/examples/figgen/plotsignal.py b/examples/figgen/plotsignal.py
index 035cb41..d3d8f29 100644
--- a/examples/figgen/plotsignal.py
+++ b/examples/figgen/plotsignal.py
@@ -9,7 +9,7 @@ def plotsignal():
     plt.plot(data["t"], data["y0"], color="black", label=r"Signal ($\mu$)")
     plt.plot(
         data["t"],
-        data["sigma_t"] * 30,
+        data["noise_amp"] * 30,
         color="red",
         linestyle="dashed",
         label=r"noise $(30 \sigma$)",
diff --git a/examples/simulation-example-2.ipynb b/examples/simulation-example-2.ipynb
index 055c11d..82e3238 100644
--- a/examples/simulation-example-2.ipynb
+++ b/examples/simulation-example-2.ipynb
@@ -11,8 +11,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 1,
    "source": [
     "from matplotlib import pyplot as plt\n",
     "from matplotlib.figure import figaspect\n",
@@ -31,8 +30,13 @@
     "from thztools.thztools import _costfuntls"
    ],
    "metadata": {
-    "collapsed": false
-   }
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T08:29:35.818460Z",
+     "start_time": "2024-01-09T08:29:35.038130Z"
+    }
+   },
+   "outputs": []
   },
   {
    "cell_type": "markdown",
@@ -50,711 +54,162 @@
    }
   },
   {
+   "metadata": {},
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
-   "outputs": [],
    "source": [
     "n = 256  # Number of samples\n",
-    "ts = 0.05  # Sampling time [ps]\n",
-    "t0 = 2.5  # Peak pulse time [ps]\n",
-    "fs = 1 / ts  # Sampling frequency [THz]\n",
+    "m = 10  # Number of waveforms\n",
+    "dt = 0.05  # Sampling time [ps]\n",
     "\n",
-    "sigma_alpha = 1e-4  # Additive noise amplitude [signal units]\n",
+    "sigma_alpha = 1e-5  # Additive noise amplitude [signal units]\n",
     "sigma_beta = 1e-2  # Multiplicative noise amplitude [dimensionless]\n",
     "sigma_tau = 1e-3  # Time base noise amplitude [ps]\n",
-    "sigma_parms = np.array([sigma_alpha, sigma_beta, sigma_tau])\n",
+    "noise_model = thz.NoiseModel(sigma_alpha, sigma_beta, sigma_tau)\n",
     "\n",
     "seed = 0  # RNG seed\n",
     "rng = default_rng(seed)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "mu, t = thz.thzgen(n, ts, t0)\n",
-    "sigma = thz.noiseamp(sigma_parms, mu, ts)\n",
-    "\n",
-    "# Use broadcasting to generate multiple simulated noisy measurements of y\n",
-    "n_sim = 11\n",
-    "x = mu + sigma * rng.standard_normal((n_sim, n))\n",
-    "\n",
-    "# Compute power spectral density\n",
-    "f, psd = signal.periodogram(x, fs=fs, detrend=False)\n",
-    "n_f = f.size"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "w, h = figaspect(0.45)\n",
-    "_, axs = plt.subplots(1, 2, figsize=(w, h))\n",
-    "\n",
-    "axs[0].plot(t, mu, c='k', label=r'Signal ($\\mu$)')\n",
-    "axs[0].plot(t, sigma * 30, ls='--', c='grey', label=r'noise $(30 \\sigma$)')\n",
-    "\n",
-    "axs[0].legend(frameon=False)\n",
-    "axs[0].text(0.025, 0.95, '(a)', transform=axs[0].transAxes)\n",
-    "\n",
-    "axs[0].set_xlabel('Time (ps)')\n",
-    "axs[0].set_ylabel(r'Amplitude (units of $\\mu_{p})$')\n",
-    "\n",
-    "axs[1].plot(f, 10 * np.log10(psd[1:, :].T / max(psd.flatten())), c='grey')\n",
-    "axs[1].plot(f, 10 * np.log10(psd[0, :].T / max(psd.flatten())), c='red')\n",
-    "\n",
-    "axs[1].text(0.025, 0.95, '(b)', transform=axs[1].transAxes)\n",
-    "\n",
-    "axs[1].set_xlabel('Frequency (THz)')\n",
-    "axs[1].set_ylabel('Relative Power (dB)')\n",
-    "\n",
-    "axs[1].set_xticks(np.arange(0, 11, 5))\n",
-    "axs[1].set_xlim(0, 10)\n",
-    "axs[1].set_yticks(np.arange(-60, 10, 20))\n",
-    "axs[1].set_ylim(-70, 10)\n",
-    "\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "source": [
-    "## Empirical transfer function estimates\n",
-    "\n",
-    "The first cell below simulates `n_sim = 250` pairs of noisy waveforms, `x` and `y`, generated from the same ideal waveforms `mu` and `psi`, respectively, where `psi` is obtained from `mu` by multiplying its amplitude by `0.5` and delaying it by `1.0`. The cell then computes the empirical transfer function estimate `etfe`.\n",
-    "\n",
-    "The next cell generates a figure that shows the real (a,c) and the imaginary (b,d) parts of `etfe` as gray dots. One estimate is highlighted in red, with the dots connected by a thin line. Thick black lines indicate the average over all `n_sim = 250` simulations. Panels (a,b) show the full bandwidth and (c,d) show the same data over the primary signal bandwidth."
    ],
-   "metadata": {
-    "collapsed": false
-   }
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
-   "source": [
-    "# Clear x and n_sim for reuse\n",
-    "del x, n_sim\n",
-    "rng = default_rng(seed)\n",
-    "\n",
-    "n_sim = 250\n",
-    "p_opt = np.array([0.5, 1.0])\n",
-    "n_p = p_opt.size\n",
-    "\n",
-    "f = rfftfreq(n, ts)\n",
-    "\n",
-    "\n",
-    "def tf_model(_p, _w):\n",
-    "    return _p[0] * np.exp(-1j * _w * _p[1])\n",
-    "\n",
-    "\n",
-    "def td_model(_p, _x):\n",
-    "    _y = irfft(rfft(_x) * tf_model(_p, 2 * pi * f), n=n)\n",
-    "    return _y\n",
-    "\n",
-    "\n",
-    "psi = td_model(p_opt, mu)\n",
-    "\n",
-    "sigma_mu = thz.noiseamp(sigma_parms, mu, ts)\n",
-    "sigma_psi = thz.noiseamp(sigma_parms, psi, ts)\n",
-    "\n",
-    "\n",
-    "x_noise = sigma_mu * rng.standard_normal((n_sim, n))\n",
-    "y_noise = sigma_psi * rng.standard_normal((n_sim, n))\n",
-    "\n",
-    "x = mu - x_noise\n",
-    "y = psi - y_noise\n",
-    "\n",
-    "x_f = np.fft.rfft(x)\n",
-    "y_f = np.fft.rfft(y)\n",
-    "etfe = y_f / x_f\n",
-    "\n",
-    "etfe_mean = np.mean(etfe, 0)"
-   ],
-   "metadata": {
-    "collapsed": false
-   }
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
-   "source": [
-    "marker_size = 1\n",
-    "marker_size_highlight = 4\n",
-    "line_width_highlight = 0.5\n",
-    "_, axs = plt.subplots(2, 2, sharex='col')\n",
-    "\n",
-    "axs[0, 0].plot(f, etfe.real[1:, :].T, '.', ms=marker_size, color='gray')\n",
-    "axs[0, 0].plot(f, etfe_mean.real, 'k-')\n",
-    "axs[0, 0].plot(\n",
-    "    f,\n",
-    "    etfe.real[0, :],\n",
-    "    'r.-',\n",
-    "    ms=marker_size_highlight,\n",
-    "    lw=line_width_highlight,\n",
-    ")\n",
-    "\n",
-    "axs[0, 0].text(0.025, 0.9, '(a)', transform=axs[0, 0].transAxes)\n",
-    "\n",
-    "axs[1, 0].plot(f, etfe.imag[1:, :].T, '.', ms=marker_size, color='gray')\n",
-    "axs[1, 0].plot(f, etfe_mean.imag, 'k-')\n",
-    "axs[1, 0].plot(\n",
-    "    f,\n",
-    "    etfe.imag[0, :],\n",
-    "    'r.-',\n",
-    "    ms=marker_size_highlight,\n",
-    "    lw=line_width_highlight,\n",
-    ")\n",
-    "\n",
-    "axs[1, 0].text(0.025, 0.9, '(b)', transform=axs[1, 0].transAxes)\n",
-    "\n",
-    "axs[0, 1].plot(f, etfe.real[1:, :].T, '.', ms=marker_size, color='gray')\n",
-    "axs[0, 1].plot(f, etfe_mean.real, 'k-')\n",
-    "axs[0, 1].plot(\n",
-    "    f,\n",
-    "    etfe.real[0, :],\n",
-    "    'r.-',\n",
-    "    ms=marker_size_highlight,\n",
-    "    lw=line_width_highlight,\n",
-    ")\n",
-    "\n",
-    "axs[0, 1].text(0.025, 0.9, '(c)', transform=axs[0, 1].transAxes)\n",
-    "\n",
-    "axs[1, 1].plot(f, etfe.imag[1:, :].T, '.', ms=marker_size, color='gray')\n",
-    "axs[1, 1].plot(f, etfe_mean.imag, 'k-')\n",
-    "axs[1, 1].plot(\n",
-    "    f,\n",
-    "    etfe.imag[0, :],\n",
-    "    'r.-',\n",
-    "    ms=marker_size_highlight,\n",
-    "    lw=line_width_highlight,\n",
-    ")\n",
-    "\n",
-    "axs[1, 1].text(0.025, 0.9, '(d)', transform=axs[1, 1].transAxes)\n",
-    "\n",
-    "axs[0, 0].set_ylim(-1.75, 2.75)\n",
-    "axs[0, 1].set_ylim(-0.6, 0.6)\n",
-    "axs[1, 0].set_ylim(-2.25, 2.25)\n",
-    "axs[1, 1].set_ylim(-0.6, 0.6)\n",
-    "\n",
-    "axs[0, 0].set_yticks(np.arange(-1, 3))\n",
-    "axs[0, 1].set_yticks(np.linspace(-0.5, 0.5, 3, endpoint=True))\n",
-    "axs[1, 0].set_yticks(np.arange(-2, 3))\n",
-    "axs[1, 1].set_yticks(np.linspace(-0.5, 0.5, 3, endpoint=True))\n",
-    "\n",
-    "axs[0, 0].set_ylabel(r'$\\mathrm{Re}\\{\\hat{H}_\\mathrm{ETFE}\\}$')\n",
-    "axs[1, 0].set_ylabel(r'$\\mathrm{Im}\\{\\hat{H}_\\mathrm{ETFE}\\}$')\n",
-    "\n",
-    "axs[1, 0].set_xlim(0, 10)\n",
-    "axs[1, 1].set_xlim(0, 4)\n",
-    "\n",
-    "axs[1, 0].set_xlabel('Frequency (THz)')\n",
-    "axs[1, 1].set_xlabel('Frequency (THz)')\n",
-    "\n",
-    "plt.show()"
-   ],
-   "metadata": {
-    "collapsed": false
-   }
-  },
-  {
-   "cell_type": "markdown",
-   "source": [
-    "## Fits to the empirical transfer function estimate\n",
-    "\n",
-    "The first cell below defines a transfer function model function `tf_model` and a residual function `res`, then computes the standard deviation of `etfe` as a function of frequency `f`. The transfer function model is\n",
-    " $$ H(\\mathbf{\\theta}; \\omega) = \\theta_0\\exp(i\\omega\\theta_1),$$\n",
-    " which rescales the input by $\\theta_0$ and shifts it by $\\theta_1$ (using the $-i\\omega t$ sign convention). In the simulations, $\\theta_0 = 0.5$ and $\\theta_1 = 1$.\n",
-    "\n",
-    "The second cell computes the parameter estimates `p_etfe` by performing a weighted least-squares fit to `etfe` with the transfer function model given above. The complex, frequency-domain fit residuals are stored in `residual_etfe` and the sum of the norm-squared residuals over all discrete frequencies is stored in `resnorm_etfe`. Note that `result.cost` returned by the SciPy `least_squares` function is equal to *half* of the sum of the squared residuals over the *given* frequencies, and the given frequencies include just those up to the Nyquist frequency, so we multiply `result.cost` by 4 and adjust for double-counting when `n` is even. The parameter covariance matrix is stored in `cv_etfe`. We use the SciPy `cumfreq` function to obtain the empirical cumulative distribution function `ecdf_etfe` of `resnorm_etfe` at the goodness-of-fit values `gof_etfe`.\n",
-    "\n",
-    "The third cell generates a figure that (a) compares the empirical cumulative distribution function `ecdf_etfe` to the $\\chi^2$ cumulative distribution function for `n - 2` degrees of freedom and (b) shows the real and imaginary parts of the normalized residuals for a fit that yields `resnorm_etfe` close to the median value. The inset of (a) shows these normalized residuals on a normal probability plot; the gray dashed line shows the standard normal distribution.\n",
-    "\n",
-    "The fourth cell compares the observed parameter distribution with the uncertainty estimates obtained from the parameter covariance matrix, and the fifth cell generates histograms of the parameter estimates."
-   ],
-   "metadata": {
-    "collapsed": false
-   }
+   "outputs": []
   },
   {
+   "metadata": {},
    "cell_type": "code",
    "execution_count": null,
-   "outputs": [],
    "source": [
-    "def etfe_res(_p, _f, _etfe, _sigma):\n",
-    "    _r = (_etfe - tf_model(_p, 2 * pi * _f)) / _sigma\n",
-    "    return np.concatenate((_r.real, _r.imag))\n",
-    "\n",
+    "thz.global_options.sampling_time = dt\n",
+    "t = thz.timebase(n)\n",
+    "mu = thz.wave(n)\n",
+    "sigma = noise_model.sigma_t(mu)\n",
     "\n",
-    "sigma_etfe = np.std(etfe, axis=0)"
+    "# Use broadcasting to generate multiple simulated noisy measurements of mu\n",
+    "n_sim = 1000"
    ],
-   "metadata": {
-    "collapsed": false
-   }
+   "outputs": []
   },
   {
+   "metadata": {},
    "cell_type": "code",
    "execution_count": null,
-   "outputs": [],
    "source": [
     "%%time\n",
     "\n",
-    "p_etfe = np.empty((2, n_sim))\n",
-    "resnorm_etfe = np.empty((n_sim,))\n",
-    "residual_etfe = np.empty((n_f, n_sim), dtype=np.complex128)\n",
-    "cv_etfe = np.empty((2, 2, n_sim))\n",
-    "for i in range(n_sim):\n",
-    "    result = opt.least_squares(\n",
-    "        lambda _p: etfe_res(_p, f, etfe[i, :], sigma_etfe), p_opt\n",
-    "    )\n",
-    "    p_etfe[:, i] = result.x\n",
-    "    residual_etfe[:, i] = result.fun[:n_f] + 1j * result.fun[n_f:]\n",
-    "    resnorm_etfe[i] = 4 * result.cost\n",
-    "    if n % 2 == 0:\n",
-    "        resnorm_etfe[i] -= np.abs(residual_etfe[-1, i]) ** 2\n",
-    "    r = np.linalg.qr(result.jac, mode='r')\n",
-    "    r_inv = np.linalg.inv(r)\n",
-    "    cv_etfe[:, :, i] = r_inv @ r_inv.T\n",
-    "    if not result.success:\n",
-    "        print(f\"Fit did not converge at i = {i}\")\n",
-    "\n",
-    "cf = stats.cumfreq(resnorm_etfe, 1000, defaultreallimits=(0, 1000))\n",
-    "gof_etfe = cf.lowerlimit + np.linspace(\n",
-    "    0, cf.binsize * cf.cumcount.size, cf.cumcount.size\n",
-    ")\n",
-    "ecdf_etfe = cf.cumcount / n_sim"
-   ],
-   "metadata": {
-    "collapsed": false
-   }
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
-   "source": [
-    "# Find the index of the median value\n",
-    "i_show = np.argsort(resnorm_etfe)[len(resnorm_etfe) // 2]\n",
-    "osm, osr = stats.probplot(\n",
-    "    np.concatenate(\n",
-    "        (residual_etfe[:, i_show].real, residual_etfe[:, i_show].imag)\n",
-    "    ),\n",
-    "    fit=False,\n",
-    ")\n",
-    "\n",
-    "w, h = figaspect(0.45)\n",
-    "_, axs = plt.subplots(1, 2, figsize=(w, h))\n",
-    "\n",
-    "# Panel (a)\n",
-    "axs[0].plot(\n",
-    "    gof_etfe, ecdf_etfe, 'k-', label=r'$S_\\mathrm{ETFE}$', clip_on=False\n",
-    ")\n",
-    "axs[0].plot(\n",
-    "    gof_etfe,\n",
-    "    stats.chi2.cdf(gof_etfe, n - 2),\n",
-    "    '--',\n",
-    "    c='gray',\n",
-    "    label=r'$\\chi^2$',\n",
-    "    clip_on=False,\n",
-    ")\n",
-    "axs[0].plot(np.median(resnorm_etfe), 0.5, 'rx')\n",
-    "\n",
-    "axs[0].legend(loc='upper right', frameon=False)\n",
-    "axs[0].text(0.025, 0.95, '(a)', transform=axs[0].transAxes)\n",
-    "\n",
-    "axs[0].set_xlim(0, 1000)\n",
-    "axs[0].set_ylim(0, 1)\n",
-    "\n",
-    "axs[0].set_xticks([0, 500, 1000])\n",
-    "axs[0].set_yticks(np.linspace(0, 1, 6, endpoint=True))\n",
-    "\n",
-    "axs[0].set_xlabel('GOF statistic')\n",
-    "axs[0].set_ylabel('Cumulative probabilty')\n",
-    "\n",
-    "# Inset to (a)\n",
-    "ins = axs[0].inset_axes([0.45, 0.2, 0.5, 0.5])\n",
-    "ins.plot(osr, osm, 'k.', ms=2)\n",
-    "ins.plot([-3, 3], [-3, 3], '--', c='gray')\n",
-    "ins.grid()\n",
-    "\n",
-    "ins.set_xlim(-3, 3)\n",
-    "ins.set_ylim(stats.norm.ppf([0.001, 0.999]))\n",
-    "\n",
-    "ins.set_xticks([-2, 0, 2])\n",
-    "ins.set_yticks(stats.norm.ppf([0.01, 0.1, 0.5, 0.9, 0.99]))\n",
-    "\n",
-    "ins.set_yticklabels(['0.01', '0.1', '0.5', '0.9', '0.99'])\n",
-    "\n",
-    "ins.set_xlabel('Normed residual')\n",
-    "ins.set_ylabel('Probability')\n",
-    "\n",
-    "# Panel (b)\n",
-    "markerline, stemlines, baseline = axs[1].stem(\n",
-    "    f,\n",
-    "    residual_etfe[:, i_show].real,\n",
-    "    linefmt='C0-',\n",
-    "    markerfmt='o',\n",
-    "    label='Real',\n",
-    ")\n",
-    "markerline.set_markerfacecolor('none')\n",
-    "markerline.set_markersize(4)\n",
-    "stemlines.set_linewidth(0.5)\n",
-    "baseline.set_linewidth(0)\n",
-    "\n",
-    "markerline, stemlines, baseline = axs[1].stem(\n",
-    "    f,\n",
-    "    residual_etfe[:, i_show].imag,\n",
-    "    linefmt='C1-',\n",
-    "    markerfmt='x',\n",
-    "    label='Imag',\n",
-    ")\n",
-    "markerline.set_markersize(4)\n",
-    "markerline.set_linewidth(0.5)\n",
-    "stemlines.set_linewidth(0.5)\n",
-    "baseline.set_linewidth(1)\n",
-    "\n",
-    "axs[1].legend(loc='upper right', frameon=False)\n",
-    "axs[1].text(0.025, 0.95, '(b)', transform=axs[1].transAxes)\n",
-    "\n",
-    "axs[1].set_xlim(0, 10)\n",
-    "axs[1].set_ylim(-3, 3)\n",
-    "\n",
-    "axs[1].set_xlabel('Frequency (THz)')\n",
-    "axs[1].set_ylabel('Normed residual')\n",
-    "\n",
-    "plt.show()"
-   ],
-   "metadata": {
-    "collapsed": false
-   }
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
-   "source": [
-    "print(\n",
-    "    f\"Amplitude parameter estimates: \"\n",
-    "    f\"{np.mean(p_etfe[0, :]):.4f} ± {np.std(p_etfe[0, :]):.4f}\"\n",
-    ")\n",
-    "print(\n",
-    "    f\"Delay parameter estimates:     \"\n",
-    "    f\"({np.mean(1e3 * p_etfe[1, :]):.2f} ± {1e3 * np.std(p_etfe[1, :]):.2f}) fs\"\n",
-    ")\n",
-    "print(\n",
-    "    f\"Mean uncertainty estimate for amplitude parameter: \"\n",
-    "    f\"{np.mean(np.sqrt(cv_etfe[0, 0, :])):.4f}\"\n",
-    ")\n",
-    "print(\n",
-    "    f\"Mean uncertainty estimate for delay parameter:     \"\n",
-    "    f\"{1e3 * np.mean(np.sqrt(cv_etfe[1, 1, :])):.2f} fs\"\n",
-    ")\n",
-    "print(\n",
-    "    f\"Standard deviation of the amplitude parameter uncertainty estimate: \"\n",
-    "    f\"{np.std(np.sqrt(cv_etfe[0, 0, :])):.2g}\"\n",
-    ")\n",
-    "print(\n",
-    "    f\"Standard deviation of the delay parameter uncertainty estimate:     \"\n",
-    "    f\"{1e3 * np.std(np.sqrt(cv_etfe[1, 1, :])):.2g} fs\"\n",
-    ")"
-   ],
-   "metadata": {
-    "collapsed": false
-   }
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
-   "source": [
-    "w, h = figaspect(0.45)\n",
-    "_, axs = plt.subplots(1, 2, figsize=(w, h), sharey=True)\n",
-    "\n",
-    "axs[0].hist(p_etfe[0, :], bins=20, range=(0.48, 0.52), ec='k')\n",
-    "\n",
-    "axs[0].set_xticks(np.arange(0.48, 0.521, 0.01))\n",
-    "axs[0].set_xlim(0.48, 0.52)\n",
-    "\n",
-    "axs[0].set_xlabel(r'$\\hat{\\theta}_0$')\n",
-    "axs[0].set_ylabel('Frequency')\n",
-    "\n",
-    "axs[1].hist(\n",
-    "    1e3 * (p_etfe[1, :] - p_opt[1]), bins=20, range=(-3, 3), ec='k', lw=0.5\n",
-    ")\n",
-    "\n",
-    "axs[1].set_xticks(np.arange(-3, 3.1, 1))\n",
-    "axs[1].set_xlim(-3, 3)\n",
-    "\n",
-    "axs[1].set_xlabel(r'$\\hat{\\theta}_1 - \\theta_1^*$ (fs)')\n",
-    "\n",
-    "plt.show()"
-   ],
-   "metadata": {
-    "collapsed": false
-   }
-  },
-  {
-   "cell_type": "markdown",
-   "source": [
-    "## Maximum-likelihood transfer function fits\n",
-    "\n",
-    "Functions from `thztools`: `costfunlsq`, `fftfreq`\n",
-    "\n",
-    "The first cell below fits the transfer function model `tf_model` to the same `x` and `y` data as above, but by optimizing the total-least-squares cost function `costfunlsq` for the time domain, which yields the maximum-likelihood parameter estimates `p_tls`. The time-domain fit residuals are stored in `residual_tls` and the sum of the squared residuals over all time points is stored in `resnorm_tls`. As with the fits to `etfe`, the `result.cost` returned by the SciPy `least_squares` function is equal to *half* of the sum of the squared residuals, so we multiply `result.cost` by 2 to obtain `resnorm_tls`. The parameter covariance matrix is stored in `cv_tls` and `ecdf_tls` is the empirical cumulative distribution function of `resnorm_tls` at the goodness-of-fit values `gof_etfe`.\n",
-    "\n",
-    "The second cell generates a figure that (a) compares the empirical cumulative distribution function `ecdf_tls` to the $\\chi^2$ cumulative distribution function for `n - 2` degrees of freedom and (b) shows the normalized residuals for a fit that yields `resnorm_tls` close to the median value. The inset of (a) shows these normalized residuals on a normal probability plot; the gray dashed line shows the standard normal distribution.\n",
-    "\n",
-    "The third cell compares the observed parameter distribution with the uncertainty estimates obtained from the parameter covariance matrix, and the fourth cell generates histograms of the parameter estimates."
-   ],
-   "metadata": {
-    "collapsed": false
-   }
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
-   "source": [
-    "p_tls = np.empty((n_p, n_sim))\n",
-    "resnorm_tls = np.empty((n_sim,))\n",
-    "delta_tls = np.empty((n, n_sim))\n",
-    "epsilon_tls = np.empty((n, n_sim))\n",
-    "cv_tls = np.empty((n_p, n_p, n_sim))\n",
+    "sigma_alpha_est = np.zeros(n_sim)\n",
+    "sigma_alpha_err = np.zeros(n_sim)\n",
+    "sigma_beta_est = np.zeros(n_sim)\n",
+    "sigma_beta_err = np.zeros(n_sim)\n",
+    "sigma_tau_est = np.zeros(n_sim)\n",
+    "sigma_tau_err = np.zeros(n_sim)\n",
+    "mu_est = np.zeros((n, n_sim))\n",
+    "mu_err = np.zeros((n, n_sim))\n",
+    "a_est = np.zeros((m, n_sim))\n",
+    "a_err = np.zeros((m, n_sim))\n",
+    "eta_est = np.zeros((m, n_sim))\n",
+    "eta_err = np.zeros((m, n_sim))\n",
     "for i in range(n_sim):\n",
-    "    sigma_x = thz.noiseamp(sigma_parms, x[i, :], ts)\n",
-    "    sigma_y = thz.noiseamp(sigma_parms, y[i, :], ts)\n",
-    "\n",
-    "    def jac_fun(_x):\n",
-    "        p_est = _x[:n_p]\n",
-    "        mu_est = x[i, :] + _x[n_p:]\n",
-    "        jac_tl = np.zeros((n, n_p))\n",
-    "        jac_tr = np.diag(-1 / sigma_x)\n",
-    "        jac_bl = -(\n",
-    "            np.stack(\n",
-    "                [\n",
-    "                    td_model([1, p_est[1]], mu_est),\n",
-    "                    irfft(\n",
-    "                        rfft(mu_est)\n",
-    "                        * (-1j * 2 * pi * f * tf_model(p_est, 2 * pi * f)),\n",
-    "                        n=n,\n",
-    "                    ),\n",
-    "                ],\n",
-    "                axis=0,\n",
-    "            )\n",
-    "            / sigma_y\n",
-    "        ).T\n",
-    "        jac_br = -(\n",
-    "            la.circulant(td_model(p_est, signal.unit_impulse(n))).T / sigma_y\n",
-    "        ).T\n",
-    "        jac_tot = np.block([[jac_tl, jac_tr], [jac_bl, jac_br]])\n",
-    "        return jac_tot\n",
-    "\n",
-    "    result = opt.least_squares(\n",
-    "        lambda _p: _costfuntls(\n",
-    "            tf_model,\n",
-    "            _p[:n_p],\n",
-    "            x[i, :] + _p[n_p:],\n",
-    "            x[i, :],\n",
-    "            y[i, :],\n",
-    "            sigma_x,\n",
-    "            sigma_y,\n",
-    "            ts,\n",
-    "        ),\n",
-    "        np.concatenate((p_opt, np.zeros(n))),\n",
-    "        jac=jac_fun,\n",
-    "        method='lm',\n",
-    "        x_scale=np.concatenate((np.array([1, 1]), sigma_x)),\n",
-    "    )\n",
-    "    p_tls[:, i] = result.x[:n_p]\n",
-    "    delta_tls[:, i] = result.fun[:n]\n",
-    "    epsilon_tls[:, i] = result.fun[n:]\n",
-    "    resnorm_tls[i] = 2 * result.cost\n",
-    "    r = np.linalg.qr(result.jac, mode='r')\n",
-    "    r_inv = np.linalg.inv(r)\n",
-    "    cv_tls[:, :, i] = (r_inv @ r_inv.T)[:n_p, :n_p]\n",
-    "    if not result.success:\n",
-    "        print(f\"Fit did not converge at i = {i}\")\n",
-    "\n",
-    "cf_tls = stats.cumfreq(resnorm_tls, 1000, defaultreallimits=(0, 1000))\n",
-    "gof_tls = cf_tls.lowerlimit + np.linspace(\n",
-    "    0, cf_tls.binsize * cf_tls.cumcount.size, cf_tls.cumcount.size\n",
-    ")\n",
-    "ecdf_tls = cf_tls.cumcount / n_sim"
-   ],
-   "metadata": {
-    "collapsed": false
-   }
+    "    x = mu + noise_model.noise_sim(np.tile(mu, (m, 1)), seed=i)\n",
+    "    fit_result = thz.noisefit(x.T)\n",
+    "    sigma_alpha_est[i] = fit_result.noise_model.sigma_alpha\n",
+    "    sigma_alpha_err[i] = fit_result.err_sigma_alpha\n",
+    "    sigma_beta_est[i] = fit_result.noise_model.sigma_beta\n",
+    "    sigma_beta_err[i] = fit_result.err_sigma_beta\n",
+    "    sigma_tau_est[i] = fit_result.noise_model.sigma_tau\n",
+    "    sigma_tau_err[i] = fit_result.err_sigma_tau\n",
+    "    mu_est[:, i] = fit_result.mu\n",
+    "    mu_err[:, i] = fit_result.err_mu\n",
+    "    a_est[:, i] = fit_result.a\n",
+    "    a_err[:, i] = fit_result.err_a\n",
+    "    eta_est[:, i] = fit_result.eta\n",
+    "    eta_err[:, i] = fit_result.err_eta\n",
+    "\n",
+    "sigma_correct = np.sqrt(m / (m - 1))"
+   ],
+   "outputs": []
   },
   {
+   "metadata": {},
    "cell_type": "code",
    "execution_count": null,
-   "outputs": [],
    "source": [
-    "# Find the index of the median value\n",
-    "i_show = np.argsort(resnorm_tls)[len(resnorm_tls) // 2]\n",
-    "norm_res_eps = np.sqrt(512 / 254) * epsilon_tls[:, i_show]\n",
-    "norm_res_delta = np.sqrt(512 / 254) * delta_tls[:, i_show]\n",
-    "osm, osr = stats.probplot(\n",
-    "    np.concatenate((norm_res_eps, norm_res_delta)), fit=False\n",
-    ")\n",
+    "sigma_alpha_est_mean = np.mean(sigma_alpha_est)\n",
+    "sigma_alpha_est_std = np.std(sigma_alpha_est, ddof=1)\n",
+    "sigma_alpha_err_mean = np.mean(sigma_alpha_err)\n",
     "\n",
-    "w, h = figaspect(0.45)\n",
-    "_, axd = plt.subplot_mosaic(\n",
-    "    [[\"CDF\", \"x-res\"], [\"CDF\", \"y-res\"]], figsize=(w, h)\n",
-    ")\n",
+    "sigma_beta_est_mean = np.mean(sigma_beta_est)\n",
+    "sigma_beta_est_std = np.std(sigma_beta_est, ddof=1)\n",
+    "sigma_beta_err_mean = np.mean(sigma_beta_err)\n",
     "\n",
-    "# Panel (a)\n",
-    "axd[\"CDF\"].plot(\n",
-    "    gof_tls, ecdf_tls, 'k-', label=r'$S_\\mathrm{TLS}$', clip_on=False\n",
-    ")\n",
-    "axd[\"CDF\"].plot(\n",
-    "    gof_tls,\n",
-    "    stats.chi2.cdf(gof_tls, n - n_p),\n",
-    "    '--',\n",
-    "    c='gray',\n",
-    "    label=r'$\\chi^2$',\n",
-    "    clip_on=False,\n",
-    ")\n",
-    "axd[\"CDF\"].plot(np.median(resnorm_tls), 0.5, 'rx')\n",
-    "\n",
-    "axd[\"CDF\"].legend(loc='upper right', frameon=False)\n",
-    "axd[\"CDF\"].text(0.025, 0.95, '(a)', transform=axd[\"CDF\"].transAxes)\n",
-    "\n",
-    "axd[\"CDF\"].set_xlim(0, 1000)\n",
-    "axd[\"CDF\"].set_ylim(0, 1)\n",
-    "\n",
-    "axd[\"CDF\"].set_xticks([0, 500, 1000])\n",
-    "axd[\"CDF\"].set_yticks(np.linspace(0, 1, 6, endpoint=True))\n",
-    "\n",
-    "axd[\"CDF\"].set_xlabel('GOF statistic')\n",
-    "axd[\"CDF\"].set_ylabel('Cumulative probabilty')\n",
-    "\n",
-    "# Inset to (a)\n",
-    "ins = axd[\"CDF\"].inset_axes([0.5, 0.2, 0.45, 0.45])\n",
-    "ins.plot(osr, osm, 'k.', ms=2)\n",
-    "ins.plot([-3, 3], [-3, 3], '--', c='gray')\n",
-    "ins.grid()\n",
-    "\n",
-    "ins.set_xlim(-5, 5)\n",
-    "ins.set_ylim(stats.norm.ppf([0.0005, 0.9995]))\n",
-    "\n",
-    "ins.set_xticks(np.arange(-4, 4.5, 2))\n",
-    "ins.set_yticks(stats.norm.ppf([0.005, 0.1, 0.5, 0.9, 0.995]))\n",
+    "sigma_tau_est_mean = np.mean(sigma_tau_est)\n",
+    "sigma_tau_est_std = np.std(sigma_tau_est, ddof=1)\n",
+    "sigma_tau_err_mean = np.mean(sigma_tau_err)\n",
     "\n",
-    "ins.set_yticklabels(['0.005', '0.1', '0.5', '0.9', '0.995'])\n",
-    "\n",
-    "ins.set_xlabel('Normed residuals\\n(x and y)')\n",
-    "ins.set_ylabel('Probability')\n",
-    "\n",
-    "# Panel (b)\n",
-    "markerline, stemlines, baseline = axd[\"x-res\"].stem(\n",
-    "    t, norm_res_delta, linefmt='k-', markerfmt='.'\n",
-    ")\n",
-    "markerline.set_markersize(2)\n",
-    "stemlines.set_linewidth(0.5)\n",
-    "baseline.set_linewidth(1)\n",
-    "\n",
-    "axd[\"x-res\"].text(0.025, 0.9, '(b)', transform=axd[\"x-res\"].transAxes)\n",
-    "\n",
-    "axd[\"x-res\"].set_xlim(0, 10)\n",
-    "axd[\"x-res\"].set_ylim(-3.5, 3.5)\n",
-    "\n",
-    "axd[\"x-res\"].set_xticklabels([])\n",
-    "axd[\"x-res\"].set_ylabel('Normed x-residual')\n",
-    "\n",
-    "# Panel (c)\n",
-    "markerline, stemlines, baseline = axd[\"y-res\"].stem(\n",
-    "    t, norm_res_eps, linefmt='k-', markerfmt='.'\n",
-    ")\n",
-    "markerline.set_markersize(2)\n",
-    "stemlines.set_linewidth(0.5)\n",
-    "baseline.set_linewidth(1)\n",
+    "mu_est_mean = np.mean(mu_est, axis=1)\n",
+    "mu_est_std = np.std(mu_est, axis=1, ddof=1)\n",
+    "mu_err_mean = np.mean(mu_err, axis=1)\n",
     "\n",
-    "axd[\"y-res\"].text(0.025, 0.9, '(c)', transform=axd[\"y-res\"].transAxes)\n",
+    "a_est_mean = np.mean(a_est, axis=1)\n",
+    "a_est_std = np.std(a_est, axis=1)\n",
+    "a_err_mean = np.mean(a_err, axis=1)\n",
     "\n",
-    "axd[\"y-res\"].set_xlim(0, 10)\n",
-    "axd[\"y-res\"].set_ylim(-3.5, 3.5)\n",
+    "eta_est_mean = np.mean(eta_est, axis=1)\n",
+    "eta_est_std = np.std(eta_est, axis=1, ddof=1)\n",
+    "eta_err_mean = np.mean(eta_err, axis=1)\n",
     "\n",
-    "axd[\"y-res\"].set_xlabel('Time (ps)')\n",
-    "axd[\"y-res\"].set_ylabel('Normed y-residual')\n",
-    "\n",
-    "plt.show()"
-   ],
-   "metadata": {
-    "collapsed": false
-   }
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
-   "source": [
-    "print(\n",
-    "    f\"Amplitude parameter estimates: \"\n",
-    "    f\"{np.mean(p_tls[0, :]):.4f} ± {np.std(p_tls[0, :]):.4f}\"\n",
-    ")\n",
-    "print(\n",
-    "    f\"Delay parameter estimates:     \"\n",
-    "    f\"({np.mean(1e3 * p_tls[1, :]):.2f} ± {1e3 * np.std(p_tls[1, :]):.2f}) fs\"\n",
-    ")\n",
-    "print(\n",
-    "    f\"Mean uncertainty estimate for amplitude parameter: \"\n",
-    "    f\"{np.mean(np.sqrt(cv_tls[0, 0, :])):.4f}\"\n",
-    ")\n",
     "print(\n",
-    "    f\"Mean uncertainty estimate for delay parameter:     \"\n",
-    "    f\"{1e3 * np.mean(np.sqrt(cv_tls[1, 1, :])):.2f} fs\"\n",
+    "    f\"sigma_alpha: {sigma_alpha_est_mean / sigma_alpha:.4f} \"\n",
+    "    f\"± {sigma_alpha_est_std / sigma_alpha:.4f}\"\n",
     ")\n",
     "print(\n",
-    "    f\"Standard deviation of the amplitude parameter uncertainty estimate: \"\n",
-    "    f\"{np.std(np.sqrt(cv_tls[0, 0, :])):.2g}\"\n",
+    "    f\"sigma_beta: {sigma_beta_est_mean / sigma_beta:.4f} \"\n",
+    "    f\"± {sigma_beta_est_std / sigma_beta:.4f}\"\n",
     ")\n",
     "print(\n",
-    "    f\"Standard deviation of the delay parameter uncertainty estimate:     \"\n",
-    "    f\"{1e3 * np.std(np.sqrt(cv_tls[1, 1, :])):.2g} fs\"\n",
+    "    f\"sigma_tau: {sigma_tau_est_mean / sigma_tau:.4f} \"\n",
+    "    f\"± {sigma_tau_est_std / sigma_tau:.4f}\"\n",
+    ")\n",
+    "\n",
+    "np.savez(\n",
+    "    \"sim_dump\",\n",
+    "    n=n,\n",
+    "    m=m,\n",
+    "    dt=dt,\n",
+    "    sigma_alpha=sigma_alpha,\n",
+    "    sigma_beta=sigma_beta,\n",
+    "    sigma_tau=sigma_tau,\n",
+    "    noise_model=noise_model,\n",
+    "    t=t,\n",
+    "    mu=mu,\n",
+    "    sigma=sigma,\n",
+    "    n_sim=n_sim,\n",
+    "    sigma_alpha_est=sigma_alpha_est,\n",
+    "    sigma_alpha_err=sigma_alpha_err,\n",
+    "    sigma_beta_est=sigma_beta_est,\n",
+    "    sigma_beta_err=sigma_beta_err,\n",
+    "    sigma_tau_est=sigma_tau_est,\n",
+    "    sigma_tau_err=sigma_tau_err,\n",
+    "    mu_est=mu_est,\n",
+    "    mu_err=mu_err,\n",
+    "    a_est=a_est,\n",
+    "    a_err=a_err,\n",
+    "    eta_est=eta_est,\n",
+    "    eta_err=eta_err,\n",
     ")"
    ],
-   "metadata": {
-    "collapsed": false
-   }
+   "outputs": []
   },
   {
+   "metadata": {},
    "cell_type": "code",
    "execution_count": null,
-   "outputs": [],
    "source": [
-    "w, h = figaspect(0.45)\n",
-    "_, axs = plt.subplots(1, 2, figsize=(w, h), sharey=True)\n",
-    "\n",
-    "axs[0].hist(p_tls[0, :], bins=20, range=(0.48, 0.52), ec='k')\n",
-    "\n",
-    "axs[0].set_xticks(np.arange(0.48, 0.521, 0.01))\n",
-    "axs[0].set_xlim(0.48, 0.52)\n",
-    "\n",
-    "axs[0].set_xlabel(r'$\\hat{\\theta}_0$')\n",
-    "axs[0].set_ylabel('Frequency')\n",
-    "\n",
-    "axs[1].hist(\n",
-    "    1e3 * (p_tls[1, :] - p_opt[1]), bins=20, range=(-3, 3), ec='k', lw=0.5\n",
-    ")\n",
-    "\n",
-    "axs[1].set_xticks(np.arange(-3, 3.1, 1))\n",
-    "axs[1].set_xlim(-3, 3)\n",
-    "\n",
-    "axs[1].set_xlabel(r'$\\hat{\\theta}_1 - \\theta_1^*$ (fs)')\n",
-    "\n",
+    "w, h = figaspect(3.0)\n",
+    "_, axs = plt.subplots(3, 1, figsize=(w, h))\n",
+    "axs[0].hist(sigma_correct * sigma_alpha_est / sigma_alpha)\n",
+    "axs[1].hist(sigma_correct * sigma_beta_est / sigma_beta)\n",
+    "axs[2].hist(sigma_correct * sigma_tau_est / sigma_tau)\n",
     "plt.show()"
    ],
-   "metadata": {
-    "collapsed": false
-   }
+   "outputs": []
   }
  ],
  "metadata": {
diff --git a/examples/simulation-example.ipynb b/examples/simulation-example.ipynb
index d301a05..551a681 100644
--- a/examples/simulation-example.ipynb
+++ b/examples/simulation-example.ipynb
@@ -11,8 +11,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 1,
    "source": [
     "from matplotlib import pyplot as plt\n",
     "from matplotlib.figure import figaspect\n",
@@ -31,8 +30,13 @@
     "from thztools.thztools import _costfuntls"
    ],
    "metadata": {
-    "collapsed": false
-   }
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T05:14:34.941285Z",
+     "start_time": "2024-01-09T05:14:34.290563Z"
+    }
+   },
+   "outputs": []
   },
   {
    "cell_type": "markdown",
@@ -51,47 +55,62 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
+   "execution_count": 2,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T05:14:34.946800Z",
+     "start_time": "2024-01-09T05:14:34.943133Z"
+    }
+   },
    "source": [
     "n = 256  # Number of samples\n",
-    "ts = 0.05  # Sampling time [ps]\n",
-    "t0 = 2.5  # Peak pulse time [ps]\n",
-    "fs = 1 / ts  # Sampling frequency [THz]\n",
+    "dt = 0.05  # Sampling time [ps]\n",
+    "fs = 1 / dt  # Sampling frequency [THz]\n",
     "\n",
     "sigma_alpha = 1e-4  # Additive noise amplitude [signal units]\n",
     "sigma_beta = 1e-2  # Multiplicative noise amplitude [dimensionless]\n",
     "sigma_tau = 1e-3  # Time base noise amplitude [ps]\n",
-    "sigma_parms = np.array([sigma_alpha, sigma_beta, sigma_tau])\n",
+    "noise_model = thz.NoiseModel(sigma_alpha, sigma_beta, sigma_tau)\n",
     "\n",
     "seed = 0  # RNG seed\n",
     "rng = default_rng(seed)"
-   ]
+   ],
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
+   "execution_count": 3,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T05:14:34.951307Z",
+     "start_time": "2024-01-09T05:14:34.946137Z"
+    }
+   },
    "source": [
-    "mu, t = thz.thzgen(n, ts, t0)\n",
-    "sigma = thz.noiseamp(sigma_parms, mu, ts)\n",
+    "thz.global_options.sampling_time = dt\n",
+    "t = thz.timebase(n)\n",
+    "mu = thz.wave(n)\n",
+    "sigma = noise_model.amplitude(mu)\n",
     "\n",
-    "# Use broadcasting to generate multiple simulated noisy measurements of y\n",
+    "# Use broadcasting to generate multiple simulated noisy measurements of mu\n",
     "n_sim = 11\n",
-    "x = mu + sigma * rng.standard_normal((n_sim, n))\n",
+    "x = mu + noise_model.noise(np.tile(mu, (n_sim, 1)), seed=1234)\n",
     "\n",
     "# Compute power spectral density\n",
     "f, psd = signal.periodogram(x, fs=fs, detrend=False)\n",
     "n_f = f.size"
-   ]
+   ],
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
+   "execution_count": 4,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T05:14:35.133633Z",
+     "start_time": "2024-01-09T05:14:34.950966Z"
+    }
+   },
    "source": [
     "w, h = figaspect(0.45)\n",
     "_, axs = plt.subplots(1, 2, figsize=(w, h))\n",
@@ -119,7 +138,8 @@
     "axs[1].set_ylim(-70, 10)\n",
     "\n",
     "plt.show()"
-   ]
+   ],
+   "outputs": []
   },
   {
    "cell_type": "markdown",
@@ -136,16 +156,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 5,
    "source": [
     "# Clear x and n_sim for reuse\n",
     "del x, n_sim\n",
-    "rng = default_rng(seed)\n",
     "\n",
     "n_sim = 250\n",
-    "x = mu + sigma * rng.standard_normal((n_sim, n))\n",
-    "y = mu + sigma * rng.standard_normal((n_sim, n))\n",
+    "x = mu + noise_model.noise(np.tile(mu, (n_sim, 1)), seed=1234)\n",
+    "y = mu + noise_model.noise(np.tile(mu, (n_sim, 1)), seed=5678)\n",
     "\n",
     "x_f = np.fft.rfft(x)\n",
     "y_f = np.fft.rfft(y)\n",
@@ -154,13 +172,17 @@
     "etfe_mean = np.mean(etfe, 0)"
    ],
    "metadata": {
-    "collapsed": false
-   }
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T05:14:35.143236Z",
+     "start_time": "2024-01-09T05:14:35.136427Z"
+    }
+   },
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 6,
    "source": [
     "marker_size = 1\n",
     "marker_size_highlight = 4\n",
@@ -237,8 +259,13 @@
     "plt.show()"
    ],
    "metadata": {
-    "collapsed": false
-   }
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T05:14:35.491921Z",
+     "start_time": "2024-01-09T05:14:35.147099Z"
+    }
+   },
+   "outputs": []
   },
   {
    "cell_type": "markdown",
@@ -261,8 +288,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 7,
    "source": [
     "def tf_model(_p, _w):\n",
     "    return _p[0] * np.exp(-1j * _w * _p[1])\n",
@@ -276,13 +302,17 @@
     "sigma_etfe = np.std(etfe, axis=0)"
    ],
    "metadata": {
-    "collapsed": false
-   }
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T05:14:35.496728Z",
+     "start_time": "2024-01-09T05:14:35.491713Z"
+    }
+   },
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 8,
    "source": [
     "%%time\n",
     "\n",
@@ -314,13 +344,17 @@
     "ecdf_etfe = cf.cumcount / n_sim"
    ],
    "metadata": {
-    "collapsed": false
-   }
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T05:14:35.611513Z",
+     "start_time": "2024-01-09T05:14:35.498113Z"
+    }
+   },
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 9,
    "source": [
     "# Find the index of the median value\n",
     "i_show = np.argsort(resnorm_etfe)[len(resnorm_etfe) // 2]\n",
@@ -414,13 +448,17 @@
     "plt.show()"
    ],
    "metadata": {
-    "collapsed": false
-   }
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T05:14:35.735426Z",
+     "start_time": "2024-01-09T05:14:35.622249Z"
+    }
+   },
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 10,
    "source": [
     "print(\n",
     "    f\"Amplitude parameter estimates: \"\n",
@@ -448,13 +486,17 @@
     ")"
    ],
    "metadata": {
-    "collapsed": false
-   }
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T05:14:35.740297Z",
+     "start_time": "2024-01-09T05:14:35.735198Z"
+    }
+   },
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 11,
    "source": [
     "w, h = figaspect(0.45)\n",
     "_, axs = plt.subplots(1, 2, figsize=(w, h), sharey=True)\n",
@@ -477,8 +519,13 @@
     "plt.show()"
    ],
    "metadata": {
-    "collapsed": false
-   }
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T05:14:35.835719Z",
+     "start_time": "2024-01-09T05:14:35.741776Z"
+    }
+   },
+   "outputs": []
   },
   {
    "cell_type": "markdown",
@@ -499,8 +546,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 12,
    "source": [
     "%%time\n",
     "\n",
@@ -516,8 +562,8 @@
     "epsilon_tls = np.empty((n, n_sim))\n",
     "cv_tls = np.empty((n_p, n_p, n_sim))\n",
     "for i in range(n_sim):\n",
-    "    sigma_x = thz.noiseamp(sigma_parms, x[i, :], ts)\n",
-    "    sigma_y = thz.noiseamp(sigma_parms, y[i, :], ts)\n",
+    "    sigma_x = noise_model.amplitude(x[i, :])\n",
+    "    sigma_y = noise_model.amplitude(y[i, :])\n",
     "\n",
     "    def jac_fun(_x):\n",
     "        p_est = _x[:n_p]\n",
@@ -553,7 +599,7 @@
     "            y[i, :],\n",
     "            sigma_x,\n",
     "            sigma_y,\n",
-    "            ts,\n",
+    "            dt,\n",
     "        ),\n",
     "        np.concatenate((p_opt, np.zeros(n))),\n",
     "        jac=jac_fun,\n",
@@ -577,13 +623,17 @@
     "ecdf_tls = cf_tls.cumcount / n_sim"
    ],
    "metadata": {
-    "collapsed": false
-   }
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T05:15:28.413774Z",
+     "start_time": "2024-01-09T05:14:35.834404Z"
+    }
+   },
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 13,
    "source": [
     "# Find the index of the median value\n",
     "i_show = np.argsort(resnorm_tls)[len(resnorm_tls) // 2]\n",
@@ -676,13 +726,17 @@
     "plt.show()"
    ],
    "metadata": {
-    "collapsed": false
-   }
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T05:15:28.609093Z",
+     "start_time": "2024-01-09T05:15:28.414312Z"
+    }
+   },
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 14,
    "source": [
     "print(\n",
     "    f\"Amplitude parameter estimates: \"\n",
@@ -710,13 +764,17 @@
     ")"
    ],
    "metadata": {
-    "collapsed": false
-   }
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T05:15:28.616727Z",
+     "start_time": "2024-01-09T05:15:28.611144Z"
+    }
+   },
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 15,
    "source": [
     "w, h = figaspect(0.45)\n",
     "_, axs = plt.subplots(1, 2, figsize=(w, h), sharey=True)\n",
@@ -739,8 +797,13 @@
     "plt.show()"
    ],
    "metadata": {
-    "collapsed": false
-   }
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T05:15:28.713957Z",
+     "start_time": "2024-01-09T05:15:28.642717Z"
+    }
+   },
+   "outputs": []
   }
  ],
  "metadata": {
diff --git a/examples/tdnll-scaled-dev.ipynb b/examples/tdnll-scaled-dev.ipynb
index 9c8a3ff..64d0a9e 100644
--- a/examples/tdnll-scaled-dev.ipynb
+++ b/examples/tdnll-scaled-dev.ipynb
@@ -12,8 +12,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 1,
    "source": [
     "from __future__ import annotations\n",
     "\n",
@@ -29,12 +28,17 @@
     "from numpy.random import default_rng\n",
     "\n",
     "import thztools as thz\n",
-    "from thztools.thztools import _tdnll_scaled as tdnll_scaled"
+    "from thztools.thztools import _costfun_noisefit as costfun"
    ],
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T04:30:48.329605Z",
+     "start_time": "2024-01-09T04:30:47.174755Z"
+    }
    },
-   "id": "416f6fcfaf869f9d"
+   "id": "416f6fcfaf869f9d",
+   "outputs": []
   },
   {
    "cell_type": "markdown",
@@ -48,33 +52,42 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 2,
    "source": [
     "rng = np.random.default_rng(0)\n",
     "n = 256\n",
     "m = 64\n",
     "ts = 0.05\n",
-    "t = np.arange(n) * ts\n",
-    "mu, _ = thz.thzgen(n, ts=ts, t0=n * ts / 3)\n",
+    "thz.global_options.sampling_time = ts\n",
+    "t = thz.timebase(n)\n",
+    "mu = thz.wave(n)\n",
     "sigma = np.array([1e-5, 1e-2, 1e-3])\n",
-    "noise = thz.noiseamp(sigma, mu, ts) * rng.standard_normal((m, n))\n",
+    "noise_model = thz.NoiseModel(\n",
+    "    sigma_alpha=sigma[0], sigma_beta=sigma[1], sigma_tau=sigma[2]\n",
+    ")\n",
+    "noise = noise_model.noise((np.ones((m, 1)) * mu), seed=0)\n",
     "x = np.array(mu + noise)\n",
-    "delta = np.zeros(n)\n",
-    "alpha = np.zeros(m - 1)\n",
+    "delta_mu = np.zeros(n)\n",
+    "delta_a = np.zeros(m - 1)\n",
     "eta = np.zeros(m - 1)\n",
     "\n",
     "logv = np.log(sigma**2)\n",
-    "scale_logv = 1e-1 * np.ones(3)\n",
-    "scale_delta = 1e-1 * thz.noiseamp(sigma, mu, ts)\n",
-    "scale_alpha = 1e-1 * np.ones(m - 1)\n",
-    "scale_eta = 1e-2 * np.ones(m - 1)\n",
-    "scale_v = 1.0e-2"
+    "scale_sigma_alpha = noise_model.sigma_alpha\n",
+    "scale_sigma_beta = noise_model.sigma_beta\n",
+    "scale_sigma_tau = noise_model.sigma_tau\n",
+    "scale_delta_mu = 1e-0 * noise_model.amplitude(mu)\n",
+    "scale_delta_a = 1e-4 * np.ones(m - 1)\n",
+    "scale_eta = 1e-3 * np.ones(m - 1)"
    ],
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T04:30:48.333891Z",
+     "start_time": "2024-01-09T04:30:48.329647Z"
+    }
    },
-   "id": "6bc9be6f87caaf9f"
+   "id": "6bc9be6f87caaf9f",
+   "outputs": []
   },
   {
    "cell_type": "markdown",
@@ -88,290 +101,354 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 3,
    "source": [
-    "_, grad_delta_tdnll = tdnll_scaled(\n",
+    "_, grad_delta_mu_tdnll = costfun(\n",
     "    x,\n",
-    "    logv,\n",
-    "    delta,\n",
-    "    alpha,\n",
-    "    eta,\n",
-    "    ts,\n",
-    "    fix_logv=True,\n",
-    "    fix_delta=False,\n",
-    "    fix_alpha=True,\n",
+    "    logv[0],\n",
+    "    logv[1],\n",
+    "    logv[2],\n",
+    "    delta_mu,\n",
+    "    delta_a,\n",
+    "    eta / ts,\n",
+    "    fix_logv_alpha=True,\n",
+    "    fix_logv_beta=True,\n",
+    "    fix_logv_tau=True,\n",
+    "    fix_delta_mu=False,\n",
+    "    fix_delta_a=True,\n",
     "    fix_eta=True,\n",
-    "    scale_logv=np.ones(3),\n",
-    "    scale_delta=np.ones(n),\n",
-    "    scale_alpha=np.ones(m - 1),\n",
+    "    scale_sigma_alpha=1.0,\n",
+    "    scale_sigma_beta=1.0,\n",
+    "    scale_sigma_tau=1.0,\n",
+    "    scale_delta_mu=np.ones(n),\n",
+    "    scale_delta_a=np.ones(m - 1),\n",
     "    scale_eta=np.ones(m - 1),\n",
-    "    scale_v=1.0,\n",
     ")\n",
     "\n",
-    "_, grad_delta_tdnll_scaled = tdnll_scaled(\n",
+    "_, grad_delta_mu_costfun = costfun(\n",
     "    x,\n",
-    "    np.log(np.exp(logv) / scale_v) / scale_logv,\n",
-    "    delta / scale_delta,\n",
-    "    alpha,\n",
-    "    eta,\n",
-    "    ts,\n",
-    "    fix_logv=True,\n",
-    "    fix_delta=False,\n",
-    "    fix_alpha=True,\n",
+    "    np.log(np.exp(logv[0]) / scale_sigma_alpha**2),\n",
+    "    np.log(np.exp(logv[1]) / scale_sigma_beta**2),\n",
+    "    np.log(np.exp(logv[2]) / scale_sigma_tau**2),\n",
+    "    delta_mu / scale_delta_mu,\n",
+    "    delta_a,\n",
+    "    eta / ts,\n",
+    "    fix_logv_alpha=True,\n",
+    "    fix_logv_beta=True,\n",
+    "    fix_logv_tau=True,\n",
+    "    fix_delta_mu=False,\n",
+    "    fix_delta_a=True,\n",
     "    fix_eta=True,\n",
-    "    scale_logv=scale_logv,\n",
-    "    scale_delta=scale_delta,\n",
-    "    scale_alpha=np.ones(m - 1),\n",
+    "    scale_sigma_alpha=scale_sigma_alpha,\n",
+    "    scale_sigma_beta=scale_sigma_beta,\n",
+    "    scale_sigma_tau=scale_sigma_tau,\n",
+    "    scale_delta_mu=scale_delta_mu,\n",
+    "    scale_delta_a=np.ones(m - 1),\n",
     "    scale_eta=np.ones(m - 1),\n",
-    "    scale_v=scale_v,\n",
     ")\n",
     "\n",
-    "grad_delta_nd = nd.Gradient(\n",
-    "    lambda _delta: tdnll_scaled(\n",
+    "grad_delta_mu_nd = nd.Gradient(\n",
+    "    lambda _delta_mu: costfun(\n",
     "        x,\n",
-    "        np.log(np.exp(logv) / scale_v),\n",
-    "        _delta,\n",
-    "        alpha,\n",
-    "        eta,\n",
-    "        ts,\n",
-    "        fix_logv=True,\n",
-    "        fix_delta=True,\n",
-    "        fix_alpha=True,\n",
+    "        logv[0],\n",
+    "        logv[1],\n",
+    "        logv[2],\n",
+    "        _delta_mu,\n",
+    "        delta_a,\n",
+    "        eta / ts,\n",
+    "        fix_logv_alpha=True,\n",
+    "        fix_logv_beta=True,\n",
+    "        fix_logv_tau=True,\n",
+    "        fix_delta_mu=True,\n",
+    "        fix_delta_a=True,\n",
     "        fix_eta=True,\n",
-    "        scale_logv=np.ones(3),\n",
-    "        scale_delta=np.ones(n),\n",
-    "        scale_alpha=np.ones(m - 1),\n",
+    "        scale_sigma_alpha=1.0,\n",
+    "        scale_sigma_beta=1.0,\n",
+    "        scale_sigma_tau=1.0,\n",
+    "        scale_delta_mu=np.ones(n),\n",
+    "        scale_delta_a=np.ones(m - 1),\n",
     "        scale_eta=np.ones(m - 1),\n",
-    "        scale_v=scale_v,\n",
     "    )[0],\n",
     "    step=1e-6,\n",
-    ")(delta)\n",
+    ")(delta_mu)\n",
     "\n",
     "np.stack(\n",
     "    (\n",
-    "        grad_delta_tdnll,\n",
-    "        grad_delta_tdnll_scaled / scale_v / scale_delta,\n",
-    "        grad_delta_nd / scale_v,\n",
+    "        grad_delta_mu_tdnll,\n",
+    "        grad_delta_mu_costfun / scale_delta_mu,\n",
+    "        grad_delta_mu_nd,\n",
     "    )\n",
     ").T"
    ],
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T04:30:48.525146Z",
+     "start_time": "2024-01-09T04:30:48.336677Z"
+    }
    },
-   "id": "e626769313b54ca5"
+   "id": "e626769313b54ca5",
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 4,
    "source": [
-    "plt.plot(t, grad_delta_tdnll)\n",
-    "plt.plot(t, grad_delta_tdnll_scaled / scale_v / scale_delta)\n",
+    "plt.plot(t, grad_delta_mu_tdnll)\n",
+    "plt.plot(t, grad_delta_mu_costfun / scale_delta_mu)\n",
     "plt.show()\n",
-    "plt.plot(t, grad_delta_tdnll - grad_delta_tdnll_scaled / scale_v / scale_delta)\n",
+    "plt.plot(t, grad_delta_mu_tdnll - grad_delta_mu_costfun / scale_delta_mu)\n",
     "plt.show()"
    ],
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T04:30:48.630961Z",
+     "start_time": "2024-01-09T04:30:48.525289Z"
+    }
    },
-   "id": "b4550f1fd812049"
+   "id": "b4550f1fd812049",
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 5,
    "source": [
-    "val, grad_logv_tdnll = tdnll_scaled(\n",
+    "val, grad_logv_tdnll = costfun(\n",
     "    x,\n",
-    "    logv,\n",
-    "    delta,\n",
-    "    alpha,\n",
-    "    eta,\n",
-    "    ts,\n",
-    "    fix_logv=False,\n",
-    "    fix_delta=True,\n",
-    "    fix_alpha=True,\n",
+    "    logv[0],\n",
+    "    logv[1],\n",
+    "    logv[2],\n",
+    "    delta_mu,\n",
+    "    delta_a,\n",
+    "    eta / ts,\n",
+    "    fix_logv_alpha=False,\n",
+    "    fix_logv_beta=False,\n",
+    "    fix_logv_tau=False,\n",
+    "    fix_delta_mu=True,\n",
+    "    fix_delta_a=True,\n",
     "    fix_eta=True,\n",
-    "    scale_logv=np.ones(3),\n",
-    "    scale_delta=scale_delta,\n",
-    "    scale_alpha=np.ones(m - 1),\n",
+    "    scale_sigma_alpha=1.0,\n",
+    "    scale_sigma_beta=1.0,\n",
+    "    scale_sigma_tau=1.0,\n",
+    "    scale_delta_mu=scale_delta_mu,\n",
+    "    scale_delta_a=np.ones(m - 1),\n",
     "    scale_eta=np.ones(m - 1),\n",
-    "    scale_v=1.0,\n",
     ")\n",
     "\n",
-    "val_scaled_v, grad_logv_tdnll_scaled_v = tdnll_scaled(\n",
+    "val_scaled_v, grad_logv_costfun_v = costfun(\n",
     "    x,\n",
-    "    np.log(np.exp(logv) / scale_v),\n",
-    "    delta,\n",
-    "    alpha,\n",
-    "    eta,\n",
-    "    ts,\n",
-    "    fix_logv=False,\n",
-    "    fix_delta=True,\n",
-    "    fix_alpha=True,\n",
+    "    logv[0],\n",
+    "    logv[1],\n",
+    "    logv[2],\n",
+    "    delta_mu,\n",
+    "    delta_a,\n",
+    "    eta / ts,\n",
+    "    fix_logv_alpha=False,\n",
+    "    fix_logv_beta=False,\n",
+    "    fix_logv_tau=False,\n",
+    "    fix_delta_mu=True,\n",
+    "    fix_delta_a=True,\n",
     "    fix_eta=True,\n",
-    "    scale_logv=np.ones(3),\n",
-    "    scale_delta=scale_delta,\n",
-    "    scale_alpha=np.ones(m - 1),\n",
+    "    scale_sigma_alpha=1.0,\n",
+    "    scale_sigma_beta=1.0,\n",
+    "    scale_sigma_tau=1.0,\n",
+    "    scale_delta_mu=scale_delta_mu,\n",
+    "    scale_delta_a=np.ones(m - 1),\n",
     "    scale_eta=np.ones(m - 1),\n",
-    "    scale_v=scale_v,\n",
     ")\n",
     "\n",
-    "val_scaled, grad_logv_tdnll_scaled = tdnll_scaled(\n",
+    "val_scaled, grad_logv_costfun = costfun(\n",
     "    x,\n",
-    "    np.log(np.exp(logv) / scale_v) / scale_logv,\n",
-    "    delta,\n",
-    "    alpha,\n",
-    "    eta,\n",
-    "    ts,\n",
-    "    fix_logv=False,\n",
-    "    fix_delta=True,\n",
-    "    fix_alpha=True,\n",
+    "    np.log(np.exp(logv[0]) / scale_sigma_alpha**2),\n",
+    "    np.log(np.exp(logv[1]) / scale_sigma_beta**2),\n",
+    "    np.log(np.exp(logv[2]) / scale_sigma_tau**2),\n",
+    "    delta_mu,\n",
+    "    delta_a,\n",
+    "    eta / ts,\n",
+    "    fix_logv_alpha=False,\n",
+    "    fix_logv_beta=False,\n",
+    "    fix_logv_tau=False,\n",
+    "    fix_delta_mu=True,\n",
+    "    fix_delta_a=True,\n",
     "    fix_eta=True,\n",
-    "    scale_logv=scale_logv,\n",
-    "    scale_delta=scale_delta,\n",
-    "    scale_alpha=np.ones(m - 1),\n",
+    "    scale_sigma_alpha=scale_sigma_alpha,\n",
+    "    scale_sigma_beta=scale_sigma_beta,\n",
+    "    scale_sigma_tau=scale_sigma_tau,\n",
+    "    scale_delta_mu=scale_delta_mu,\n",
+    "    scale_delta_a=np.ones(m - 1),\n",
     "    scale_eta=np.ones(m - 1),\n",
-    "    scale_v=scale_v,\n",
     ")\n",
     "\n",
     "grad_logv_nd = nd.Gradient(\n",
-    "    lambda _logv: tdnll_scaled(\n",
+    "    lambda _logv: costfun(\n",
     "        x,\n",
-    "        np.log(np.exp(_logv) / scale_v) / scale_logv,\n",
-    "        delta,\n",
-    "        alpha,\n",
-    "        eta,\n",
-    "        ts,\n",
-    "        fix_logv=True,\n",
-    "        fix_delta=True,\n",
-    "        fix_alpha=True,\n",
+    "        np.log(np.exp(_logv[0]) / scale_sigma_alpha**2),\n",
+    "        np.log(np.exp(_logv[1]) / scale_sigma_beta**2),\n",
+    "        np.log(np.exp(_logv[2]) / scale_sigma_tau**2),\n",
+    "        delta_mu,\n",
+    "        delta_a,\n",
+    "        eta / ts,\n",
+    "        fix_logv_alpha=True,\n",
+    "        fix_logv_beta=True,\n",
+    "        fix_logv_tau=True,\n",
+    "        fix_delta_mu=True,\n",
+    "        fix_delta_a=True,\n",
     "        fix_eta=True,\n",
-    "        scale_logv=scale_logv,\n",
-    "        scale_delta=scale_delta,\n",
-    "        scale_alpha=np.ones(m - 1),\n",
+    "        scale_sigma_alpha=scale_sigma_alpha,\n",
+    "        scale_sigma_beta=scale_sigma_beta,\n",
+    "        scale_sigma_tau=scale_sigma_tau,\n",
+    "        scale_delta_mu=scale_delta_mu,\n",
+    "        scale_delta_a=scale_delta_a,\n",
     "        scale_eta=np.ones(m - 1),\n",
-    "        scale_v=scale_v,\n",
     "    )[0]\n",
     ")(logv)\n",
     "\n",
     "print(f\"{val =}\")\n",
-    "print(f\"{val_scaled_v / scale_v =}\")\n",
-    "print(f\"{val_scaled / scale_v =}\")\n",
+    "print(f\"{val_scaled_v =}\")\n",
+    "print(f\"{val_scaled =}\")\n",
     "np.stack(\n",
     "    (\n",
     "        grad_logv_tdnll,\n",
-    "        grad_logv_tdnll_scaled_v,\n",
-    "        grad_logv_tdnll_scaled / scale_logv,\n",
+    "        grad_logv_costfun_v,\n",
+    "        grad_logv_costfun,\n",
     "        grad_logv_nd,\n",
     "    )\n",
     ").T"
    ],
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T04:30:48.700114Z",
+     "start_time": "2024-01-09T04:30:48.635896Z"
+    }
    },
-   "id": "26b4830fe79ebc8"
+   "id": "26b4830fe79ebc8",
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 6,
    "source": [
-    "_, grad_alpha_tdnll = tdnll_scaled(\n",
+    "_, grad_delta_a_tdnll = costfun(\n",
     "    x,\n",
-    "    logv,\n",
-    "    delta,\n",
-    "    alpha,\n",
-    "    eta,\n",
-    "    ts,\n",
-    "    fix_logv=True,\n",
-    "    fix_delta=True,\n",
-    "    fix_alpha=False,\n",
+    "    logv[0],\n",
+    "    logv[1],\n",
+    "    logv[2],\n",
+    "    delta_mu,\n",
+    "    delta_a,\n",
+    "    eta / ts,\n",
+    "    fix_logv_alpha=True,\n",
+    "    fix_logv_beta=True,\n",
+    "    fix_logv_tau=True,\n",
+    "    fix_delta_mu=True,\n",
+    "    fix_delta_a=False,\n",
     "    fix_eta=True,\n",
-    "    scale_logv=np.ones(3),\n",
-    "    scale_delta=scale_delta,\n",
-    "    scale_alpha=np.ones(m - 1),\n",
+    "    scale_sigma_alpha=1.0,\n",
+    "    scale_sigma_beta=1.0,\n",
+    "    scale_sigma_tau=1.0,\n",
+    "    scale_delta_mu=scale_delta_mu,\n",
+    "    scale_delta_a=np.ones(m - 1),\n",
     "    scale_eta=np.ones(m - 1),\n",
-    "    scale_v=scale_v,\n",
     ")\n",
     "\n",
-    "grad_alpha_nd = nd.Gradient(\n",
-    "    lambda _alpha: tdnll_scaled(\n",
+    "grad_delta_a_nd = nd.Gradient(\n",
+    "    lambda _delta_a: costfun(\n",
     "        x,\n",
-    "        logv,\n",
-    "        delta,\n",
-    "        _alpha,\n",
-    "        eta,\n",
-    "        ts,\n",
-    "        fix_logv=True,\n",
-    "        fix_delta=True,\n",
-    "        fix_alpha=True,\n",
+    "        logv[0],\n",
+    "        logv[1],\n",
+    "        logv[2],\n",
+    "        delta_mu,\n",
+    "        _delta_a,\n",
+    "        eta / ts,\n",
+    "        fix_logv_alpha=True,\n",
+    "        fix_logv_beta=True,\n",
+    "        fix_logv_tau=True,\n",
+    "        fix_delta_mu=True,\n",
+    "        fix_delta_a=True,\n",
     "        fix_eta=True,\n",
-    "        scale_logv=np.ones(3),\n",
-    "        scale_delta=scale_delta,\n",
-    "        scale_alpha=np.ones(m - 1),\n",
+    "        scale_sigma_alpha=1.0,\n",
+    "        scale_sigma_beta=1.0,\n",
+    "        scale_sigma_tau=1.0,\n",
+    "        scale_delta_mu=scale_delta_mu,\n",
+    "        scale_delta_a=np.ones(m - 1),\n",
     "        scale_eta=np.ones(m - 1),\n",
-    "        scale_v=scale_v,\n",
     "    )[0]\n",
-    ")(alpha)\n",
+    ")(delta_a)\n",
     "\n",
-    "np.stack((grad_alpha_tdnll, grad_alpha_nd)).T"
+    "np.stack((grad_delta_a_tdnll, grad_delta_a_nd)).T"
    ],
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T04:30:49.276220Z",
+     "start_time": "2024-01-09T04:30:48.672931Z"
+    }
    },
-   "id": "5136dc2abb3fcd50"
+   "id": "5136dc2abb3fcd50",
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 7,
    "source": [
-    "_, grad_eta_tdnll = tdnll_scaled(\n",
+    "_, grad_eta_tdnll = costfun(\n",
     "    x,\n",
-    "    logv,\n",
-    "    delta,\n",
-    "    alpha,\n",
-    "    eta,\n",
-    "    ts,\n",
-    "    fix_logv=True,\n",
-    "    fix_delta=True,\n",
-    "    fix_alpha=True,\n",
+    "    logv[0],\n",
+    "    logv[1],\n",
+    "    logv[2],\n",
+    "    delta_mu,\n",
+    "    delta_a,\n",
+    "    eta / ts,\n",
+    "    fix_logv_alpha=True,\n",
+    "    fix_logv_beta=True,\n",
+    "    fix_logv_tau=True,\n",
+    "    fix_delta_mu=True,\n",
+    "    fix_delta_a=True,\n",
     "    fix_eta=False,\n",
-    "    scale_logv=np.ones(3),\n",
-    "    scale_delta=scale_delta,\n",
-    "    scale_alpha=np.ones(m - 1),\n",
+    "    scale_sigma_alpha=1.0,\n",
+    "    scale_sigma_beta=1.0,\n",
+    "    scale_sigma_tau=1.0,\n",
+    "    scale_delta_mu=scale_delta_mu,\n",
+    "    scale_delta_a=np.ones(m - 1),\n",
     "    scale_eta=np.ones(m - 1),\n",
-    "    scale_v=scale_v,\n",
     ")\n",
     "\n",
     "grad_eta_nd = nd.Gradient(\n",
-    "    lambda _eta: tdnll_scaled(\n",
+    "    lambda _eta_scaled: costfun(\n",
     "        x,\n",
-    "        logv,\n",
-    "        delta,\n",
-    "        alpha,\n",
-    "        _eta,\n",
-    "        ts,\n",
-    "        fix_logv=True,\n",
-    "        fix_delta=True,\n",
-    "        fix_alpha=True,\n",
+    "        logv[0],\n",
+    "        logv[1],\n",
+    "        logv[2],\n",
+    "        delta_mu,\n",
+    "        delta_a,\n",
+    "        _eta_scaled,\n",
+    "        fix_logv_alpha=True,\n",
+    "        fix_logv_beta=True,\n",
+    "        fix_logv_tau=True,\n",
+    "        fix_delta_mu=True,\n",
+    "        fix_delta_a=True,\n",
     "        fix_eta=True,\n",
-    "        scale_logv=np.ones(3),\n",
-    "        scale_delta=scale_delta,\n",
-    "        scale_alpha=np.ones(m - 1),\n",
+    "        scale_sigma_alpha=1.0,\n",
+    "        scale_sigma_beta=1.0,\n",
+    "        scale_sigma_tau=1.0,\n",
+    "        scale_delta_mu=scale_delta_mu,\n",
+    "        scale_delta_a=np.ones(m - 1),\n",
     "        scale_eta=np.ones(m - 1),\n",
-    "        scale_v=scale_v,\n",
     "    )[0]\n",
     ")(eta)\n",
     "\n",
     "np.stack((grad_eta_tdnll, grad_eta_nd)).T"
    ],
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T04:30:49.887798Z",
+     "start_time": "2024-01-09T04:30:49.275409Z"
+    }
    },
-   "id": "340e8a4843e9d0a5"
+   "id": "340e8a4843e9d0a5",
+   "outputs": []
   },
   {
    "cell_type": "markdown",
@@ -385,108 +462,167 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 8,
    "source": [
-    "result = thz.tdnoisefit(x.T, v0=sigma**2, ts=ts, fix_a=False, fix_eta=False)"
+    "result = thz.noisefit(\n",
+    "    x.T,\n",
+    "    sigma_alpha0=sigma[0],\n",
+    "    sigma_beta0=sigma[1],\n",
+    "    sigma_tau0=sigma[2],\n",
+    "    dt=ts,\n",
+    "    fix_a=False,\n",
+    "    fix_eta=False,\n",
+    ")"
    ],
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T04:30:51.300096Z",
+     "start_time": "2024-01-09T04:30:49.887238Z"
+    }
    },
-   "id": "dc0d1efd00bb35d7"
+   "id": "dc0d1efd00bb35d7",
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 9,
    "source": [
-    "print(result[2][\"message\"])"
+    "print(result.diagnostic[\"message\"])"
    ],
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T04:30:51.304279Z",
+     "start_time": "2024-01-09T04:30:51.300273Z"
+    }
    },
-   "id": "445f8f0c1a609ac3"
+   "id": "445f8f0c1a609ac3",
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 10,
    "source": [
-    "var_out = result[0][\"var\"] * m / (m - 1)\n",
-    "var_err = result[2][\"err\"][\"var\"] * m / (m - 1)\n",
-    "for val_in, val_out, err in zip(sigma**2, var_out, var_err):\n",
-    "    print(f\"Input: {val_in:6.4g}\\t Output: {val_out:6.4g} ± {err:6.4g}\")"
+    "sigma_out = np.array(\n",
+    "    [\n",
+    "        result.noise_model.sigma_alpha,\n",
+    "        result.noise_model.sigma_beta,\n",
+    "        result.noise_model.sigma_tau,\n",
+    "    ]\n",
+    ") * np.sqrt(m / (m - 1))\n",
+    "sigma_err = np.sqrt(\n",
+    "    np.array(\n",
+    "        [result.err_sigma_alpha, result.err_sigma_beta, result.err_sigma_tau]\n",
+    "    )\n",
+    "    * m\n",
+    "    / (m - 1)\n",
+    ")\n",
+    "for sigma_in, sigma_out, err in zip(sigma, sigma_out, sigma_err):\n",
+    "    print(f\"Input: {sigma_in:6.4g}\\t Output: {sigma_out:6.4g} ± {err:6.4g}\")"
    ],
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T04:30:51.305622Z",
+     "start_time": "2024-01-09T04:30:51.303501Z"
+    }
    },
-   "id": "24c887e71ac696f6"
+   "id": "24c887e71ac696f6",
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 11,
    "source": [
-    "plt.plot(result[2][\"grad_scaled\"])\n",
+    "plt.plot(result.diagnostic[\"jac\"])\n",
     "plt.show()"
    ],
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T04:30:51.358621Z",
+     "start_time": "2024-01-09T04:30:51.312097Z"
+    }
    },
-   "id": "2461337648eae5fa"
+   "id": "2461337648eae5fa",
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 12,
    "source": [
-    "plt.semilogy(np.diag(result[2][\"hess_inv_scaled\"]))\n",
+    "plt.semilogy(np.diag(result.diagnostic[\"hess_inv\"]))\n",
     "plt.show()"
    ],
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T04:30:51.493972Z",
+     "start_time": "2024-01-09T04:30:51.368942Z"
+    }
    },
-   "id": "9fbe795d084c20db"
+   "id": "9fbe795d084c20db",
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 13,
    "source": [
-    "np.diag(result[2][\"hess_inv_scaled\"][:3, :3])"
+    "np.diag(result.diagnostic[\"hess_inv\"][:3, :3])"
    ],
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T04:30:51.498098Z",
+     "start_time": "2024-01-09T04:30:51.494797Z"
+    }
    },
-   "id": "d69c3719ad1e9a4a"
+   "id": "d69c3719ad1e9a4a",
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 14,
    "source": [
-    "plt.plot(t, np.log10(np.diag(result[2][\"hess_inv_scaled\"])[3 : 3 + n]))\n",
+    "plt.plot(t, np.log10(np.diag(result.diagnostic[\"hess_inv\"])[3 : 3 + n]))\n",
     "plt.plot(t, mu)\n",
     "plt.show()"
    ],
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T04:30:51.551680Z",
+     "start_time": "2024-01-09T04:30:51.502409Z"
+    }
    },
-   "id": "f75beb5c4868ea14"
+   "id": "f75beb5c4868ea14",
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 15,
    "source": [
-    "plt.plot(t, result[2][\"err\"][\"delta\"] * np.sqrt(m))\n",
-    "plt.plot(t, thz.noiseamp(np.sqrt(result[0][\"var\"]), result[0][\"mu\"], ts))\n",
+    "noise_model_fit = thz.NoiseModel(\n",
+    "    sigma_alpha=result.noise_model.sigma_alpha,\n",
+    "    sigma_beta=result.noise_model.sigma_beta,\n",
+    "    sigma_tau=result.noise_model.sigma_tau,\n",
+    ")\n",
+    "\n",
+    "plt.plot(t, result.err_mu)\n",
+    "plt.plot(t, noise_model_fit.amplitude(result.mu))\n",
+    "plt.plot(t, np.std(x, axis=0))\n",
     "plt.show()"
    ],
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T04:30:51.615113Z",
+     "start_time": "2024-01-09T04:30:51.549414Z"
+    }
    },
-   "id": "f8565beffa655791"
+   "id": "f8565beffa655791",
+   "outputs": []
   },
   {
    "cell_type": "markdown",
@@ -500,46 +636,78 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 16,
    "source": [
-    "result = thz.tdnoisefit(x.T, v0=sigma**2, ts=ts, fix_a=True, fix_eta=True)\n",
-    "print(result[2][\"message\"])"
+    "result = thz.noisefit(\n",
+    "    x.T,\n",
+    "    sigma_alpha0=sigma[0],\n",
+    "    sigma_beta0=sigma[1],\n",
+    "    sigma_tau0=sigma[2],\n",
+    "    dt=ts,\n",
+    "    fix_a=True,\n",
+    "    fix_eta=True,\n",
+    ")\n",
+    "print(result.diagnostic[\"message\"])"
    ],
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T04:30:52.006086Z",
+     "start_time": "2024-01-09T04:30:51.589023Z"
+    }
    },
-   "id": "460ada6a325cdef8"
+   "id": "460ada6a325cdef8",
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 17,
    "source": [
-    "var_out = result[0][\"var\"] * m / (m - 1)\n",
-    "var_err = result[2][\"err\"][\"var\"] * m / (m - 1)\n",
-    "for val_in, val_out, err in zip(sigma**2, var_out, var_err):\n",
-    "    print(f\"Input: {val_in:6.4g}\\t Output: {val_out:6.4g} ± {err:6.4g}\")"
+    "sigma_out = np.array(\n",
+    "    [\n",
+    "        result.noise_model.sigma_alpha,\n",
+    "        result.noise_model.sigma_beta,\n",
+    "        result.noise_model.sigma_tau,\n",
+    "    ]\n",
+    ") * np.sqrt(m / (m - 1))\n",
+    "sigma_err = np.sqrt(\n",
+    "    np.array(\n",
+    "        [result.err_sigma_alpha, result.err_sigma_beta, result.err_sigma_tau]\n",
+    "    )\n",
+    "    * m\n",
+    "    / (m - 1)\n",
+    ")\n",
+    "for sigma_in, sigma_out, err in zip(sigma, sigma_out, sigma_err):\n",
+    "    print(f\"Input: {sigma_in:6.4g}\\t Output: {sigma_out:6.4g} ± {err:6.4g}\")"
    ],
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T04:30:52.007733Z",
+     "start_time": "2024-01-09T04:30:52.005553Z"
+    }
    },
-   "id": "b209ca09ceff8cd2"
+   "id": "b209ca09ceff8cd2",
+   "outputs": []
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "outputs": [],
+   "execution_count": 18,
    "source": [
-    "plt.plot(t, result[2][\"err\"][\"delta\"] * np.sqrt(m))\n",
-    "plt.plot(t, thz.noiseamp(np.sqrt(result[0][\"var\"]), result[0][\"mu\"], ts))\n",
+    "plt.plot(t, result.err_mu)\n",
+    "plt.plot(t, noise_model_fit.amplitude(result.mu))\n",
     "plt.plot(t, np.std(x, axis=0))\n",
     "plt.show()"
    ],
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T04:30:52.051917Z",
+     "start_time": "2024-01-09T04:30:52.009585Z"
+    }
    },
-   "id": "e8e7a5403d19035b"
+   "id": "e8e7a5403d19035b",
+   "outputs": []
   }
  ],
  "metadata": {
diff --git a/paper/paper.bib b/paper/paper.bib
new file mode 100644
index 0000000..1f49b97
--- /dev/null
+++ b/paper/paper.bib
@@ -0,0 +1,47 @@
+@article{lee2023,
+  title = {The {{dotTHz Project}}: {{A Standard Data Format}} for {{Terahertz Time-Domain Data}}},
+  shorttitle = {The {{dotTHz Project}}},
+  author = {Lee, Jongmin and Leung, Chi Ki and Ma, Mingrui and {Ward-Berry}, Jasper and Santitewagun, Supawan and Zeitler, J. Axel},
+  year = {2023},
+  journal = {Journal of Infrared, Millimeter, and Terahertz Waves},
+  volume = {44},
+  number = {11},
+  pages = {795--813},
+  doi = {10.1007/s10762-023-00947-w}
+}
+
+@article{mohtashemi2021,
+  title = {Maximum-likelihood parameter estimation in terahertz time-domain spectroscopy},
+  author = {Mohtashemi, Laleh and Westlund, Paul and Sahota, Derek G. and Lea, Graham B. and Bushfield, Ian and Mousavi, Payam and Dodge, J. Steven},
+  year = {2021},
+  journal = {Optics Express},
+  volume = {29},
+  number = {4},
+  pages = {4912-4926},
+  publisher = {Optical Society of America},
+  doi = {10.1364/OE.417724}
+}
+
+@article{peretti2019,
+  title = {{{THz-TDS Time-Trace Analysis}} for the {{Extraction}} of {{Material}} and {{Metamaterial Parameters}}},
+  author = {Peretti, Romain and Mitryukovskiy, Sergey and Froberger, Kevin and Mebarki, Mohamed Aniss and Eliet, Sophie and Vanwolleghem, Mathias and Lampin, Jean-Francois},
+  year = {2019},
+  journal = {IEEE Transactions on Terahertz Science and Technology},
+  volume = {9},
+  number = {2},
+  pages = {136--149},
+  doi = {10.1109/TTHZ.2018.2889227}
+}
+
+@article{tayvah2021,
+  title = {Nelly: {{A User-Friendly}} and {{Open-Source Implementation}} of {{Tree-Based Complex Refractive Index Analysis}} for {{Terahertz Spectroscopy}}},
+  shorttitle = {Nelly},
+  author = {Tayvah, Uriel and Spies, Jacob A. and Neu, Jens and Schmuttenmaer, Charles A.},
+  year = {2021},
+  journal = {Analytical Chemistry},
+  volume = {93},
+  number = {32},
+  pages = {11243--11250},
+  publisher = {American Chemical Society},
+  doi = {10.1021/acs.analchem.1c02132}
+}
diff --git a/paper/paper.md b/paper/paper.md
new file mode 100644
index 0000000..b9c94c6
--- /dev/null
+++ b/paper/paper.md
@@ -0,0 +1,53 @@
+---
+title: 'THzTools: data analysis software for terahertz time-domain spectroscopy'
+
+tags:
+  - Python
+  - terahertz spectroscopy
+  - time-domain spectroscopy
+  - data analysis
+  - parameter estimation
+  - maximum-likelihood analysis
+
+authors:
+  - name: Jonathan Posada Loaiza
+    orcid: 0009-0004-4443-6161
+    equal-contrib: true
+    affiliation: 1
+  - name: Santiago Higuera-Quintero
+    orcid: 0000-0001-6905-4662
+    equal-contrib: true
+    affiliation: 1
+  - name: Alireza Noori
+    orcid: 0009-0009-1091-0802
+    equal-contrib: true
+    affiliation: 1
+  - name: Laleh Mohtashemi
+    orcid: 0000-0002-2706-857X
+    equal-contrib: true
+    affiliation: 1
+  - name: R. P. Hall
+    orcid: 0009-0004-1137-3651
+    affiliation: 1
+  - name: Naod Ayalew Yimam
+    orcid: 0009-0007-7522-5027
+    affiliation: 1
+  - name: J. Steven Dodge
+    orcid: 0000-0003-3219-7314
+    corresponding: true
+    affiliation: 1
+
+affiliations:
+ - name: Department of Physics, Simon Fraser University, Canada
+   index: 1
+
+date: 4 May 2024
+
+bibliography: paper.bib
+---
+# Summary
+Example text with citations [@lee2023; @mohtashemi2021; @peretti2019; @tayvah2021]
+
+# Statement of need
+
+# References
\ No newline at end of file
diff --git a/pyproject.toml b/pyproject.toml
index 08e1528..a77b9a1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -5,16 +5,21 @@ build-backend = "hatchling.build"
 [project]
 name = "thztools"
 dynamic = ["version"]
-description = "Tools for terahertz time-domain spectroscopy (THz-TDS)"
+description = """
+Data analysis software tools for terahertz time-domain spectroscopy (THz-TDS)
+"""
 readme = "README.md"
 requires-python = ">=3.8"
 license = "MIT"
 keywords = ["terahertz", "THz", "time-domain spectroscopy", "data analysis"]
 authors = [
-    { name="Steve Dodge", email="jsdodge@sfu.ca" },
-    { name="Santiago Higuera Quintero", email="s.higuera@uniandes.edu.co" },
-    { name="Jonathan Posada", email="jonathan.posada1@udea.edu.co" },
-    { name="Alireza Noori", email="alireza_noori@sfu.ca"}
+    { name="Jonathan Posada Loaiza", email="jonathan.posada1@udea.edu.co" },
+    { name="Santiago Higuera-Quintero", email="s.higuera@uniandes.edu.co" },
+    { name="Alireza Noori", email="alireza_noori@sfu.ca" },
+    { name="Laleh Mohtashemi", email="lalehmo@gmail.com" },
+    { name="R. P. Hall", email="rph4@sfu.ca" },
+    { name="Naod Ayalew Yimam", email="nayimam@gmail.com"},
+    { name="J. Steven Dodge", email="jsdodge@sfu.ca" },
 ]
 maintainers = [
     { name="Steve Dodge", email="jsdodge@sfu.ca" },
@@ -34,14 +39,14 @@ classifiers = [
     "Topic :: Scientific/Engineering",
 ]
 dependencies = [
-    "numpy",
+    "numpy>=1.21",
     "scipy",
-    "numdifftools",
 ]
 
 [project.urls]
 Documentation = "https://dodge-research-group.github.io/thztools/"
 Issues = "https://github.com/dodge-research-group/thztools/issues"
+Discussion = "https://github.com/dodge-research-group/thztools/discussions"
 Source = "https://github.com/dodge-research-group/thztools"
 
 [tool.hatch.version]
@@ -52,6 +57,7 @@ dependencies = [
     "coverage[toml]>=6.5",
     "jupyter",
     "matplotlib",
+    "numdifftools",
     "pytest",
 ]
 [tool.hatch.envs.default.scripts]
@@ -69,7 +75,7 @@ cov = [
 make-logo = ["python logo/make_logo.py; mv thztools_logo.svg docs/source/_static/;mv thztools_logo_dark.svg docs/source/_static/"]
 
 [[tool.hatch.envs.all.matrix]]
-python = ["3.8", "3.9", "3.10", "3.11"]
+python = ["3.8", "3.9", "3.10", "3.11", "3.12",]
 
 [tool.hatch.envs.lint]
 detached = true
@@ -188,14 +194,20 @@ addopts = [
 
 [tool.hatch.envs.docs]
 dependencies = [
+    "ipykernel",
+    "IPython",
     "matplotlib",
+    "nbsphinx",
     "numpydoc",
     "pydata-sphinx-theme",
     "sphinx",
+    "sphinx-design",
+    "sphinx-gallery",
     "tomli",
 ]
 [tool.hatch.envs.docs.scripts]
 build = "sphinx-build -v -b html docs/source docs/build/html"
+test = "sphinx-build -v -b doctest docs/source docs/doctest"
 
 [tool.mypy]
 plugins = "numpy.typing.mypy_plugin"
diff --git a/src/thztools/__about__.py b/src/thztools/__about__.py
index d7b30e1..6a9beea 100644
--- a/src/thztools/__about__.py
+++ b/src/thztools/__about__.py
@@ -1,1 +1,1 @@
-__version__ = "0.3.6"
+__version__ = "0.4.0"
diff --git a/src/thztools/__init__.py b/src/thztools/__init__.py
index c0f6c2f..cfdad3e 100644
--- a/src/thztools/__init__.py
+++ b/src/thztools/__init__.py
@@ -4,22 +4,30 @@
     NoiseModel,
     NoiseResult,
     fit,
-    global_options,
+    get_option,
+    noisefit,
+    options,
     scaleshift,
-    tdnoisefit,
-    transfer_out,
+    set_option,
+    timebase,
+    transfer,
     wave,
 )
 
+set_option("sampling_time", None)
+
 __all__ = [
     "GlobalOptions",
     "NoiseModel",
     "NoiseResult",
     "fit",
-    "global_options",
+    "get_option",
+    "noisefit",
+    "options",
     "scaleshift",
-    "tdnoisefit",
-    "transfer_out",
+    "set_option",
+    "timebase",
+    "transfer",
     "wave",
     "__version__",
 ]
diff --git a/src/thztools/thztools.py b/src/thztools/thztools.py
index a8af2a8..41c1637 100644
--- a/src/thztools/thztools.py
+++ b/src/thztools/thztools.py
@@ -2,17 +2,17 @@
 
 import warnings
 from dataclasses import dataclass
-from typing import Callable
+from typing import Any, Callable
 
 import numpy as np
 import scipy.linalg as la
 import scipy.optimize as opt
 from numpy.fft import irfft, rfft, rfftfreq
 from numpy.random import default_rng
-from numpy.typing import ArrayLike
+from numpy.typing import ArrayLike, NDArray
 from scipy import signal
+from scipy.optimize import OptimizeResult, minimize
 from scipy.optimize import approx_fprime as fprime
-from scipy.optimize import minimize
 
 NUM_NOISE_PARAMETERS = 3
 NUM_NOISE_DATA_DIMENSIONS = 2
@@ -23,9 +23,7 @@ class GlobalOptions:
     r"""
     Class for global options.
 
-    Currently, this just includes the ``sampling_time`` parameter.
-
-    Parameters
+    Attributes
     ----------
     sampling_time : float | None, optional
         Global sampling time, normally in picoseconds. When set to None, the
@@ -35,25 +33,59 @@ class GlobalOptions:
     sampling_time: float | None
 
 
-global_options = GlobalOptions(
-    sampling_time=None,
-)
+options = GlobalOptions
+
+
+def set_option(key: str, value: Any) -> None:
+    r"""
+    Set global option.
+
+    Parameters
+    ----------
+    key : str
+        Option name.
+    value : object
+        Option value.
+
+    Returns
+    -------
+    None
+        No return value.
+    """
+    setattr(options, key, value)
+
+
+def get_option(key: str) -> Any:
+    r"""
+    Get global option.
+
+    Parameters
+    ----------
+    key : str
+        Option name.
+
+    Returns
+    -------
+    val : object
+        Option value.
+    """
+    return getattr(options, key)
 
 
-def _validate_sampling_time(dt: float | None) -> float:
-    if dt is None and global_options.sampling_time is None:
-        dt_out = 1.0
-    if dt is None and global_options.sampling_time is not None:
-        dt_out = global_options.sampling_time
-    if dt is not None and global_options.sampling_time is None:
+def _assign_sampling_time(dt: float | None) -> float:
+    dt_out = 1.0
+    if dt is None and get_option("sampling_time") is not None:
+        dt_out = get_option("sampling_time")
+    elif dt is not None and get_option("sampling_time") is None:
         dt_out = dt
-    if dt is not None and global_options.sampling_time is not None:
-        if np.isclose(dt, global_options.sampling_time):
+    elif dt is not None and get_option("sampling_time") is not None:
+        if np.isclose(dt, get_option("sampling_time")):
             dt_out = dt
         else:
+            opt_sampling_time = get_option("sampling_time")
             msg = (
                 f"Input sampling time {dt=} conflicts with "
-                f"{global_options.sampling_time=}, using {dt=}"
+                f"{opt_sampling_time=}, using {dt=}"
             )
             warnings.warn(msg, category=UserWarning, stacklevel=2)
             dt_out = dt
@@ -63,11 +95,11 @@ def _validate_sampling_time(dt: float | None) -> float:
 @dataclass
 class NoiseModel:
     r"""
-    Terahertz noise model class.
+    Noise model class.
 
     For noise parameters :math:`\sigma_\alpha`, :math:`\sigma_\beta`,
     :math:`\sigma_\tau` and signal vector :math:`\boldsymbol{\mu}`, the
-    :math:`k`-th element of the time-domain noise variance
+    :math:`k`-th element of the time-domain noise variance vector
     :math:`\boldsymbol{\sigma}^2` is given by [1]_
 
     .. math:: \sigma_k^2 = \sigma_\alpha^2 + \sigma_\beta^2\mu_k^2 \
@@ -77,15 +109,29 @@ class NoiseModel:
 
     Parameters
     ----------
-    alpha : float
+    sigma_alpha : float
         Additive noise amplitude.
-    beta : float
+    sigma_beta : float
         Multiplicative noise amplitude.
-    tau : float
+    sigma_tau : float
         Timebase noise amplitude.
     dt : float or None, optional
-        Sampling time, normally in picoseconds. Default set to None, for which
-        the ``tau`` parameter is given in units of the sampling time.
+        Sampling time, normally in picoseconds. Default is None, which sets
+        the sampling time to ``thztools.options.sampling_time``. If both
+        :attr:`dt` and ``thztools.get_option("sampling_time")`` are ``None``,
+        the :attr:`sigma_tau` parameter is given in units of the sampling time.
+
+    Warns
+    -----
+    UserWarning
+        If ``thztools.options.sampling_time`` and the :attr:`dt`
+        parameter are both not ``None`` and are set to different ``float``
+        values, the function will set the sampling time to :attr:`dt` and raise
+        a :class:`UserWarning`.
+
+    See Also
+    --------
+    noisefit : Estimate noise model parameters.
 
     References
     ----------
@@ -104,33 +150,34 @@ class NoiseModel:
     its peak magnitude, :math:`\mu_0`. The noise variance is normalized to
     :math:`(\sigma_\beta\mu_0)^2`.
 
-    .. plot::
-       :include-source: True
-
-        >>> import matplotlib.pyplot as plt
-        >>> import thztools as thz
-        >>> n, dt, t0 = 256, 0.05, 2.5
-        >>> mu, t = thz.wave(n, dt, t0)
-        >>> alpha, beta, tau = 1e-4, 1e-2, 1e-3
-        >>> noise_mod = thz.NoiseModel(alpha=alpha, beta=beta, tau=tau,
-        ... dt=dt)
-        >>> var_t = noise_mod.variance(mu)
-
-        >>> _, axs = plt.subplots(2, 1, sharex=True, layout="constrained")
-        >>> axs[0].plot(t, var_t / beta**2)
-        >>> axs[0].set_ylabel(r"$\sigma^2/(\sigma_\beta\mu_0)^2$")
-        >>> axs[1].plot(t, mu)
-        >>> axs[1].set_ylabel(r"$\mu/\mu_0$")
-        >>> axs[1].set_xlabel("t (ps)")
-        >>> plt.show()
+    >>> import thztools as thz
+    >>> from matplotlib import pyplot as plt
+    >>> n, dt = 256, 0.05
+    >>> thz.set_option("sampling_time", dt)
+    >>> t = thz.timebase(n)
+    >>> mu = thz.wave(n)
+    >>> alpha, beta, tau = 1e-4, 1e-2, 1e-3
+    >>> noise_model = thz.NoiseModel(sigma_alpha=alpha, sigma_beta=beta,
+    ...  sigma_tau=tau)
+    >>> noise_variance = noise_model.noise_var(mu)
+
+    >>> _, axs = plt.subplots(2, 1, sharex=True, layout="constrained")
+    >>> axs[0].plot(t, noise_variance / beta**2)
+    >>> axs[0].set_ylabel(r"$\sigma^2/(\sigma_\beta\mu_0)^2$")
+    >>> axs[1].plot(t, mu)
+    >>> axs[1].set_ylabel(r"$\mu/\mu_0$")
+    >>> axs[1].set_xlabel("t (ps)")
+    >>> plt.show()
     """
-    alpha: float
-    beta: float
-    tau: float
+    sigma_alpha: float
+    sigma_beta: float
+    sigma_tau: float
     dt: float | None = None
 
     # noinspection PyShadowingNames
-    def variance(self, x: ArrayLike, *, axis: int = -1) -> np.ndarray:
+    def noise_var(
+        self, x: ArrayLike, *, axis: int = -1
+    ) -> NDArray[np.float64]:
         r"""
         Compute the time-domain noise variance.
 
@@ -144,7 +191,7 @@ def variance(self, x: ArrayLike, *, axis: int = -1) -> np.ndarray:
 
         Returns
         -------
-        var_t : ndarray
+        noise_variance : ndarray
             Time-domain noise variance, in signal units (squared).
 
         Notes
@@ -176,30 +223,27 @@ def variance(self, x: ArrayLike, *, axis: int = -1) -> np.ndarray:
         its peak magnitude, :math:`\mu_0`. The noise variance is normalized to
         :math:`(\sigma_\beta\mu_0)^2`.
 
-        .. plot::
-           :include-source: True
-
-            >>> import matplotlib.pyplot as plt
-            >>> import thztools as thz
-            >>> n, dt, t0 = 256, 0.05, 2.5
-            >>> mu, t = thz.wave(n, dt, t0)
-            >>> alpha, beta, tau = 1e-4, 1e-2, 1e-3
-            >>> noise_mod = thz.NoiseModel(alpha=alpha, beta=beta, tau=tau,
-            ... dt=dt)
-            >>> var_t = noise_mod.variance(mu)
-
-            >>> _, axs = plt.subplots(2, 1, sharex=True, layout="constrained")
-            >>> axs[0].plot(t, var_t / beta**2)
-            >>> axs[0].set_ylabel(r"$\sigma^2/(\sigma_\beta\mu_0)^2$")
-            >>> axs[1].plot(t, mu)
-            >>> axs[1].set_ylabel(r"$\mu/\mu_0$")
-            >>> axs[1].set_xlabel("t (ps)")
-            >>> plt.show()
+        >>> import thztools as thz
+        >>> from matplotlib import pyplot as plt
+        >>> n, dt = 256, 0.05
+        >>> thz.set_option("sampling_time", dt)
+        >>> t = thz.timebase(n)
+        >>> mu = thz.wave(n)
+        >>> alpha, beta, tau = 1e-4, 1e-2, 1e-3
+        >>> noise_model = thz.NoiseModel(sigma_alpha=alpha, sigma_beta=beta,
+        ... sigma_tau=tau)
+        >>> noise_variance = noise_model.noise_var(mu)
+
+        >>> _, axs = plt.subplots(2, 1, sharex=True, layout="constrained")
+        >>> axs[0].plot(t, noise_variance / beta**2)
+        >>> axs[0].set_ylabel(r"$\sigma^2/(\sigma_\beta\mu_0)^2$")
+        >>> axs[1].plot(t, mu)
+        >>> axs[1].set_ylabel(r"$\mu/\mu_0$")
+        >>> axs[1].set_xlabel("t (ps)")
+        >>> plt.show()
         """
-        dt = self.dt
-        if dt is None:
-            dt = 1.0
-        x = np.asarray(x)
+        dt = _assign_sampling_time(self.dt)
+        x = np.asarray(x, dtype=np.float64)
         axis = int(axis)
         if x.ndim > 1:
             if axis != -1:
@@ -209,16 +253,22 @@ def variance(self, x: ArrayLike, *, axis: int = -1) -> np.ndarray:
         w_scaled = 2 * np.pi * rfftfreq(n)
         xdot = irfft(1j * w_scaled * rfft(x), n=n) / dt
 
-        var_t = self.alpha**2 + (self.beta * x) ** 2 + (self.tau * xdot) ** 2
+        noise_variance = (
+            self.sigma_alpha**2
+            + (self.sigma_beta * x) ** 2
+            + (self.sigma_tau * xdot) ** 2
+        )
 
         if x.ndim > 1:
             if axis != -1:
-                var_t = np.moveaxis(var_t, -1, axis)
+                noise_variance = np.moveaxis(noise_variance, -1, axis)
 
-        return var_t
+        return noise_variance
 
     # noinspection PyShadowingNames
-    def amplitude(self, x: ArrayLike, *, axis: int = -1) -> np.ndarray:
+    def noise_amp(
+        self, x: ArrayLike, *, axis: int = -1
+    ) -> NDArray[np.float64]:
         r"""
         Compute the time-domain noise amplitude.
 
@@ -232,7 +282,7 @@ def amplitude(self, x: ArrayLike, *, axis: int = -1) -> np.ndarray:
 
         Returns
         -------
-        sigma_t : ndarray
+        noise_amplitude : ndarray
             Time-domain noise amplitude, in signal units.
 
         Notes
@@ -264,39 +314,37 @@ def amplitude(self, x: ArrayLike, *, axis: int = -1) -> np.ndarray:
         its peak magnitude, :math:`\mu_0`. The noise amplitude is normalized to
         :math:`\sigma_\beta\mu_0`.
 
-        .. plot::
-           :include-source: True
-
-            >>> import matplotlib.pyplot as plt
-            >>> import thztools as thz
-
-            >>> n, dt, t0 = 256, 0.05, 2.5
-            >>> mu, t = thz.wave(n, dt, t0)
-            >>> alpha, beta, tau = 1e-4, 1e-2, 1e-3
-            >>> noise_mod = thz.NoiseModel(alpha=alpha, beta=beta, tau=tau,
-            ... dt=dt)
-            >>> sigma_t = noise_mod.amplitude(mu)
-
-            >>> _, axs = plt.subplots(2, 1, sharex=True, layout="constrained")
-            >>> axs[0].plot(t, sigma_t / beta)
-            >>> axs[0].set_ylabel(r"$\sigma/(\sigma_\beta\mu_0)$")
-            >>> axs[1].plot(t, mu)
-            >>> axs[1].set_ylabel(r"$\mu/\mu_0$")
-            >>> axs[1].set_xlabel("t (ps)")
-            >>> plt.show()
+        >>> import thztools as thz
+        >>> from matplotlib import pyplot as plt
+        >>> n, dt = 256, 0.05
+        >>> thz.set_option("sampling_time", dt)
+        >>> t = thz.timebase(n)
+        >>> mu = thz.wave(n)
+        >>> alpha, beta, tau = 1e-4, 1e-2, 1e-3
+        >>> noise_model = thz.NoiseModel(sigma_alpha=alpha, sigma_beta=beta,
+        ... sigma_tau=tau)
+        >>> noise_amplitude = noise_model.noise_amp(mu)
+
+        >>> _, axs = plt.subplots(2, 1, sharex=True, layout="constrained")
+        >>> axs[0].plot(t, noise_amplitude / beta)
+        >>> axs[0].set_ylabel(r"$\sigma/(\sigma_\beta\mu_0)$")
+        >>> axs[1].plot(t, mu)
+        >>> axs[1].set_ylabel(r"$\mu/\mu_0$")
+        >>> axs[1].set_xlabel("t (ps)")
+        >>> plt.show()
         """
-        return np.sqrt(self.variance(x, axis=axis))
+        return np.sqrt(self.noise_var(x, axis=axis))
 
     # noinspection PyShadowingNames
-    def noise(
+    def noise_sim(
         self,
         x: ArrayLike,
         *,
         axis: int = -1,
         seed: int | None = None,
-    ) -> np.ndarray:
+    ) -> NDArray[np.float64]:
         r"""
-        Compute a time-domain noise array.
+        Simulate time-domain noise.
 
         Parameters
         ----------
@@ -310,8 +358,8 @@ def noise(
 
         Returns
         -------
-        noise : ndarray
-            Time-domain noise, in signal units.
+        noise_simulation : ndarray
+            Simulated time-domain noise, in signal units.
 
         Notes
         -----
@@ -340,34 +388,32 @@ def noise(
         :math:`\sigma_\alpha = 10^{-4}`, :math:`\sigma_\beta = 10^{-2}`,
         :math:`\sigma_\tau = 10^{-3}` and the simulated signal :math:`\mu(t)`.
 
-        .. plot::
-           :include-source: True
-
-            >>> import matplotlib.pyplot as plt
-            >>> import thztools as thz
-
-            >>> n, dt, t0 = 256, 0.05, 2.5
-            >>> mu, t = thz.wave(n, dt, t0)
-            >>> alpha, beta, tau = 1e-4, 1e-2, 1e-3
-            >>> noise_mod = thz.NoiseModel(alpha=alpha, beta=beta, tau=tau,
-            ... dt=dt)
-            >>> noise = noise_mod.noise(mu)
-
-            >>> _, axs = plt.subplots(2, 1, sharex=True, layout="constrained")
-            >>> axs[0].plot(t, noise / beta)
-            >>> axs[0].set_ylabel(r"$\sigma_\mu\epsilon/(\sigma_\beta\mu_0)$")
-            >>> axs[1].plot(t, mu)
-            >>> axs[1].set_ylabel(r"$\mu/\mu_0$")
-            >>> axs[1].set_xlabel("t (ps)")
-            >>> plt.show()
+        >>> import thztools as thz
+        >>> from matplotlib import pyplot as plt
+        >>> n, dt = 256, 0.05
+        >>> thz.set_option("sampling_time", dt)
+        >>> t = thz.timebase(n)
+        >>> mu = thz.wave(n)
+        >>> alpha, beta, tau = 1e-4, 1e-2, 1e-3
+        >>> noise_model = thz.NoiseModel(sigma_alpha=alpha, sigma_beta=beta,
+        ... sigma_tau=tau)
+        >>> noise = noise_model.noise_sim(mu, seed=1234)
+
+        >>> _, axs = plt.subplots(2, 1, sharex=True, layout="constrained")
+        >>> axs[0].plot(t, noise / beta)
+        >>> axs[0].set_ylabel(r"$\sigma_\mu\epsilon/(\sigma_\beta\mu_0)$")
+        >>> axs[1].plot(t, mu)
+        >>> axs[1].set_ylabel(r"$\mu/\mu_0$")
+        >>> axs[1].set_xlabel("t (ps)")
+        >>> plt.show()
         """
-        x = np.asarray(x)
+        x = np.asarray(x, dtype=np.float64)
         axis = int(axis)
         if x.ndim > 1:
             if axis != -1:
                 x = np.moveaxis(x, axis, -1)
 
-        amp = self.amplitude(x)
+        amp = self.noise_amp(x)
         rng = default_rng(seed)
         noise = amp * rng.standard_normal(size=x.shape)
         if x.ndim > 1:
@@ -381,63 +427,66 @@ def noise(
 @dataclass
 class NoiseResult:
     r"""
-    Represents the noise parameter estimate output.
+    Dataclass for the output of :func:`noisefit`.
 
     Parameters
     ----------
     noise_model : NoiseModel
-        Noise parameters, represented as a NoiseModel object.
+        Noise parameters, represented as a :class:`NoiseModel` object.
     mu : ndarray, shape (n,)
         Signal vector.
     a : ndarray, shape (m,)
-        Amplitude vector.
+        Signal amplitude drift vector.
     eta : ndarray, shape (m,)
-        Delay vector.
+        Signal delay drift vector.
     fval : float
-        Value of NLL cost function from FMINUNC.
-    diagnostic : dict
-        Dictionary containing diagnostic information:
-
-            grad_scaled : ndarray
-                Gradient of the scaled negative log-likelihood function with
-                respect to the scaled fit parameters.
-            hess_inv_scaled : ndarray
-                Inverse of the Hessian obtained from scipy.optimize.minimize
-                using the BFGS method, which is determined for the scaled
-                negative log-likelihood function with respect to the scaled
-                fit parameters.
-            err : dict
-                Dictionary containing  error of the parameters. Uses the same
-                keys as ``p``.
-            success : bool
-                Whether the fit terminated successfully.
-            status : int
-                Termination status of fit.
-            message : str
-                Description of the termination condition.
-            nfev, njev : int
-                Number of evaluations of the objective function and the
-                Jacobian.
-            nit : int
-                Number of iterations performed by `scipy.optimize.minimize`.
+        Value of optimized NLL cost function.
+    hess_inv : ndarray
+        Inverse Hessian matrix of optimized NLL cost function.
+    err_sigma_alpha, err_sigma_beta, err_sigma_tau : float
+        Estimated uncertainty in the noise model parameters. Set equal to 0.0
+        when the parameter is fixed.
+    err_mu : ndarray
+        Estimated uncertainty in ``mu``.
+    err_a : ndarray
+        Estimated uncertainty in ``a``.
+    err_eta : ndarray
+        Estimated uncertainty in ``eta``.
+    diagnostic : scipy.optimize.OptimizeResult
+        Instance of :class:`scipy.optimize.OptimizeResult` returned by
+        :func:`scipy.optimize.minimize`. Note that the attributes ``fun``,
+        ``jac``, and ``hess_inv`` represent functions over the internally
+        scaled parameters.
+
+    See Also
+    --------
+    NoiseModel : Noise model class.
+    noisefit : Estimate noise model parameters.
     """
     noise_model: NoiseModel
-    mu: ArrayLike
-    a: ArrayLike
-    eta: ArrayLike
+    mu: NDArray[np.float64]
+    a: NDArray[np.float64]
+    eta: NDArray[np.float64]
     fval: float
-    diagnostic: dict
+    hess_inv: NDArray[np.float64]
+    err_sigma_alpha: float
+    err_sigma_beta: float
+    err_sigma_tau: float
+    err_mu: NDArray[np.float64]
+    err_a: NDArray[np.float64]
+    err_eta: NDArray[np.float64]
+    diagnostic: OptimizeResult
 
 
 # noinspection PyShadowingNames
-def transfer_out(
+def transfer(
     tfun: Callable,
     x: ArrayLike,
-    dt: float | None,
     *,
-    fft_sign: bool = True,
+    dt: float | None = None,
+    numpy_sign_convention: bool = True,
     args: tuple = (),
-) -> np.ndarray:
+) -> NDArray[np.float64]:
     r"""
     Apply a transfer function to a waveform.
 
@@ -446,7 +495,7 @@ def transfer_out(
     tfun : callable
         Transfer function.
 
-            ``tfun(omega, *args) -> np.ndarray``
+            ``tfun(omega, *args) -> ndarray``
 
         where ``omega`` is an array of angular frequencies and ``args`` is a
         tuple of the fixed parameters needed to completely specify
@@ -454,18 +503,40 @@ def transfer_out(
         of ``dt``, such as radians/picosecond.
     x : array_like
         Data array.
-    dt : float or None
-        Sampling time, normally in picoseconds.
-    fft_sign : bool, optional
-        Complex exponential sign convention for harmonic time dependence.
+    dt : float or None, optional
+        Sampling time, normally in picoseconds. Default is None, which sets
+        the sampling time to ``thztools.options.sampling_time``. If both
+        ``dt`` and ``thztools.options.sampling_time`` are ``None``, the
+        sampling time is set to ``1.0``. In this case, the angular frequency
+        ``omega`` must be given in units of radians per sampling time, and any
+        parameters in ``args`` must be expressed with the sampling time as the
+        unit of time.
+    numpy_sign_convention : bool, optional
+        Adopt NumPy sign convention for harmonic time dependence, e.g, express
+        a harmonic function with frequency :math:`\omega` as
+        :math:`x(t) = a e^{i\omega t}`. Default is ``True``. Uses the
+        convention more common in physics, :math:`x(t) = a e^{-i\omega t}`,
+        when set to ``False``.
     args : tuple, optional
         Extra arguments passed to the transfer function.
 
     Returns
     -------
-    y : np.ndarray
+    y : ndarray
         Result of applying the transfer function to ``x``.
 
+    Warns
+    -----
+    UserWarning
+        If ``thztools.options.sampling_time`` and the ``dt`` parameter
+        are both not ``None`` and are set to different ``float`` values, the
+        function will set the sampling time to ``dt`` and raise a
+        :class:`UserWarning`.
+
+    See Also
+    --------
+    fit : Fit a transfer function to time-domain data.
+
     Notes
     -----
     The output waveform is computed by transforming :math:`x[n]` into the
@@ -476,66 +547,65 @@ def transfer_out(
 
     Examples
     --------
-
-    .. plot::
-       :include-source: True
-
-        Apply a transfer function that rescales and shifts the input.
-
-            .. math:: H(\omega) = a\exp(-i\omega\tau).
-
-        Note that this form assumes the :math:`e^{+i\omega t}` representation
-        of harmonic time dependence, which corresponds to the default setting
-        ``fft_sign=True``.
-
-        >>> import matplotlib.pyplot as plt
-        >>> import thztools as thz
-
-        >>> n, dt, t0 = 256, 0.05, 2.5
-        >>> x, t = thz.wave(n, dt, t0)
-
-        >>> def shiftscale(_w, _a, _tau):
-        >>>     return _a * np.exp(-1j * _w * _tau)
-        >>>
-        >>> y = thz.transfer_out(shiftscale, x, dt=dt, fft_sign=True,
-        ...                      args=(0.5, 1))
-
-        >>> _, ax = plt.subplots()
-        >>>
-        >>> ax.plot(t, x, label='x')
-        >>> ax.plot(t, y, label='y')
-        >>>
-        >>> ax.legend()
-        >>> ax.set_xlabel('t (ps)')
-        >>> ax.set_ylabel('Amplitude (arb. units)')
-        >>>
-        >>> plt.show()
-
-        If the transfer function is expressed using the :math:`e^{-i\omega t}`
-        representation that is more common in physics,
-
-            .. math:: H(\omega) = a\exp(i\omega\tau),
-
-        set ``fft_sign=False``.
-
-        >>> def shiftscale_phys(_w, _a, _tau):
-        >>>     return _a * np.exp(1j * _w * _tau)
-        >>>
-        >>> y_p = thz.transfer_out(shiftscale_phys, x, dt=dt, fft_sign=False,
-        ...                        args=(0.5, 1))
-
-        >>> _, ax = plt.subplots()
-        >>>
-        >>> ax.plot(t, x, label='x')
-        >>> ax.plot(t, y_p, label='y')
-        >>>
-        >>> ax.legend()
-        >>> ax.set_xlabel('t (ps)')
-        >>> ax.set_ylabel('Amplitude (arb. units)')
-        >>>
-        >>> plt.show()
+    Apply a transfer function that rescales the input by :math:`a` and shifts
+    it by :math:`\tau`.
+
+    .. math:: H(\omega) = a\exp(-i\omega\tau).
+
+    Note that this form assumes the :math:`e^{+i\omega t}` representation
+    of harmonic time dependence, which corresponds to the default setting
+    ``fft_sign=True``.
+
+    If the transfer function is expressed using the :math:`e^{-i\omega t}`
+    representation, more common in physics,
+
+    .. math:: H(\omega) = a\exp(i\omega\tau),
+
+    set ``fft_sign=False``.
+
+    >>> import thztools as thz
+    >>> from matplotlib import pyplot as plt
+    >>> n, dt = 256, 0.05
+    >>> thz.set_option("sampling_time", dt)
+    >>> t = thz.timebase(n)
+    >>> x = thz.wave(n)
+    >>> def shiftscale(_w, _a, _tau):
+    ...     return _a * np.exp(-1j * _w * _tau)
+    >>>
+    >>> y = thz.transfer(shiftscale, x,
+    ...                  numpy_sign_convention=True,
+    ...                  args=(0.5, 1))
+
+    >>> _, ax = plt.subplots()
+    >>>
+    >>> ax.plot(t, x, label='x')
+    >>> ax.plot(t, y, label='y')
+    >>>
+    >>> ax.legend()
+    >>> ax.set_xlabel('t (ps)')
+    >>> ax.set_ylabel('Amplitude (arb. units)')
+    >>>
+    >>> plt.show()
+
+    >>> def shiftscale_phys(_w, _a, _tau):
+    ...     return _a * np.exp(1j * _w * _tau)
+    >>>
+    >>> y_p = thz.transfer(shiftscale_phys, x,
+    ...                    numpy_sign_convention=False,
+    ...                    args=(0.5, 1))
+
+    >>> _, ax = plt.subplots()
+    >>>
+    >>> ax.plot(t, x, label='x');
+    >>> ax.plot(t, y_p, label='y');
+    >>>
+    >>> ax.legend();
+    >>> ax.set_xlabel('t (ps)');
+    >>> ax.set_ylabel('Amplitude (arb. units)');
+    >>>
+    >>> plt.show()
     """
-    x = np.asarray(x)
+    x = np.asarray(x, dtype=np.float64)
     if x.ndim != 1:
         msg = "x must be a one-dimensional array"
         raise ValueError(msg)
@@ -543,12 +613,11 @@ def transfer_out(
     if not isinstance(args, tuple):
         args = (args,)
 
-    dt = _validate_sampling_time(dt)
+    dt = _assign_sampling_time(dt)
     n = x.size
-    f_scaled = rfftfreq(n)
-    w = 2 * np.pi * f_scaled / dt
-    h = tfun(w, *args)
-    if not fft_sign:
+    w_scaled = 2 * np.pi * rfftfreq(n)
+    h = tfun(w_scaled / dt, *args)
+    if not numpy_sign_convention:
         h = np.conj(h)
 
     y = np.fft.irfft(np.fft.rfft(x) * h, n=n)
@@ -557,16 +626,78 @@ def transfer_out(
 
 
 # noinspection PyShadowingNames
+def timebase(
+    n: int, *, dt: float | None = None, t_init: float = 0.0
+) -> NDArray[np.float64]:
+    r"""
+    Timebase for time-domain waveforms.
+
+    Parameters
+    ----------
+    n : int
+        Number of samples.
+    dt : float or None, optional
+        Sampling time, normally in picoseconds. Default is None, which sets
+        the sampling time to ``thztools.options.sampling_time``. If both
+        ``dt`` and ``thztools.options.sampling_time`` are ``None``, the
+        sampling time is set to 1.0.
+    t_init : float, optional
+        Value of the initial time point. Default is ``0.0``.
+
+    Returns
+    -------
+    t : ndarray
+        Array of time samples.
+
+    Warns
+    -----
+    UserWarning
+        If ``thztools.options.sampling_time`` and the ``dt`` parameter
+        are both not ``None`` and are set to different ``float`` values, the
+        function will set the sampling time to ``dt`` and raise a
+        :class:`UserWarning`.
+
+    Notes
+    -----
+    This is a helper function that computes
+
+        ``t = t_init + dt * numpy.arange(n)``.
+
+    Examples
+    --------
+    The following example computes the timebase with different methods for
+    assigning the sampling time.
+
+    >>> import thztools as thz
+    >>> n, dt, t_init = 256, 0.05, 2.5
+    >>> t_1 = thz.timebase(n)
+    >>> t_2 = thz.timebase(n, dt=dt)
+    >>> thz.set_option("sampling_time", dt)
+    >>> t_3 = thz.timebase(n)
+    >>> t_4 = thz.timebase(n, t_init=t_init)
+    >>> print(t_1[:3])
+    [0. 1. 2.]
+    >>> print(t_2[:3])
+    [0.   0.05 0.1 ]
+    >>> print(t_3[:3])
+    [0.   0.05 0.1 ]
+    >>> print(t_4[:3])
+    [2.5  2.55 2.6 ]
+    """
+    dt = _assign_sampling_time(dt)
+    return t_init + dt * np.arange(n)
+
+
 def wave(
     n: int,
-    dt: float | None,
-    t0: float,
     *,
+    dt: float | None = None,
+    t0: float | None = None,
     a: float = 1.0,
-    taur: float = 0.3,
-    tauc: float = 0.1,
-    fwhm: float = 0.05,
-) -> tuple[np.ndarray, np.ndarray]:
+    taur: float | None = None,
+    tauc: float | None = None,
+    fwhm: float | None = None,
+) -> NDArray[np.float64]:
     r"""
     Simulate a terahertz waveform.
 
@@ -575,25 +706,33 @@ def wave(
 
     n : int
         Number of samples.
-    dt : float or None
-        Sampling time, normally in picoseconds.
-    t0 : float
-        Pulse center, normally in picoseconds.
+    dt : float or None, optional
+        Sampling time, normally in picoseconds. Default is None, which sets
+        the sampling time to ``thztools.options.sampling_time``. If both
+        ``dt`` and ``thztools.options.sampling_time`` are ``None``, the
+        sampling time is set to 1.0.
+    t0 : float or None, optional
+        Pulse location, normally in picoseconds. Default is ``0.3 * n * dt``.
     a : float, optional
         Peak amplitude. The default is one.
-    taur, tauc, fwhm : float, optional
+    taur, tauc, fwhm : float or None, optional
         Current pulse rise time, current pulse decay time, and laser pulse
-        FWHM, respectively. The defaults are 0.3 ps, 0.1 ps, and 0.05 ps,
-        respectively, and assume that ``ts`` and ``t0`` are also given in
-        picoseconds.
+        FWHM, respectively. The defaults are ``6.0 * dt``, ``2.0 * dt``, and
+        ``1.0 * dt``, respectively.
 
     Returns
     -------
 
     x : ndarray
         Signal array.
-    t : ndarray
-        Array of time samples.
+
+    Warns
+    -----
+    UserWarning
+        If ``thztools.options.sampling_time`` and the ``dt`` parameter
+        are both not ``None`` and are set to different ``float`` values, the
+        function will set the sampling time to ``dt`` and raise a
+        :class:`UserWarning`.
 
     Notes
     -----
@@ -620,21 +759,31 @@ def wave(
     The following example shows the simulated signal :math:`\mu(t)` normalized
     to its peak magnitude, :math:`\mu_0`.
 
-    .. plot::
-       :include-source: True
+    >>> import thztools as thz
+    >>> from matplotlib import pyplot as plt
+    >>> n, dt = 256, 0.05
+    >>> thz.set_option("sampling_time", dt)
+    >>> t = thz.timebase(n)
+    >>> mu = thz.wave(n)
+
+    >>> _, ax = plt.subplots(layout="constrained")
+    >>> ax.plot(t, mu);
+    >>> ax.set_xlabel("t (ps)");
+    >>> ax.set_ylabel(r"$\mu/\mu_0$");
+    >>> plt.show()
+    """
+    dt = _assign_sampling_time(dt)
+    if t0 is None:
+        t0 = 0.3 * n * dt
 
-        >>> import matplotlib.pyplot as plt
-        >>> import thztools as thz
-        >>> n, dt, t0 = 256, 0.05, 2.5
-        >>> mu, t = thz.wave(n, dt, t0)
+    if taur is None:
+        taur = 6.0 * dt
 
-        >>> _, ax = plt.subplots(layout="constrained")
-        >>> ax.plot(t, mu)
-        >>> ax.set_xlabel("t (ps)")
-        >>> ax.set_ylabel(r"$\mu/\mu_0$")
-        >>> plt.show()
-    """
-    dt = _validate_sampling_time(dt)
+    if tauc is None:
+        tauc = 2.0 * dt
+
+    if fwhm is None:
+        fwhm = dt
 
     taul = fwhm / np.sqrt(2 * np.log(2))
 
@@ -645,23 +794,21 @@ def wave(
     r = 1 / (1 / taur - 1j * w) - 1 / (1 / taur + 1 / tauc - 1j * w)
     s = -1j * w * (ell * r) ** 2 * np.exp(1j * w * t0)
 
-    t = dt * np.arange(n)
-
     x = irfft(np.conj(s), n=n)
     x = a * x / np.max(x)
 
-    return x, t
+    return x
 
 
 # noinspection PyShadowingNames
 def scaleshift(
     x: ArrayLike,
     *,
+    dt: float | None = None,
     a: ArrayLike | None = None,
     eta: ArrayLike | None = None,
-    dt: float | None = 1.0,
     axis: int = -1,
-) -> np.ndarray:
+) -> NDArray[np.float64]:
     r"""
     Rescale and shift waveforms.
 
@@ -669,54 +816,63 @@ def scaleshift(
     ----------
     x : array_like
         Data array.
+    dt : float or None, optional
+        Sampling time, normally in picoseconds. Default is None, which sets
+        the sampling time to ``thztools.options.sampling_time``. If both
+        ``dt`` and ``thztools.options.sampling_time`` are ``None``, the
+        sampling time is set to ``1.0``. In this case, ``eta`` must be given in
+        units of the sampling time.
     a : array_like, optional
         Scale array.
     eta : array_like, optional
         Shift array.
-    dt : float or None, optional
-        Sampling time. Default is 1.0.
     axis : int, optional
         Axis over which to apply the correction. If not given, applies over the
         last axis in ``x``.
 
     Returns
     -------
-    x_adjusted : ndarray
+    x_adj : ndarray
         Adjusted data array.
 
+    Warns
+    -----
+    UserWarning
+        If ``thztools.options.sampling_time`` and the ``dt`` parameter
+        are both not ``None`` and are set to different ``float`` values, the
+        function will set the sampling time to ``dt`` and raise a
+        :class:`UserWarning`.
+
     Examples
     --------
     The following example makes an array with 4 identical copies of the
-    signal ``mu`` returned by ``thztools.thzgen``. It then uses
-    ``thztools.scaleshift`` to rescale each copy by
-    ``a = [1.0, 0.5, 0.25, 0.125]`` and shift it by
+    signal ``mu`` returned by :func:`wave`. It then uses :func:`scaleshift` to
+    rescale each copy by ``a = [1.0, 0.5, 0.25, 0.125]`` and shift it by
     ``eta = [0.0, 1.0, 2.0, 3.0]``.
 
-    .. plot::
-       :include-source: True
-
-        >>> import matplotlib.pyplot as plt
-        >>> import thztools as thz
-        >>> n, dt, t0 = 256, 0.05, 2.5
-        >>> mu, t = thz.wave(n, dt, t0)
-        >>> m = 4
-        >>> x = np.repeat(np.atleast_2d(mu), m, axis=0)
-        >>> a = 0.5**np.arange(m)
-        >>> eta = np.arange(m)
-        >>> x_adj = thz.scaleshift(x, a=a, eta=eta, dt=dt)
-
-        >>> _, ax = plt.subplots(layout="constrained")
-        >>> ax.plot(t, x_adj.T, label=[f"{k=}" for k in range(4)])
-        >>> ax.legend()
-        >>> ax.set_xlabel("t (ps)")
-        >>> ax.set_ylabel(r"$x_{\mathrm{adj}, k}$")
-        >>> plt.show()
+    >>> import thztools as thz
+    >>> from matplotlib import pyplot as plt
+    >>> n, dt = 256, 0.05
+    >>> thz.set_option("sampling_time", dt)
+    >>> t = thz.timebase(n)
+    >>> mu = thz.wave(n)
+    >>> m = 4
+    >>> x = np.repeat(np.atleast_2d(mu), m, axis=0)
+    >>> a = 0.5**np.arange(m)
+    >>> eta = np.arange(m)
+    >>> x_adj = thz.scaleshift(x, a=a, eta=eta, dt=dt)
+
+    >>> plt.plot(t, x_adj.T, label=[f"{k=}" for k in range(4)]);
+    >>> plt.legend();
+    >>> plt.xlabel("t (ps)");
+    >>> plt.ylabel(r"$x_{\mathrm{adj}, k}$");
+    >>> plt.show()
     """
-    x = np.asarray(x)
+    x = np.asarray(x, dtype=np.float64)
     if x.size == 0:
         return np.empty(x.shape)
 
-    dt = _validate_sampling_time(dt)
+    dt = _assign_sampling_time(dt)
 
     axis = int(axis)
     if x.ndim > 1:
@@ -729,7 +885,7 @@ def scaleshift(
     if a is None:
         a = np.ones(m)
     else:
-        a = np.asarray(a)
+        a = np.asarray(a, dtype=np.float64)
         if a.shape != m:
             msg = (
                 f"Scale correction with shape {a.shape} can not be applied "
@@ -740,7 +896,7 @@ def scaleshift(
     if eta is None:
         eta = np.zeros(m)
     else:
-        eta = np.asarray(eta)
+        eta = np.asarray(eta, dtype=np.float64)
         if eta.shape != m:
             msg = (
                 f"Shift correction with shape {a.shape} can not be applied "
@@ -772,7 +928,7 @@ def _costfuntls(
     sigma_x: ArrayLike,
     sigma_y: ArrayLike,
     dt: float | None = 1.0,
-) -> np.ndarray:
+) -> NDArray[np.float64]:
     r"""Computes the residual vector for the total least squares cost function.
 
     Parameters
@@ -780,45 +936,51 @@ def _costfuntls(
     fun : callable
         Transfer function.
 
-            ``fun(p, w, *args, **kwargs) -> np.ndarray``
+            ``fun(p, w, *args, **kwargs) -> ndarray``
 
         Assumes the :math:`+i\omega t` convention for harmonic time dependence.
-
     theta : array_like
         Input parameters for the function.
-
     mu : array_like
         Estimated input signal.
-
     x : array_like
         Measured input signal.
-
     y : array_like
         Measured output signal.
-
     sigma_x : array_like
         Noise vector of the input signal.
-
     sigma_y : array_like
         Noise vector of the output signal.
-
     dt : float or None, optional
-        Sampling time. Default set to 1.0.
+        Sampling time, normally in picoseconds. Default is None, which sets
+        the sampling time to ``thztools.options.sampling_time``. If both
+        ``dt`` and ``thztools.options.sampling_time`` are ``None``, the
+        sampling time is set to ``1.0``. In this case, the angular frequency
+        ``omega`` must be given in units of radians per sampling time, and any
+        parameters in ``args`` must be expressed with the sampling time as the
+        unit of time.
 
     Returns
     -------
-    res : array_like
-
+    res : ndarray
+        Residual array.
 
+    Warns
+    -----
+    UserWarning
+        If ``thztools.options.sampling_time`` and the ``dt`` parameter
+        are both not ``None`` and are set to different ``float`` values, the
+        function will set the sampling time to ``dt`` and raise a
+        :class:`UserWarning`.
     """
-    theta = np.asarray(theta)
-    mu = np.asarray(mu)
-    x = np.asarray(x)
-    y = np.asarray(y)
-    sigma_x = np.asarray(sigma_x)
-    sigma_y = np.asarray(sigma_y)
+    theta = np.asarray(theta, dtype=np.float64)
+    mu = np.asarray(mu, dtype=np.float64)
+    x = np.asarray(x, dtype=np.float64)
+    y = np.asarray(y, dtype=np.float64)
+    sigma_x = np.asarray(sigma_x, dtype=np.float64)
+    sigma_y = np.asarray(sigma_y, dtype=np.float64)
 
-    dt = _validate_sampling_time(dt)
+    dt = _assign_sampling_time(dt)
 
     n = x.shape[-1]
     delta_norm = (x - mu) / sigma_x
@@ -832,98 +994,103 @@ def _costfuntls(
     return res
 
 
-def _tdnll_scaled(
-    x: ArrayLike,
-    logv: ArrayLike,
-    delta: ArrayLike,
-    alpha: ArrayLike,
-    eta_on_dt: ArrayLike,
+def _costfun_noisefit(
+    x: NDArray[np.float64],
+    logv_alpha: float,
+    logv_beta: float,
+    logv_tau: float,
+    delta_mu: NDArray[np.float64],
+    delta_a: NDArray[np.float64],
+    eta: NDArray[np.float64],
     *,
-    fix_logv: bool,
-    fix_delta: bool,
-    fix_alpha: bool,
+    fix_logv_alpha: bool,
+    fix_logv_beta: bool,
+    fix_logv_tau: bool,
+    fix_delta_mu: bool,
+    fix_delta_a: bool,
     fix_eta: bool,
-    scale_logv: ArrayLike,
-    scale_delta: ArrayLike,
-    scale_alpha: ArrayLike,
-    scale_eta: ArrayLike,
-    scale_v: float,
-) -> tuple[float, np.ndarray]:
+    scale_sigma_alpha: float,
+    scale_sigma_beta: float,
+    scale_sigma_tau: float,
+    scale_delta_mu: NDArray[np.float64],
+    scale_delta_a: NDArray[np.float64],
+    scale_eta: NDArray[np.float64],
+) -> tuple[float, NDArray[np.float64]]:
     r"""
-    Compute the scaled negative log-likelihood for the time-domain noise model.
+    Compute the scaled cost function for the time-domain noise model.
 
-    Computes the scaled negative log-likelihood function for obtaining the
+    Computes the scaled maximum-likelihood cost function for obtaining the
     data matrix ``x`` given ``logv``, ``delta``, ``alpha``, ``eta``,
-    ``scale_logv``, ``scale_delta``, ``scale_alpha``, ``scale_eta``, and
-    ``scale_v``.
+    ``scale_sigma_alpha``, ``scale_sigma_beta``, ``scale_sigma_tau``,
+    ``scale_delta_mu``, ``scale_delta_a``, and ``scale_eta``.
 
     Parameters
     ----------
-    x : array_like
+    x : ndarray
         Data matrix with shape (m, n), row-oriented.
-    logv : array_like
-        Array of three noise parameters.
-    delta : array_like
+    logv_alpha, logv_beta, logv_tau : float
+        Logarithm of the associated noise variance parameter, with sigma_tau
+        given in units of sampling time.
+    delta_mu : ndarray
         Signal deviation vector with shape (n,).
-    alpha: array_like
+    delta_a: ndarray
         Amplitude deviation vector with shape (m - 1,).
-    eta_on_dt : array_like
-        Normalized delay deviation vector with shape (m - 1,), equal to
-        ``eta``/``dt``.
-    fix_logv : bool
-        Exclude noise parameters from gradiate calculation when ``True``.
-    fix_delta : bool
+    eta : ndarray
+        Delay deviation vector with shape (m - 1,), given in units of sampling
+        time.
+    fix_logv_alpha, fix_logv_beta, fix_logv_tau : bool
+        Exclude noise parameter from gradiate calculation when ``True``.
+    fix_delta_mu : bool
         Exclude signal deviation vector from gradiate calculation when
         ``True``.
-    fix_alpha : bool
-        Exclude amplitude deviation vector from gradiate calculation when
-        ``True``.
+    fix_delta_a : bool
+        Exclude signal amplitude deviation vector from gradiate calculation
+        when ``True``.
     fix_eta : bool
-        Exclude delay deviation vector from gradiate calculation when ``True``.
-    scale_logv : array_like
-        Array of three scale parameters for ``logv``.
-    scale_delta : array_like
+        Exclude signal delay deviation vector from gradiate calculation when
+        ``True``.
+    scale_sigma_alpha : float
+        Scale parameter for ``sigma_alpha``.
+    scale_sigma_beta : float
+        Scale parameter for ``sigma_beta``.
+    scale_sigma_tau : float
+        Scale parameter for ``sigma_tau``.
+    scale_delta_mu : ndarray
         Array of scale parameters for ``delta`` with shape (n,).
-    scale_alpha : array_like
+    scale_delta_a : ndarray
         Array of scale parameters for ``alpha`` with shape (m - 1,).
-    scale_eta : array_like
+    scale_eta : ndarray
         Array of scale parameters for ``eta`` with shape (m - 1,).
-    scale_v : float
-        Scale parameter for overall variance.
 
     Returns
     -------
-    nll_scaled : float
+    costfun_scaled : float
         Scaled negative log-likelihood.
-    gradnll_scaled : array_like
+    gradnll_scaled : ndarray
         Gradient of the negative scaled log-likelihood function with respect to
         free parameters.
     """
-    x = np.asarray(x)
-    logv = np.asarray(logv)
-    delta = np.asarray(delta)
-    alpha = np.asarray(alpha)
-    eta_on_dt = np.asarray(eta_on_dt)
-
-    scale_logv = np.asarray(scale_logv)
-    scale_delta = np.asarray(scale_delta)
-    scale_alpha = np.asarray(scale_alpha)
-    scale_eta = np.asarray(scale_eta)
-
     m, n = x.shape
 
-    # Compute scaled variance, mu, a, and eta
-    v_scaled = np.exp(logv * scale_logv)
-    mu = x[0, :] - delta * scale_delta
-    a = np.insert(1.0 + alpha * scale_alpha, 0, 1.0)
-    eta_on_dt = np.insert(eta_on_dt * scale_eta, 0, 0.0)
+    logv = np.asarray([logv_alpha, logv_beta, logv_tau], dtype=np.float64)
+
+    # Compute scaled noise model parameters, mu, a, and eta
+    scale_sigma = np.array(
+        [scale_sigma_alpha, scale_sigma_beta, scale_sigma_tau],
+        dtype=np.float64,
+    )
+    scale_v = scale_sigma**2
+    v_scaled = np.exp(logv) * scale_v
+    mu = x[0, :] - delta_mu * scale_delta_mu
+    a = np.insert(1.0 + delta_a * scale_delta_a, 0, 1.0)
+    eta = np.insert(eta * scale_eta, 0, 0.0)
 
     # Compute frequency vector and Fourier coefficients of mu
     f_scaled = rfftfreq(n)
     w_scaled = 2 * np.pi * f_scaled
     mu_f = rfft(mu)
 
-    exp_iweta = np.exp(1j * np.outer(eta_on_dt, w_scaled))
+    exp_iweta = np.exp(1j * np.outer(eta, w_scaled))
     zeta_f = ((np.conj(exp_iweta) * mu_f).T * a).T
     zeta = irfft(zeta_f, n=n)
 
@@ -939,51 +1106,57 @@ def _tdnll_scaled(
     valpha = v_scaled[0]
     vbeta = v_scaled[1] * zeta**2
     vtau = v_scaled[2] * dzeta**2
-    vtot_scaled = valpha + vbeta + vtau
+    vtot = valpha + vbeta + vtau
 
-    resnormsq_scaled = ressq / vtot_scaled
-    nll_scaled = (
-        scale_v * m * n * (np.log(2 * np.pi) + np.log(scale_v)) / 2
-        + scale_v * np.sum(np.log(vtot_scaled)) / 2
-        + np.sum(resnormsq_scaled) / 2
-    )
+    resnormsq_scaled = ressq / vtot
+    costfun_scaled = np.sum(np.log(vtot)) + np.sum(resnormsq_scaled)
 
     # Compute gradient
-    gradnll_scaled = np.array([])
-    if not (fix_logv & fix_delta & fix_alpha & fix_eta):
-        reswt = res / vtot_scaled
-        dvar = (scale_v * vtot_scaled - ressq) / vtot_scaled**2
-        if not fix_logv:
-            # Gradient wrt logv
+    gradnll_scaled = np.array([], dtype=np.float64)
+    if not (
+        fix_logv_alpha
+        & fix_logv_beta
+        & fix_logv_tau
+        & fix_delta_mu
+        & fix_delta_a
+        & fix_eta
+    ):
+        reswt = res / vtot
+        dvar = (vtot - ressq) / vtot**2
+        # Gradient wrt logv
+        if not fix_logv_alpha:
             gradnll_scaled = np.append(
                 gradnll_scaled,
-                0.5 * np.sum(dvar) * v_scaled[0] * scale_logv[0],
+                np.sum(dvar) * v_scaled[0],
             )
+        if not fix_logv_beta:
             gradnll_scaled = np.append(
                 gradnll_scaled,
-                0.5 * np.sum(zeta**2 * dvar) * v_scaled[1] * scale_logv[1],
+                np.sum(zeta**2 * dvar) * v_scaled[1],
             )
+        if not fix_logv_tau:
             gradnll_scaled = np.append(
                 gradnll_scaled,
-                0.5 * np.sum(dzeta**2 * dvar) * v_scaled[2] * scale_logv[2],
+                np.sum(dzeta**2 * dvar) * v_scaled[2],
             )
-        if not fix_delta:
+        if not fix_delta_mu:
             # Gradient wrt delta
             p = rfft(v_scaled[1] * dvar * zeta - reswt) - 1j * v_scaled[
                 2
             ] * w_scaled * rfft(dvar * dzeta)
             gradnll_scaled = np.append(
                 gradnll_scaled,
-                -np.sum((irfft(exp_iweta * p, n=n).T * a).T, axis=0)
-                * scale_delta,
+                -2
+                * np.sum((irfft(exp_iweta * p, n=n).T * a).T, axis=0)
+                * scale_delta_mu,
             )
-        if not fix_alpha:
+        if not fix_delta_a:
             # Gradient wrt alpha
-            term = (vtot_scaled - valpha) * dvar - reswt * zeta
+            term = (vtot - valpha) * dvar - reswt * zeta
             dnllda = np.sum(term, axis=1).T / a
             # Exclude first term, which is held fixed
             gradnll_scaled = np.append(
-                gradnll_scaled, dnllda[1:] * scale_alpha
+                gradnll_scaled, 2 * dnllda[1:] * scale_delta_a
             )
         if not fix_eta:
             # Gradient wrt eta
@@ -996,272 +1169,635 @@ def _tdnll_scaled(
             )
             # Exclude first term, which is held fixed
             gradnll_scaled = np.append(
-                gradnll_scaled, dnlldeta[1:] * scale_eta
+                gradnll_scaled, 2 * dnlldeta[1:] * scale_eta
             )
 
-    return nll_scaled, gradnll_scaled
+    return costfun_scaled / (m * n), gradnll_scaled / (m * n)
 
 
-def tdnoisefit(
+def noisefit(
     x: ArrayLike,
     *,
-    dt: float | None = 1.0,
-    v0: ArrayLike | None = None,
+    dt: float | None = None,
+    sigma_alpha0: float | None = None,
+    sigma_beta0: float | None = None,
+    sigma_tau0: float | None = None,
     mu0: ArrayLike | None = None,
     a0: ArrayLike | None = None,
     eta0: ArrayLike | None = None,
-    fix_v: bool = False,
+    fix_sigma_alpha: bool = False,
+    fix_sigma_beta: bool = False,
+    fix_sigma_tau: bool = False,
     fix_mu: bool = False,
-    fix_a: bool = True,
-    fix_eta: bool = True,
+    fix_a: bool = False,
+    fix_eta: bool = False,
+    scale_sigma_alpha: float | None = None,
+    scale_sigma_beta: float | None = None,
+    scale_sigma_tau: float | None = None,
+    scale_delta_mu: ArrayLike | None = None,
+    scale_delta_a: ArrayLike | None = None,
+    scale_eta: ArrayLike | None = None,
 ) -> NoiseResult:
     r"""
-    Estimate noise model parameters.
+    Estimate noise model from a set of nominally identical waveforms.
 
-    Computes the noise parameters sigma and the underlying signal vector ``mu``
-    for the data matrix ``x``, where the columns of ``x`` are each noisy
-    measurements of ``mu``.
+    The data array ``x`` should include `m` nominally identical waveform
+    measurements, where each waveform comprises `n` samples that are spaced
+    equally in time. The ``noise_model`` attribute of the return object
+    is an instance of the :class:`NoiseModel` class that represents the
+    estimated noise parameters. See `Notes` for details.
 
     Parameters
     ----------
-    x : ndarray, shape(n, m)
-        Data array with `m` waveforms, each composed of `n` points.
+    x : array_like with shape (n, m)
+        Data array composed of ``m`` waveforms, each of which is sampled at
+        ``n`` points.
     dt : float or None, optional
-        Sampling time. Default is 1.0.
-    v0 : ndarray, shape (3,), optional
-        Initial guess, noise model parameters with size (3,), expressed as
-        variance amplitudes.
-    mu0 : ndarray, shape(n,), optional
-        Initial guess, signal vector with size (n,). Default is first column of
-        ``x``.
-    a0 : ndarray, shape(m,), optional
-        Initial guess, amplitude vector with size (m,). Default is one for all
-        entries.
-    eta0 : ndarray, shape(m,), optional
-        Initial guess, delay vector with size (m,). Default is zero for all
-        entries.
-    fix_v : bool, optional
-        Fix noise variance parameters. Default is False.
+        Sampling time, normally in picoseconds. Default is None, which sets
+        the sampling time to ``thztools.options.sampling_time``. If both
+        ``dt`` and ``thztools.options.sampling_time`` are ``None``, the
+        sampling time is set to ``1.0``. In this case, the unit of time is the
+        sampling time for the noise model parameter ``sigma_tau``, the delay
+        parameter array ``eta``, and the initial guesses for both quantities.
+    sigma_alpha0, sigma_beta0, sigma_tau0 : float, optional
+        Initial values for noise parameters. When set to
+        ``None``, the default, use a linear least-squares fit of the noise
+        model to ``np.var(x, 1, ddof=1)``.
+    mu0 : array_like with shape(n,), optional
+        Initial guess, signal vector with shape (n,). Default is first column
+        of ``x``.
+    a0 : array_like with shape(m,), optional
+        Initial guess, signal amplitude drift vector with shape (m,). Default
+        is ``np.ones(m)``.
+    eta0 : array_like with shape(m,), optional
+        Initial guess, signal delay drift vector with shape (m,). Default is
+        ``np.zeros(m)``.
+    fix_sigma_alpha, fix_sigma_beta, fix_sigma_tau : bool, optional
+        Fix the associated noise parameter. Default is False.
     fix_mu : bool, optional
         Fix signal vector. Default is False.
     fix_a : bool, optional
-        Fix amplitude vector. Default is True.
+        Fix signal amplitude drift vector. Default is False.
     fix_eta : bool, optional
-        Fix delay vector. Default is True.
+        Fix signal delay drift vector. Default is False.
 
     Returns
     -------
     res : NoiseResult
-        Fit result, represented as a ``NoiseResult`` object.
+        Fit result represented as a ``NoiseResult`` object. Important
+        attributes are: ``noise_model``, an instance of :class:`NoiseModel`
+        with the estimated noise parameters; ``mu``, the estimated signal
+        vector; ``a``, the estimated signal amplitude drift vector; ``eta``,
+        the estimated signal delay drift vector; ``fval``, the value of the
+        optimized NLL cost function; and ``diagnostic``, an instance of
+        :class:`scipy.optimize.OptimizeResult` returned by
+        :func:`scipy.optimize.minimize`. Note that the fit parameters are
+        scaled internally to improve convergence (see `Other Parameters`
+        below), which affects the attributes ``fun``, ``jac``, and ``hess_inv``
+        of ``diagnostic``. See :class:`NoiseResult` for more details and for a
+        description of other attributes.
+
+    Other Parameters
+    ----------------
+    scale_sigma_alpha, scale_sigma_beta, scale_sigma_tau : float, optional
+        Scale for varying noise parameters. Default is ``1.0``.
+    scale_delta_mu : array_like with shape(n,), optional
+        Scale for varying signal vector. When set to ``None``, the default,
+        use ``NoiseModel(sigma_alpha0, sigma_beta0,
+        sigma_tau0).amplitude(mu0)``.
+    scale_delta_a : array_like with shape(n,), optional
+            Scale for varying signal amplitude drift vector. Default is
+            ``1.0e-2`` for all entries.
+    scale_eta : array_like with shape(n,), optional
+            Scale for varying signal delay drift vector. Default is
+            ``1.0e-3 / dt`` for all entries.
+
+    Raises
+    ------
+    ValueError
+        If all parameters are held fixed or if the input arrays have
+        incompatible shapes.
+
+    Warns
+    -----
+    UserWarning
+        If ``thztools.options.sampling_time`` and the ``dt`` parameter
+        are both not ``None`` and are set to different ``float`` values, the
+        function will set the sampling time to ``dt`` and raise a
+        :class:`UserWarning`.
+
+    See Also
+    --------
+    NoiseModel : Noise model class.
+
+    Notes
+    -----
+    Given an :math:`N\times M` data array :math:`\mathbf{X}`, the function uses
+    the BFGS method of :func:`scipy.optimize.minimize` to minimize the
+    maximum-likelihood cost function [1]_
+
+    .. math:: \begin{split}\
+        Q_\text{ML}\
+        (\sigma_\alpha,\sigma_\beta,\sigma_\tau,\boldsymbol{\mu},\
+        \mathbf{A},\boldsymbol{\eta};\mathbf{X})\
+        = \frac{1}{MN}\sum_{k=0}^{N-1}\sum_{l=0}^{M-1}& \
+        \left[\ln\sigma_{kl}^2 + \frac{(X_{kl} \
+        - Z_{kl})^2}{\sigma_{kl}^2}\right]\
+        \end{split}
+
+    with respect to the unknown model parameters :math:`\sigma_\alpha`,
+    :math:`\sigma_\beta`, :math:`\sigma_\tau`, :math:`\boldsymbol{\mu}`,
+    :math:`\mathbf{A}` and :math:`\boldsymbol{\eta}`. The model assumes that
+    the elements of :math:`\mathbf{X}` are normally distributed around the
+    corresponding elements of :math:`\mathbf{Z}` with variances given by
+    :math:`\boldsymbol{\sigma}^2`, which are in turn given by
+
+    .. math:: Z_{kl} = A_l\mu(t_k - \eta_l)
+
+    and
+
+    .. math:: \sigma_{kl}^2 = \sigma_\alpha^2 + \sigma_\beta^2 Z_{kl}^2 \
+        + \sigma_\tau^2(\mathbf{D}\mathbf{Z})_{kl}^2,
+
+    where :math:`\mathbf{D}` is the time-domain derivative operator. The
+    function :math:`\mu(t)` represents an ideal primary waveform, which drifts
+    in amplitude and temporal location between each of the :math:`M` waveform
+    measurements. Each waveform comprises :math:`N` samples at the nominal
+    times :math:`t_k`, and the drift in the amplitude and the temporal location
+    of the :math:`l` waveform are given by :math:`A_l` and :math:`\eta_l`,
+    respectively, with initial values fixed at :math:`A_0 = 1.0` and
+    :math:`\eta_0 = 0.0` by definition.
+
+    References
+    ----------
+    .. [1] Laleh Mohtashemi, Paul Westlund, Derek G. Sahota, Graham B. Lea,
+        Ian Bushfield, Payam Mousavi, and J. Steven Dodge, "Maximum-
+        likelihood parameter estimation in terahertz time-domain
+        spectroscopy," Opt. Express **29**, 4912-4926 (2021),
+        `<https://doi.org/10.1364/OE.417724>`_.
+
+    Examples
+    --------
+    In basic usage, it should be possible to call ``noisefit`` with just
+    the data array. In the code below, we simulate ``m = 50`` noisy waveforms
+    with ``n = 256`` time points each, then verify that the fit results are
+    consistent with the true noise parameters.
+
+    >>> import thztools as thz
+    >>> from matplotlib import pyplot as plt
+    >>> rng = np.random.default_rng(0)
+    >>> n, m, dt = 256, 50, 0.05
+    >>> thz.set_option("sampling_time", dt)
+    >>> t = thz.timebase(n)
+    >>> mu = thz.wave(n)
+    >>> alpha, beta, tau = 1e-5, 1e-2, 1e-3
+    >>> noise_model = thz.NoiseModel(sigma_alpha=alpha, sigma_beta=beta,
+    ...  sigma_tau=tau)
+    >>> a = 1.0 + 1e-2 * np.concatenate(([0.0],
+    ...                                 rng.standard_normal(m - 1)))
+    >>> eta = 1e-3 * np.concatenate(([0.0], rng.standard_normal(m - 1)))
+    >>> z = thz.scaleshift(np.repeat(np.atleast_2d(mu), m, axis=0),
+    ...                    a=a, eta=eta).T  # Orient the array columnwise
+    >>> x = z + noise_model.noise_sim(z, axis=1, seed=1234)
+    >>> noise_res = thz.noisefit(x)
+    >>> noise_res.noise_model  #doctest: +NORMALIZE_WHITESPACE
+    NoiseModel(sigma_alpha=1.000...e-05, sigma_beta=0.009609...,
+    sigma_tau=0.0009243..., dt=0.05)
+
+    >>> plt.plot(t, np.std(thz.scaleshift(x, a=1 / noise_res.a,
+    ... eta=-noise_res.eta, axis=0), axis=1), "-",
+    ... label="Data")
+    >>> plt.plot(t, noise_res.noise_model.noise_amp(noise_res.mu)
+    ...  * np.sqrt(m / (m + 1)), "--", label="Fit")
+    >>> plt.legend()
+    >>> plt.xlabel("t (ps)")
+    >>> plt.ylabel(r"$\sigma(t)$")
+    >>> plt.show()
     """
-    if fix_v and fix_mu and fix_a and fix_eta:
-        msg = "All variables are fixed"
-        raise ValueError(msg)
-    # Parse and validate function inputs
-    x = np.asarray(x)
+    x = np.asarray(x, dtype=np.float64)
+    dt = _assign_sampling_time(dt)
+
+    parsed = _parse_noisefit_input(
+        x,
+        dt=dt,
+        sigma_alpha0=sigma_alpha0,
+        sigma_beta0=sigma_beta0,
+        sigma_tau0=sigma_tau0,
+        mu0=mu0,
+        a0=a0,
+        eta0=eta0,
+        fix_sigma_alpha=fix_sigma_alpha,
+        fix_sigma_beta=fix_sigma_beta,
+        fix_sigma_tau=fix_sigma_tau,
+        fix_mu=fix_mu,
+        fix_a=fix_a,
+        fix_eta=fix_eta,
+        scale_sigma_alpha=scale_sigma_alpha,
+        scale_sigma_beta=scale_sigma_beta,
+        scale_sigma_tau=scale_sigma_tau,
+        scale_delta_mu=scale_delta_mu,
+        scale_delta_a=scale_delta_a,
+        scale_eta=scale_eta,
+    )
+
+    objective, x0, input_parsed = parsed
+
+    # Minimize cost function with respect to free parameters
+    out = minimize(
+        objective,
+        x0,
+        method="BFGS",
+        jac=True,
+    )
+    fit_result = _parse_noisefit_output(out, x, dt=dt, **input_parsed)
+    return fit_result
+
+
+def _parse_noisefit_input(
+    x: NDArray[np.float64],
+    dt: float,
+    *,
+    sigma_alpha0: float | None,
+    sigma_beta0: float | None,
+    sigma_tau0: float | None,
+    mu0: ArrayLike | None,
+    a0: ArrayLike | None,
+    eta0: ArrayLike | None,
+    fix_sigma_alpha: bool,
+    fix_sigma_beta: bool,
+    fix_sigma_tau: bool,
+    fix_mu: bool,
+    fix_a: bool,
+    fix_eta: bool,
+    scale_sigma_alpha: float | None,
+    scale_sigma_beta: float | None,
+    scale_sigma_tau: float | None,
+    scale_delta_mu: ArrayLike | None,
+    scale_delta_a: ArrayLike | None,
+    scale_eta: ArrayLike | None,
+) -> tuple[Callable, NDArray[np.float64], dict]:
+    """Parse noisefit inputs"""
     if x.ndim != NUM_NOISE_DATA_DIMENSIONS:
         msg = "Data array x must be 2D"
         raise ValueError(msg)
     n, m = x.shape
 
-    dt = _validate_sampling_time(dt)
-
-    if v0 is None:
-        v0 = np.mean(np.var(x, 1)) * np.ones(NUM_NOISE_PARAMETERS)
-    else:
-        v0 = np.asarray(v0)
-        if v0.size != NUM_NOISE_PARAMETERS:
-            msg = (
-                "Noise parameter array logv must have "
-                f"{NUM_NOISE_PARAMETERS} elements."
-            )
-            raise ValueError(msg)
+    if (
+        fix_sigma_alpha
+        and fix_sigma_beta
+        and fix_sigma_tau
+        and fix_mu
+        and fix_a
+        and fix_eta
+    ):
+        msg = "All variables are fixed"
+        raise ValueError(msg)
 
+    # Use the first column of x as a default
     if mu0 is None:
         mu0 = x[:, 0]
     else:
-        mu0 = np.asarray(mu0)
+        mu0 = np.asarray(mu0, dtype=np.float64)
         if mu0.size != n:
             msg = "Size of mu0 is incompatible with data array x."
             raise ValueError(msg)
 
-    scale_logv = 1e0 * np.ones(3)
-    alpha, beta, tau = np.sqrt(v0)
-    noise_model = NoiseModel(alpha, beta, tau, dt)
-    scale_delta = 1e-0 * noise_model.amplitude(x[:, 0])
-    scale_alpha = 1e-2 * np.ones(m - 1)
-    scale_eta = 1e-3 * np.ones(m - 1)
-    scale_v = 1.0e-6
+    # Compute time-dependent variance of signal array and
+    # the minimum value of the time-dependent noise amplitude
+    v_t = np.var(x, 1, ddof=1)
+    sigma_min = np.sqrt(np.min(v_t))
+
+    # If any initial guess for the noise parameters is unspecified,
+    # estimate all noise parameters with a linear least-squares fit
+    # to the time-dependent variance
+    if None in [sigma_alpha0, sigma_beta0, sigma_tau0]:
+        mu0_f = np.fft.rfft(mu0)
+        w = 2 * np.pi * np.fft.rfftfreq(n, dt)
+        dmu0_dt = np.fft.irfft(1j * w * mu0_f, n=n)
+        a_matrix = np.stack([np.ones(n), mu0**2, dmu0_dt**2], axis=1)
+        sol = np.linalg.lstsq(a_matrix, v_t, rcond=None)
+        sigma_est = np.ma.sqrt(sol[0]).filled(sigma_min)
+
+    if sigma_alpha0 is None:
+        sigma_alpha0 = sigma_est[0]
+
+    if sigma_beta0 is None:
+        sigma_beta0 = sigma_est[1]
+
+    if sigma_tau0 is None:
+        sigma_tau0 = sigma_est[2]
+
+    noise_model = NoiseModel(sigma_alpha0, sigma_beta0, sigma_tau0, dt=dt)
+
+    if scale_sigma_alpha is None:
+        if np.isclose(sigma_alpha0, 0.0):
+            scale_sigma_alpha = 1.0
+        else:
+            scale_sigma_alpha = sigma_alpha0
+
+    if scale_sigma_beta is None:
+        if np.isclose(sigma_beta0, 0.0):
+            scale_sigma_beta = 1.0
+        else:
+            scale_sigma_beta = sigma_beta0
+
+    if scale_sigma_tau is None:
+        if np.isclose(sigma_tau0, 0.0):
+            scale_sigma_tau = 1.0
+        else:
+            scale_sigma_tau = sigma_tau0 / dt
+
+    if scale_delta_mu is None:
+        scale_delta_mu = noise_model.noise_amp(mu0)
+        scale_delta_mu[np.isclose(scale_delta_mu, 0.0)] = np.sqrt(
+            np.finfo(float).eps
+        )
+
+    if scale_delta_a is None:
+        scale_delta_a_amp = np.max((sigma_min, sigma_beta0))
+        scale_delta_a = scale_delta_a_amp * np.ones(m - 1)
 
-    # Replace log(x) with -inf when x <= 0
-    v0_scaled = np.asarray(v0, dtype=float) / scale_v
-    logv0_scaled = np.ma.log(v0_scaled).filled(-np.inf) / scale_logv
-    delta0 = (x[:, 0] - mu0) / scale_delta
+    scale_delta_a = np.asarray(scale_delta_a, dtype=np.float64)
+
+    if scale_eta is None:
+        scale_eta_amp = np.max((sigma_min, sigma_tau0))
+        scale_eta = scale_eta_amp * np.ones(m - 1) / dt
+
+    scale_eta = np.asarray(scale_eta, dtype=np.float64)
+
+    scale_sigma = np.array(
+        [scale_sigma_alpha, scale_sigma_beta, scale_sigma_tau],
+        dtype=np.float64,
+    )
+
+    # Replace log(x) with -1e2 when x <= 0
+    v0_scaled = (
+        np.asarray(
+            [sigma_alpha0, sigma_beta0, sigma_tau0 / dt], dtype=np.float64
+        )
+        / scale_sigma
+    ) ** 2
+    logv0_scaled = np.ma.log(v0_scaled).filled(-1.0e2)
+    delta0 = (x[:, 0] - mu0) / scale_delta_mu
 
     if a0 is None:
         a0 = np.ones(m)
     else:
-        a0 = np.asarray(a0)
+        a0 = np.asarray(a0, dtype=np.float64)
         if a0.size != m:
             msg = "Size of a0 is incompatible with data array x."
             raise ValueError(msg)
 
-    alpha0 = (a0[1:] - 1.0) / (a0[0] * scale_alpha)
+    epsilon0 = (a0[1:] - 1.0) / (a0[0] * scale_delta_a)
 
     if eta0 is None:
         eta0 = np.zeros(m)
     else:
-        eta0 = np.asarray(eta0)
+        eta0 = np.asarray(eta0, dtype=np.float64)
         if eta0.size != m:
             msg = "Size of eta0 is incompatible with data array x."
             raise ValueError(msg)
 
-    eta0 = eta0[1:] / scale_eta
+    eta_on_dt0 = eta0[1:] / dt / scale_eta
 
     # Set initial guesses for all free parameters
-    x0 = np.array([])
-    if not fix_v:
-        x0 = np.concatenate((x0, logv0_scaled))
+    x0 = np.array([], dtype=np.float64)
+    if not fix_sigma_alpha:
+        x0 = np.append(x0, logv0_scaled[0])
+    if not fix_sigma_beta:
+        x0 = np.append(x0, logv0_scaled[1])
+    if not fix_sigma_tau:
+        x0 = np.append(x0, logv0_scaled[2])
     if not fix_mu:
         x0 = np.concatenate((x0, delta0))
     if not fix_a:
-        x0 = np.concatenate((x0, alpha0))
+        x0 = np.concatenate((x0, epsilon0))
     if not fix_eta:
-        x0 = np.concatenate((x0, eta0))
+        x0 = np.concatenate((x0, eta_on_dt0))
 
     # Bundle free parameters together into objective function
     def objective(_p):
-        if fix_v:
-            _logv = logv0_scaled
+        if fix_sigma_alpha:
+            _logv_alpha = logv0_scaled[0]
+        else:
+            _logv_alpha = _p[0]
+            _p = _p[1:]
+        if fix_sigma_beta:
+            _logv_beta = logv0_scaled[1]
         else:
-            _logv = _p[:3]
-            _p = _p[3:]
+            _logv_beta = _p[0]
+            _p = _p[1:]
+        if fix_sigma_tau:
+            _logv_tau = logv0_scaled[2]
+        else:
+            _logv_tau = _p[0]
+            _p = _p[1:]
         if fix_mu:
             _delta = delta0
         else:
             _delta = _p[:n]
             _p = _p[n:]
         if fix_a:
-            _alpha = alpha0
+            _epsilon = epsilon0
         else:
-            _alpha = _p[: m - 1]
+            _epsilon = _p[: m - 1]
             _p = _p[m - 1 :]
         if fix_eta:
-            _eta = eta0
+            _eta_on_dt = eta_on_dt0
         else:
-            _eta = _p[: m - 1]
-        return _tdnll_scaled(
+            _eta_on_dt = _p[: m - 1]
+        return _costfun_noisefit(
             x.T,
-            _logv,
+            _logv_alpha,
+            _logv_beta,
+            _logv_tau,
             _delta,
-            _alpha,
-            _eta / dt,
-            fix_logv=fix_v,
-            fix_delta=fix_mu,
-            fix_alpha=fix_a,
+            _epsilon,
+            _eta_on_dt,
+            fix_logv_alpha=fix_sigma_alpha,
+            fix_logv_beta=fix_sigma_beta,
+            fix_logv_tau=fix_sigma_tau,
+            fix_delta_mu=fix_mu,
+            fix_delta_a=fix_a,
             fix_eta=fix_eta,
-            scale_logv=scale_logv,
-            scale_delta=scale_delta,
-            scale_alpha=scale_alpha,
+            scale_sigma_alpha=scale_sigma_alpha,
+            scale_sigma_beta=scale_sigma_beta,
+            scale_sigma_tau=scale_sigma_tau,
+            scale_delta_mu=scale_delta_mu,
+            scale_delta_a=scale_delta_a,
             scale_eta=scale_eta,
-            scale_v=scale_v,
         )
 
-    # Minimize cost function with respect to free parameters
-    out = minimize(
-        objective,
-        x0,
-        method="BFGS",
-        jac=True,
-    )
+    input_parsed = {
+        "sigma_alpha0": sigma_alpha0,
+        "sigma_beta0": sigma_beta0,
+        "sigma_tau0": sigma_tau0,
+        "mu0": mu0,
+        "a0": a0,
+        "eta0": eta0,
+        "fix_sigma_alpha": fix_sigma_alpha,
+        "fix_sigma_beta": fix_sigma_beta,
+        "fix_sigma_tau": fix_sigma_tau,
+        "fix_mu": fix_mu,
+        "fix_a": fix_a,
+        "fix_eta": fix_eta,
+        "scale_sigma_alpha": scale_sigma_alpha,
+        "scale_sigma_beta": scale_sigma_beta,
+        "scale_sigma_tau": scale_sigma_tau,
+        "scale_delta_mu": scale_delta_mu,
+        "scale_delta_a": scale_delta_a,
+        "scale_eta": scale_eta,
+    }
+    return objective, x0, input_parsed
+
 
+def _parse_noisefit_output(
+    out: OptimizeResult,
+    x: NDArray[np.float64],
+    dt: float,
+    *,
+    sigma_alpha0: float,
+    sigma_beta0: float,
+    sigma_tau0: float,
+    mu0: NDArray[np.float64],
+    a0: NDArray[np.float64],
+    eta0: NDArray[np.float64],
+    fix_sigma_alpha: bool,
+    fix_sigma_beta: bool,
+    fix_sigma_tau: bool,
+    fix_mu: bool,
+    fix_a: bool,
+    fix_eta: bool,
+    scale_sigma_alpha: float,
+    scale_sigma_beta: float,
+    scale_sigma_tau: float,
+    scale_delta_mu: NDArray[np.float64],
+    scale_delta_a: NDArray[np.float64],
+    scale_eta: NDArray[np.float64],
+) -> NoiseResult:
+    """Parse noisefit output"""
     # Parse output
+    n, m = x.shape
+
     x_out = out.x
-    if fix_v:
-        v_out = v0_scaled * scale_v
+    if fix_sigma_alpha:
+        alpha = sigma_alpha0
     else:
-        v_out = np.exp(x_out[:3] * scale_logv) * scale_v
-        x_out = x_out[3:]
-    alpha, beta, tau = np.sqrt(v_out)
-    noise_model = NoiseModel(alpha, beta, tau, dt)
+        alpha = np.exp(x_out[0] / 2) * scale_sigma_alpha
+        x_out = x_out[1:]
+    if fix_sigma_beta:
+        beta = sigma_beta0
+    else:
+        beta = np.exp(x_out[0] / 2) * scale_sigma_beta
+        x_out = x_out[1:]
+    if fix_sigma_tau:
+        tau = sigma_tau0
+    else:
+        tau = np.exp(x_out[0] / 2) * scale_sigma_tau * dt
+        x_out = x_out[1:]
+    # noinspection PyArgumentList
+    noise_model = NoiseModel(alpha, beta, tau, dt=dt)
 
     if fix_mu:
         mu_out = mu0
     else:
-        mu_out = x[:, 0] - x_out[:n] * scale_delta
+        mu_out = x[:, 0] - x_out[:n] * scale_delta_mu
         x_out = x_out[n:]
 
     if fix_a:
         a_out = a0
     else:
-        a_out = np.concatenate(([1.0], 1.0 + x_out[: m - 1] * scale_alpha))
+        a_out = np.concatenate(([1.0], 1.0 + x_out[: m - 1] * scale_delta_a))
         x_out = x_out[m - 1 :]
 
     if fix_eta:
         eta_out = eta0
     else:
-        eta_out = np.concatenate(([0.0], x_out[: m - 1] * scale_eta))
-
-    diagnostic = {
-        "grad_scaled": out.jac,
-        "hess_inv_scaled": out.hess_inv,
-        "err": {
-            "var": np.array([]),
-            "delta": np.array([]),
-            "a": np.array([]),
-            "eta": np.array([]),
-        },
-        "success": out.success,
-        "status": out.status,
-        "message": out.message,
-        "nfev": out.nfev,
-        "njev": out.njev,
-        "nit": out.nit,
-    }
+        eta_out = np.concatenate(([0.0], x_out[: m - 1] * dt * scale_eta))
+
+    diagnostic = out
+    fun = out.fun
 
     # Concatenate scaling vectors for all sets of free parameters
     scale_hess_inv = np.concatenate(
         [
             val
             for tf, val in zip(
-                [fix_v, fix_mu, fix_a, fix_eta],
-                [scale_logv, scale_delta, scale_alpha, scale_eta],
+                [
+                    fix_sigma_alpha,
+                    fix_sigma_beta,
+                    fix_sigma_tau,
+                    fix_mu,
+                    fix_a,
+                    fix_eta,
+                ],
+                [
+                    [scale_sigma_alpha],
+                    [scale_sigma_beta],
+                    [scale_sigma_tau],
+                    scale_delta_mu,
+                    scale_delta_a,
+                    scale_eta,
+                ],
             )
             if not tf
         ]
     )
 
     # Convert inverse Hessian into unscaled parameters
-    hess_inv = scale_v * (
-        np.diag(scale_hess_inv)
-        @ diagnostic["hess_inv_scaled"]
-        @ np.diag(scale_hess_inv)
-    )
+    hess_inv = np.diag(scale_hess_inv) @ out.hess_inv @ np.diag(scale_hess_inv)
 
     # Determine parameter uncertainty vector from diagonal entries
     err = np.sqrt(np.diag(hess_inv))
+    err_mu = np.array([])
+    err_a = np.array([])
+    err_eta = np.array([])
 
     # Parse error vector
-    if not fix_v:
-        # Propagate error from log(V) to V
-        diagnostic["err"]["var"] = np.sqrt(
-            np.diag(np.diag(v_out) @ hess_inv[0:3, 0:3]) @ np.diag(v_out)
-        )
-        err = err[3:]
+    # Propagate error from log(V) to sigma
+    if not fix_sigma_alpha:
+        err_sigma_alpha = np.sqrt(0.5 * alpha * err[0])
+        err = err[1:]
+    else:
+        err_sigma_alpha = 0.0
+
+    if not fix_sigma_beta:
+        err_sigma_beta = np.sqrt(0.5 * beta * err[0])
+        err = err[1:]
+    else:
+        err_sigma_beta = 0.0
+
+    if not fix_sigma_tau:
+        err_sigma_tau = np.sqrt(0.5 * tau * err[0])
+        err = err[1:]
+    else:
+        err_sigma_tau = 0.0
 
     if not fix_mu:
-        diagnostic["err"]["delta"] = err[:n]
+        err_mu = err[:n]
         err = err[n:]
 
     if not fix_a:
-        diagnostic["err"]["alpha"] = np.concatenate(([0], err[: m - 1]))
+        err_a = np.concatenate(([0], err[: m - 1]))
         err = err[m - 1 :]
 
     if not fix_eta:
-        diagnostic["err"]["eta"] = np.concatenate(([0], err[: m - 1]))
+        err_eta = np.concatenate(([0], err[: m - 1]))
 
     return NoiseResult(
-        noise_model, mu_out, a_out, eta_out, out.fun / scale_v, diagnostic
+        noise_model,
+        mu_out,
+        a_out,
+        eta_out,
+        fun,
+        hess_inv,
+        err_sigma_alpha,
+        err_sigma_beta,
+        err_sigma_tau,
+        err_mu,
+        err_a,
+        err_eta,
+        diagnostic,
     )
 
 
@@ -1271,7 +1807,7 @@ def fit(
     x: ArrayLike,
     y: ArrayLike,
     *,
-    dt: float = 1.0,
+    dt: float | None = None,
     sigma_parms: ArrayLike = (1.0, 0.0, 0.0),
     f_bounds: ArrayLike = (0.0, np.inf),
     p_bounds: ArrayLike | None = None,
@@ -1291,7 +1827,7 @@ def fit(
     fun : callable
         Transfer function.
 
-            ``fun(p, w, *args, **kwargs) -> np.ndarray``
+            ``fun(p, w, *args, **kwargs) -> ndarray``
 
         Assumes the :math:`+i\omega t` convention for harmonic time dependence.
     p0 : array_like
@@ -1300,8 +1836,14 @@ def fit(
         Measured input signal.
     y : array_like
         Measured output signal.
-    dt : float, optional
-        Sampling time.
+    dt : float or None, optional
+        Sampling time, normally in picoseconds. Default is None, which sets
+        the sampling time to ``thztools.options.sampling_time``. If both
+        ``dt`` and ``thztools.options.sampling_time`` are ``None``, the
+        sampling time is set to ``1.0``. In this case, the angular frequency
+        ``omega`` must be given in units of radians per sampling time, and any
+        parameters in ``args`` must be expressed with the sampling time as the
+        unit of time.
     sigma_parms : None or array_like, optional
         Noise parameters with size (3,), expressed as standard deviation
         amplitudes.
@@ -1338,15 +1880,23 @@ def fit(
                 Resiuals of the output waveform ``y``.
             success : bool
                 True if one of the convergence criteria is satisfied.
+
+    Warns
+    -----
+    UserWarning
+        If ``thztools.options.sampling_time`` and the ``dt`` parameter
+        are both not ``None`` and are set to different ``float`` values, the
+        function will set the sampling time to ``dt`` and raise a
+        :class:`UserWarning`.
     """
     fit_method = "trf"
 
-    p0 = np.asarray(p0)
-    x = np.asarray(x)
-    y = np.asarray(y)
-    sigma_parms = np.asarray(sigma_parms)
+    p0 = np.asarray(p0, dtype=np.float64)
+    x = np.asarray(x, dtype=np.float64)
+    y = np.asarray(y, dtype=np.float64)
+    sigma_parms = np.asarray(sigma_parms, dtype=np.float64)
 
-    dt = _validate_sampling_time(dt)
+    dt = _assign_sampling_time(dt)
 
     n = y.shape[-1]
     n_p = len(p0)
@@ -1355,7 +1905,7 @@ def fit(
         p_bounds = (-np.inf, np.inf)
         fit_method = "lm"
     else:
-        p_bounds = np.asarray(p_bounds)
+        p_bounds = np.asarray(p_bounds, dtype=np.float64)
         if p_bounds.shape[0] == 2:  # noqa: PLR2004
             p_bounds = (
                 np.concatenate((p_bounds[0], np.full((n,), -np.inf))),
@@ -1369,16 +1919,17 @@ def fit(
         kwargs = {}
 
     w = 2 * np.pi * rfftfreq(n, dt)
-    w_bounds = 2 * np.pi * np.asarray(f_bounds)
+    w_bounds = 2 * np.pi * np.asarray(f_bounds, dtype=np.float64)
     w_below_idx = w <= w_bounds[0]
     w_above_idx = w > w_bounds[1]
     w_in_idx = np.invert(w_below_idx) * np.invert(w_above_idx)
     n_f = len(w)
 
     alpha, beta, tau = sigma_parms
-    noise_model = NoiseModel(alpha, beta, tau, dt)
-    sigma_x = noise_model.amplitude(x)
-    sigma_y = noise_model.amplitude(y)
+    # noinspection PyArgumentList
+    noise_model = NoiseModel(alpha, beta, tau, dt=dt)
+    sigma_x = noise_model.noise_amp(x)
+    sigma_y = noise_model.noise_amp(y)
     p0_est = np.concatenate((p0, np.zeros(n)))
 
     def etfe(_x, _y):


Test patch:
diff --git a/examples/tls-test.ipynb b/examples/tls-test.ipynb
index c38588d..5e0d6b3 100644
--- a/examples/tls-test.ipynb
+++ b/examples/tls-test.ipynb
@@ -19,7 +19,6 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "outputs": [],
    "source": [
     "import numpy as np\n",
     "from numpy import pi\n",
@@ -39,7 +38,8 @@
    ],
    "metadata": {
     "collapsed": false
-   }
+   },
+   "outputs": []
   },
   {
    "cell_type": "markdown",
@@ -53,7 +53,6 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "outputs": [],
    "source": [
     "n = 256  # Number of samples\n",
     "ts = 0.05  # Sampling time [ps]\n",
@@ -70,7 +69,8 @@
    ],
    "metadata": {
     "collapsed": false
-   }
+   },
+   "outputs": []
   },
   {
    "cell_type": "markdown",
@@ -86,7 +86,6 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "outputs": [],
    "source": [
     "def tf_model(_p, _w):\n",
     "    return _p[0] * np.exp(-1j * _w * _p[1])\n",
@@ -101,12 +100,12 @@
    ],
    "metadata": {
     "collapsed": false
-   }
+   },
+   "outputs": []
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "outputs": [],
    "source": [
     "n_sim = 250\n",
     "p_opt = np.array([0.5, 1])\n",
@@ -125,12 +124,12 @@
    ],
    "metadata": {
     "collapsed": false
-   }
+   },
+   "outputs": []
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "outputs": [],
    "source": [
     "plt.plot(t, mu)\n",
     "plt.plot(t, td_model(p_opt, mu))\n",
@@ -138,7 +137,8 @@
    ],
    "metadata": {
     "collapsed": false
-   }
+   },
+   "outputs": []
   },
   {
    "cell_type": "markdown",
@@ -152,7 +152,6 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "outputs": [],
    "source": [
     "%%time\n",
     "\n",
@@ -228,12 +227,12 @@
    ],
    "metadata": {
     "collapsed": false
-   }
+   },
+   "outputs": []
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "outputs": [],
    "source": [
     "# Find the index of the median value\n",
     "i_show = np.argsort(resnorm_tls)[len(resnorm_tls) // 2]\n",
@@ -325,12 +324,12 @@
    ],
    "metadata": {
     "collapsed": false
-   }
+   },
+   "outputs": []
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "outputs": [],
    "source": [
     "print(\n",
     "    f\"Amplitude parameter estimates: \"\n",
@@ -359,7 +358,8 @@
    ],
    "metadata": {
     "collapsed": false
-   }
+   },
+   "outputs": []
   },
   {
    "cell_type": "markdown",
@@ -375,19 +375,18 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "outputs": [],
    "source": [
     "osm_a_tls, osr_a_tls = probplot(p_tls[0, :], fit=False)\n",
     "osm_eta_tls, osr_eta_tls = probplot(p_tls[1, :], fit=False)"
    ],
    "metadata": {
     "collapsed": false
-   }
+   },
+   "outputs": []
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "outputs": [],
    "source": [
     "w, h = figaspect(0.45)\n",
     "\n",
@@ -427,7 +426,8 @@
    ],
    "metadata": {
     "collapsed": false
-   }
+   },
+   "outputs": []
   }
  ],
  "metadata": {
diff --git a/tests/test_thztools.py b/tests/test_thztools.py
index 2be29b5..6e42653 100644
--- a/tests/test_thztools.py
+++ b/tests/test_thztools.py
@@ -2,7 +2,6 @@
 
 import numpy as np
 import pytest
-import scipy
 from numpy import pi
 from numpy.testing import assert_allclose
 from numpy.typing import ArrayLike
@@ -10,13 +9,14 @@
 import thztools
 from thztools.thztools import (
     NoiseModel,
+    _assign_sampling_time,
+    _costfun_noisefit,
     _costfuntls,
-    _tdnll_scaled,
-    _validate_sampling_time,
+    _parse_noisefit_input,
     fit,
+    noisefit,
     scaleshift,
-    tdnoisefit,
-    transfer_out,
+    transfer,
     wave,
 )
 
@@ -41,39 +41,32 @@ def jac_fun(p, w):
     return np.stack((exp_ipw, 1j * w * p[0] * exp_ipw)).T
 
 
-class TestGlobalOptions:
-    thztools.global_options.sampling_time = None
-
+class TestOptions:
     def test_sampling_time(self):
-        assert thztools.global_options.sampling_time is None
-        thztools.global_options.sampling_time = 0.1
-        assert_allclose(thztools.global_options.sampling_time, 0.1)
+        assert thztools.get_option("sampling_time") is None
 
     @pytest.mark.parametrize("global_sampling_time", [None, 0.1])
     @pytest.mark.parametrize("dt", [None, 0.1, 1.0])
-    def test_validation(self, global_sampling_time, dt):
-        thztools.global_options.sampling_time = global_sampling_time
+    def test_assignment(self, global_sampling_time, dt, monkeypatch):
+        monkeypatch.setattr(
+            thztools.options, "sampling_time", global_sampling_time
+        )
         if global_sampling_time is None and dt is None:
-            assert np.isclose(_validate_sampling_time(dt), 1.0)
+            assert np.isclose(_assign_sampling_time(dt), 1.0)
         elif global_sampling_time is None and dt is not None:
-            assert np.isclose(_validate_sampling_time(dt), dt)
+            assert np.isclose(_assign_sampling_time(dt), dt)
         elif global_sampling_time is not None and dt is None:
-            assert np.isclose(
-                _validate_sampling_time(dt), global_sampling_time
-            )
+            assert np.isclose(_assign_sampling_time(dt), global_sampling_time)
         elif global_sampling_time is not None and np.isclose(
             dt, global_sampling_time
         ):
-            assert np.isclose(
-                _validate_sampling_time(dt), global_sampling_time
-            )
+            assert np.isclose(_assign_sampling_time(dt), global_sampling_time)
         else:
             with pytest.warns(UserWarning):
-                assert np.isclose(_validate_sampling_time(dt), dt)
+                assert np.isclose(_assign_sampling_time(dt), dt)
 
 
 class TestNoiseModel:
-    thztools.global_options.sampling_time = None
     n = 16
     dt = 1.0 / n
     t = np.arange(n) * dt
@@ -103,10 +96,10 @@ def test_var_definition(
     ) -> None:
         if dt is None:
             noise_model = NoiseModel(alpha, beta, tau)
-            result = noise_model.variance(mu, axis=axis)
+            result = noise_model.noise_var(mu, axis=axis)
         else:
-            noise_model = NoiseModel(alpha, beta, tau, dt)
-            result = noise_model.variance(mu, axis=axis)
+            noise_model = NoiseModel(alpha, beta, tau, dt=dt)
+            result = noise_model.noise_var(mu, axis=axis)
         assert_allclose(result, expected, atol=atol, rtol=rtol)  # type: ignore
 
     @pytest.mark.parametrize(
@@ -132,10 +125,10 @@ def test_amp_definition(
     ) -> None:
         if dt is None:
             noise_model = NoiseModel(alpha, beta, tau)
-            result = noise_model.amplitude(mu, axis=axis)
+            result = noise_model.noise_amp(mu, axis=axis)
         else:
-            noise_model = NoiseModel(alpha, beta, tau, dt)
-            result = noise_model.amplitude(mu, axis=axis)
+            noise_model = NoiseModel(alpha, beta, tau, dt=dt)
+            result = noise_model.noise_amp(mu, axis=axis)
         assert_allclose(result, expected, atol=atol, rtol=rtol)  # type: ignore
 
     @pytest.mark.parametrize(
@@ -159,10 +152,10 @@ def test_noise_definition(
     ) -> None:
         if dt is None:
             noise_model = NoiseModel(alpha, beta, tau)
-            result = noise_model.noise(mu, axis=axis)
+            result = noise_model.noise_sim(mu, axis=axis)
         else:
-            noise_model = NoiseModel(alpha, beta, tau, dt)
-            result = noise_model.noise(mu, axis=axis)
+            noise_model = NoiseModel(alpha, beta, tau, dt=dt)
+            result = noise_model.noise_sim(mu, axis=axis)
         assert result.shape == expected
 
 
@@ -182,9 +175,8 @@ class TestTransferOut:
     )
     def test_inputs(self, t_fun, x, fft_sign, p, expected):
         ts = self.dt
-        thztools.global_options.sampling_time = None
         assert_allclose(
-            transfer_out(t_fun, x, dt=ts, fft_sign=fft_sign, args=p),
+            transfer(t_fun, x, dt=ts, numpy_sign_convention=fft_sign, args=p),
             expected,
         )
 
@@ -192,49 +184,77 @@ def test_inputs(self, t_fun, x, fft_sign, p, expected):
     def test_error(self, x):
         dt = self.dt
         with pytest.raises(ValueError):
-            _ = transfer_out(x, tfun1, dt=dt, args=[1.0, 0.0])
+            _ = transfer(tfun1, x, dt=dt, args=[1.0, 0.0])
 
 
-class TestTHzGen:
-    thztools.global_options.sampling_time = None
+class TestTimebase:
+    @pytest.mark.parametrize(
+        "dt",
+        [
+            None,
+            1.0,
+            2.0,
+        ],
+    )
+    @pytest.mark.parametrize(
+        "t_init",
+        [
+            None,
+            1.0,
+            2.0,
+        ],
+    )
+    def test_timebase(self, dt, t_init):
+        n = 8
+        if t_init is None:
+            t = thztools.timebase(n, dt=dt)
+            t_init = 0.0
+        else:
+            t = thztools.timebase(n, dt=dt, t_init=t_init)
+        dt = _assign_sampling_time(dt)
+        t_expected = t_init + np.arange(n) * dt
+        assert_allclose(t, t_expected, rtol=rtol, atol=atol)
+
 
+class TestWave:
+    @pytest.mark.parametrize(
+        "t0",
+        [2.4, None],
+    )
     @pytest.mark.parametrize(
         "kwargs",
         [
             {},
             {"a": 1.0},
-            {"taur": 0.3},
-            {"tauc": 0.1},
-            {"fwhm": 0.05},
+            {"taur": 6.0},
+            {"tauc": 2.0},
+            {"fwhm": 1.0},
         ],
     )
-    def test_inputs(self, kwargs: dict) -> None:
+    def test_inputs(self, t0: float | None, kwargs: dict) -> None:
         n = 8
         dt = 1.0
-        t0 = 2.0
         y_expected = np.array(
             [
-                0.05651348,
-                -0.13522073,
+                0.07767792,
+                -0.63041598,
+                -1.03807638,
+                -0.81199085,
+                -0.03955531,
+                0.72418661,
                 1.0,
-                -0.65804181,
-                -0.28067975,
-                0.05182924,
-                -0.01837401,
-                -0.01602642,
+                0.71817398,
             ]
         )
-        t_expected = np.arange(n)
         assert_allclose(
-            wave(n, dt, t0, **kwargs),  # type: ignore
-            (y_expected, t_expected),  # type: ignore
+            wave(n, dt=dt, t0=t0, **kwargs),  # type: ignore
+            y_expected,  # type: ignore
             atol=atol,
             rtol=rtol,
         )
 
 
 class TestScaleShift:
-    thztools.global_options.sampling_time = None
     n = 16
     dt = 1.0 / n
     t = np.arange(n) * dt
@@ -302,12 +322,11 @@ def test_inputs(
         "x, kwargs", [[x, {"a": [2, 0.5]}], [x, {"eta": [1, -1]}]]
     )
     def test_errors(self, x: ArrayLike, kwargs: dict) -> None:
-        with pytest.raises(ValueError):
+        with pytest.raises(ValueError, match="correction with shape"):
             scaleshift(x, **kwargs)
 
 
 class TestCostFunTLS:
-    thztools.global_options.sampling_time = None
     theta = (1, 0)
     mu = np.arange(8)
     xx = mu
@@ -323,21 +342,48 @@ class TestCostFunTLS:
 
 
 class TestTDNLLScaled:
-    thztools.global_options.sampling_time = None
     m = 2
     n = 16
     dt = 1.0 / n
     t = np.arange(n) * dt
     mu = np.cos(2 * pi * t)
     x = np.tile(mu, [m, 1])
-    logv = (0, -np.inf, -np.inf)
-    delta = np.zeros(n)
-    alpha = np.zeros(m - 1)
+    logv_alpha = 0
+    logv_beta = -np.inf
+    logv_tau = -np.inf
+    delta_mu = np.zeros(n)
+    delta_a = np.zeros(m - 1)
     eta = np.zeros(m - 1)
     desired_nll = x.size * np.log(2 * pi) / 2
 
     @pytest.mark.parametrize(
-        "fix_logv, desired_gradnll_logv",
+        "fix_logv_alpha, desired_gradnll_logv_alpha",
+        [
+            [
+                True,
+                [],
+            ],
+            [
+                False,
+                [1.0],
+            ],
+        ],
+    )
+    @pytest.mark.parametrize(
+        "fix_logv_beta, desired_gradnll_logv_beta",
+        [
+            [
+                True,
+                [],
+            ],
+            [
+                False,
+                [0.0],
+            ],
+        ],
+    )
+    @pytest.mark.parametrize(
+        "fix_logv_tau, desired_gradnll_logv_tau",
         [
             [
                 True,
@@ -345,12 +391,12 @@ class TestTDNLLScaled:
             ],
             [
                 False,
-                [16.0, 0.0, 0.0],
+                [0.0],
             ],
         ],
     )
     @pytest.mark.parametrize(
-        "fix_delta, desired_gradnll_delta",
+        "fix_delta_mu, desired_gradnll_delta_mu",
         [
             [
                 True,
@@ -363,7 +409,7 @@ class TestTDNLLScaled:
         ],
     )
     @pytest.mark.parametrize(
-        "fix_alpha, desired_gradnll_alpha",
+        "fix_delta_a, desired_gradnll_delta_a",
         [
             [
                 True,
@@ -390,117 +436,226 @@ class TestTDNLLScaled:
     )
     def test_gradnll_calc(
         self,
-        fix_logv,
-        fix_delta,
-        fix_alpha,
+        fix_logv_alpha,
+        fix_logv_beta,
+        fix_logv_tau,
+        fix_delta_mu,
+        fix_delta_a,
         fix_eta,
-        desired_gradnll_logv,
-        desired_gradnll_delta,
-        desired_gradnll_alpha,
+        desired_gradnll_logv_alpha,
+        desired_gradnll_logv_beta,
+        desired_gradnll_logv_tau,
+        desired_gradnll_delta_mu,
+        desired_gradnll_delta_a,
         desired_gradnll_eta,
     ):
         n = self.n
         m = self.m
         x = self.x
-        logv = self.logv
-        delta = self.delta
-        alpha = self.alpha
+        logv_alpha = self.logv_alpha
+        logv_beta = self.logv_beta
+        logv_tau = self.logv_tau
+        delta_mu = self.delta_mu
+        delta_a = self.delta_a
         eta_on_dt = self.eta / self.dt
         desired_gradnll = np.concatenate(
             (
-                desired_gradnll_logv,
-                desired_gradnll_delta,
-                desired_gradnll_alpha,
+                desired_gradnll_logv_alpha,
+                desired_gradnll_logv_beta,
+                desired_gradnll_logv_tau,
+                desired_gradnll_delta_mu,
+                desired_gradnll_delta_a,
                 desired_gradnll_eta,
             )
         )
-        _, gradnll = _tdnll_scaled(
+        _, gradnll = _costfun_noisefit(
             x,
-            logv,
-            delta,
-            alpha,
+            logv_alpha,
+            logv_beta,
+            logv_tau,
+            delta_mu,
+            delta_a,
             eta_on_dt,
-            fix_logv=fix_logv,
-            fix_delta=fix_delta,
-            fix_alpha=fix_alpha,
+            fix_logv_alpha=fix_logv_alpha,
+            fix_logv_beta=fix_logv_beta,
+            fix_logv_tau=fix_logv_tau,
+            fix_delta_mu=fix_delta_mu,
+            fix_delta_a=fix_delta_a,
             fix_eta=fix_eta,
-            scale_logv=np.ones(3),
-            scale_delta=np.ones(n),
-            scale_alpha=np.ones(m - 1),
+            scale_sigma_alpha=1.0,
+            scale_sigma_beta=1.0,
+            scale_sigma_tau=1.0,
+            scale_delta_mu=np.ones(n),
+            scale_delta_a=np.ones(m - 1),
             scale_eta=np.ones(m - 1),
-            scale_v=1.0,
         )
         assert_allclose(
             gradnll, desired_gradnll, atol=10 * np.finfo(float).eps
         )
 
 
-class TestTDNoiseFit:
-    thztools.global_options.sampling_time = None
+class TestNoiseFit:
     rng = np.random.default_rng(0)
     n = 256
     m = 64
     dt = 0.05
     t = np.arange(n) * dt
-    mu, _ = wave(n, dt=dt, t0=n * dt / 3)
-    alpha, beta, tau = 1e-5, 1e-2, 1e-3
+    mu = wave(n, dt=dt, t0=n * dt / 3)
+    alpha, beta, tau = 1e-5, 1e-3, 1e-3
     sigma = np.array([alpha, beta, tau])
     noise_model = NoiseModel(alpha, beta, tau, dt=dt)
-    noise = noise_model.amplitude(mu) * rng.standard_normal((m, n))
+    noise = noise_model.noise_sim(np.ones((m, 1)) * mu, seed=0)
+    noise_amp = noise_model.noise_amp(mu)
     x = np.array(mu + noise)
     a = np.ones(m)
     eta = np.zeros(m)
+    scale_delta_a = 1e-2 * np.ones(m - 1)
+    scale_eta = 1e-3 * np.ones(m - 1) / dt
 
-    @pytest.mark.parametrize("x", [x, x[:, 0]])
-    @pytest.mark.parametrize("v0", [None, sigma**2, []])
-    @pytest.mark.parametrize("mu0", [None, mu, []])
-    @pytest.mark.parametrize("a0", [None, a, []])
-    @pytest.mark.parametrize("eta0", [None, eta, []])
-    @pytest.mark.parametrize("fix_v", [True, False])
-    @pytest.mark.parametrize("fix_mu", [True, False])
-    @pytest.mark.parametrize("fix_a", [True, False])
-    @pytest.mark.parametrize("fix_eta", [True, False])
-    def test_inputs(self, x, v0, mu0, a0, eta0, fix_v, fix_mu, fix_a, fix_eta):
-        print(f"{scipy.__version__=}")
-        n = self.n
+    @pytest.mark.parametrize(
+        "x, mu0, a0, eta0, fix_sigma_alpha, fix_sigma_beta, fix_sigma_tau, "
+        "fix_mu, fix_a, fix_eta",
+        [
+            [x[:, 0], mu, a, eta, False, False, False, False, False, False],
+            [x, [], a, eta, False, False, False, False, False, False],
+            [x, mu, [], eta, False, False, False, False, False, False],
+            [x, mu, a, [], False, False, False, False, False, False],
+            [x, mu, a, eta, True, True, True, True, True, True],
+        ],
+    )
+    def test_exceptions(
+        self,
+        x,
+        mu0,
+        a0,
+        eta0,
+        fix_sigma_alpha,
+        fix_sigma_beta,
+        fix_sigma_tau,
+        fix_mu,
+        fix_a,
+        fix_eta,
+    ):
         m = self.m
-        if (
-            x.ndim < 2
-            or (v0 is not None and len(v0) != 3)
-            or (mu0 is not None and len(mu0) != n)
-            or (a0 is not None and len(a0) != m)
-            or (eta0 is not None and len(eta0) != m)
-            or (fix_v and fix_mu and fix_a and fix_eta)
-        ):
-            with pytest.raises(ValueError):
-                _ = tdnoisefit(
-                    x.T,
-                    v0=v0,
-                    mu0=mu0,
-                    a0=a0,
-                    eta0=eta0,
-                    fix_v=fix_v,
-                    fix_mu=fix_mu,
-                    fix_a=fix_a,
-                    fix_eta=fix_eta,
-                )
-        else:
-            res = tdnoisefit(
+        n = self.n
+        dt = self.dt
+        sigma_alpha0 = self.alpha
+        sigma_beta0 = self.beta
+        sigma_tau0 = self.tau
+        with pytest.raises(ValueError):
+            _ = _parse_noisefit_input(
                 x.T,
-                v0=v0,
+                dt=dt,
+                sigma_alpha0=sigma_alpha0,
+                sigma_beta0=sigma_beta0,
+                sigma_tau0=sigma_tau0,
                 mu0=mu0,
                 a0=a0,
                 eta0=eta0,
-                fix_v=fix_v,
+                fix_sigma_alpha=fix_sigma_alpha,
+                fix_sigma_beta=fix_sigma_beta,
+                fix_sigma_tau=fix_sigma_tau,
                 fix_mu=fix_mu,
                 fix_a=fix_a,
                 fix_eta=fix_eta,
+                scale_sigma_alpha=1.0,
+                scale_sigma_beta=1.0,
+                scale_sigma_tau=1.0,
+                scale_delta_mu=np.ones(n),
+                scale_delta_a=np.ones(m - 1),
+                scale_eta=np.ones(m - 1),
             )
-            assert res.diagnostic["status"] == 0
+
+    @pytest.mark.parametrize("sigma_alpha0", [None, alpha, 0.0])
+    @pytest.mark.parametrize("sigma_beta0", [None, beta, 0.0])
+    @pytest.mark.parametrize("sigma_tau0", [None, tau, 0.0])
+    @pytest.mark.parametrize("mu0", [None, mu])
+    @pytest.mark.parametrize("a0", [None, a])
+    @pytest.mark.parametrize("eta0", [None, eta])
+    @pytest.mark.parametrize(
+        "scale_sigma_alpha, scale_sigma_beta, scale_sigma_tau, "
+        "scale_delta_mu, scale_delta_a, scale_eta",
+        [
+            [None, None, None, None, None, None],
+            [alpha, beta, tau / dt, noise_amp, scale_delta_a, scale_eta],
+        ],
+    )
+    def test_inputs(
+        self,
+        sigma_alpha0,
+        sigma_beta0,
+        sigma_tau0,
+        mu0,
+        a0,
+        eta0,
+        scale_sigma_alpha,
+        scale_sigma_beta,
+        scale_sigma_tau,
+        scale_delta_mu,
+        scale_delta_a,
+        scale_eta,
+    ):
+        x = self.x
+        dt = self.dt
+        _ = _parse_noisefit_input(
+            x.T,
+            dt=dt,
+            sigma_alpha0=sigma_alpha0,
+            sigma_beta0=sigma_beta0,
+            sigma_tau0=sigma_tau0,
+            mu0=mu0,
+            a0=a0,
+            eta0=eta0,
+            fix_sigma_alpha=False,
+            fix_sigma_beta=False,
+            fix_sigma_tau=False,
+            fix_mu=False,
+            fix_a=False,
+            fix_eta=False,
+            scale_sigma_alpha=scale_sigma_alpha,
+            scale_sigma_beta=scale_sigma_beta,
+            scale_sigma_tau=scale_sigma_tau,
+            scale_delta_mu=scale_delta_mu,
+            scale_delta_a=scale_delta_a,
+            scale_eta=scale_eta,
+        )
+
+    @pytest.mark.parametrize(
+        "fix_sigma_alpha, fix_sigma_beta, fix_sigma_tau, fix_mu, fix_a, "
+        "fix_eta",
+        [
+            [False, False, False, False, False, False],
+            [True, False, False, False, False, False],
+            [False, True, False, False, False, False],
+            [False, False, True, False, False, False],
+            [False, False, False, True, False, False],
+            [False, False, False, False, True, False],
+            [False, False, False, False, False, True],
+        ],
+    )
+    def test_noisefit(
+        self,
+        fix_sigma_alpha,
+        fix_sigma_beta,
+        fix_sigma_tau,
+        fix_mu,
+        fix_a,
+        fix_eta,
+    ):
+        x = self.x
+        result = noisefit(
+            x.T,
+            fix_sigma_alpha=fix_sigma_alpha,
+            fix_sigma_beta=fix_sigma_beta,
+            fix_sigma_tau=fix_sigma_tau,
+            fix_mu=fix_mu,
+            fix_a=fix_a,
+            fix_eta=fix_eta,
+        )
+        assert result.diagnostic["status"] == 0
 
 
 class TestFit:
-    thztools.global_options.sampling_time = None
     rng = np.random.default_rng(0)
     n = 16
     dt = 1.0 / n
@@ -511,7 +666,7 @@ class TestFit:
     alpha, beta, tau = 1e-5, 0, 0
     sigma = np.array([alpha, beta, tau])
     noise_model = NoiseModel(alpha, beta, tau, dt=dt)
-    noise_amp = noise_model.amplitude(mu)
+    noise_amp = noise_model.noise_amp(mu)
     x = mu + noise_amp * rng.standard_normal(n)
     y = psi + noise_amp * rng.standard_normal(n)