Custom Baseline Fitting

This notebook shows how to customize the baseline fitting method. dysh includes a baseline fitting method for its Spectrum objects, the baseline function. For simplicity, baseline uses astropy.modeling.polynomial as models. However, it is possible to fit any astropy.modeling.Model with some extra steps. Here we describe one possible approach.

Under the hood, the baseline method uses specutils.fitting.fit_continuum, which can take as input any astropy.modeling.Model. We can leverage this functionality to fit a custom model.

To illustrate this, we will create some fake data and then fit it. The notebook also shows another way of overplotting on a spectrum.

Dysh Commands

The following dysh commands are introduced (leaving out all the function arguments):

  sp = Spectrum.fake_spectrum()
  p = sp.plot()
  p.axis.plot()

Loading Modules

We start by loading the modules we will use to create the fake data and the baseline fit.

import numpy as np

from dysh.spectra import Spectrum
from dysh.log import init_logging
from specutils.fitting import fit_continuum
from astropy import modeling
from astropy import units as u

# These modules are used for file I/O
from dysh.util.files import dysh_data
from pathlib import Path

Setup

We start the dysh logging, so we get more information about what is happening. This is only needed if working on a notebook. If using the CLI through the dysh command, then logging is setup for you.

init_logging(2)

# also create a local "output" directory where temporary notebook files can be stored.
output_dir = Path.cwd() / "output"
output_dir.mkdir(exist_ok=True)

Make Fake Data

We create a Spectrum object. This could be the result of time averaging some calibrated data. Note for reasons of reproducibility, we fix the random number seed.

spectrum = Spectrum.fake_spectrum(seed=123)

Next, we add a power law to it, to represent a type of model that is not one of astropy.modeling.polynomial models. We use a power law with an exponent, alpha, of -12.7 and an amplitude of 10 K at the first channel of the spectrum, x_0.

amplitude = 10 * u.K
alpha = -12.7
x_0 = spectrum.spectral_axis.quantity[0]
power_law_values = amplitude * np.power(spectrum.spectral_axis.quantity / x_0, alpha)

We add the power law to our spectrum. This will be our starting point for the fit.

spectrum += power_law_values

See how it looks.

plot = spectrum.plot(xaxis_unit="GHz")

The power law is not obvious, as our spectrum only covers a small frequency range. But we know it is there, because the spectrum is offset from zero and it has a slope.

Custom Model Fitting

Now we define our model using astropy.modeling. When we define our model we must specify the starting guess for the model parameters. We will use starting values not too far from the answer.

model = modeling.powerlaws.PowerLaw1D(amplitude=5*u.K, 
                                      x_0=spectrum.spectral_axis.quantity[0],
                                      alpha=-10.)

Now, we use this model to fit our spectrum.

model_fit = fit_continuum(spectrum, model=model)

The return from fit_continuum is our model with the best fit parameters.

model_fit

<PowerLaw1D(amplitude=10.94015241 K, x_0=1.39456181e+09 Hz, alpha=13.3753621)>

Now, we evaluate the model over the same frequency range as our spectrum and plot it on top of the data. To add it to the plot we leverage the SpectrumPlot.show_baseline function that takes as input the y-axis values, and any of the arguments available to plot, except for gid. We call plot.show at the end of the cell to show the figure.

model_fit_values = model_fit(spectrum.spectral_axis)
plot.show_baseline(model_fit_values, label="Best fit model")
plot.show()

So, even if the best fit alpha was far from the input, this model still describes the fake data reasonably well.

Model Removal

We can remove this baseline model from the spectrum using arithmetic:

baseline_removed_spectrum = spectrum - model_fit_values
baseline_removed_spectrum.plot();

The baseline model has been removed from baseline_removed_spectrum, as the spectrum is now centered around zero.

Please check astropy.modeling for more models and additional options that can be used during the fitting.

Final Stats

Finally, at the end we compute some statistics over a spectrum, merely as a checksum if the notebook is reproducible.

baseline_removed_spectrum.check_stats(0.1008555 * u.K)

00:09:44.923 I rms is OK