Aligning Spectra#

This guide shows how to align spectra before they can be averaged.

We use an observation with a mixed observing strategy; position switching and frequency switching. In this case the on and off source observations also switch the signal in the frequency domain between a signal and a reference state, so there are four switching states: signal with the noise diode, signal without the noise diode, reference with the noise diode and reference without the noise diode.

We will calibrate the signal and reference states independently, using position switching. Then, we’ll look at the calibrated spectra as a function of channel and vizualize the shift between the signal and reference states. Finally we will shift the reference state and average it with the signal state.

from pathlib import Path
import matplotlib.pyplot as plt
from dysh.util.download import from_url
from dysh.fits.gbtfitsload import GBTFITSLoad
from dysh.spectra import average_spectra

Data Retrieval#

Download the example SDFITS data, if necessary.

url = "http://www.gb.nrao.edu/dysh/example_data/mixed-fs-ps/data/TGBT24B_613_04.raw.vegas.trim.fits"
savepath = Path.cwd() / "data"
filename = from_url(url, savepath)


Data Loading#

Next, we use GBTFITSLoad to load the data, and then its summary method to inspect its contents.

sdf = GBTFITSLoad(filename)
sdf.summary()
SCAN OBJECT VELOCITY PROC PROCSEQN RESTFREQ DOPFREQ # IF # POL # INT # FEED AZIMUTH ELEVATIO
0 35 MESSIER32 -200.0 OnOff 1 1.420406 1.420406 1 1 5 1 73.115543 62.295527
1 36 MESSIER32 -200.0 OnOff 2 1.420406 1.420406 1 1 5 1 72.308676 59.322465

Data Reduction#

Position Switched Calibration#

We calibrate the signal and reference states independently.

# Signal.
ps_sig = sdf.getps(scan=35, sig="T").timeaverage()
# Reference.
ps_ref = sdf.getps(scan=35, sig="F").timeaverage()

Plot the signal and reference spectra as a function of channel.

plt.figure()
plt.plot(ps_sig.flux)
plt.plot(ps_ref.flux)
plt.ylabel(f"Antenna temperature ({ps_sig.flux.unit})")
plt.xlabel("Channel");
../../_images/1d3fecc63ace480849ed36ee5d24bf7a26134afe966b07dee67d1ec1b916cec1.png
print(f"Noise on the signal state {ps_sig[5000:10000].stats()['rms']:.2f}")
print(f"Noise on the reference state {ps_ref[5000:10000].stats()['rms']:.2f}")
Noise on the signal state 0.41 K
Noise on the reference state 0.40 K

If we were to average the data, then the result would have the spectral line misaligned, since the averaging is done per channel. We can see this if we calibrate the data without separating the signal and reference states. This is less than ideal because the signal gets reduced and the line shows up in two different places even though it is the same line.

ps = sdf.getps(scan=35).timeaverage()
ps.plot(xaxis_unit="GHz")
../../_images/24da3761d6144968798a27f9ad61c5bda34393fe272caae4cf71126700663fa1.png

Align Spectra#

Instead, we must align the spectra before averaging. In this case we will shift the reference state spectrum to match the signal state spectrum.

ps_ref_aligned = ps_ref.align_to(ps_sig)

plt.figure()
plt.plot(ps_sig.flux)
plt.plot(ps_ref_aligned.flux)
plt.ylabel(f"Antenna temperature ({ps_sig.flux.unit})")
plt.xlabel("Channel");
../../_images/67675f4e930bfd951889520e5f4ab9546b2ad5aacc5d756c64a25778272d064c.png

Now the spectra are aligned and can be averaged together.

ps_avg = average_spectra((ps_sig, ps_ref_aligned))
ps_avg.plot(xaxis_unit="chan")
../../_images/7151f8e286f36f43d37fcf605a193d7f0fe6630683809d9bb35b74222aff5c95.png
print(f"Noise on the average {ps_avg[5000:10000].stats()['rms']:.2f}")
Noise on the average 0.26 K

The resulting spectrum shows only one peak, at the frequency of the line, and the noise is smaller than in the individual states where the spectra overlapped (channels below ~18000). The noise improvement is better than \(\sqrt{2}\) since shifting the reference spectrum for the alignment artificially lowers its noise.