{ "cells": [ { "cell_type": "markdown", "id": "9359fc63-3229-4fd3-8631-09c181e6a27b", "metadata": {}, "source": [ "# Aligning Spectra\n", "\n", "This notebook shows how to align spectra before they can be averaged.\n", "\n", "We use an observation with a mixed observing strategy; position switching and frequency switching.\n", "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.\n", "\n", "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.\n", "\n", "## Dysh commands\n", "\n", "The following dysh commands are introduced (leaving out all the function arguments):\n", "\n", " filename = dysh_data()\n", " sdf = GBTFITSLoad()\n", " sb = sdf.getps()\n", " ta = sb.timeaverage()\n", " ta1 = ta.align_to()\n", " ta.average()\n", " p1 = ta.plot()\n", " p1.oshow(ta1)\n", "\n", "## Loading Modules\n", "We start by loading the modules we will use for the data reduction. " ] }, { "cell_type": "code", "execution_count": 1, "id": "e34b86e9-6438-4e84-8503-d6dd1699e243", "metadata": {}, "outputs": [], "source": [ "# These modules are required for the data reduction.\n", "from dysh.fits.gbtfitsload import GBTFITSLoad\n", "from dysh.log import init_logging\n", "from astropy import units as u\n", "\n", "# This module is used for custom plotting.\n", "import matplotlib.pyplot as plt\n", "\n", "# These modules are used for file I/O\n", "from dysh.util.files import dysh_data\n", "from pathlib import Path" ] }, { "cell_type": "markdown", "id": "333aa096-b7dc-4812-a77c-88b07dcc116d", "metadata": {}, "source": [ "## Setup\n", "dysh uses a logger to communicate. If you are working in the command\n", "line, then the logging is setup for you. If you are working in a\n", "jupyter lab instance, then you need to set it up. You can do so using\n", "the init_logging function imported above. As an argument, init_logging\n", "takes a number, the verbosity level. level 0 is for error messages\n", "only, 1 for warning, 2 for info and 3 for debug. Here we set it to\n", "level 2.\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "4b612a07-7193-41dd-b7ed-5a4f695d0b50", "metadata": {}, "outputs": [], "source": [ "init_logging(2)\n", "\n", "# also create a local \"output\" directory where temporary notebook files can be stored.\n", "output_dir = Path.cwd() / \"output\"\n", "output_dir.mkdir(exist_ok=True)" ] }, { "cell_type": "markdown", "id": "74f02c88-5632-4f64-abda-6e70f63ff0fd", "metadata": {}, "source": [ "## Data Retrieval\n", "\n", "Download the example SDFITS data, if necessary." ] }, { "cell_type": "code", "execution_count": 3, "id": "cf801734-2ba0-48cd-a664-3ca15ce58890", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "06:24:09.652 I Resolving example=align -> mixed-fs-ps/data/TGBT24B_613_04.raw.vegas.trim.fits\n", "06:24:09.653 I url: http://www.gb.nrao.edu/dysh//example_data/mixed-fs-ps/data/TGBT24B_613_04.raw.vegas.trim.fits\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Odd-1, did not find /home/teuben/GBT/dysh_data/example_data/mixed-fs-ps/data/TGBT24B_613_04.raw.vegas.trim.fits\n", "TGBT24B_613_04.raw.vegas.trim.fits already downloaded\n" ] } ], "source": [ "filename = dysh_data(example=\"align\")" ] }, { "cell_type": "markdown", "id": "66d69520-58ac-46ac-be99-071c5812217c", "metadata": {}, "source": [ "## Data Loading\n", "\n", "Next, we use `GBTFITSLoad` to load the data, and then its `summary` method to inspect its contents." ] }, { "cell_type": "code", "execution_count": 4, "id": "fb10b18f-4510-4707-ad2a-70a1ec2b4d79", "metadata": {}, "outputs": [], "source": [ "sdf = GBTFITSLoad(filename)" ] }, { "cell_type": "code", "execution_count": 5, "id": "62f47e96-7553-48c9-9dfb-032d2af60897", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
SCANOBJECTVELOCITYPROCPROCSEQNRESTFREQDOPFREQ# IF# POL# INT# FEEDAZIMUTHELEVATION
35MESSIER32-200.0OnOff11.4204061.420406115173.115562.2955
36MESSIER32-200.0OnOff21.4204061.420406115172.308759.3225
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sdf.summary()" ] }, { "cell_type": "markdown", "id": "17ac95ab-15f2-4f45-80cd-cc9d6ca0f001", "metadata": {}, "source": [ "## Data Reduction\n", "\n", "### Position Switched Calibration\n", "\n", "We calibrate the signal and reference states independently." ] }, { "cell_type": "code", "execution_count": 6, "id": "d076291b-e69e-4018-ac95-718a7fa3c350", "metadata": {}, "outputs": [], "source": [ "# Signal.\n", "ps_sig = sdf.getps(scan=35, sig=\"T\", ifnum=0, plnum=0, fdnum=0).timeaverage()\n", "# Reference.\n", "ps_ref = sdf.getps(scan=35, sig=\"F\", ifnum=0, plnum=0, fdnum=0).timeaverage()" ] }, { "cell_type": "markdown", "id": "5d0d337f-68b5-4efb-a0f4-cc16eec85a74", "metadata": {}, "source": [ "Plot the signal and reference spectra as a function of channel.\n", "Here we use `matplotlib.pyplot` (`plt`) directly, since `dysh` is still not able to show more than one `Spectrum` in a plot.\n", "\n", "We now have the oshow()" ] }, { "cell_type": "code", "execution_count": 7, "id": "d5fdbfce-c703-44d8-b813-b4f6022be3b2", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "a4b5496d1e0847418004b53825d62409", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(HBox(children=(Button(description='Clear All Regions', style=ButtonStyle(), tooltip='Clear all …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p1 = ps_sig.plot(xaxis_unit=\"chan\")\n", "p1.oshow(ps_ref)" ] }, { "cell_type": "markdown", "id": "da57d514-95d2-4868-9f95-901387028983", "metadata": {}, "source": [ "```\n", "plt.plot(ps_sig.flux, label=\"sig\")\n", "plt.plot(ps_ref.flux, label=\"ref\")\n", "plt.legend() \\# Add a legend to display the labels.\n", "\\# Set the axis labels.\n", "plt.ylabel(f\"Antenna temperature ({ps_sig.flux.unit})\")\n", "plt.xlabel(\"Channel\")\n", "plt.show()\n", "```" ] }, { "cell_type": "markdown", "id": "13fb0551-7661-4d00-ab14-940b8547fbd1", "metadata": {}, "source": [ "Compute the rms in a line-free portion of the spectra." ] }, { "cell_type": "code", "execution_count": 8, "id": "037dc86e-05f4-4b44-9018-bfc15bc0949e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Noise on the signal state 0.41 K\n", "Noise on the reference state 0.40 K\n" ] } ], "source": [ "print(f\"Noise on the signal state {ps_sig[5000:10000].stats()['rms']:.2f}\")\n", "print(f\"Noise on the reference state {ps_ref[5000:10000].stats()['rms']:.2f}\")" ] }, { "cell_type": "markdown", "id": "1d1b09e3-969d-40db-89bc-f933756423e2", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 9, "id": "95b851e3-75cb-4b7b-b943-f9dbce5033ad", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "50fee8cc8f6e403eac136bc6f51b8470", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(HBox(children=(Button(description='Clear All Regions', style=ButtonStyle(), tooltip='Clear all …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ps = sdf.getps(scan=35, fdnum=0, plnum=0, ifnum=0).timeaverage()\n", "ps.plot(xaxis_unit=\"GHz\");" ] }, { "cell_type": "markdown", "id": "bc32e74f-de9c-4175-a60f-dc2014a4cac1", "metadata": {}, "source": [ "### Align Spectra\n", "\n", "Instead, we must align the spectra before averaging. In this case we will shift the reference state spectrum to match the signal state spectrum.\n", "\n", "We now have the oshow()" ] }, { "cell_type": "code", "execution_count": 10, "id": "b45ca858-fd11-4ce0-ae98-f8c835db81a4", "metadata": {}, "outputs": [], "source": [ "ps_ref_aligned = ps_ref.align_to(ps_sig)" ] }, { "cell_type": "markdown", "id": "9a70a61a-7e51-45e7-b017-0a2b213d5cb3", "metadata": {}, "source": [ "```\n", "plt.figure()\n", "plt.plot(ps_sig.flux)\n", "plt.plot(ps_ref_aligned.flux)\n", "plt.ylabel(f\"Antenna temperature ({ps_sig.flux.unit})\")\n", "plt.xlabel(\"Channel\")\n", "plt.show()\n", "```" ] }, { "cell_type": "code", "execution_count": 11, "id": "e9d6fcfc-e65d-4567-b1fe-0c2675c37944", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "5481b67796364357b01fc070c34ca3b6", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(HBox(children=(Button(description='Clear All Regions', style=ButtonStyle(), tooltip='Clear all …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p1 = ps_sig.plot(xaxis_unit=\"chan\")\n", "p1.oshow(ps_ref_aligned)" ] }, { "cell_type": "markdown", "id": "1961a051-3c7a-4435-ab0a-b608b06c4145", "metadata": {}, "source": [ "Now the spectra are aligned and can be averaged together." ] }, { "cell_type": "code", "execution_count": 12, "id": "cdf33158-62c0-4f66-a27d-f35fd22b12ea", "metadata": {}, "outputs": [], "source": [ "ps_avg = ps_sig.average(ps_ref_aligned)" ] }, { "cell_type": "code", "execution_count": 13, "id": "1b0c1152-0072-4bd9-9a11-0d5c5cd09fec", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "a3817fd390554aeb88cdf0089b518c72", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(HBox(children=(Button(description='Clear All Regions', style=ButtonStyle(), tooltip='Clear all …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ps_avg.plot(xaxis_unit=\"chan\");" ] }, { "cell_type": "code", "execution_count": 14, "id": "8ef254f2-4323-453b-825d-3ba8e9d1938e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Noise on the average 0.26 K\n" ] } ], "source": [ "print(f\"Noise on the average {ps_avg[5000:10000].stats()['rms']:.2f}\")" ] }, { "cell_type": "markdown", "id": "320a0e99-9042-44a9-b3f1-2a43bff985b8", "metadata": {}, "source": [ "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." ] }, { "cell_type": "markdown", "id": "2db3450b-0b8f-4869-94be-02fa4da20242", "metadata": {}, "source": [ "Finally, at the end we compute a fairly meaningless statistic over the whole spectrum, merely as a checksum if the notebook produces the same answer as when it was designed." ] }, { "cell_type": "markdown", "id": "1ff39aed-d620-4bfd-8d2f-b981d5248639", "metadata": {}, "source": [ "## Final Stats\n", "\n", "Finally, at the end we compute some statistics over a spectrum, merely as a checksum if the notebook is reproducible.\n" ] }, { "cell_type": "code", "execution_count": 15, "id": "722c2490-f653-412a-a200-c64adc30095b", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "06:24:11.233 I rms is OK \n" ] } ], "source": [ "ps_avg.check_stats(1.03823425 * u.K)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.9" } }, "nbformat": 4, "nbformat_minor": 5 }