Aligning Spectra

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.

You can find a copy of this tutorial as a Jupyter notebook here or download it by right clicking here and selecting “Save Link As”.

Loading Modules#

We start by loading the modules we will use for the data reduction.

For display purposes, we use the static (non-interactive) matplotlib backend in this tutorial. However, you can tell matplotlib to use the ipympl backend to enable interactive plots. This is only needed if working on jupyter lab or notebook.

# Set interactive plots in jupyter.
#%matplotlib ipympl

# These modules are required for the data reduction.
from dysh.fits.gbtfitsload import GBTFITSLoad

# This module is used for custom plotting.
import matplotlib.pyplot as plt

# These modules are only used to download the data.
from pathlib import Path
from dysh.util.download import from_url

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"
savepath.mkdir(exist_ok=True) # Create the data directory if it does not exist.
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
      ELEVATION
    

  
  
    

      35
      MESSIER32
      -200.0
      OnOff
      1
      1.420406
      1.420406
      1
      1
      5
      1
      73.1155
      62.2955
    

    

      36
      MESSIER32
      -200.0
      OnOff
      2
      1.420406
      1.420406
      1
      1
      5
      1
      72.3087
      59.3225
    

  









Data Reduction#




Position Switched Calibration#


We calibrate the signal and reference states independently.






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









Plot the signal and reference spectra as a function of channel.
Here we use matplotlib.pyplot (plt) directly, since dysh is still not able to show more than one Spectrum in a plot.






plt.figure()
plt.plot(ps_sig.flux, label="sig")
plt.plot(ps_ref.flux, label="ref")
plt.legend() # Add a legend to display the labels.
# Set the axis labels.
plt.ylabel(f"Antenna temperature ({ps_sig.flux.unit})")
plt.xlabel("Channel");








../../_images/87a789b633143ca29cc4ba973a4c3dc766a260f37e324cba5a3f32dba5cf6ff4.png





Compute the rms in a line-free portion of the spectra.






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, fdnum=0, plnum=0, ifnum=0).timeaverage()
ps.plot(xaxis_unit="GHz")









<dysh.plot.specplot.SpectrumPlot at 0x712131da3700>




../../_images/ed98d8454ee12f7a3ae2b08e5762e38a1ad47e723d450eaaead338f74172907e.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/118d6fdef5cc4e3578027dab28440ab6eeed927696ebf00b00eb2e9b1062bd55.png





Now the spectra are aligned and can be averaged together.






ps_avg = ps_sig.average(ps_ref_aligned)













ps_avg.plot(xaxis_unit="chan")









<dysh.plot.specplot.SpectrumPlot at 0x712133f231f0>




../../_images/d821e4772c306198621fff8ad83adcb6f9411bd5f5e48c0c771c4075dc2a9202.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.









                

              

              
              
              
              
                
              
            

            
            
              
                



  

  

     Contents