Skip to content

## 🐞 Bug: .set_spectral_resolution() not working as expected #36

Description

@brackham

Describe the issue
When using Spectrum.set_spectral_resolution() to convolve a high-resolution model spectrum to a lower resolving power, the output does not match the expected behavior:

  • It seems to apply the kernel incorrectly, resulting in an over-smoothed or incorrectly broadened spectrum.
  • When comparing to a manually implemented log(λ) convolution with scipy.ndimage.gaussian_filter1d, the manual method correctly preserves the resolving power across the wavelength range, while .set_spectral_resolution() does not.

Expected behavior

  • For a target constant resolving power R, the convolved spectrum should have a Gaussian LSF with FWHM = λ / R at each wavelength.
  • The broadening should appear consistent with how a real instrument works: constant Δλ/λ (constant in log(λ)).

Steps to reproduce

  1. Load a high-resolution Spectrum from grid:

    spec = Spectrum.from_grid(teff=4800, logg=4.5, feh=0.0)
  2. Call:

    spec_lowres = spec.set_spectral_resolution(500 * u.AA)  # or target_R
  3. Plot a sharp feature: compare the width to what’s expected for the desired R.

  4. Compare to this working log(λ) Gaussian convolution:

    import numpy as np
    from scipy.ndimage import gaussian_filter1d
    from speclib import Spectrum
    
    def convolve_spectrum_to_R(spec, target_R):
        """
        Convolve a speclib.Spectrum object to a target resolving power R.
        """
        wavelength = spec.wavelength.to("Angstrom").value
        flux = spec.flux.value
    
        log_lambda = np.log(wavelength)
        fwhm_log_lambda = 1.0 / target_R
        sigma_log_lambda = fwhm_log_lambda / (2.0 * np.sqrt(2.0 * np.log(2)))
        dlog_lambda = np.median(np.diff(log_lambda))
        sigma_pixels = sigma_log_lambda / dlog_lambda
    
        flux_smoothed = gaussian_filter1d(flux, sigma_pixels, mode="nearest")
    
        spec_smoothed = Spectrum(
            spectral_axis=spec.wavelength,
            flux=flux_smoothed * spec.flux.unit
        )
    
        return spec_smoothed

Metadata

Metadata

Assignees

Labels

bugSomething isn't working

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions