Merge ef14dde into a5139a9

docs/analysis.rst

@@ -159,7 +159,7 @@ Each of the width analysis functions are applied to this spectrum below:
    >>> gaussian_fwhm(noisy_gaussian) # doctest: +FLOAT_CMP
    <Quantity 1.81144683 GHz>
    >>> fwhm(noisy_gaussian) # doctest: +FLOAT_CMP
-   <Quantity 1.90954774 GHz>
+   <Quantity 1.91686888 GHz>
    >>> fwzi(noisy_gaussian) # doctest: +FLOAT_CMP
    <Quantity 89.99999455 GHz>
 

docs/fitting.rst

@@ -158,7 +158,7 @@ Then estimate the line  parameters it it for a Gaussian line profile::
               amplitude            mean             stddev
                   Jy                um                um
           ------------------ ---------------- ------------------
-          1.1845669151078486 4.57517271067525 0.3015075376884422
+          1.1845669151078486 4.57517271067525 0.3199324375933905
 
 
 If an `~astropy.modeling.Model` is used that does not have the predefined
@@ -206,7 +206,7 @@ Model (Line) Fitting
 
 The generic model fitting machinery is well-suited to fitting spectral lines.
 The first step is to create a set of models with initial guesses as the
-parameters. To acheive better fits it may be wise to include a set of bounds for
+parameters. To achieve better fits it may be wise to include a set of bounds for
 each parameter, but that is optional.
 
 .. note::

specutils/analysis/width.py

@@ -213,26 +213,51 @@ def _compute_gaussian_sigma_width(spectrum, regions=None):
 
 
 def _compute_single_fwhm(flux, spectral_axis):
+    """
+    This is a helper function for the above `fwhm()` method.
+    """
 
-    argmax = np.argmax(flux)
-    halfval = flux[argmax] / 2
+    # The .value attribute is used here as the following algorithm does not
+    # use any array operations and would otherwise introduce a relatively
+    # significant overhead factor.  Two-point linear interpolation is used to
+    # achieve sub-pixel precision.
+    flux_value = flux.value
+    spectral_value = spectral_axis.value
 
-    left = flux[:argmax] <= halfval
-    right = flux[argmax+1:] <= halfval
+    argmax = flux_value.argmax()
+    halfval = flux_value[argmax] / 2
+    left = flux_value[:argmax] < halfval
+    right = flux_value[argmax + 1:] < halfval
 
     # Highest signal at the first point
-    if np.sum(left) == 0:
-        l_idx = 0
+    i0 = np.nonzero(left)[0]
+    if i0.size == 0:
+        left_value = spectral_value[0]
     else:
-        l_idx = np.where(left == True)[0][-1]
+        i0 = i0[-1]
+        i1 = i0 + 1
+        left_flux = flux_value[i0]
+        left_spectral = spectral_value[i0]
+        left_value = ((halfval - left_flux)
+                      * (spectral_value[i1] - left_spectral)
+                      / (flux_value[i1] - left_flux)
+                      + left_spectral)
 
     # Highest signal at the last point
-    if np.sum(right) == 0:
-        r_idx = len(flux)-1
+    i1 = np.nonzero(right)[0]
+    if i1.size == 0:
+        right_value = spectral_value[-1]
     else:
-        r_idx = np.where(right == True)[0][0] + argmax
-
-    return spectral_axis[r_idx] - spectral_axis[l_idx]
+        i1 = i1[0] + argmax + 1
+        i0 = i1 - 1
+        left_flux = flux_value[i0]
+        left_spectral = spectral_value[i0]
+        right_value = ((halfval - left_flux)
+                       * (spectral_value[i1] - left_spectral)
+                       / (flux_value[i1] - left_flux)
+                       + left_spectral)
+
+    return spectral_axis.unit * (right_value - left_value)
 
 
 def _compute_fwhm(spectrum, regions=None):

specutils/tests/test_analysis.py

@@ -476,14 +476,27 @@ def test_fwhm():
     expected = stddev * gaussian_sigma_to_fwhm
     assert quantity_allclose(result, expected, atol=0.01*u.GHz)
 
-
     # Highest point at the first point
     wavelengths = np.linspace(1, 10, 100) * u.um
-    flux = (1.0 / wavelengths.value)*u.Jy # highest point first.
+    flux = (1.0 / wavelengths.value) * u.Jy  # highest point first.
+
+    spectrum = Spectrum1D(spectral_axis=wavelengths, flux=flux)
+    result = fwhm(spectrum)
+    # Note that this makes a little more sense than the previous version;
+    # since the maximum value occurs at wavelength=1, and the half-value of
+    # flux (0.5) occurs at exactly wavelength=2, the result should be
+    # exactly 1 (2 - 1).
+    assert result == 1.0 * u.um
+
+    # Test the interpolation used in FWHM for wavelength values that are not
+    # on the grid
+    wavelengths = np.linspace(1, 10, 31) * u.um
+    flux = (1.0 / wavelengths.value) * u.Jy  # highest point first.
 
     spectrum = Spectrum1D(spectral_axis=wavelengths, flux=flux)
     result = fwhm(spectrum)
-    assert result == 0.9090909090909092*u.um
+
+    assert quantity_allclose(result, 1.01 * u.um)
 
     # Highest point at the last point
     wavelengths = np.linspace(1, 10, 100) * u.um

specutils/tests/test_fitting.py

@@ -107,7 +107,7 @@ def test_single_peak_estimate():
 
     assert np.isclose(g_init.amplitude.value, 3.354169257846847)
     assert np.isclose(g_init.mean.value, 6.218588636687762)
-    assert np.isclose(g_init.stddev.value, 1.608040201005025)
+    assert np.isclose(g_init.stddev.value, 1.6339001193853715)
 
     assert g_init.amplitude.unit == u.Jy
     assert g_init.mean.unit == u.um
@@ -121,7 +121,7 @@ def test_single_peak_estimate():
 
     assert np.isclose(g_init.amplitude.value, 3.354169257846847)
     assert np.isclose(g_init.x_0.value, 6.218588636687762)
-    assert np.isclose(g_init.fwhm.value, 1.608040201005025)
+    assert np.isclose(g_init.fwhm.value, 1.6339001193853715)
 
     assert g_init.amplitude.unit == u.Jy
     assert g_init.x_0.unit == u.um
@@ -135,8 +135,8 @@ def test_single_peak_estimate():
 
     assert np.isclose(g_init.amplitude_L.value, 3.354169257846847)
     assert np.isclose(g_init.x_0.value, 6.218588636687762)
-    assert np.isclose(g_init.fwhm_L.value, 1.1370561305512321)
-    assert np.isclose(g_init.fwhm_G.value, 1.1370561305512321)
+    assert np.isclose(g_init.fwhm_L.value, 1.1553418541989058)
+    assert np.isclose(g_init.fwhm_G.value, 1.1553418541989058)
 
     assert g_init.amplitude_L.unit == u.Jy
     assert g_init.x_0.unit == u.um
@@ -159,7 +159,7 @@ def test_single_peak_estimate():
 
     assert np.isclose(g_init.amplitude.value, 3.354169257846847)
     assert np.isclose(g_init.x_0.value, 6.218588636687762)
-    assert np.isclose(g_init.stddev.value, 1.608040201005025)
+    assert np.isclose(g_init.stddev.value, 1.6339001193853715)
 
     assert g_init.amplitude.unit == u.Jy
     assert g_init.x_0.unit == u.um