diff --git a/notebooks/single_baseline_postprocessing_and_pspec.ipynb b/notebooks/single_baseline_postprocessing_and_pspec.ipynb index a348073..b8bcfb1 100644 --- a/notebooks/single_baseline_postprocessing_and_pspec.ipynb +++ b/notebooks/single_baseline_postprocessing_and_pspec.ipynb @@ -7,7 +7,7 @@ "source": [ "# Single Baseline Filtering and Power Spectrum Estimation\n", "\n", - "**by Josh Dillon, Bobby Pascua, and Mike Wilensky**, last updated October 28, 2024\n", + "**by Josh Dillon, Bobby Pascua, and Mike Wilensky**, last updated November 5, 2024\n", "\n", "This notebook is designed to take a single redundantly-averaged unique baseline (typically after LST-binning) and push it through all the way to the power spectrum. It operates on single files that contain a single baseline for all LSTs and both `'ee'` and `'nn'` polarizations. It then can:\n", "* Throw out highly flagged times and/or channels\n", @@ -151,6 +151,8 @@ "FR_EIGENVAL_CUTOFF = float(os.environ.get(\"FR_EIGENVAL_CUTOFF\", 1e-12))\n", "TARGET_AVERAGING_TIME = float(os.environ.get(\"TARGET_AVERAGING_TIME\", 300.0)) # coherent integration time in seconds. Actual time might be less to so that all interleaves have the same number of samples averaged\n", "USE_CORR_MATRIX = os.environ.get(\"USE_CORR_MATRIX\", \"TRUE\").upper() == \"TRUE\" # Whether or not to use the correlation matrix to do the noise statistics correction factor\n", + "CORR_MATRIX_FREQ_DECIMATION = int(os.environ.get(\"CORR_MATRIX_FREQ_DECIMATION\", 10)) # How much to decimate the frequency axis when getting the noise statistic correction factor\n", + "CORR_MATRIX_NOTCH_CUTOFF = float(os.environ.get(\"CORR_MATRIX_NOTCH_CUTOFF\", 30.)) # Maximum EW-projection in meters for including the notch filter in error covariance calculation\n", "\n", "# Power spectrum settings\n", "EFIELD_HEALPIX_BEAM_FILE = os.environ.get(\"EFIELD_HEALPIX_BEAM_FILE\", \"/lustre/aoc/projects/hera/h6c-analysis/IDR2/beams/NF_HERA_Vivaldi_efield_beam_healpix.fits\")\n", @@ -1417,21 +1419,26 @@ "outputs": [], "source": [ "def get_frop_wrapper(mode=\"main_lobe\", pol=\"nn\", stream_ind=0, band_ind=0, t_avg=AVERAGING_TIME, \n", - " rephase=True, wgt_tavg_by_nsample=True, nsamples=None, bl_vec=None,\n", - " dlst=None, coherent_avg=True, times=None):\n", + " rephase=True, wgt_tavg_by_nsample=True, bl_vec=None,\n", + " dlst=None, coherent_avg=True, times=None,\n", + " freq_decimation=CORR_MATRIX_FREQ_DECIMATION):\n", " \"\"\"\n", " This wraps hera_cal.frf.get_frop_for_noise using specific information from this notebook so that we can do\n", " post FRF time-time visibility covariance calculation. It returns an operator that can be used on a dynamic \n", " spectrum to filter it down the time axis.\n", " \"\"\"\n", " bl=(ANTPAIR[0], ANTPAIR[1], pol)\n", + " \n", " band = bands[band_ind]\n", " band_slice = band_slices[band_ind]\n", - " Nfreqs = band_slice.stop - band_slice.start\n", + " if freq_decimation > 1:\n", + " band_slice = slice(band_slice.start, band_slice.stop, freq_decimation)\n", "\n", " weights=deint_wgts[stream_ind][bl][:, band_slice]\n", + " nsamples = deint_nsamples[stream_ind][bl][:, band_slice]\n", + " \n", " if times is None:\n", - " times = np.modf(single_bl_times[tslice][stream_ind::4])[0] * 24 * 3600\n", + " times = np.modf(single_bl_times[tslice][stream_ind::NINTERLEAVE])[0] * 24 * 3600\n", " times = times - times[0]\n", "\n", " if mode == \"main_lobe\":\n", @@ -1461,8 +1468,12 @@ "metadata": {}, "outputs": [], "source": [ + "\n", "coherent_avg_correction_factors = []\n", "for spw, band in enumerate(bands):\n", + " freq_slice = band_slices[spw]\n", + " decimated_freq_slice = slice(freq_slice.start, freq_slice.stop, CORR_MATRIX_FREQ_DECIMATION)\n", + " \n", " if USE_CORR_MATRIX:\n", " corr_factors = []\n", " for stream_ind in range(NINTERLEAVE):\n", @@ -1470,18 +1481,36 @@ " cov_here = 0\n", " for pol in (\"ee\", \"nn\"):\n", " cross_antpairpol = ANTPAIR + (pol,)\n", - " freq_slice = band_slices[spw]\n", + "\n", " var = frf.prep_var_for_frop(deint_filt_data[stream_ind],\n", " deint_nsamples[stream_ind],\n", " deint_wgts[stream_ind],\n", " cross_antpairpol,\n", - " freq_slice,\n", + " decimated_freq_slice,\n", " auto_ant=0)\n", - " frop = get_frop_wrapper(pol=pol, stream_ind=stream_ind, band_ind=spw,\n", - " nsamples=deint_nsamples[stream_ind][cross_antpairpol][:, freq_slice],\n", - " dlst=dlst, bl_vec=bl_vec[cross_antpairpol],\n", - " times=times)\n", - "\n", + " main_lobe_frop = get_frop_wrapper(pol=pol, stream_ind=stream_ind, band_ind=spw,\n", + " dlst=dlst, bl_vec=bl_vec[cross_antpairpol],\n", + " times=times)\n", + " \n", + " # Calculate whether the EW-baseline projection is short enough to care about the notch\n", + " EW_proj = data.antpos[ANTPAIR[0]][0] - data.antpos[ANTPAIR[1]][0]\n", + " if np.abs(EW_proj) < CORR_MATRIX_NOTCH_CUTOFF:\n", + " xtalk_frop = get_frop_wrapper(pol=pol, stream_ind=stream_ind, band_ind=spw,\n", + " dlst=dlst, bl_vec=bl_vec[cross_antpairpol],\n", + " times=times, rephase=False, mode=\"xtalk\", coherent_avg=True,\n", + " t_avg=times[1] - times[0])\n", + "\n", + "\n", + " frop = np.zeros_like(main_lobe_frop)\n", + " Nfreqs_band = freq_slice.stop - freq_slice.start\n", + " Nfreqs_calc = int(np.ceil(Nfreqs_band / CORR_MATRIX_FREQ_DECIMATION))\n", + "\n", + " for freq_ind in range(Nfreqs_calc):\n", + " frop[:, :, freq_ind] = np.tensordot(main_lobe_frop[:, :, freq_ind],\n", + " xtalk_frop[:, :, freq_ind],\n", + " axes=1)\n", + " else:\n", + " frop = main_lobe_frop\n", " cov_here += frf.get_FRF_cov(frop, var) # computes covariance for pI by adding ee and nn\n", " if hd.pol_convention == \"avg\":\n", " cov_here /= 4\n", @@ -1494,7 +1523,8 @@ "\n", " else:\n", " # use the mode-counting method, which is simpler but less accurate\n", - " coherent_avg_correction_factors.append(dpss_coherent_avg_correction(spw))" + " coherent_avg_correction_factors.append(dpss_coherent_avg_correction(spw))\n", + " " ] }, { @@ -2068,7 +2098,9 @@ "cell_type": "code", "execution_count": null, "id": "93240fe7", - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ "plot_Pk_SNR_vs_LST(func=np.real)\n", @@ -2222,7 +2254,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.0" + "version": "3.9.6" }, "toc": { "base_numbering": 1,