{ "cells": [ { "cell_type": "markdown", "id": "8c8138ef-e687-444e-95b1-9797654f3f62", "metadata": {}, "source": [ "# Smoothing\n", "\n", "This notebook shows how to use dysh to smooth a spectrum. \n", "\n", "For the example below we will use data from the Position-Switch example again.\n", "The following dysh commands are the simplest to get and smooth a spectrum (leaving out all the function arguments):\n", "\n", " sdf = GBTFITSLoad()\n", " sb = sdf.getps()\n", " ta = sb.timeaverage()\n", " tb = ta.smooth()\n", " tb.plot()\n", "\n", "or if you wish to make use of the Python object chaining you can do this in one line:\n", "\n", " GBTFITSLoad().getps().timeaverage().smooth().plot()\n", "\n", "## Loading Modules\n", "We start by loading the modules we will use for the data reduction. \n", "\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", "from astropy import units as u\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": "bea96434-ab5d-4c13-9f4e-c20258cad627", "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": "2fb86f56-0d81-4347-916f-f1e3a7893c06", "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": "6bc88bc5-986d-4eae-b1c7-6398cc9ddd5a", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "wget" ] }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "11:40:07.485 I Resolving test=getps -> AGBT05B_047_01/AGBT05B_047_01.raw.acs/\n" ] } ], "source": [ "filename = dysh_data(test=\"getps\")" ] }, { "cell_type": "markdown", "id": "01c33131-69a2-4021-b184-6f80bb31bb20", "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": [], "source": [ "sdfits = GBTFITSLoad(filename)" ] }, { "cell_type": "code", "execution_count": 5, "id": "6bf49756-6af2-41ff-8bc0-01d6a973837a", "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", " \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
51NGC52914386.0OnOff11.4204051.42040512111198.343118.6427
52NGC52914386.0OnOff21.4204051.42040512111198.930618.7872
53NGC52914386.0OnOff11.4204051.42040512111199.330518.3561
54NGC52914386.0OnOff21.4204051.42040512111199.915718.4927
55NGC52914386.0OnOff11.4204051.42040512111200.304218.0575
56NGC52914386.0OnOff21.4204051.42040512111200.890618.1860
57NGC52914386.0OnOff11.4204051.42040512111202.327517.3853
58NGC52914386.0OnOff21.4204051.42040512111202.919217.4949
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sdfits.summary()" ] }, { "cell_type": "markdown", "id": "86227d72-cdff-42d0-b98f-b4dff9fe13dd", "metadata": {}, "source": [ "## Data Reduction\n", "\n", "### Calibration and Time Averaging\n", "\n", "This test data has 8 scans making up 4 pairs of OnOff observations. We will calibrate all of the OnOff scan pairs, time average them and then compare the results with and without smoothing. To calibrate all of the scans, we can omit the `scan=` keyword when using the `sdfits.getps` method.\n", "\n", "*Technical note*: ``getps`` returns a **ScanBlock** with in this case four `PSScan`s, since there are four pairs of OnOff observations.\n" ] }, { "cell_type": "markdown", "id": "df12d3c8-ec1b-4ea2-acfb-729a06d0ce68", "metadata": {}, "source": [ "The noise in the time averaged spectrum is 36 mK. However, this has not been baseline subtracted and there is a slope in the continuum.\n", "In addition, the edges look horrible, and there is a strong \"spike\" (galactic) emission near channel 3000. We like to avoid using those data for the baseline fit." ] }, { "cell_type": "code", "execution_count": 6, "id": "e0a33c19-576b-46e8-9a57-f34b8f85f6aa", "metadata": { "editable": true, "scrolled": true, "slideshow": { "slide_type": "" }, "tags": [ "timeaverage", "test" ] }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "f6e7967739834b5b8cd12d000d8ee09e", "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": [ "Stats : mean: 0.2708 K median: 0.2711 K rms: 0.0359 K min: 0.14 K max: 0.42 K\n", "Expect: mean: 0.2708 K median: 0.2711 K rms: 0.0359 K min: 0.14 K max: 0.42 K\n", "Stats : mean: 0.2916 K median: 0.2916 K rms: 0.0366 K min: 0.14 K max: 0.41 K\n", "Expect: mean: 0.2916 K median: 0.2916 K rms: 0.0366 K min: 0.14 K max: 0.41 K\n" ] } ], "source": [ "scan_block = sdfits.getps(ifnum=0, plnum=0, fdnum=0)\n", "ta = scan_block.timeaverage()\n", "ta.plot(xaxis_unit=\"chan\", yaxis_unit=\"mK\", ymin=100, ymax=600, grid=True)\n", "\n", "# Define a string for printing the spectrum statistics.\n", "fmt_str = \"mean: {mean:.4f} median: {median:.4f} rms: {rms:.4f} min: {min:.2f} max: {max:.2f}\"\n", "print(f\"Stats : {fmt_str}\".format(**ta[5000:14000].stats()))\n", "print(\"Expect: mean: 0.2708 K median: 0.2711 K rms: 0.0359 K min: 0.14 K max: 0.42 K\")\n", "print(f\"Stats : {fmt_str}\".format(**ta[20000:30000].stats()))\n", "print(\"Expect: mean: 0.2916 K median: 0.2916 K rms: 0.0366 K min: 0.14 K max: 0.41 K\")" ] }, { "cell_type": "markdown", "id": "d19573c5-ef7e-4395-a88d-91c3587b1c96", "metadata": {}, "source": [ "### Baseline Removal\n", "\n", "We remove the slope by baseline fitting an order 1 polynomial from the data while ignoring the channels that contain a signal. This is done using the `Spectrum.baseline()` method." ] }, { "cell_type": "code", "execution_count": 7, "id": "652e5500-dfda-420d-8ca3-eed516967415", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "11:40:09.183 I EXCLUDING [Spectral Region, 1 sub-regions:\n", " (1374818364.0 Hz, 1379040470.9335938 Hz) \n", ", Spectral Region, 1 sub-regions:\n", " (1397351017.8085938 Hz, 1403454533.4335938 Hz) \n", ", Spectral Region, 1 sub-regions:\n", " (1419476261.9492188 Hz, 1424816838.1210938 Hz) \n", "]\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "ac31088a32f64122ae70abcc52f1787e", "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": [ "Stats : mean: -0.0007 K median: -0.0007 K rms: 0.0358 K min: -0.13 K max: 0.15 K\n", "Expect: mean: -0.0007 K median: -0.0007 K rms: 0.0358 K min: -0.13 K max: 0.15 K\n" ] } ], "source": [ "ta.baseline(model=\"poly\", degree=1, include=[[3500,14000],[18000,30000]], remove=True)\n", "ta.plot(xaxis_unit=\"chan\", yaxis_unit=\"mK\", ymin=-200, ymax=300, grid=True)\n", "print(f\"Stats : {fmt_str}\".format(**ta[5000:14000].stats()))\n", "print(\"Expect: mean: -0.0007 K median: -0.0007 K rms: 0.0358 K min: -0.13 K max: 0.15 K\")" ] }, { "cell_type": "markdown", "id": "571a36b9-57b0-4448-b7a0-4aea4e09edb3", "metadata": {}, "source": [ "### Smoothing in a Few Ways\n", "\n", "By default smoothing will also decimate the signal, to (roughly) make each channel independant from the next. This assuming the input signal had independant channels. If your input was oversampled by a factor of 2, the smoothed signal will be as well, although you can manually decimate by a different value too, for example by using ``decimate=8`` .\n", "\n", "Since we smooth to a gauss of FWHM 16 channels, the noise should go down by a factor of 4 (36 mK to 9 mK)." ] }, { "cell_type": "code", "execution_count": 8, "id": "54db3d01-d3a8-4fad-bc2a-6138e1cb7d58", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "2c685df8ae264d5ab2e6d7dbe7be4852", "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": [ "Stats : mean: -0.0007 K median: -0.0007 K rms: 0.0075 K min: -0.02 K max: 0.02 K\n", "Expect: mean: 0.0005 K median: 0.0004 K rms: 0.0075 K min: -0.02 K max: 0.03 K\n" ] } ], "source": [ "ts1 = ta.smooth('gaussian', 16)\n", "ts1.plot(xaxis_unit=\"chan\", yaxis_unit=\"mK\", ymin=-200, ymax=300, grid=True)\n", "print(f\"Stats : {fmt_str}\".format(**ts1[5000//16:14000//16].stats()))\n", "print(\"Expect: mean: 0.0005 K median: 0.0004 K rms: 0.0075 K min: -0.02 K max: 0.03 K\")" ] }, { "cell_type": "markdown", "id": "224c299c-e9b4-48e7-9089-468c11f0670d", "metadata": {}, "source": [ "The noise is now 8 mK, very close to the improvement we expected.\n", "\n", "Now smoothing by 320 channels should result in a noise of $55/\\sqrt{320}$ or 2 mK. In this case we will use a box as the convolving kernel.\n", "\n", "We will also estimate the correlation between channels by comparing the noise in the spectra with the noise in the spectra when we subtract adjacent channels (using `roll=1` in `Spectrum.stats`)." ] }, { "cell_type": "code", "execution_count": 9, "id": "ed472d60-d23a-41c4-aa2d-a3044da52c22", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "91a52327fbb349b68779b04f2c800c01", "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": [ "Stats : mean: -0.0007 K median: -0.0010 K rms: 0.0023 K min: -0.01 K max: 0.00 K\n", "Rolled rms/rms ratio: 0.7845596461442604\n" ] } ], "source": [ "ts2 = ta.smooth('box', 320)\n", "ts2.plot(xaxis_unit=\"chan\", yaxis_unit=\"mK\", ymin=-200, ymax=300, grid=True)\n", "\n", "print(f\"Stats : {fmt_str}\".format(**ts2[5000//320:14000//320].stats()))\n", "rratio = ts2[5000//320:14000//320].stats(roll=1)[\"rms\"]/ts2[5000//320:14000//320].stats(roll=0)[\"rms\"]\n", "print(f\"Rolled rms/rms ratio: {rratio}\")" ] }, { "cell_type": "markdown", "id": "11a1f74b-ab1b-4cdc-9aca-8aa86f637b55", "metadata": {}, "source": [ "The noise in the smoothed spectrum is ~2 mK, as expected.\n", "\n", "The rolled RMS ratio is very close to 1, so neighboring channels are not related. If you would decimate by 160, you would see this ratio drop. Be sure to adjust the range of channels for any new computations of the spectrum statistics (`Spectrum.stats`)." ] }, { "cell_type": "markdown", "id": "15afb100-99da-4eca-aed1-1b79f3987930", "metadata": {}, "source": [ "### Smoothing the Reference (\"OFF\") Scans\n", "\n", "Under certain circumstances it can be beneficial to (boxcar) smooth the reference (OFF) signal before the usual\n", "(ON-OFF)/OFF calibration. \n", "\n", "*Technical note*: if you want to achieve identical results to GBTIDL, the width of the boxcar needs to be odd.\n" ] }, { "cell_type": "code", "execution_count": 10, "id": "ea6820db-4bd6-4d47-bc04-12d77578651d", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "836149dd9f504a98a100d58af36d4fc3", "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": [ "Stats : mean: 0.2667 K median: 0.2666 K rms: 0.0254 K min: 0.18 K max: 0.35 K\n", "Expect: mean: 0.2667 K median: 0.2666 K rms: 0.02536 K min: 0.18 K max: 0.35 K\n" ] } ], "source": [ "scan_block2 = sdfits.getps(ifnum=0, plnum=0, fdnum=0, smoothref=31)\n", "ta2 = scan_block2.timeaverage()\n", "ta2.plot(xaxis_unit=\"chan\", yaxis_unit=\"mK\", ymin=100, ymax=600, grid=True)\n", "print(f\"Stats : {fmt_str}\".format(**ta2[5000:10000].stats())) # avoid spikes\n", "print(\"Expect: mean: 0.2667 K median: 0.2666 K rms: 0.02536 K min: 0.18 K max: 0.35 K\")" ] }, { "cell_type": "markdown", "id": "3c76eb98-a5d0-4abc-93f2-f87038179687", "metadata": {}, "source": [ "When smoothing the reference spectra, the noise in the time average of the calibrated spectrum is 25 mK. For a given \n", "value of the `smoothref=`$N$ the improvement is \n", "\n", "$$\n", "\\sqrt{\\frac{N+1}{2N}}\n", "$$\n", "\n", "Given that the RMS was 36 mK we would expect $\\sqrt(32/62) \\approx 0.72$ times, or 26 mK. At best this will give about 30% improvement.\n", "See our [sdmath](https://dysh.readthedocs.io/en/latest/reference/sdmath.html) for more details.\n", "\n", "\n", "We repeat the baseline subtraction using the same model." ] }, { "cell_type": "code", "execution_count": 11, "id": "41f61fa9-05d8-4e9f-afbd-87121eb71b96", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "11:40:10.559 I EXCLUDING [Spectral Region, 1 sub-regions:\n", " (1397351017.8085938 Hz, 1403454533.4335938 Hz) \n", "]\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "6c870a6b01324d71a3014397bd25dcaa", "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": [ "Stats : mean: 0.0009 K median: 0.0007 K rms: 0.0270 K min: -0.19 K max: 0.50 K\n", "Expect: mean: 0.0009 K median: 0.0007 K rms: 0.0270 K min: -0.19 K max: 0.50 K\n" ] } ], "source": [ "ta2.baseline(model=\"poly\", degree=1, exclude=[(14000,18000)], remove=True)\n", "ta2.plot(xaxis_unit=\"chan\", yaxis_unit=\"mK\", ymin=-200, ymax=300, grid=True)\n", "print(f\"Stats : {fmt_str}\".format(**ta2[5000:14000].stats()))\n", "print(\"Expect: mean: 0.0009 K median: 0.0007 K rms: 0.0270 K min: -0.19 K max: 0.50 K\")" ] }, { "cell_type": "markdown", "id": "5e695829-ad09-4956-870d-86f564742b6e", "metadata": {}, "source": [ "We could smooth this spectrum the normal way, as was done a few cells ago, and not much difference is visible, except for the noise." ] }, { "cell_type": "code", "execution_count": 12, "id": "26d5240f-c475-48bc-84ea-e817dc38363d", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "c6a4c00c54b94839a4d77d68df183429", "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": [ "Stats : mean: 0.0009 K median: 0.0007 K rms: 0.0023 K min: -0.00 K max: 0.01 K\n", "Rolled rms/rms ratio: 0.7802633164999722\n" ] } ], "source": [ "ts3 = ta2.smooth('box', 320)\n", "ts3.plot(xaxis_unit=\"chan\", yaxis_unit=\"mK\", ymin=-200, ymax=300, grid=True)\n", "print(f\"Stats : {fmt_str}\".format(**ts3[5000//320:14000//320].stats()))\n", "rratio = ts3[5000//320:14000//320].stats(roll=1)[\"rms\"]/ts3[5000//320:14000//320].stats(roll=0)[\"rms\"]\n", "print(f\"Rolled rms/rms ratio: {rratio}\")" ] }, { "cell_type": "markdown", "id": "bcd7355a-8cd6-4ea8-88bf-128486406089", "metadata": {}, "source": [ "The RMS has gone down (2.4 mK to 2.3 mK), but the signal correlation has also increased from 1.06 to 1.10. This increase is due to the added correlation of the reference smoothing." ] }, { "cell_type": "markdown", "id": "d950a619-e889-411e-82a2-6eb8d7d54604", "metadata": {}, "source": [ "## Successive Smoothing\n", "\n", "Smoothing with a Gaussian twice in a row should be the same as smoothing with a single Gaussian of the square root of the sum of their squares. Note that the **width** is the final width (FWHM) of the smoothing operation. We need to ignore decimation here, which is the third argument (`decimate=-1`), otherwise the spectra will have a different number of channels.\n", "\n", "We plot the difference between the two smoothing routes, which is perhaps surprisingly still below 1e-8. It is not closer to 0 due to the finite range the Gaussian is sampled (to 4 sigma)." ] }, { "cell_type": "code", "execution_count": 13, "id": "252f27f1-d63f-4b1e-bda2-5493b85aec65", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "759618eaac1649b8894d3f585bf2910f", "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": [ "Stats : mean: -2.945e-10 K median: 4.285e-09 K rms: 7.291e-07 K min: -2.8e-06 K max: 2.7e-06 K\n" ] } ], "source": [ "ts5a = ta.smooth('gaussian', width=3, decimate=-1).smooth('gaussian', width=5, decimate=-1)\n", "ts5b = ta.smooth('gaussian', width=5, decimate=-1)\n", "diff = ts5a-ts5b\n", "diff.plot(ymin=-0.00001, ymax=0.00001)\n", "\n", "fmt_str = \"mean: {mean:.4g} median: {median:.4g} rms: {rms:.4g} min: {min:.2g} max: {max:.2g}\"\n", "print(f\"Stats : {fmt_str}\".format(**diff[5000:14000].stats()))" ] }, { "cell_type": "markdown", "id": "c1993d96-6194-419e-bc2f-3dd751df3aef", "metadata": {}, "source": [ "## Final Stats\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": 14, "id": "d6f11d87-b230-45f4-973a-ccc90988e93c", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "11:40:11.314 I rms is OK \n" ] } ], "source": [ "ts5b.check_stats(0.04618235 * 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" }, "toc": { "base_numbering": 0 } }, "nbformat": 4, "nbformat_minor": 5 }