From 21ed4ecc5aade02f0d1c77fc20e188028ac3d9d1 Mon Sep 17 00:00:00 2001 From: Daniel Perera Date: Sat, 13 Jul 2019 17:02:07 -0500 Subject: [PATCH 1/2] issue 378 --- specutils/analysis/width.py | 49 ++++++++++++++++++++++++++++--------- 1 file changed, 37 insertions(+), 12 deletions(-) diff --git a/specutils/analysis/width.py b/specutils/analysis/width.py index cb6ca1021..9956a34f9 100644 --- a/specutils/analysis/width.py +++ b/specutils/analysis/width.py @@ -141,26 +141,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 + # acheive 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 + 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): From 2efba82d252fb007554d009daeb07638b2e56cde Mon Sep 17 00:00:00 2001 From: Daniel Perera Date: Sun, 14 Jul 2019 11:03:36 -0500 Subject: [PATCH 2/2] Sub-pixel precision for calculate_fwhm. --- docs/analysis.rst | 2 +- docs/fitting.rst | 2 +- specutils/analysis/width.py | 6 +++--- specutils/tests/test_analysis.py | 6 +++++- specutils/tests/test_fitting.py | 10 +++++----- 5 files changed, 15 insertions(+), 11 deletions(-) diff --git a/docs/analysis.rst b/docs/analysis.rst index e6fca59d9..2c03e305f 100644 --- a/docs/analysis.rst +++ b/docs/analysis.rst @@ -152,7 +152,7 @@ Each of the width analysis functions are applied to this spectrum below: >>> gaussian_fwhm(noisy_gaussian) # doctest: +FLOAT_CMP >>> fwhm(noisy_gaussian) # doctest: +FLOAT_CMP - + Reference/API diff --git a/docs/fitting.rst b/docs/fitting.rst index 7b0d369ba..56cac2bc3 100644 --- a/docs/fitting.rst +++ b/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 diff --git a/specutils/analysis/width.py b/specutils/analysis/width.py index 9956a34f9..bca44f18a 100644 --- a/specutils/analysis/width.py +++ b/specutils/analysis/width.py @@ -154,8 +154,8 @@ def _compute_single_fwhm(flux, spectral_axis): argmax = flux_value.argmax() halfval = flux_value[argmax] / 2 - left = flux_value[:argmax] <= halfval - right = flux_value[argmax + 1:] <= halfval + left = flux_value[:argmax] < halfval + right = flux_value[argmax + 1:] < halfval # Highest signal at the first point i0 = np.nonzero(left)[0] @@ -176,7 +176,7 @@ def _compute_single_fwhm(flux, spectral_axis): if i1.size == 0: right_value = spectral_value[-1] else: - i1 = i1[0] + argmax + i1 = i1[0] + argmax + 1 i0 = i1 - 1 left_flux = flux_value[i0] left_spectral = spectral_value[i0] diff --git a/specutils/tests/test_analysis.py b/specutils/tests/test_analysis.py index 943c9719b..3281f6080 100644 --- a/specutils/tests/test_analysis.py +++ b/specutils/tests/test_analysis.py @@ -483,7 +483,11 @@ def test_fwhm(): spectrum = Spectrum1D(spectral_axis=wavelengths, flux=flux) result = fwhm(spectrum) - assert result == 0.9090909090909092*u.um + # 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 # Highest point at the last point wavelengths = np.linspace(1, 10, 100) * u.um diff --git a/specutils/tests/test_fitting.py b/specutils/tests/test_fitting.py index 1cf700f95..abe3c7689 100644 --- a/specutils/tests/test_fitting.py +++ b/specutils/tests/test_fitting.py @@ -106,7 +106,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 @@ -120,7 +120,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 @@ -134,8 +134,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 @@ -158,7 +158,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