Merge 157b63c into b5720f6

docs/analysis.rst

@@ -135,24 +135,33 @@ computing a second-moment-based approximation of the standard deviation.
 The `~gaussian_fwhm` function estimates the width of the spectrum at half max,
 again by computing an approximation of the standard deviation.
 
-Both of these functions assume that spectrum is approximately gaussian.
+Both of these functions assume that the spectrum is approximately gaussian.
 
 The function `~fwhm` provides an estimate of the full width of the spectrum at
 half max that does not assume the spectrum is gaussian. It locates the maximum,
 and then locates the value closest to half of the maximum on either side, and
 measures the distance between them.
 
+A function to calculate the full width at zero intensity (i.e. the width of a
+spectral feature at the continuum) is provided as `~fwzi`. Like the `~fwhm`
+calculation, it does not make assumptions about the shape of the feature
+and calculates the width by finding the points at either side of maximum
+that reach the continuum value. In this case, it assumes the provided
+spectrum has been continuum subtracted.
+
 Each of the width analysis functions are applied to this spectrum below:
 
 .. code-block:: python
 
-   >>> from specutils.analysis import gaussian_sigma_width, gaussian_fwhm, fwhm
+   >>> from specutils.analysis import gaussian_sigma_width, gaussian_fwhm, fwhm, fwzi
    >>> gaussian_sigma_width(noisy_gaussian) # doctest: +FLOAT_CMP
    <Quantity 0.76925064 GHz>
    >>> gaussian_fwhm(noisy_gaussian) # doctest: +FLOAT_CMP
    <Quantity 1.81144683 GHz>
    >>> fwhm(noisy_gaussian) # doctest: +FLOAT_CMP
    <Quantity 1.90954774 GHz>
+   >>> fwzi(noisy_gaussian) # doctest: +FLOAT_CMP
+   <Quantity 89.99999455 GHz>
 
 
 Reference/API

specutils/analysis/width.py

@@ -5,12 +5,14 @@
 
 import numpy as np
 from astropy.stats.funcs import gaussian_sigma_to_fwhm
+from astropy.modeling.models import Gaussian1D
 from ..manipulation import extract_region
 from . import centroid
 from .utils import computation_wrapper
+from scipy.signal import chirp, find_peaks, peak_widths
 
 
-__all__ = ['gaussian_sigma_width', 'gaussian_fwhm', 'fwhm']
+__all__ = ['gaussian_sigma_width', 'gaussian_fwhm', 'fwhm', 'fwzi']
 
 
 def gaussian_sigma_width(spectrum, regions=None):
@@ -105,6 +107,76 @@ def fwhm(spectrum, regions=None):
     return computation_wrapper(_compute_fwhm, spectrum, regions)
 
 
+def fwzi(spectrum, regions=None):
+    """
+    Compute the true full width at zero intensity (i.e. the continuum level)
+    of the spectrum.
+
+    This makes no assumptions about the shape of the spectrum. It uses the
+    scipy peak-finding to determine the index of the highest flux value, and
+    then calculates width at that base of the feature.
+
+    Parameters
+    ----------
+    spectrum : `~specutils.spectra.spectrum1d.Spectrum1D`
+        The spectrum object over which the width will be calculated.
+
+    regions: `~specutils.utils.SpectralRegion` or list of `~specutils.utils.SpectralRegion`
+        Region within the spectrum to calculate the FWZI value. If regions is
+        `None`, computation is performed over entire spectrum.
+
+    Returns
+    -------
+    `~astropy.units.Quantity` or list (based on region input)
+        Full width of the signal at zero intensity.
+
+    Notes
+    -----
+    The spectrum must be continuum subtracted before being passed to this
+    function.
+    """
+    return computation_wrapper(_compute_fwzi, spectrum, regions)
+
+
+def _compute_fwzi(spectrum, regions=None):
+    if regions is not None:
+        calc_spectrum = extract_region(spectrum, regions)
+    else:
+        calc_spectrum = spectrum
+
+    # Create a copy of the flux array to ensure the value on the spectrum
+    # object is not altered.
+    disp, flux = calc_spectrum.spectral_axis, calc_spectrum.flux.copy()
+
+    # For noisy data, ensure that the search from the centroid stops on
+    # either side once the flux value reaches zero.
+    flux[flux < 0] = 0
+
+    def find_width(data):
+        # Find the peaks in the flux data
+        peaks, _ = find_peaks(data)
+
+        # Find the index of the maximum peak value in the found peak list
+        peak_ind = [peaks[np.argmin(np.abs(
+            np.array(peaks) - np.argmin(np.abs(data - np.max(data)))))]]
+
+        # Calculate the width for the given feature
+        widths, _, _, _ = \
+            peak_widths(data, peak_ind, rel_height=1-1e-7)
+
+        return widths[0] * disp.unit
+
+    if flux.ndim > 1:
+        tot_widths = []
+
+        for i in range(flux.shape[0]):
+            tot_widths.append(find_width(flux[i]))
+
+        return tot_widths
+
+    return find_width(flux)
+
+
 def _compute_gaussian_fwhm(spectrum, regions=None):
     """
     This is a helper function for the above `gaussian_fwhm()` method.

specutils/tests/test_analysis.py

@@ -10,7 +10,7 @@
 from ..spectra import Spectrum1D, SpectralRegion
 from ..analysis import (line_flux, equivalent_width, snr, centroid,
                         gaussian_sigma_width, gaussian_fwhm, fwhm,
-                        snr_derived)
+                        snr_derived, fwzi)
 from ..tests.spectral_examples import simulated_spectra
 
 
@@ -524,3 +524,44 @@ def test_fwhm_multi_spectrum():
 
     expected = stddevs * gaussian_sigma_to_fwhm
     assert quantity_allclose(results, expected, atol=0.01*u.GHz)
+
+
+def test_fwzi():
+    np.random.seed(42)
+
+    disp = np.linspace(0, 100, 1000) * u.AA
+
+    g = models.Gaussian1D(mean=np.mean(disp),
+                          amplitude=1 * u.Jy,
+                          stddev=2 * u.AA)
+
+    flux = g(disp)
+    flux_with_noise = g(disp) + ((np.random.sample(disp.size) - 0.5) * 0.1) * u.Jy
+
+    spec = Spectrum1D(spectral_axis=disp, flux=flux)
+    spec_with_noise = Spectrum1D(spectral_axis=disp, flux=flux_with_noise)
+
+    assert quantity_allclose(fwzi(spec), 226.89732509 * u.AA)
+    assert quantity_allclose(fwzi(spec_with_noise), 106.99998944 * u.AA)
+
+
+def test_fwzi_multi_spectrum():
+    np.random.seed(42)
+
+    disp = np.linspace(0, 100, 1000) * u.AA
+
+    amplitudes = [0.1, 1, 10] * u.Jy
+    means = [25, 50, 75] * u.AA
+    stddevs = [1, 5, 10] * u.AA
+    params = list(zip(amplitudes, means, stddevs))
+
+    flux = np.zeros(shape=(3, len(disp)))
+
+    for i in range(3):
+        flux[i] = models.Gaussian1D(*params[i])(disp)
+
+    spec = Spectrum1D(spectral_axis=disp, flux=flux * u.Jy)
+
+    expected = [113.51706001 * u.AA, 567.21252727 * u.AA, 499.5024546 * u.AA]
+
+    assert quantity_allclose(fwzi(spec), expected)