{ "cells": [ { "cell_type": "markdown", "id": "8c8138ef-e687-444e-95b1-9797654f3f62", "metadata": {}, "source": [ "# Frequency-switched Data Reduction\n", "\n", "This notebook shows how to use dysh to calibrate a frequency switched observation. \n", "The idea is similar to an OnOff observation, except the telescope does not move to an Off in position on the sky, but moves its intermediate frequency in frequency space. Here we call the On and Off the Sig and Ref, but both will have the signal, just shifted in the band. Since the telescope is always tracking the target, combining the Sig and a shifted (folded) Ref, a $\\sqrt{2}$ improvement in signal-to-noise can be achieved.\n", "\n", "The retrieval and calibration of frequency-switched observations uses `GBTFITSLoad.getfs`, which returns a `ScanBlock` object. \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", " sdf.select()\n", " sb = sdf.getfs()\n", " ta = sb.timeaverage()\n", " ta.baseline()\n", " ta.average()\n", " ta.plot()\n", "\n", "## Loading Modules\n", "We start by loading the modules we will use for this example. \n" ] }, { "cell_type": "code", "execution_count": 1, "id": "b4967550-2ca1-4931-b53b-6f9868718490", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "import" ] }, "outputs": [], "source": [ "# These modules are required for working with the data.\n", "from dysh.fits.gbtfitsload import GBTFITSLoad\n", "from dysh.log import init_logging\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": "5a180707-0731-448f-92dc-8635eff146ee", "metadata": {}, "source": [ "## Setup\n", "We start the dysh logging, so we get more information about what is happening.\n", "This is only needed if working on a notebook.\n", "If using the CLI through the dysh command, then logging is setup for you." ] }, { "cell_type": "code", "execution_count": 2, "id": "fe33a62b-e224-4615-875f-468a21cb327c", "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": "87669763-8d96-4521-9e81-f69d7213e133", "metadata": {}, "source": [ "## Data Retrieval\n", "\n", "Download the example SDFITS data, if necessary." ] }, { "cell_type": "code", "execution_count": 3, "id": "20611559-a142-4a8e-ae64-6690c93834b8", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "16:31:37.665 I Resolving example=getfs2 -> frequencyswitch/data/TREG_050627/TREG_050627.raw.acs/TREG_050627.raw.acs.fits\n" ] } ], "source": [ "filename = dysh_data(example=\"getfs2\")" ] }, { "cell_type": "markdown", "id": "8017d5a9-7e4d-47db-bd95-082fa56e40f4", "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": "93a62e3a-c95d-475b-8602-b5b8b7934733", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "load" ] }, "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", " \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
90W3OH0.0Track11.6673591.667359126122.228318.1145
91W3OH0.0Track11.6673591.667359126122.352118.2098
92W3OH0.0Track11.6673591.667359126122.473918.3043
93W3OH0.0Track11.6673591.667359126122.595318.3993
94W3OH0.0Track11.6673591.667359126122.716318.4949
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sdfits = GBTFITSLoad(filename)\n", "sdfits.summary()" ] }, { "cell_type": "markdown", "id": "80b45884-dd0c-41da-87af-60f4759fdc27", "metadata": {}, "source": [ "This data set contains 5 scans, 90 through 94. Each scan observed a single spectral window in two polarizations." ] }, { "cell_type": "markdown", "id": "0fef320a-3c1f-4a25-97fc-b8a870ae7205", "metadata": {}, "source": [ "## Data Reduction\n", "\n", "### Single Scan\n", "\n", "The default is to fold the Sig and Ref to create the final spectrum. Use `fold=False` to not fold them and `use_sig=False` to reverse the role of Sig and Ref. We will start calibrating `scan=90` (there are 5), `ifnum=0` (there is 1) and `plnum=1` (there are 2)." ] }, { "cell_type": "code", "execution_count": 5, "id": "defa4a6f-9da8-4488-846a-7dceccfa9446", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "getps", "test" ] }, "outputs": [], "source": [ "fs_scan_block = sdfits.getfs(scan=90, ifnum=0, plnum=1, fdnum=0)" ] }, { "cell_type": "markdown", "id": "85b26a7e-7e40-4877-88ae-e186272e429d", "metadata": {}, "source": [ "The return of `sdfits.getfs` is a `ScanBlock` which contains a single `FSScan`. The `FSScan` contains the calibrated data for all of the integrations. Now we will time average the integrations.\n", "\n", "#### Time Averaging\n", "To time average the contents of a `ScanBlock` use its `timeaverage` method. Be aware that time averging will not check if the source is the same. \n", "\n", "By default time averaging uses the following weights: \n", "$$\n", "\\frac{T^{2}_{sys}}{\\Delta\\nu\\Delta t}\n", "$$\n", "with $T_{sys}$ the system temperature, $\\Delta\\nu$ the channel width and $\\Delta t$ the integration time. In `dysh` these are set using `weights='tsys'` (the default)." ] }, { "cell_type": "code", "execution_count": 6, "id": "e0a33c19-576b-46e8-9a57-f34b8f85f6aa", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "timeaverage", "test" ] }, "outputs": [], "source": [ "ta = fs_scan_block.timeaverage(weights='tsys')" ] }, { "cell_type": "markdown", "id": "23229196-579b-4ca1-a2fa-e71e3a83a2b2", "metadata": {}, "source": [ "#### Plotting\n", "Plot the data and use different units for the spectral axis." ] }, { "cell_type": "code", "execution_count": 7, "id": "40747f0b-d63d-45fa-8910-fef56fe5c751", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "plot1", "test" ] }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "84a0c34922bc4e34b2fb138928fa82eb", "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": [ "ta.plot();" ] }, { "cell_type": "markdown", "id": "a94d57f3-2364-4779-9a66-452884ad8777", "metadata": {}, "source": [ "Now change the x-axis units to km s$^{-1}$ and restrict the range being shown." ] }, { "cell_type": "code", "execution_count": 8, "id": "6c07cf19-3a5d-47a9-9a10-ccd9f305efb9", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "plot2", "test" ] }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "9859bad55fc44e07a1d143e1c905a686", "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": [ "ta.plot(xaxis_unit=\"km/s\", yaxis_unit=\"K\", ymin=-5, ymax=5, xmin=-2000, xmax=2500);" ] }, { "cell_type": "markdown", "id": "dbde7ba7-0a34-49ac-bde7-ecabbb2264c3", "metadata": {}, "source": [ "#### Baseline Subtraction\n", "\n", "Next we remove a baseline from the data.\n", "\n", "We zoom in into the baseline to determine which model the use." ] }, { "cell_type": "code", "execution_count": 9, "id": "51fdd166-57d8-43b2-b292-61b706095a46", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "077ffcd7996341e78a35f4720f628534", "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": [ "ta.plot(xaxis_unit=\"chan\", yaxis_unit=\"K\", ymin=-2, ymax=2, title=\"before baseline subtraction\");" ] }, { "cell_type": "markdown", "id": "f2091fd1-2ce6-4cae-9854-4e5194e8e066", "metadata": {}, "source": [ "The baseline looks quasi-periodic, so a Chebyshev (`model='chebyshev'`) polynomial model may be a good model to use. Other available alternatives are: **hermite**, **legendre** and classic **polynomial**\n", "\n", "For baseline subtraction it is possible to specify the range of channels to be included in the fit (using the `include` argument) or excluded (using the `exclude` argument). Only one channel selection can be used at a time.\n", "\n", "We will also look at the mean, rms, min and max of the data before and after baseline subtraction." ] }, { "cell_type": "code", "execution_count": 10, "id": "603ae347-568a-41ba-a02f-7b6824c1af79", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "16:31:38.565 I EXCLUDING [Spectral Region, 1 sub-regions:\n", " (1642454441.0 Hz, 1643624790.1210938 Hz) \n", ", Spectral Region, 1 sub-regions:\n", " (1658883579.1835938 Hz, 1677194126.0585938 Hz) \n", ", Spectral Region, 1 sub-regions:\n", " (1686044223.7148438 Hz, 1692452915.1210938 Hz) \n", "]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Before baseline subtraction -- mean: 0.1370 K median: 0.1368 K rms: 0.140 K min: -0.95 K max: 0.56 K\n", "After baseline subtraction -- mean: -0.0089 K median: -0.0093 K rms: 0.145 K min: -1.13 K max: 0.45 K\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "84462444016f44199e118c92d5b5dc83", "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": [ "# Define a string.\n", "fmt_str = \"mean: {mean:.4f} median: {median:.4f} rms: {rms:.3f} min: {min:.2f} max: {max:.2f}\"\n", "# Print the statistics before baseline subtraction.\n", "print(f\"Before baseline subtraction -- {fmt_str}\".format(**ta[4200:5300].stats()))\n", "\n", "# Subtract the baseline.\n", "ta.baseline(model=\"chebyshev\", degree=5, include=[(4200,10000),(22000,32000)], remove=True)\n", "\n", "# Print the statistics after baseline subtraction.\n", "print(f\"After baseline subtraction -- {fmt_str}\".format(**ta[4200:5300].stats()))\n", "\n", "# Now plot the baseline subtracted spectrum.\n", "ta.plot(xaxis_unit=\"chan\", yaxis_unit=\"K\", ymin=-1, ymax=1, title=\"after baseline subtraction\");" ] }, { "cell_type": "markdown", "id": "55462dfb-cdcd-496c-9fa3-ebbe566dd888", "metadata": {}, "source": [ "Now we will undo this baseline fit, and plot it again to see if we get the same spectrum back. Also adding the statistics on the first section of the baseline to confirm the statistics are all the same as before we started fitting." ] }, { "cell_type": "code", "execution_count": 11, "id": "35423008-bd1a-450d-8686-adae3df6a880", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "967707e74db2410797889258331a1b68", "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" }, { "name": "stdout", "output_type": "stream", "text": [ "After undoing the baseline subtraction -- mean: 0.1370 K median: 0.1368 K rms: 0.140 K min: -0.95 K max: 0.56 K\n" ] } ], "source": [ "ta.undo_baseline()\n", "ta_plt = ta.plot(xaxis_unit=\"chan\", yaxis_unit=\"K\", ymin=-1, ymax=1, title='undo the baseline subtraction')\n", "print(f\"After undoing the baseline subtraction -- {fmt_str}\".format(**ta[4200:5300].stats()))" ] }, { "cell_type": "code", "execution_count": 12, "id": "3a7e4887-ba55-447a-9274-8275e9896d39", "metadata": {}, "outputs": [], "source": [ "ta_plt.savefig(output_dir / \"baselined_removed.png\")" ] }, { "cell_type": "markdown", "id": "835afcaa-bdb5-4d37-a65b-0d25d510b1a1", "metadata": {}, "source": [ "#### Using Selection\n", "\n", "We will repeat the calibration of `scan=90` using selection. To do this we pre-select the data using the `sdfits.select()` method." ] }, { "cell_type": "code", "execution_count": 13, "id": "dc3445e6-3b0a-47c5-8478-50453258a468", "metadata": {}, "outputs": [], "source": [ "sdfits.select(scan=90)" ] }, { "cell_type": "code", "execution_count": 14, "id": "7ca6e18f-ff84-4303-adad-59cc4614c0a1", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "a7da81d1b0064cf19f3f143853e4a7d3", "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" }, { "name": "stdout", "output_type": "stream", "text": [ "Using selection polarization 1 -- mean: 0.1370 K median: 0.1368 K rms: 0.140 K min: -0.95 K max: 0.56 K\n" ] } ], "source": [ "fs_scan_block2 = sdfits.getfs(plnum=1, ifnum=0, fdnum=0)\n", "ta2 = fs_scan_block2.timeaverage()\n", "ta2.plot(title='Time averaged plnum=1')\n", "print(f\"Using selection polarization 1 -- {fmt_str}\".format(**ta2[4200:5300].stats()))" ] }, { "cell_type": "markdown", "id": "54ecff3e-bffd-4236-808f-94dc5b754d71", "metadata": {}, "source": [ "#### Polarization Average\n", "\n", "Now we will calibrate the other polarization and average the two polarizations together.\n", "\n", "First we calibrate the second polarization, then time average it and inspect the time-averaged calibrated spectrum." ] }, { "cell_type": "code", "execution_count": 15, "id": "7be3d731-7da4-41da-ad9d-2f32281765b6", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "146fe61b9a3d4504b8a6e33ecfd99dac", "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" }, { "name": "stdout", "output_type": "stream", "text": [ "Using selection polarization 0 -- mean: 0.3070 K median: 0.3110 K rms: 0.134 K min: -0.45 K max: 0.79 K\n" ] } ], "source": [ "fs_scan_block3 = sdfits.getfs(plnum=0, ifnum=0,fdnum=0)\n", "ta3 = fs_scan_block3.timeaverage()\n", "ta3.plot(title='Time averaged plnum=0')\n", "print(f\"Using selection polarization 0 -- {fmt_str}\".format(**ta3[4200:5300].stats()))" ] }, { "cell_type": "markdown", "id": "2b2c3744-f312-41ca-9a90-8c0fda4266e9", "metadata": {}, "source": [ "Now we average the polarizations." ] }, { "cell_type": "code", "execution_count": 16, "id": "cb4c0c01-2f9d-4be2-954c-e75272dfd51f", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "88d6bfc5e0fe40b0a566d9381b851ce0", "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" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "d2ecc3787b114f7988645cd58324d2d6", "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" }, { "name": "stdout", "output_type": "stream", "text": [ "Polarization average -- mean: 0.2179 K median: 0.2208 K rms: 0.105 K min: -0.71 K max: 0.52 K\n" ] } ], "source": [ "# avg = 0.5*(ta2 + ta3) also works but the average() function does proper tsys weighting\n", "avg = ta2.average(ta3)\n", "avg.plot(ymin=-1,ymax=1, title='Averaged spectrum')\n", "avg[4200:5300].plot()\n", "print(f\"Polarization average -- {fmt_str}\".format(**avg[4200:5300].stats()))" ] }, { "cell_type": "markdown", "id": "f6473e81-a2ba-4125-8b1a-747e96a29da3", "metadata": {}, "source": [ "The noise is reduced by 1.34. Not exactly a factor of $\\sqrt{2}$, likely because of the baseline." ] }, { "cell_type": "markdown", "id": "f2fd1002-cf8a-486f-be9e-2f072ce219dd", "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." ] }, { "cell_type": "code", "execution_count": 17, "id": "cff1ba9c-4eff-413f-aa0b-d8660be37417", "metadata": { "scrolled": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "16:31:39.985 I rms is OK (no unit was given, assumed K)\n" ] } ], "source": [ "avg.check_stats(4.72584433) " ] }, { "cell_type": "code", "execution_count": 18, "id": "06160ebd-7ae6-430b-a58c-99d3e41eda72", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "np.float64(0.9421295947531141)" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "avg[23000:27000].radiometer() # 0.9542747726517753" ] } ], "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" }, "toc": { "base_numbering": 0 } }, "nbformat": 4, "nbformat_minor": 5 }