High cadence kurtosis based RFI excision forCHIMEbyArash MirhosseiniB.Sc., Leipzig University, 2016M.Sc., Paris-Saclay University, 2018A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Astronomy)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)October 2020© Arash Mirhosseini 2020The following individuals certify that they have read, and recommend to theFaculty of Graduate and Postdoctoral Studies for acceptance, the thesis entitledHigh cadence kurtosis based RFI excision for CHIME submitted by ArashMirhosseini in partial fulfillment of the requirements for the degree of Master ofScience in Astronomy:• Mark Halpern, Physics and Astronomy (Supervisor)• Richard Shaw, Physics and Astronomy (Examining committee member)iiAbstractThis document describes the real-time Radio Frequency Interference detection sys-tem for the Canadian Hydrogen Intensity Mapping Experiment (CHIME). CHIMEis a transit radio telescope located at Dominion Radio Astrophysical Observatory(DRAO) in Penticton, BC, and it is originally designed to map the large-scalestructures in the redshift range 0.8 < 푧 < 2.5 by observing the 21-cm emissionline of the neutral hydrogen atom. One of the common problems for astrophysicalradio observations is Radio Frequency interference (RFI) from terrestrial sourcessuch as TV stations, airplanes, cellphones, etc. RFI detection and mitigation is anessential part of any research in radio astronomy, because RFI contaminates theastrophysical data and reduces the sensitivity of the telescope to the faint sources.Since most of the RFI is non-Gaussian and lasts less than one second, we developeda real-time high cadence RFI excision system for CHIME which uses the fourthstatistical moment (kurtosis) to detect non-gaussianity in the signal.In this thesis, I introduce the algorithms for kurtosis based RFI excision that wehave used in CHIME. The algorithms were tested and the results were comparedwith each other. I also discuss the effect of truncation of the samples in CHIMEcorrelator on the spectral kurtosis estimates. I show that truncation bias causes theRFI system to flag bright celestial sources. I derive a correction for the truncationbias with a polynomial fitting and a cubic spline interpolation. Moreover, I evaluatethe quality of the CHIME data taken between May 2019 and September 2020. Ifind that the RFI excision system can mitigate many types of RFI by excising lessthan 20% of the data (on average), from intermittent RFI caused by satellites orairplanes to static RFI, especially between 400 MHz and 500 MHz.iiiLay SummaryThe Canadian Hydrogen Intensity Mapping Experiment (CHIME) is a radio tele-scope in Penticton, BC, that is originally designed to probe the dark energy andmake the largest 3D map of the large-scale structures in the observable universe.One of the common problems in radio astronomy research is the Radio FrequencyInterference (RFI). RFI is a human made signal that interferes with the cosmic sig-nal and degrades the quality of astronomical data. Since CHIME aims to measurevery weak signals coming from deep space, we need to detect and mitigate RFI toincrease the sensitivity of the telescope. In this thesis, I introduce the RFI excisionsystem for CHIME and I evaluate its performance from May 2019 to September2020. I find that the system is able to detect and mitigate different types of RFI andincrease the sensitivity of the system.ivPrefaceThis thesis is original, unpublished work by the author, Arash Mirhosseini, con-ducted as part of the CHIME collaboration, under supervision of Mark Halpernwith input and ideas from Richard Shaw, as well as other members of the CHIMEteam. The RFI excision system for CHIME was previously designed and explainedin Taylor et al. (2019) [1]. The contribution I made was to evaluate the performanceof different excision algorithms with different thresholds, find and correct the bugsin the spectral kurtosis estimator (for example 4+4-bits truncation bias) and assessthe data quality after RFI excision.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Theoretical background . . . . . . . . . . . . . . . . . . . . . . 11.1.1 Cosmological models . . . . . . . . . . . . . . . . . . . 21.2 Dark Energy probes . . . . . . . . . . . . . . . . . . . . . . . . 31.2.1 Supernovae Type Ia . . . . . . . . . . . . . . . . . . . . 31.2.2 Baryon Acoustic Oscillations . . . . . . . . . . . . . . . 51.2.3 Weak Gravitational lensing . . . . . . . . . . . . . . . . 81.3 CHIME . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101.3.1 21-cm intensity mapping . . . . . . . . . . . . . . . . . 101.3.2 Radio Frequency Interference . . . . . . . . . . . . . . . 111.3.3 CHIME instrument . . . . . . . . . . . . . . . . . . . . 121.3.4 RFI environment at DRAO . . . . . . . . . . . . . . . . 142 Kurtosis-based RFI excision . . . . . . . . . . . . . . . . . . . . . . 162.1 Statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.1.1 Moments of a probability distribution . . . . . . . . . . . 162.1.2 Cumulant . . . . . . . . . . . . . . . . . . . . . . . . . . 172.2 Spectral kurtosis . . . . . . . . . . . . . . . . . . . . . . . . . . 182.3 Spectral kurtosis estimator . . . . . . . . . . . . . . . . . . . . . 192.4 Spectral kurtosis for a compact array . . . . . . . . . . . . . . . 21viTable of Contents2.5 Implementation for CHIME . . . . . . . . . . . . . . . . . . . . 233 High Cadence Excision Algorithms . . . . . . . . . . . . . . . . . . 253.1 Excision algorithms . . . . . . . . . . . . . . . . . . . . . . . . 253.1.1 Single Stage . . . . . . . . . . . . . . . . . . . . . . . . 253.1.2 Excision on 30 ms frames (EOF) . . . . . . . . . . . . . 263.1.3 Two stage . . . . . . . . . . . . . . . . . . . . . . . . . . 273.2 Evaluation of excision algorithms . . . . . . . . . . . . . . . . . 283.2.1 Gaussianity test . . . . . . . . . . . . . . . . . . . . . . 283.2.2 Offline test . . . . . . . . . . . . . . . . . . . . . . . . . 293.2.3 Comparison of different excision algorithms . . . . . . . 354 SK estimator biases . . . . . . . . . . . . . . . . . . . . . . . . . . . 364.1 4+4 bits truncation . . . . . . . . . . . . . . . . . . . . . . . . . 364.1.1 Effect of truncation bias on CHIME data . . . . . . . . . 394.1.2 Correction for truncation bias . . . . . . . . . . . . . . . 414.1.3 Results of truncation bias correction on CHIME data . . . 424.2 Common mode signal . . . . . . . . . . . . . . . . . . . . . . . 465 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 495.1 May to November 2019: Single stage algorithm . . . . . . . . . . 495.1.1 Moving blob . . . . . . . . . . . . . . . . . . . . . . . . 515.2 December 2019 to February 2020: Excision on 30 ms frames withmultiple thresholds . . . . . . . . . . . . . . . . . . . . . . . . . 525.3 March to September 2020: Two stage algorithm . . . . . . . . . . 555.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 575.4.1 Reliability of Gaussianity test . . . . . . . . . . . . . . . 575.4.2 Solar transit issue . . . . . . . . . . . . . . . . . . . . . 595.4.3 RFI polarization . . . . . . . . . . . . . . . . . . . . . . 606 Summary and future work . . . . . . . . . . . . . . . . . . . . . . . 61Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62viiList of Tables3.1 the results of single stage excision for different values of parameter푛 on August 29 between 5:00 and 6:00 PM. Second column showsthe number of frequency channels having 퐺푇 < 2 after excision.Third and fourth columns show the average fraction of exciseddata for the frequency channels that passed the test and for allthe frequency channels, respectively. As expected, we recovermore frequency channels by decreasing the threshold value at theexpense of loosing more data. Note that without any RFI excision288 frequency channels have passed the Gaussianity test. . . . . . 333.2 Results of excision with EOF algorithm. Third column shows thenumber of frequency channels having 퐺푇 < 2 after excision. Thenumber of good frequency channels (퐺푇 < 2) is smaller than thesingle stage algorithm. This is because the single stage algorithmis more sensitive to the very short RFI events where only a fewsamples exceed the detection threshold. EOF algorithm is sensitiveto the longer RFI, when more than a few percent of the samples in30 ms exceed the threshold, but it is insensitive to the very shortRFI. The meaning of the third and fourth columns is similar to theones in the table 3.1. . . . . . . . . . . . . . . . . . . . . . . . . 333.3 The results of gaussianity test for the two stage algorithm. I used5휎 threshold for the first stage, and different values of (푛, 푓 ) forthe second stage. The meaning of the third and fourth columns issimilar to the ones in the table 3.1. . . . . . . . . . . . . . . . . . 34viiiList of Figures1.1 Left: Hubble diagram from the JLA sample (top) and Residualsfrom the best ΛCDMfit (bottom). Right: 68% and 95% confidencecontours of 푤 and 푤푎 parameters for the flat 푤 − ΛCDM model(figures taken from [14]). . . . . . . . . . . . . . . . . . . . . . . 51.2 Figure shows the generation of BAO peak from initial overdensity.Dark matter overdensity sits where it is initially, as it has no elec-tromagnetic interaction. Photons (red) and baryons (blue) movetogether until recombination (푧푟푒푐 ∼ 1100). After recombinationphotons free-stream and baryons attract surrounding dark mattergravitationally. In the end, there will be a mass concentration at acharacteristic scale of 푟푠 ≈ 148 Mpc (figure from [16]) . . . . . . 71.3 CHIME telescope. It is an array of 4 cylinders in the east-westdirection, with axes oriented in the north-south direction. . . . . . 131.4 The autocorrelations of 30 ms samples for one feed. We can seevarious RFI features in this plot, namely LTE band from ∼ 730MHz to 755 MHz, TV station bands around 500-600 MHz, andrepeating RFI spots between 400 MHz and 500 MHz. 6 MHzwide bursts from distant broadcast TV bands appear as blob around00:15 PT in 480 MHz or 520 MHz. White lines correspond to themissing frequencies due to the malfunctioning of some of the GPUnodes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.1 Left: Histogram of 30000 realizations of SK estimator. Each SKvalue is generated from 20000 complex random variables drawnfrom a Gaussian distribution. Red region indicates 1휎 intervalwhere 휎 is the theoretical standard deviation. Right: Cumulativehistogram of the SK values. Yellow region corresponds to 1휎interval. The SK distribution is symmetric around 1, and %68 ofthe estimates lie within one standard deviation. These two confirmthe derived mean and variance of the SK estimator. . . . . . . . . 21ixList of Figures3.1 Histogram of ∼ 500000 random numbers drawn from a Gaussiandistribution with mean 1 and some variance 휎2. These valuesrepresent the RFI-free SK values. RFI events change the shapeof the RFI-free SK distribution. So, a few non-Gaussian featuresare added to the the distribution to represent the RFI. The singlestage algorithm is sensitive to the features beyond 5휎, while EOFis sensitive to lower power, but longer duration RFI. . . . . . . . . 263.2 The probability for a Gaussian 30 ms frame to pass the EOF algo-rithm without being removed for different values of (푛, 푓 ) param-eters. For example it is very probable for a 30 ms frame to not tobe masked by EOF algorithm with (푛, 푓 ) = (3, 0.4). In this casethe probability for the frame to pass the algorithm is more than 8휎.The two red lines show two thresholds that we have used used forRFI detection between December 2019 and February 2020 usingEOF algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.3 Left: Result of gaussianity test without anyRFI excision (horizontalaxis) and with the single stage algorithm (vertical axis). Right:Gaussianity test with EOF algorithm using (푛, 푓 ) = (1.5, 0.4).Each data point corresponds to one frequency bin and the colorshows the average excised rate at that frequency. In other words, thefigure shows three things for each frequency bin: GT value beforeRFI excision, GT value after applying the excision algorithm, andthe average excised rate over ∼ 1 hour. The color bar is the same forboth figures. Most of the data points are below the diagonal line.This means that Excision algorithm has improved the Gaussianityof those frequencies. . . . . . . . . . . . . . . . . . . . . . . . . 313.4 Average excision rate for 8 s samples over 1 hour as a function offrequency. 600-700 MHz band is the cleanest part of the spectrum.There are many RFI features in 400-500 MHz band that. Thekurtosis based RFI excision system detected those RFI and somefrequency channels are cleaned by excising only ∼ 20% of thesamples in the integrated 8 second sample. . . . . . . . . . . . . 323.5 A comparison between two excision algorithms with different 푓parameters. The number above data points shows the average per-centage of excised data in good frequency channels. . . . . . . . . 34xList of Figures3.6 The performance of single stage (black), EOF (green) and twostage (red) algorithms over ∼ 8 hours. In general, two stage andEOF outperform the single stage. The difference between twostage and EOF is not significant because this test was performed atnight, when there are no strong source of transient RFI. If there arepowerful, short time scale RFI, the first stage of two stage algorithmcan remove it, while EOF is not able to do so. . . . . . . . . . . . 354.1 Effect of 4-bits truncation on the shape of a Gaussian signal withan rms of ∼ 5.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . 374.2 Effect of 4+4-bits truncation on SK estimates. The SK estimatesare generated from 6000 Gaussian dataset with rms ∼ 2.9. The SKestimates tend to be smaller after truncation. . . . . . . . . . . . . 384.3 Bias of the SK estimator as a function of the rms of truncatedsamples. Each data point shows the deviation of the average SKvalue (over 6000 Gaussian dataset) from unity for a given rms. . . 394.4 Bias of SK values in terms of the nominal 휎 for CHIME. TheSK values are more negatively biased during the transit of Cyg-awhich is around ∼ 9:29 AM. Note that lower frequencies are morenegatively biased than the higher frequencies. . . . . . . . . . . . 404.5 The excised rate during Cyg-A transit . Each data point representsthe fraction of excised samples in a 2 seconds frame averagedbetween 400 MHz and 500 MHz. . . . . . . . . . . . . . . . . . . 414.6 Left: Fitting residuals in terms of nominal 휎 as a function of rmsfor different polynomials. Right: Spectral kurtosis estimated fromtruncated samples as function of the rms of the truncated samples. 424.7 The bias of SK estimates in terms of nominal휎 for CHIMEwithout(Left) and with (Right) the correction for truncation bias.Cyg-Atransit is around 4:45 PT in this plot. Each data point is the averageSK value over 30 ms and the white spots are 30 ms frames thatare masked by EOF algorithm with (푛, 푓 ) = (1.5, 0.4). The colorscale is the same in both plots. . . . . . . . . . . . . . . . . . . . 444.8 Left: The excised fraction of the samples in an 2 second frame,averaged between 400 MHz and 500 MHz during Cyg-A transit(∼ 4:45 PT) without truncation bias correction. Right: The samequantity as in the left figure, but corrected for truncation bias witha polynomial of degree 5. The bump corresponding to Cyg-Adisappeared, i.e. the RFI system no longer excise the Cyg-A. . . . 44xiList of Figures4.9 Deviation of SK values from unity after correcting the kurtosisvalues for truncation bias with a cubic spline function during Cyg-A transit around ∼ 4 : 45푃푇 . The bias is closer to zero for allfrequencies compared to the figure 4.7. . . . . . . . . . . . . . . . 454.10 Variance of the SK estimator in presence of a commonmode signal.The variance reaches to 10−2 if the signal is totally dominated bythe common mode signal, as there are only 푛 = 256 independenttime samples. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 484.11 Number of effective antenna as a function of fractional commonmode amplitude. . . . . . . . . . . . . . . . . . . . . . . . . . . 485.1 Left: Normalized autocorrelations in April 10, 2019 without anyRFI excision. Right: Normalized autocorrelations in June 7, 2019after deploying high cadence RFI excision system with single stagealgorithm. Most of the repeating RFI features disappeared. . . . . 505.2 Excised fraction of 10 s samples averaged over 1.5 hours as afunction of frequency. . . . . . . . . . . . . . . . . . . . . . . . . 505.3 The signature of the Meridian satellite is visible in the plot ofexcision rate of 10 s samples (in percent) as a function of timearound frequency ∼484 MHz. Each diagonal line corresponds to24 hours. The blob appears every day, but it slowly drifts in rightascension. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 515.4 Excised rate of 10 s samples averaged over all frequencies for 3consecutive days. The red box is the Cyg-A transit time and twogreen boxes correspond to the meridian satellite. The excision rateincreases every day during Cyg-A transit. This is due to 4+4-bittruncation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 535.5 Average excised rate over 1.5 hours as a function of frequency fora typical day in February 2020. The excision algorithm completelymasks the LTE and TV station bands. . . . . . . . . . . . . . . . 535.6 Number of good frequency channels reported every 1.5 hours for 12consecutive days. We can see a repeatable pattern every day. Notethat there are zero good frequency channels during solar transit, asthe RFI system flags the Sun. This issue is discussed in section 5.4.2 545.7 Change of the number of good frequency channels in 24 hours. Thenumber of good frequency channels is more at night and it goesdown in the morning, when there are more RFI at site. . . . . . . 545.8 Excised fraction as a function of time and frequency. The greencircles and rectangle show the blob which is fixed in RA (whichturned out to be Cyg-A). . . . . . . . . . . . . . . . . . . . . . . 55xiiList of Figures5.9 Excised rate of 10 second samples averaged between 400 and 430MHz on May 19, 2020; without truncation bias correction. Redarrow shows the Cyg-A transit time. . . . . . . . . . . . . . . . . 565.10 Excised rate of 10 s samples averaged between 400 and 430MHz onJune 16, 2020. The SK estimator is corrected for the truncation bias.The The bump corresponding to the Cyg-A transit disappeared. . 565.11 Number of good frequencies from July 9, 2020 to July 14,2020.Two stage algorithm was used in this period with 5휎 threshold onindividual 0.65 ms samples and (3휎, 0.13) on 30 ms frames. Zeroscorrespond to the solar transit, when the RFI excision system is off. 575.12 Standard deviation of the data between 3:00 and 4:00 AM for thosechannels which pass the gaussianity test. Two of the outliers areshown with black circles. . . . . . . . . . . . . . . . . . . . . . . 585.13 Left: Histogram of the autocorrelations one of the outliers at fre-quency 419.14MHz. Right: Histogram after subtracting the neigh-bouring samples. . . . . . . . . . . . . . . . . . . . . . . . . . . 585.14 Deviation of SK values from unity in terms of nominal 휎 forCHIME during solar transit (around 13:00). Left: Without trun-cation bias correction SK values have a very high negative bias.Right: Since SK values sill highly deviate from expected valueafter truncation bias correction with a fifth degree polynomial, theSun is still being flagged by the RFI system. This might be ex-plained by the fact that the variance of the estimator is increasedwhen the common mode signal dominates. . . . . . . . . . . . . 595.15 Result of the Gaussianity test for all feeds at a single frequency. TheNS and EW corresponds to two different polarizations. Note thatthe first 256 inputs on each cylinder have a different polarizationthan the next 256. . . . . . . . . . . . . . . . . . . . . . . . . . . 60xiiiChapter 1IntroductionTwenty years ago, measurements of the distance to the Type Ia supernovae bytwo independent teams confirmed the discovery of accelerating expansion of theuniverse [2, 3]. The nature of this acceleration is still unknown and it is usuallyattributed to an enigmatic component of the universe: dark energy. Dark energycan be considered as a fluid with an equation of state 푃 = 푤휌, where P and 휌are the pressure and density of the fluid, respectively, and 푤 is the equation ofstate parameter. The current observations are compatible with 푤 = −1 whichcorresponds to a static density. More accurate observations are needed to imposetighter constraints on 푤 as well as other cosmological parameters. Putting tighterconstraints on the parameter 푤 is an essential step towards understanding the natureof dark energy.In this chapter, I review somebasic concepts that are necessary for understandingthe scientific goal of CHIME; a project whose primary goal is to characterize thedark energy bymapping the large scale structures between the redshift 0.8 ≤ 푧 ≤ 2.5using the 21-cm intensity mapping technique.1.1 Theoretical backgroundAccording to the cosmological principle, theUniverse is homogeneous and isotropicat large scales. This principle is observationally supported by the quasi isotropy ofthe cosmic microwave background (CMB) and the homogeneity of the distributionof the galaxies at large scales. In the framework of General Relativity, symmetriesof the cosmological principle can be used to define the metric for a homogeneousand isotropic universe, the Friedman-Lemaitre-Robertson-Walker (FLRW) metric:푑푠2 = 푑푡2 − 푎2(푡)( 푑푟21 − 퐾푟2 + 푟2(푑휃2 + 푠푖푛2휃푑휙2)), (1.1)where 푎(푡) is the dimensionless scale factor, and 퐾 = −1, 0 or +1 is the spatialcurvature parameter corresponding to an open, flat or closed universe, respectively.Then, the evolution of such a homogeneous and isotropic universe is described bythe Friedman equations. They are derived from Einstein equations by substituting11.1. Theoretical backgroundthe Ricci tensor (computed from FLRW metric) together with the stress-energytensor for a perfect fluid into Einstein field equations:퐻2(푧) ≡ ( ¤푎(푡)푎(푡))2 = − 퐾푎2(푡)+8휋퐺3∑휌 (1.2)¥푎(푡)푎(푡)= −4휋퐺3∑(휌 + 3푝), (1.3)where 퐻(푧) is the Hubble constant at redshift 푧, and 휌 and 푝 are the density andpressure of different components of the Universe, respectively. The sum is takenover all components of the Universe, namely, matter, radiation and dark energy.The equation of state for the component 푥 of the Universe is 푝푥 = 푤휌푥 . Thecritical density is defined as the total density of a flat universe 휌푐 = 3퐻2(푧)/8휋퐺(current critical density). The density parameter of the component 푥 of the Universeis Ω푥 = 휌푥/휌푐 . The first Friedman equation (Eq. 1.2) for a flat universe can bewritten in terms of density parameters of different components of the Universe:퐻2(푧) = 퐻20(Ω푚(1 + 푧)3 + Ω푟 (1 + 푧)4 + Ω퐷퐸휌퐷퐸(푧)휌퐷퐸(푧 = 0)), (1.4)where indices 푚, 푟, 퐷퐸 correspond to the matter, radiation and a generic form ofdark energy.The simplest form of dark energy is a cosmological constant. It turns out thata cosmological constant behaves like a fluid with static density ( ¤휌 = 0). In thiscase, from continuity equation 푝 = −휌 and so equation of state parameter forthe cosmological constant would be 푤 = −1. Although current observations arecompatible with 푤 = −1, dark energy in the form of a cosmological constant is notthe unique possibility.1.1.1 Cosmological modelsThe simplest model which provides a good fit to a wide range of observationaldata is the standard ΛCDM model. It assumes a flat universe with the simplestpossible form of dark energy, i.e. a cosmological constant (Λ), and Cold DarkMatter (CDM). A true cosmological constant has an equation of state parameter푤 = −1 which does not change in time. Standard ΛCDM model comprises 6 freeparameters, including density parameters of baryons Ω푏 and dark matter Ω퐷푀 ,optical depth at the time of reionization 휏 , current value of Hubble constant 퐻0,and primordial amplitude of scalar fluctuations 퐴푠, and spectral index of scalarfluctuations 푛푠 related to the inflationary epoch of the Universe.Although this model can be well fitted to the data to find the parameters, it is notunique and extensions to the standard ΛCDMmodel are possible. Some extensions21.2. Dark Energy probesto standard Λ CDM are models with dynamical dark energy (time-varying 푤),neutrino masses and additional relativistic particles. For example, if the darkenergy is not a cosmological constant, then it is possible that the 푤 parameterchanges with time. In this case, a time-varying 푤 parameter can be considered asa Taylor expansion of 푤 at the first order in the scale factor:푤(푡) = 푤0 + [1 − 푎(푡)]푤푎 . (1.5)Note that at low redshifts (i.e. 푎(푡) ≈ 1) we have 푤(푡) ≈ 푤0 and only high redshiftmeasurements are sensitive to 푤푎 parameters.A detailed list of extensions to the standard ΛCDM model is given in [4].Because there ismore than onemodelwhich can befitted to the current observationaldata, one needs to combine different cosmological probes (CMB, BAOs, SNe Ia,weak lensing, etc.) to constrain cosmological parameters and invalidate some ofthese models.1.2 Dark Energy probes1.2.1 Supernovae Type IaObjects with known luminosity can be used as distance estimators. Luminositydistance is defined by 푑2퐿 = 퐿/4휋 푓 , where 퐿 is the luminosity and 푓 is theobserved flux of the source. Hence, for an object with known luminosity, theproblem of finding the distance is reduced to flux measurement. Such an objectwith a known luminosity is called a standard candle. Type Ia supernova (SN Ia) isan example of standard candles. Although SNe Ia show a dispersion of about 40%at their peak luminosity (which makes them not good distance estimator), it turnsout that the dispersion is reduced if one considers the so called brighter-bluer andbrighter-slower correlations. Brighter-bluer relation is the correlation between therest frame color of the supernova and its maximum luminosity: bluer supernovaeare brighter than redder ones. Brighter-slower relation, known as Phillips relation[5], is the correlation between the maximum luminosity of SNe Ia, and the decayrate of their light curve after reaching the peak luminosity: The light curve of abright supernova decays slower than fainter ones. An empirical correction of thelight curve using these two correlations reduces the peak luminosity dispersion ofSNe Ia [6] and makes type Ia SNe a good distance estimator.In a flat universe, the luminosity distance of a source observed at a redshift 푧 isrelated to the other cosmological parameters 휃푖 , through:푑퐿(푧) = (1 + 푧) ×∫ 푧0푑푧′퐻(푧′, 휃). (1.6)31.2. Dark Energy probesHence, measurement of luminosity distance 푑퐿 will constrain the Hubble param-eter at redshift 푧, and consequently, the density parameters through equation 1.4.Therefore, type Ia SNe are powerful probes of the cosmic expansion history andthey are very sensitive to the 푤 parameter in the equation of state of dark energy.After the discovery of the dark energy in the late 1990s, several second generationsurveys have collected samples of a few hundred well-measured SNe Ia up to 푧 ∼ 1.The Supernova Legacy Survey (SNLS) used data taken as part of deep componentof the five-year Canada-France-Hawaii survey (CFHT-LS). Using a rolling-searchapproach (i.e., repeatedly imaging the same sky patch), it targeted four one squaredegree fields during 5 to 7 consecutive lunations per year using four different broad-band filters g푀 ,r푀 ,i푀 ,z푀 . During its 5 years of operation (mid-2003 to mid-2008)it delivered about 500 SNe Ia [7].SDSS-II supernova survey used SDSS camera on SDSS 2.5 m telescope atthe Apache Point Observatory (APO) and it searched for supernovae during thenorthern fall season of 2005 to 2007. Images were taken in ugriz SDSS passbands[8] with a typical cadence of once every four nights. The supernovae candidateswere identified by a computing cluster at APO within 24 hours of data collection.Then, the spectroscopic follow up was performed using a dozen of telescopes.Details of the SDSS-II supernova survey are given in [9] and [10]. Out of 4607candidates, 500 SN Ia were confirmed by spectroscopic follow up at a redshift푧 < 0.5. However, only 413 of them were used to constrain the cosmologicalparameters. Assuming a ΛCDM model, acceleration (Ω푀 < 2Λ) was detected ata confidence of 3.1휎. Moreover, with a flat geometry ΩΛ > 0 at a confidence of5.7휎 is required and Ω푀 = 0.315 ± 0.093. The details of the analysis is given in[11].Joint Light Curve Analysis (hereafter JLA) was part of SNLS-SDSS collabora-tive effort that was initiated in 2010 to improve the previous analysis of SNLS andSDSS teams. The main goals of this collaborative effort were to improve the accu-racy of photometric calibration of both surveys [12], to determine more rigorouslythe uncertainties in SNe Ia light curve models [13], and to include the full SDSS-IISNe Ia spectroscopic sample in cosmology analysis. JLA sample of SNe Ia whichis used to improve the cosmological constraints includes 740 supernovae selectedfrom the data of SDSS-II, SNLS, HST and several nearby experiments.JLA sample is currently state-of-the-art collection of SNe Ia in the redshiftrange 0.01 < 푧 < 1.2 [14]. However, distances to supernovae with a redshift above0.8 are not well measured due to photometric uncertainties. The next generationof surveys, including the Large Synoptic Survey Telescope (LSST) aim to increasethe sample size of SNe Ia at higher redshifts.41.2. Dark Energy probesFigure 1.1: Left: Hubble diagram from the JLA sample (top) and Residuals fromthe best ΛCDM fit (bottom). Right: 68% and 95% confidence contours of 푤 and푤푎 parameters for the flat 푤 − ΛCDM model (figures taken from [14]).1.2.2 Baryon Acoustic OscillationsAn alternative to the standard candles to constrain the cosmological parameters isthe standard rulers. A standard ruler is an object or statistical featurewhose intrinsicsize is known and do not change with time. One example of standard rulers is theBaryon Acoustic Oscillations (BAOs). In the following, I briefly discuss the originof BAOs.The primordial universe was a hot and dense plasma of photons, baryonsand dark matter, where baryons and photons were coupled together via Comptonscattering. The initial overdensities of matter attract the surrounding baryons anddark matter. Because dark matter has no electromagnetic interaction, it only feelsthe gravitational pull from the overdensity. The baryons feel two competing forces:An inward gravitational force from the overdensity and a photon pressure outward.Baryon-photon fluid is compressed, the temperature increases, and the outwardphoton pressure increases. At some point, the outward force from photon pressurewill be stronger than the compressing force of gravity, and thus the region will startto expand in the form of a sound wave. The speed of the sound wave in the plasmain the unit of the speed of light is:퐶푠(푧) =푐√3(1 + 푅(푧))(1.7)where 푅(푧) is the baryon to photon density ratio 푅(푧) = 휌푏(푧)/휌푟 (푧). At a redshift51.2. Dark Energy probesof 푧푟푒푐 ∼ 1100 the plasma is cool enough so that the electrons combine with protonsand form neutral hydrogen. At this time photons decouple from baryon. The photonpressure is removed, leaving a spherical shell of baryons at a fixed radius from theinitial overdensity. The characteristic comoving radius of the baryonic shell at thetime of recombination is:푟푠 =∫∞푧푟푒푐퐶푆(푧)푑푧퐻(푧)=∫∞푧푟푒푐퐶푆(푧)푑푧√퐻0(Ω푚(1 + 푧)3 + Ω푟 (1 + 푧)4, (1.8)which is the comoving distance the baryonic shell has travelled until the time ofrecombination. Redshift at recombination 푧푟푒푐 is precisely determined by atomicphysics, and the cosmological parameters in the above integral are measured bythe CMB observations. These two fix the comoving radius of the baryonic shell atrecombination. Planck has measured the acoustic scale at the time of recombination푟푠 ≈ 145 Mpc [15]. Figure 1.2 shows the evolution of acoustic scale as a functionof redshift.Now that the characteristic scale of the BAO is known, it can be used asa standard ruler to estimate the distances and constrain the 푤 parameter in theequation of state of dark energy. To use the BAO scale for distance measurements,we need to measure its observational scale. The angular diameter distance for a flatuniverse is:푑퐴(푧) =11 + 푧휒(푧) =11 + 푧∫ 푧0푐푑푧′퐻(푧′), (1.9)where 휒(푧) is the comoving distance. Since the comoving size of acoustic scale 푟푠is constant after recombination, the observed angular scale of BAO 휃 is related tothe angular diameter distance via:휃 =푟푠퐷푀=푟푠⊥푑퐴(1 + 푧), (1.10)where 퐷푀 is the comoving angular diameter distance, and 푟푠⊥ is the perpendicularcomponent of 푟푠 to the line of sight. Note that while perpendicular component of푟푠 in a redshift slice gives the angular diameter distance , the radial component 푟푠 ‖along the line of sight corresponds to different redshifts and constrain the Hubbleparameter through the Hubble’s law:푟푠 ‖ =∆푧퐻(푧)(1.11)The signature of BAOs is imprinted in the two point correlation function ofmatter distribution. Two point correlation function (hereafter, 2PCF) quantifies theexcess probability of finding a pair of galaxies at some redshift 푧 that are separated61.2. Dark Energy probesFigure 1.2: Figure shows the generation of BAO peak from initial overdensity. Darkmatter overdensity sits where it is initially, as it has no electromagnetic interaction.Photons (red) and baryons (blue) move together until recombination (푧푟푒푐 ∼ 1100).After recombination photons free-stream and baryons attract surrounding dark mat-ter gravitationally. In the end, there will be a mass concentration at a characteristicscale of 푟푠 ≈ 148 Mpc (figure from [16])71.2. Dark Energy probesby a distance 푠 compared to the case if they are uniformly random distributed. Todo this, a random catalog (that is a random sample of galaxies) is constructed andthe distribution of galaxies is compared with the data. A commonly used estimatorfor 2PCF 휁(푠) is given by [17]:ˆ휁(푠) =퐷퐷(푠) − 2퐷푅(푠) + 푅푅(푠)푅푅(푠), (1.12)where 푠 is the comoving separation between two galaxies, 퐷퐷(푠) and 푅푅(푠) arethe number of galaxy pairs with separation 푠 in the real-real catalog (data) andrandom-random catalog, respectively, and 퐷푅(푠) is the number of galaxy pairswith a separation 푠 between a galaxy in the real data and a galaxy in the randomcatalog. The BAO signal appears as a bump in the 2PCF.The first detection of the BAO signal in large scale was made in 2005 bymeasuring the large-scale correlation function from a spectroscopic sample ofluminous red galaxies from SDSS data release 3 in the redshift range of 0.16 to0.47 by Eisenstein et al. [18]. More recent measurements of the BAO scale havebeen made by different surveys includingWiggleZ Dark Energy Survey which used∼ 130000 galaxies at redshift 푧 = 0.6 [19], and SDSS-III BOSS survey of ∼ onemillion galaxies out to a redshift of 0.7 [20]. Moreover, BOSS detected the BAOsignal using the Lyman 훼 forest spectra of ∼ 48000 high redshift quasars across theredshift range 2.1 ≤ 푧 ≥ 3.5 [21].Since the BAO scale is very large, a large volume of the universe must be probedto decrease the sample variance. Therefore, we need a large sample of galaxiesto detect the BAO signal. For example, an order of 108 galaxies is required toapproach the cosmic variance limit at redshift 푧 > 1 [16]. Galaxy-redshift surveysare very expensive for this purpose: First, obtaining the redshift of each galaxythrough high resolution spectroscopy for a large sample of galaxies is very timeconsuming. Moreover, high redshift objects are fainter and harder to detect andthe measurements are dominated by shot noise. There is an alternative for thegalaxy-redshift surveys which is discussed in the section 1.3.1.1.2.3 Weak Gravitational lensingWeak gravitational lensing is a powerful probe of total matter density in the Uni-verse, without distinguishing between dark matter and baryonic matter. The inter-vening mass distribution between the galaxy and the observer distorts the imageof the galaxy. Such small distortions contain rich information about the matterdistribution on small and large scales and their evolution over time. Particularly,the dependence of the weak lensing signal on the angular diameter distance, and onthe matter power spectrum makes it a great tool to probe the dark energy [23]. A81.2. Dark Energy probesdetailed overview of the basics of of weak gravitational lensing and its applicationsin cosmology is given in [24].Two particularly powerful aspects of the method are that it is based on a geomet-rical observable, i.e. the distorted shapes of galaxy images, and that it is sensitiveto the gravitational potential of structures, without distinguishing between baryonicand dark matter. Particularly, weak lensing signal is sensitive to the matter powerspectrum over a redshift range and so it provides a measure of the growth rate oflarge scale structures. Because the growth rate depends strongly on the 푤 param-eter in the equation of state of the dark energy, cosmic shear can be used to studydifferent dark energy models. It is also sensitive to the angular diameter distance(through 휔 factor in equation 9) and so it can also be used to determine angulardiameter distances as a function of redshift and constrain 푤 parameter with thedistance-redshift relation. These two aspects make the cosmic shear measurementsa great tool to probe the dark energy.In the following, I quickly overview some of the ongoing weak lensing surveysand their results:• Kilo Degree Survey: Kilo Degree Survey (KiDS) [25] is a wide-field opticalimaging survey which covers 1500 square degree of the southern hemispheresky started on October 2011. KiDS uses OmegaCam wide-field camera on2.6-m VLT Survey Telescope (VST) of Paranal observatory in Chile. Thissurvey was designed to measure the galaxy population out to redshift ∼ 1and to measure the effect of weak lensing by structures along the line of sighton galaxy shapes. Full cosmological analysis of 450 square degree of KiDSdata (KiDS-450) including constraints on various cosmological parametersand extensions to the standard ΛCDM model are given in [26], [27] and[28]. Particularly, they found 푤 < −0.24 With the KiDS data alone, but theconstraint from the combined probe is 푤 < −0.73 (all at %95 CL).• DarkEnergy Survey: Dark Energy Survey (DES) is an optical galaxy surveyconducted using the 570 Megapixel DECam instrument mounted on 4-mBlanco Telescope in Chile. DES is planned to map 5000 square degree of thesky and the main goal of this survey is to study the dark energy using severaltechniques, including weak lensing. DES collaboration have used year-oneshape catalogs to study the shape of 26 million galaxies within the redshiftrange 0.2 to 1.3 over 1321 square degree of the sky and they constrainedvarious cosmological parameters in ΛCDM model [29], particularly, thedark energy equation of state parameter 푤 = −0.950.33−0.39. They concludedthat there is no disagreement between DES results and CMB data, and thereis no evidence for a 푤CDM model with 푤 deviating from −1.91.3. CHIME• HSC-SSP Survey: Hyper Suprime-Cam Subaru Strategic program (HSC-SSP) is an imaging survey on the 8.2-m Subaru telescope at Mauna KeaObservatory in Hawaii. The goal of this survey is to probe the nature ofdark energy and dark matter by various techniques, including weak lensing[30]. HSC-SSP is the deepest ongoing weak lensing survey. So although it’snarrower than DES, but its S/N ratio is higher at high redshift, enabling abetter measurement of the equation-of-state of dark energy at higher redshift.The cosmological constraints from cosmic shear power spectra with the firstyear data of HSC is recently published in [31]. Data Release 1 (DR1) ofHSC survey was publicly released in February 2017 and it is based on datataken using 61.5 nights between March 2014 and November 2015 [32]. The푤 parameter in the equation of state of dark energy is not well constrainedfrom shear analysis alone, 푤 = −1.37+0.43−0.36 . However, this constraint is stillconsistent with other observations and particularly it is compatible with 푤 =−1 which corresponds to the simplest form of dark energy, a cosmologicalconstant.1.3 CHIME1.3.1 21-cm intensity mappingThis section is based on Pritchard & Loeb (2012) [33].An alternative to the galaxy-redshift surveys for the measurements of BAO scaleis the 21-cm intensity mapping. The basic idea of 21-cm intensity mapping isto use 21-cm emission line to map the hydrogen contents of the universe acrosssome redshift range. 21-cm line arises from the hyperfine transition between thetriplet and singlet levels in the ground state of neutral hydrogen atom. Whensuch a transition occurs, a photon will be emitted (or absorbed) with a frequency휈 = 1420 MHz which corresponds to a wavelength of 휆 ∼ 21-cm. The probabilityof a spontaneous transition from triplet to singlet state in neutral hydrogen is verysmall and such transition can never be seen in the labs. However, the total numberof neutral hydrogen atoms in the intergalactic medium is so large that the emission(or absorption) line can be observed by radio telescopes. Expansion of the universechanges the the observed frequency 휈 of 21-cm line. So the 21-cm emission orabsorption from an object at redshift 푧 will be observed at:휈 =14201 + 푧MHz. (1.13)Therefore, there is no need for spectroscopy to obtain the redshift in this method.Instead, the frequency bin at which the sky is mapped in 21-cm directly gives the101.3. CHIMEredshift.The intensity of photons emitted through the hyperfine transition of neutralHydrogen atom is determined by the spin temperature 푇푠 which is defined by theequation:푛1푛2= 3 exp(− 푇∗푇푠), (1.14)where 푛1 and 푛2 are the number densities of electrons in the triplet and singletstates respectively, 3 is the degeneracy ratio between the triplet and singlet levels,푇∗ ≈ 0.068 K is the temperature corresponding to the energy difference between thetriplet and singlet states, and the spin temperature 푇푠 is the excitation temperatureof 21-cm line.We want to use the 21-cm line as a probe of a hydrogen gas cloud with opticaldepth 휏 along the line of sight. The radiative transfer equation in the Rayleigh-Jeanslimit for the gas along the line of sight is:푇푏 = 푇퐶푀퐵푒−휏 + 푇푠(1 − 푒−휏), (1.15)where 푇푏 is the observed brightness temperature. I have assumed that in thebackground there are only CMB photons. The interpretation of this equation iseasy: The CMB photons are exponentially attenuated by the gas cloud (hence wemultiply 푇퐶푀퐵 by 푒−휏) and the excitation temperature 푇푠 releases 21-cm photons.Then, part of the 21-cm emission is absorbed by the gas.The observed differential brightness temperature is:훿푇푏 =푇푏 − 푇퐶푀퐵1 + 푧=(푇푠 − 푇퐶푀퐵)(1 − 푒−휏)1 + 푧≈ (푇푠 − 푇퐶푀퐵)휏1 + 푧. (1.16)Therefore, 21-cm signal is observable only if the the spin temperature deviates fromthe background temperature (otherwise 훿푇푏 = 0). We can see the 21-cm signalas emission against CMB if 푇푠 > 푇퐶푀퐵 (훿푇푏 > 0) or absorption if 푇푠 < 푇퐶푀퐵(훿푇푏 < 0). In this case the 21-cm feature is seen as a spectral distortion to thebackground CMB and the diffuse distortion of the background can be studied in asimilar way as the CMB anisotropies. Then, observations at different frequenciesprobe different redshift slices of the observable universe.1.3.2 Radio Frequency InterferenceA phenomenon which is problematic not only for 21-cm intensity mapping experi-ments, but for any research in radio astronomy is the Radio Frequency Interference(RFI) which is the subject of this thesis. RFI is an unwanted human made signalthat interferes with the astronomical signal, degrades the quality of astronomical111.3. CHIMEdata and leads to data loss. RFI is usually non-Gaussian and come from TV sta-tions, cell phones, airplanes, satellites or any other activity which produces a radiosignal in a frequency that the radio telescope is working. RFI can be consideredas radio pollution in radio astronomy, similar to the light pollution for the opticalobservations.Some radio frequencies that are very important for astronomical research areprotected by regulations. For example, 1420 MHz (21-cm HI line). However, aradio telescope like CHIME covers a very broad bandwidth and such bandwidthcannot be fully protected. Therefore, radio astronomers have to deal with theRFI. There is no universal method for RFI mitigation. In the case of CHIME, weuse fourth statistical moment (kurtosis) to detect non-gaussianity in a signal. Thetheoretical foundation of the kurtosis based RFI mitigation is given in the chapter 2.1.3.3 CHIME instrumentThe Canadian Hydrogen Intensity Mapping Experiment (CHIME) is a transit radiotelescope which scans the northern sky across 400-800 MHz band. The telescopehas no moving part and the sky is scanned over a sidereal day as it is driftingeast to west. CHIME is located at Dominion Radio Astrophysical Observatory(DRAO) near Penticton, BC, Canada. It is originally designed to constrain thedark energy equation of state by measuring the BAO scale across the redshift range0.8 ≤ 푧 ≤ 2.5. This is the redshift range where the dark energy began to dominatethe total energy density of the Universe. We use the 21-cm emission line to mapthe hydrogen contents of the Universe over the required redshift range.CHIME consists of 4 cylinders with axes oriented along the north-south di-rection, each with a length of 100 meters and a diameter of 20 meters. There are256 dual polarization feeds along the focal line on each cylinder, i.e. 2048 inputscan receive the signal at the same time. The signal from the sky is focused by thecylinders onto the feeds and it is amplified by the low noise amplifier (LNA) whichis connected to the output of the feed. The analog signal then is transmitted by a50 m coaxial cable to a receiver hut, where it is further amplified by another setof amplifiers, filtered by a 400 − 800 MHz bandpass filter and sampled every 1.25ns and quantized to 8 bits. CHIME correlator has an FX design, i.e. the Fouriertransform is done before spatial cross multiplication. The F-Engine is responsiblefor sampling the analog signal with 8 bits, Fourier transforming every 2048 timesamples, and channelizing them by dividing the 400 MHz bandwidth into 1024 fre-quency channels. This process is done for all 2048 inputs. The real and imaginaryparts of the complex-valued data from each frequency channel are then separatelyquantized to 4 bits (4-bits real and 4-bits imaginary). The resulting 2.56 µs samplesfrom 2048 inputs fed into the X-Engine. There are 256 GPU nodes on X-Engine121.3. CHIMEFigure 1.3: CHIME telescope. It is an array of 4 cylinders in the east-west direction,with axes oriented in the north-south direction.where each node processes the data 4 frequency bins. In other words, a GPU nodedoes the cross multiplication of 2048 inputs for 4 frequency channels. Before thecorrelation operations, the real-time RFI excision system mask the RFI contami-nated samples by generating the spectral kurtosis estimates every 0.65 ms for eachfrequency bin . Note that every 0.65 ms one kurtosis estimate per frequency binis generated for the whole array (not for the individual inputs). Then the samplesare integrated to 30 ms and the signal of each input is correlated with the signalsof all other inputs to make 푁2 correlation matrix. Correlation is the process ofcross multiplication of each input with the complex conjugate of all other inputs.Then a second stage RFI excision is applied on 30ms frames to excise lower power,but longer duration RFI events. The theory of kurtosis based RFI excision and thedetails of the detection algorithms are described in chapters 2 and 3. After secondstage of RFI excision, samples are integrated for ∼10 seconds. The correlationproduct is called visibility, 푉푖 푗 :푉푖 푗(푡, 휈) = 〈퐸푖퐸∗푗〉, (1.17)131.3. CHIMEwhere푉푖 푗(푡, 휈) is the cross correlation of input 푖 with input 푗 at time 푡 and frequency휈, 퐸푖 and 퐸 푗 are the signals received by feeds 푖 and 푗 , respectively, and 〈.〉 standsfor the average over 10 second. The correlation products together with someother information are written to various HDF5 files. In this thesis I evaluatethe performance of the real time RFI excision system using the files containingautocorrelations (the correlation of one input with itself). The autocorrelation filesinclude the autocorrelations of all inputs over 10 s, excision rate for 10 s samples.These files are written to an HDF5 file every ∼ 43 minutes.1.3.4 RFI environment at DRAOThe mountains around the observatory at DRAO shield the site from RFI fromnearby cities. But a significant portion of the CHIME frequency band is stillcontaminated by satellites, airplanes, wireless communication and TV broadcastingbands. This includes LTE band between 730 MHz and 755 MHz range, a few TVstation bands between 500 MHz to 580 MHz, and a lot of RFI lines between 400MHz to 500 MHz, including UHF repeaters around 450 MHz. These features arevisible in the CHIME data (figure 1.4). Besides cell phone and TV station bandsthat are static in nature, there are many sources of intermittent RFI events such assatellites and airplane. Moreover, the atmosphere can be ionized by the meteorsentering to it. The ionized region of the atmosphere acts as a reflector for thedistant ground based RFI sources, such as distant TV stations. These scatteringevents typically appear as ∼6 MHz wide bursts and last for a few seconds. One cansee this type of RFI in the figure 1.4 between 480 MHz and 520 MHz around 00:15PT.141.3. CHIMEFigure 1.4: The autocorrelations of 30 ms samples for one feed. We can seevarious RFI features in this plot, namely LTE band from ∼ 730 MHz to 755 MHz,TV station bands around 500-600 MHz, and repeating RFI spots between 400 MHzand 500 MHz. 6 MHz wide bursts from distant broadcast TV bands appear as blobaround 00:15 PT in 480 MHz or 520 MHz. White lines correspond to the missingfrequencies due to the malfunctioning of some of the GPU nodes.15Chapter 2Kurtosis-based RFI excisionIn this chapter I review some definitions from statistics. Then the spectral kurtosisconcept is introduced and spectral kurtosis estimator is derived. I will show that thestatistical properties of the estimator can be used to detect non-Gaussian featuresof a signal in a Gaussian background. Most of this chapter is based on [34], [35],[36] and [1]. Note that in this chapter, I will use the words antennas, receivers andinputs interchangeably.2.1 Statistics2.1.1 Moments of a probability distributionMoments of a probability distribution are the expectation values of a randomvariable to integer powers and they often give valuable information about thedistribution of the random variable. The 푟 푡ℎ moment of a probability distribution푃(푥) of a random variable 푋 is defined as:퐸(푥푟 ) =∫∞−∞푥푟푃(푥)푑푥, (2.1)where 푥 is the value of the random variable 푋 . The first moment of a probabilitydistribution is called the mean 휇. The moments can also be derived from theso called moment generating function. Moment generating function of a randomvariable 푋 is defined as:푀퐺퐹푋 (푡) = 퐸(푒푡 푥), (2.2)where 푡 is a real-valued number. By expanding the exponential function and takingthe expectation value, it can easily be shown that the moments of 푃(푥) are relatedto the derivatives of moment generating function at t=0:퐸(푥푟 ) =푑푟푑푡푟푀퐺퐹푋 (푡)|푡=0 (2.3)The central moment is the moment around the mean:퐸 [(푥 − 휇)푟 ] =∫∞−∞(푥 − 휇)푟푃(푥)푑푥. (2.4)162.1. StatisticsClearly, the second central moment is the variance 휎2 of 푃(푥). Lastly, one candefine the normalized moment (commonly known as standardized moment) in thefollowing way:휇˜푟 =퐸 [(푥 − 휇)푟 ]휎푟, (2.5)where 휎푟 =(√퐸 [(푥 − 휇)2])푟is the 푟 푡ℎ power of the standard deviation of 푃(푥).The standardized moment facilitates the comparison of the shape of probabilitydistributions. For example, third standardized moment (skewness) is a measureof asymmetry of a distribution about its mean. The quantity of interest for thiswork is the fourth standardized moment, known as kurtosis. Kurtosis measuresthe heaviness of the tails of a distribution. From equation 2.5 the moment-basedkurtosis is:휇˜4 =퐸 [(푥 − 휇)4]휎4=퐸 [(푥 − 휇)4](퐸(푥2) − 휇2)2 .(2.6)One can evaluate the moment-based kurtosis for a Gaussian signal with mean zero,using the following Gaussian integrals:퐸(푥4) =1√2휋휎2∫∞−∞푥4푒−푥22휎2 = 3휎4, (2.7)퐸(푥2) =1√2휋휎2∫∞−∞푥2푒−푥22휎2 = 휎2. (2.8)Therefore, the moment based kurtosis of a Gaussian signal with mean zero is 3.2.1.2 CumulantAn alternative to the moments of a probability distribution is the cumulant. Similarto the moments, cumulants can be used to characterize the statistical properties of aprobability distribution. It is defined by the cumulant generating function 퐶퐺퐹(푡)which is the natural logarithm of the moment generating function:퐶퐺퐹(푡) = ln 퐸(푒푡 푥). (2.9)The power series of this function is:퐶퐺퐹(푡) =∞∑푟=1푘푟푡푟푟!, (2.10)172.2. Spectral kurtosiswhere the coefficients 푘푟 are called cumulants. Using equation 2.10 it can beshown that cumulant of 푟 푡ℎ order can be derived by taking the 푟 푡ℎ order derivativeof 퐶퐺퐹(푡) and evaluating it at 푡 = 0:푘푟 =푑푟푑푡푟퐶퐺퐹(푡)|푡=0. (2.11)The first three cumulants are same as the first three moments. However, onecan show that the cumulant-based kurtosis for a real random variable 푥 with meanzero is:푘4 = 퐸(푥4) − 3퐸2(푥2), (2.12)which is different fromequation 2.6. Using theGaussian integrals given in equations2.7 and 2.8 it can simply be shown that the cumulant-based kurtosis for a Gaussiansignal with mean zero is zero.2.2 Spectral kurtosisSpectral kurtosis (SK) is a statistical tool which can be used to identify the non-Gaussian behaviour of a signal in frequency domain. A recent cumulant-baseddefinition of the spectral kurtosis for a circularly symmetric random variable in thefrequency bin 푘 whose real and imaginary parts of the DFT have zero mean is givenby [34]:푆퐾[푋푘 ] =푘4(푥푘 , 푥∗푘)(푘2(푥푘 , 푥∗푘))2=퐸(|푥푘 |4) − 2퐸(|푥푘 |2)2퐸(|푥2푘 |2). (2.13)Note that the expression for 푘4 for a complex random variable is different from areal random variable (equation 2.12). Also, the spectral kurtosis is normalized bythe second order cumulant 푘2, which is same as the second order moment.To evaluate the spectral kurtosis for a circularly symmetric Gaussian signal withmean zero, note that 푥푘=푋푘+푖푌푘 . So, we will have:퐸(|푥푘 |2) = 퐸(푋2푘 ) + 퐸(푌2푘 ) = 2휎2, (2.14)퐸(|푥푘 |4) = 퐸(푋4푘 ) + 퐸(푌4푘 ) + 2퐸(푋2푘푌2푘 ) = 8휎4, (2.15)where I have used the Gaussian integrals (equations 2.7 and 2.8) and the fact that thereal and imaginary parts of the 푥푘 are independent, i.e. 퐸(푋2푘푌2푘 ) = 퐸(푋2푘 )퐸(푌2푘 ).Therefore, the spectral kurtosis for a circularly symmetric Gaussian signal withmean zero is zero.Although equation 2.13 gives the general definition for the spectral kurtosis of acircularly symmetric complex random variable with zero mean, it is much better torewrite it in terms of total power 푃푘 . This is because 푃푘 is an observable quantity.182.3. Spectral kurtosis estimatorThe total power 푃푘 in a spectral frequency bin 푓푘 is proportional to the absolutesquare of the frequency domain signal |푥푘 |2. So, equation 2.13 in terms of powerspectral density 푃푘 is:푆퐾[푋푘 ] =퐸(푃2푘) − 2퐸(푃푘)2퐸(푃푘)2. (2.16)Moreover, the variance 휎2푘 and mean 휇푘 of power spectral densities are:휎2푘 = 퐸(푃2푘) − 퐸(푃푘)2, (2.17)휇푘 = 퐸(푃푘). (2.18)So, equation 2.16 can be rewritten in terms of 휎2푘 and 휇푘 parameters:푆퐾[푋푘 ] = 푉2푘 − 1, (2.19)where 푉2푘 =휎2푘휇2푘is the normalized uncertainty or the spectral variability of thesignal. Equation 2.19 shows that the spectral kurtosis is equivalent to the spectralvariability up to a constant. So, we can use the two terms interchangeably. Theonly difference between 푆퐾[푋푘 ] and 푉2푘 for a Gaussian signal is that the spectralkurtosis of such a signal is zero from equation 2.13, but its spectral variability isone (equation 2.19). Since measuring the spectral variability is easier than spectralkurtosis and they are equivalent to each other, I use the term "spectral kurtosis" asa synonym for spectral variability in the rest of the thesis.2.3 Spectral kurtosis estimatorSuppose that we have a set of 푛 complex values 푥푘 which represent the post Fourier-transform time-stream for a single frequency channel. The power spectral densityestimator (푃̂푘) is proportional to |푥푘 |2. Based on the discussion in the previoussection, the unbiased spectral kurtosis estimator is defined as:푆̂퐾 ≡ 푉̂2푘 =휎̂2푘휇̂2푘, (2.20)where 휎̂2푘 and 휇̂푘 are the unbiased variance and mean estimators of the 푃̂푘 for thefrequency bin 푘 , respectively:휇̂푘 =1푛푛∑푖=1푃̂푘,푖 (2.21)휎̂2푘 =1푛 − 1푛∑푖=1(푃̂푘,푖 − 휇̂푘)2. (2.22)192.3. Spectral kurtosis estimatorIn order to simplify the final expression for the spectral kurtosis estimator, thefollowing parameters are defined:푆1 =푛∑푖=1푃̂푘 , 푆2 =푛∑푖=1푃̂2푘 . (2.23)Therefore, mean and variance estimators in terms of 푆1 and 푆2 are:휇̂푘 =푆1푛, 휎̂2푘 =1푛(푛 − 1)(푛푆2 − 푆21). (2.24)And the spectral kurtosis estimator 푆̂퐾 would be:푆̂퐾 =( 푛푛 − 1) (푛푆2푆21− 1)(2.25)The expected value of the SK estimator for a circularly symmetric Gaussian randomvariable must be one. But as it is shown by [36] this is not the case for the estimatordefined by equation 2.25:퐸(푆̂퐾) =푛푛 + 1. (2.26)Therefore equation 2.25 is a biased estimator. To get the unbiased estimator whoseexpectation value is one, one can rescale the estimator by multiplying it by 푛+ 1/푛:푆̂퐾 =(푛 + 1푛 − 1) (푛푆2푆21− 1)(2.27)Nita & Gary (2010) [36] have shown that for 푥푖 being drawn from a circularly-symmetric complexGaussian distribution, 푆퐾 estimator has the following statisticalproperties to the first order in 푛:퐸(푆̂퐾) = 1 (2.28)Variance(푆̂퐾) =4푛+ O( 1푛2). (2.29)Since the RFI events are mostly non-Gaussian [37], one can use the aboveproperties to detect andmitigate the RFI. Themean of the spectral kurtosis estimatoris invariant. So the sensitivity of the estimator to the RFI events is determinedby its variance. To see the statistical properties visually, I generated 30000 SKestimates, each from 20000 complex random samples drawn from a circularlysymmetric Gaussian distribution. Figure 2.1 shows the result of the simulation.The cumulative plot shows that half of the SK estimates are below 1 and the other202.4. Spectral kurtosis for a compact arrayFigure 2.1: Left: Histogram of 30000 realizations of SK estimator. Each SKvalue is generated from 20000 complex random variables drawn from a Gaussiandistribution. Red region indicates 1휎 interval where 휎 is the theoretical standarddeviation. Right: Cumulative histogram of the SK values. Yellow region corre-sponds to 1휎 interval. The SK distribution is symmetric around 1, and %68 of theestimates lie within one standard deviation. These two confirm the derived meanand variance of the SK estimator.half are above 1. This confirms that the mean value of the SK estimator is unity.Moreover, 68% of the SK estimates lie within 1휎 (yellow region), where 휎 is thetheoretical variance given by the equation 2.29. A non-Gaussian RFI changes theshape of the histogram, for example the tails of the histogram will be longer. Soone can simply detect non-Gaussian part of the signal by setting a threshold, e.g.at 5휎, and mask all the samples whose spectral kurtosis go beyond the threshold.2.4 Spectral kurtosis for a compact arrayEquation 2.29 shows that the discrimination potential of the kurtosis estimatordepend on the number of samples that are being used to generate the kurtosisestimates. To decrease the variance of the estimator one needs to integrate moresamples in a frequency bin, which in turns decrease the cadence at which thekurtosis estimates are being produced. The output of CHIME data are 10-secondsamples. Our aim is to flag RFI at a very short timescale (i.e. high cadence), sothat 10-second samples become RFI free. To increase the size of the input dataset212.4. Spectral kurtosis for a compact arraywithout increasing the timescale at which the kurtosis estimates are generated,one can combine the simultaneous measurements from multiple antennas in aninterferometer. Combining measurements from different antennas requires thefollowing two conditions: First, the array must be compact; i.e. the time it takesfor a signal to go across the array must be smaller than the timescale at whichthe kurtosis estimates are generated. Second, receivers must provide independentsamples. Correlated samples reduce the effective integration length and increasethe variance of the estimator. I discuss these two conditions for CHIME in section2.5.One way to combine the measurements from different antennas is to take theaverage spectral kurtosis of all antennas. For a compact array with 푁 independentantennas in which spectral kurtosis estimates are produced every 푛 time-samplesfor each antenna, the average spectral kurtosis can be found in the following way:Using equation 2.27, the spectral kurtosis for the i푡ℎ antenna is:푆̂퐾 푖 =(푛 + 1푛 − 1) (푛푆2,푖푆21,푖− 1)(2.30)where (푆1)푖 and (푆2)푖 are defined by the equation 2.23 for the i푡ℎ antenna. Then,spectral kurtosis estimator averaged over all antennas would be:푆̂퐾푎푟푟푎푦 =1푁푁∑푖푆̂퐾 푖 (2.31)Using the fact that the mean of the constituent 푆퐾푖 is unity, the average spectralkurtosis estimator for the whole array 푆̂퐾푎푟푟푎푦 is:〈푆̂퐾푎푟푟푎푦〉=1푁푁∑푖〈푆̂퐾 푖〉=푁푁= 1, (2.32)And since the total integration length is 푛푁 , the variance of the estimator is:Variance(푆̂퐾푎푟푟푎푦) =4푛푁+ O( 1푛2푁) (2.33)where I used the the statistical properties of constituents 푆퐾푖 from equations 2.28and 2.29. Comparing equations 2.29 and 2.33, it is obvious that by combining theestimates from all antennas, the variance of the estimator is reduced by a factor of1/푁 .There is an alternative formulation for the estimator to reduce the total numberof operations. The idea is to redefine 푆1,푖 and 푆2,푖 parameters by normalizing themby the mean power of the i푡ℎ receiver 휇푖:휇푖 =∑푛푗=1 푃푘, 푗푛=푆1,푖푛(2.34)222.5. Implementation for CHIMEThen, the new parameters 푆′1 and 푆′2 are defined by summing the normalized 푆1,푖and 푆2,푖 over all the receivers:푆′1 =푁∑푖=1푆1,푖휇푖=푁∑푖=1푆1,푖푛푆1,푖= 푛푁, (2.35)푆′2 =푁∑푖=1푆2,푖휇2푖=푁∑푖=1푆2,푖푆21,푖푛2 (2.36)Using equation 2.30 one can rewrite 푆′2 in terms of individual 푆̂퐾 푖:푆′2 = 푛푁∑푖=1(푛 − 1푛 + 1푆퐾푖 + 1)=푛(푛 − 1)푛 + 1푁∑푖=1푆퐾푖 + 푛푁. (2.37)Then the spectral kurtosis estimator will be:푆̂퐾 =푛푁 + 1푛푁 − 1(푛푁푆′2푆′12− 1)(2.38)=푛푁 + 1푛푁 − 1(푛푁푆′2(푛푁)2− 1)(2.39)This is obviously a biased estimator, because its expectation value is not unity. Theunbiased estimator can be found by a simple re-scaling:푆̂퐾 =푛 + 1푛 − 1( 푆′2푛푁− 1)(2.40)One can check that the expectation value of this estimator is one and its variance issame as the one in equation 2.33. The advantage of this alternative formulation isthat it requires slightly fewer number of operations to estimate the SK value [1].2.5 Implementation for CHIMEAs mentioned in the section 1.3.3, each of 256 X-Engine GPU nodes process fourfrequencies. The spectral kurtosis is estimated on each GPU node by accumulationof 256 time-samples from all inputs 푁 ∼ 2048 separately for each of the fourfrequencies . So, we would have one kurtosis estimate for each frequency channelevery 256 × 2.56휇s= 0.65 ms for the whole array. Since the expectation value andthe variance of the SK estimator for a Gaussian population are known, one can232.5. Implementation for CHIMEset a detection threshold to identify the samples that are contaminated by the RFI.Then signals from every feed for all 256 samples are removed when the SK valuelies out of bounds. Different RFI detection algorithms are discussed in detail inthe chapter 3. Note that the nominal variance in the case of CHIME (푛 = 256 and푁 ∼ 2048) is:휎2 ≈ 4푁푛≈ 7.63 × 10−3. (2.41)As already mentioned, two conditions are needed for combining the samplesfrom various antennas to increase the integration length: First, the cadence atwhich 푆퐾 estimates are generated must be shorter than the time the signal needsto travel across the array. This condition is satisfied for CHIME: The maximumdistance between the antennas is≈ 100m, so any signal arrives at all antennas within≈ 0.3휇푠. This is much smaller than the cadence at which the kurtosis estimatesare generated (0.65 ms). So the RFI appears simultaneous at all the antennas. Thesecond condition is that the samples from different inputs (that are combined toincrease the integration length) must be independent. This is normally the caseunless there is a bright astrophysical radio source in the CHIME beam or an RFI.This issue is discussed in detail in the section 4.2.24Chapter 3High Cadence ExcisionAlgorithmsIn this chapter, different excision algorithms that are used in CHIME pipeline areintroduced, and the performance of the algorithms is evaluated by applying themto an offline dataset.The idea for the high cadence excision is the following: The mean value ofCHIME visibilities are archived with 10 s cadence. But transient RFI events aremuch shorter than 10 seconds. They usually last from a few milliseconds to 1-2seconds. This means that a tiny fraction of a 10 second sample is contaminatedby the RFI, but they are so bright that they saturate the whole sample. Therefore,performing the excision at a high cadence can detect and remove those short events.So instead of removing the whole 10-second sample which is saturated by the RFI,we remove only the part which is being contaminated. Then the output of 10 secondsample which is free from RFI will be written on the disk.3.1 Excision algorithmsAs mentioned in the section 2.5, one SK estimate per frequency bin is generatedevery 0.65 ms by combining the measurements from different antennas. In chapter2 I showed that the expectation value and the variance of the SK estimator areknown for a Gaussian signal: the mean value of the estimator is invariant 〈푆퐾〉 = 1and the variance is 휎2 ∼ (0.0027)2 for CHIME. In general, a non-Gaussian signalcauses the SK estimator to deviate largely from the expected value with respect toits variance. So a non-Gaussian RFI can be detected by setting a threshold. If theSK estimate of the signal exceed that threshold, it will be flagged as RFI. In thefollowing I will describe different excision algorithms.3.1.1 Single StageIn the single stage algorithm, we mask all 0.65 ms samples whose SK value exceed1 ± 푛휎, where 푛 is the detection threshold. We set 푛 = 5, because according to[36], if the SK estimates are generated from Gaussian signals the distribution of253.1. Excision algorithmsSK estimates is also Gaussian. This means that the probability of the SK estimateof a Gaussian signal being beyond 5휎 limit is less than 10−6. So a sample whoseSK value is beyond 5휎 is most probably RFI.3.1.2 Excision on 30 ms frames (EOF)Suppose that we have∼500000 SK values that are estimated from aGaussian signal.Now I add a few values (representing RFI) that change the shape of the Gaussian.The situation is shown in the figure 3.1. Single stage algorithm is only sensitive tothe powerful RFI events whose SK values exceed 5휎 limit. But it is insensitive tothe lower power RFI events. Lower power RFI appears as non-Gaussian featuresthat are closer to the mean value; for example the bumps around 2휎 in the figure3.1. Obviously, if we lower the threshold of the single stage algorithm to detectlower power RFI, the Gaussian part of the signal will be lost as well. So, how canwe make an algorithm which is sensitive to those features without losing too muchof data?Figure 3.1: Histogram of ∼ 500000 random numbers drawn from a Gaussiandistribution with mean 1 and some variance 휎2. These values represent the RFI-free SK values. RFI events change the shape of the RFI-free SK distribution. So,a few non-Gaussian features are added to the the distribution to represent the RFI.The single stage algorithm is sensitive to the features beyond 5휎, while EOF issensitive to lower power, but longer duration RFI.263.1. Excision algorithmsOne way is to make a frame by accumulating 푀 individual 0.65 ms SK esti-mates, and if more than a few percent of the samples in the frame exceed somethreshold, mask the whole frame. In the CHIME pipeline, M=48 individual sam-ples are collected to make a 30 ms frame (48×0.65ms. ∼ 30ms). Now, if more thana few percent of the samples in the 30 ms frame exceed some threshold 1± 푛휎, thewhole frame is masked. So, in this algorithm a set of 2 parameters (푛, 푓 ) definesthe RFI detection threshold, where 푓 is the fraction of the samples in a 30 ms framewhose SK values exceed 1 ± 푛휎.To choose appropriate values for 푛 and 푓 parameters, I generated randomnumbers from a Gaussian distribution whose statistical properties are similar tothe free-RFI SK estimator for CHIME (i.e. a normal distribution with mean 1 andvariance (0.0027)2). Then, 48 of the samples are integrated to form a 30ms frame.A set of thresholds with different (푛, 푓 ) parameters are applied to the Gaussian 30ms frames, and the probability for the frames to pass the EOF algorithm withoutbeing masked is computed. For example, for (푛 = 2, 푓 = 0.3), I compute theprobability that 30 ms frames are not flagged by the EOF algorithm, that is, lessthan 30% of the samples in a frame (whose size is 48) exceed 1 ± 2휎. We try tokeep this probability high enough, for example around 99.9999%. In this case, only1 in a million of Gaussian frames will be masked. This probability is equivalent tothe probability for a Gaussian random variable to fall inside 5휎 interval. The resultis shown in the figure 3.2 in terms of the probability in the unit of 휎 interval. Wechoose the parameters so that they are equivalent to 7휎 threshold.3.1.3 Two stageThe two stage algorithm, is nothing but the combination of the single stage andEOF algorithms. In this algorithm, all the individual 0.65 ms samples whose SKvalue exceed 1 ± 5휎 are removed. So the first stage of excision is similar to thesingle stage algorithm. Then, 48 samples are integrated to form a 30 ms frame anda second stage excision with the EOF algorithm is performed on these frames.In comparison to the EOF, two stage algorithm is more sensitive to the veryshort RFI. Suppose that there are less than 13% samples in a 30 ms exceeding 3휎,in this case EOF does not remove any of those samples, but the first stage excisionof the two stage algorithm does remove those samples. It is also more sensitiveto RFI whose SK values are closer to the mean value. In other words, two stagealgorithm is both sensitive to the powerful, short time-scale RFI and to the lowerpower, but longer lasting RFI.273.2. Evaluation of excision algorithmsFigure 3.2: The probability for a Gaussian 30 ms frame to pass the EOF algorithmwithout being removed for different values of (푛, 푓 ) parameters. For example itis very probable for a 30 ms frame to not to be masked by EOF algorithm with(푛, 푓 ) = (3, 0.4). In this case the probability for the frame to pass the algorithm ismore than 8휎. The two red lines show two thresholds that we have used used forRFI detection between December 2019 and February 2020 using EOF algorithm.3.2 Evaluation of excision algorithms3.2.1 Gaussianity testGaussianity test is our main tool to assess the quality of the data and to see how closeare data are to being Gaussian distributed. It can be used to compare the gaussianityof the data before and after the RFI excision. Note that the signal received by anantenna is mainly composed of three components: astrophysical sky signal, thermalnoise and possible RFI. By thermal noise I mean the noise generated by the receiveritself, CMB radiation, galactic background radiation, atmospheric emission, andthe receiver noise. Thermal noise and the sky signal are both Gaussian with meanzero, however, the thermal noise varies in a much shorter time scale. The RFI283.2. Evaluation of excision algorithmsis usually non-Gaussian. Since sky is slowly varying in time, one can take thedifference between the neighboring time samples in autocorrelation products toremove the astrophysical sky contribution:∆푉푖푖(푡푛, 휈) = 푉푖푖(푡푛+1, 휈) −푉푖푖(푡푛, 휈), (3.1)where 푉푖푖 is the autocorrelation of input 푖, and 푛 is an integer representing thetime sample. Suppose that there is no RFI. Then what remains in ∆푉푖푖 is thethermal noise component. As already mentioned, thermal noise is Gaussian andthe fluctuation in thermal noise can be reduced by accumulating more samples,i.e. increasing the integration time ∆푡 in a bandwidth 퐵. Then, variance of ∆푉푖푖normalized by average autocorrelation is:12푉푎푟(∆푉푖푖푉 푖푖) =1푁, (3.2)where 푉 푖푖 denotes the average autocorrelation of neighboring samples and 푁 =∆푡 × 퐵 is the total number of samples in bandwidth 퐵 accumulated over time ∆푡.The 1/2 factor in front of the variance is because we are taking the difference oftwo random variables that are drawn from the same Gaussian distribution. A largedeviation from this relation is a signature of non-Gaussian RFI. Gaussianity test isthe process of checking our data to see if the equation 3.2 is satisfied. I define theGaussianity test value 퐺푇 in the following way퐺푇 =√푉푎푟(∆푉푖푖푉 푖푖)∆푡퐵2. (3.3)Note that for 퐺푇 = 1 we get the equation 3.2 which is for a Gaussian noise. Alarge deviation from 퐺푇 = 1 shows that the data is non-Gaussian.The Gaussianity test is our tool to evaluate the RFI excision system performanceand to compare different excision algorithms. For each excision algorithm, I runthe gaussianity test on all all the frequencies (and a single input) and by definition,the frequency channels whose 퐺푇 value is smaller than 2 are considered to begood frequency channels, i.e. nearly free from RFI. Therefore, we can comparethe number of good frequency channels for different algorithms and see whichalgorithm can detect and mitigate the RFI more effectively.3.2.2 Offline testTo evaluate the performance of excision algorithm, the variance and SK values of0.65ms samples for a single input were recorded on August 29, 2019, from 17:00293.2. Evaluation of excision algorithmsto 18:00 PT, and on October 10, 2019, from 00:00 to 10:00 PT. Different excisionalgorithms were applied to the data. Since the mean value of the voltage is zero(〈푣〉=0)), variance of the voltage gives the autocorrelation of the input:푉푎푟(푥) = 〈푥2〉 − 〈푥〉2 = 〈푥2〉, (3.4)where 푥 is the measured voltage. The average of SK values for samples drawn froma circularly-symmetric complex Gaussian distribution is 1. The standard deviationof the SK values are estimated by휎 ≈√4푁푛≈ 0.0027, (3.5)where 푛 = 256 is the number of time samples (2.56 µs cadence) that are collectedto generate spectral kurtosis every 0.65 ms, and 푁 ∼ 2048 is the number of inputsthat are used to estimate the SK values. Note that SK values are generated bycombining the signals from all the healthy inputs, and the variance files contain therms of the voltage for a single input. The files are processed in the following way:1. The excision algorithm is applied to the SK files. For example in the singlestage scenario, an upper and lower limit at 1±5휎푆퐾 is used to mask the SKoutliers.2. This mask is applied to the variance values to remove the samples exceedingthe upper and lower limits.3. The resulting variances are averaged over ∼ 8 seconds.4. Fraction of the excised data (those that exceed the upper and lower limits) arecomputed for every 8 seconds samples. I call this quantity the excision rate.5. Those samples that are affected by malfunctioning of the GPU nodes aredetected and removed.6. Variances and excised rate of 8 s time samples are reported.Unflagged-variances (i.e. no excision) are reported in the same way by doingsteps 3 to 6. To see the effect of excision algorithm I run the Gaussianity test onflagged and unflagged data (∼1 hour, 8 second samples) and compare the results.Figure 3.3 show the results of Gaussianity test on the data with the single stagealgorithm with a threshold at 5휎 (left) and with EOF algorithm with (푛, 푓 ) =(1.5, 0.4), and compare both with the case with no RFI excision. Each data pointcorresponds to one frequency channel and the color shows the fraction of excised303.2. Evaluation of excision algorithmsFigure 3.3: Left: Result of gaussianity test without any RFI excision (horizontalaxis) and with the single stage algorithm (vertical axis). Right: Gaussianity testwith EOF algorithm using (푛, 푓 ) = (1.5, 0.4). Each data point corresponds to onefrequency bin and the color shows the average excised rate at that frequency. Inother words, the figure shows three things for each frequency bin: GT value beforeRFI excision, GT value after applying the excision algorithm, and the averageexcised rate over ∼ 1 hour. The color bar is the same for both figures. Most of thedata points are below the diagonal line. This means that Excision algorithm hasimproved the Gaussianity of those frequencies.data in an 8 s sample. In figure 3.3, anything below the diagonal line has becomemore Gaussian after the excision. Note that there is no data point above the diagonalline, except a few that correspond to LTE and TV station bands (mostly blue andcyan points, see figure 3.4)). Before excision, ∼ 15% of the data points had a GTvalue less than 2. After excision, this increases to ∼ 34% of frequency channels(those that are below the horizontal gray line). In other words, we have recovered∼ 19% of frequency channels that were initially contaminated by RFI. In most ofthe frequency channels we can increase the quality of the data by excising less than0.5% of time samples (red points).Figure 3.4 shows the average excised rate as a function of frequency, and colorshows the difference in GT value between unflagged and flagged data . Largepositive difference shows an improvement on gaussianity of the data. Cyan andgreen points are UHF repeaters around 450 MHz. These frequencies are cleanedby excising < 25% of the samples and their GT values are decreased by a factorof 100 or more. These two figures together show the potential of high cadence313.2. Evaluation of excision algorithmsFigure 3.4: Average excision rate for 8 s samples over 1 hour as a function offrequency. 600-700 MHz band is the cleanest part of the spectrum. There aremany RFI features in 400-500 MHz band that. The kurtosis based RFI excisionsystem detected those RFI and some frequency channels are cleaned by excisingonly ∼ 20% of the samples in the integrated 8 second sample.kurtosis based RFI excision.To further evaluate the excision algorithms, the parameters of the excisionalgorithm are changed and the results are compared. The results for three algorithmsare shown in the tables 3.1, 3.2 and 3.3. It is clear that two-stage scenario canrecover more frequency channels without excising too much data. The results aresummarized in figure 3.5. For comparison, note that 288 frequency channels passedthe gaussianity test before excision and 387 passed the test using the single stagescenario (at 5휎) while 0.15% of the data are excised. Using two-stage algorithmwith 푝 = 0.4, 푛 = 1.4, 22 more frequency channels pass the gaussianity test withexcision of another 0.25%323.2. Evaluation of excision algorithmsn # of frequencies Excised fraction (%), Excised fraction (%),with 퐺푇 < 2 퐺푇 < 2 all frequencies5 387 0.15 11.634 394 0.19 12.353.5 398 0.31 12.833 400 0.81 13.692.5 408 2.58 15.652 408 7.38 20.27Table 3.1: the results of single stage excision for different values of parameter 푛on August 29 between 5:00 and 6:00 PM. Second column shows the number offrequency channels having 퐺푇 < 2 after excision. Third and fourth columns showthe average fraction of excised data for the frequency channels that passed the testand for all the frequency channels, respectively. As expected, we recover morefrequency channels by decreasing the threshold value at the expense of loosingmore data. Note that without any RFI excision 288 frequency channels have passedthe Gaussianity test.푛 푓 # of frequencies Excised fraction (%), Excised fraction (%),with 퐺푇 < 2 퐺푇 < 2 all frequencies3 0.13 348 0.35 15.372 0.25 341 0.19 15.751.5 0.4 327 0.14 15.83Table 3.2: Results of excisionwith EOF algorithm. Third column shows the numberof frequency channels having퐺푇 < 2 after excision. The number of good frequencychannels (퐺푇 < 2) is smaller than the single stage algorithm. This is because thesingle stage algorithm is more sensitive to the very short RFI events where onlya few samples exceed the detection threshold. EOF algorithm is sensitive to thelonger RFI, when more than a few percent of the samples in 30 ms exceed thethreshold, but it is insensitive to the very short RFI. The meaning of the third andfourth columns is similar to the ones in the table 3.1.333.2. Evaluation of excision algorithms푓 푛 # of frequencies, Excised fraction(%), Excised fraction(%),퐺푇 < 2 퐺푇 < 2 all frequencies0.41.5 406 0.23 15.691.4 409 0.39 16.241.3 414 1.18 17.411.2 417 4.16 20.541.1 419 12.75 28.590.51.5 401 0.17 14.851.4 403 0.18 15.051.3 405 0.20 15.331.2 410 0.29 15.771.1 413 0.85 16.72Table 3.3: The results of gaussianity test for the two stage algorithm. I used 5휎threshold for the first stage, and different values of (푛, 푓 ) for the second stage. Themeaning of the third and fourth columns is similar to the ones in the table 3.1.Figure 3.5: A comparison between two excision algorithms with different 푓 param-eters. The number above data points shows the average percentage of excised datain good frequency channels.343.2. Evaluation of excision algorithms3.2.3 Comparison of different excision algorithmsIn this section, single-stage, two-stage and EOM algorithms are applied to the 10hour dataset recorded between 12:00 and 10:00 AM on October 11. Note that inthe single stage scenario 푛 = 5 is used, and in the two-stage case I used the fixedvalues 푛2 = 1.4 and 푓 = 0.4 as the excision parameters. Then, the gaussianitytest is performed on 8-s samples every ∼ 1.5 hour and those frequencies whose 퐺푇value is less than 2 are labeled as good frequencies. The fraction of excised dataaveraged over all frequency channels for the single-stage and two-stage excisionalgorithms are ∼ 12% and 17%, respectively. The average excised fraction for goodfrequencies is ∼ 0.2% for the single-stage and ∼ 0.5% for the two-stage algorithm.Average excised fraction in all frequencies for EOM algorithm is ∼ 17% and ∼0.5% for good frequency channels (similar to two-stage).Figure 3.6: The performance of single stage (black), EOF (green) and two stage(red) algorithms over ∼ 8 hours. In general, two stage and EOF outperform thesingle stage. The difference between two stage and EOF is not significant becausethis test was performed at night, when there are no strong source of transient RFI.If there are powerful, short time scale RFI, the first stage of two stage algorithmcan remove it, while EOF is not able to do so.35Chapter 4SK estimator biasesThe SK estimator and its statistical properties discussed in chapter 2 are derivedfrom floating point Gaussian samples. Moreover, the value of the variance ofthe estimator for the CHIME array (equation 2.41) is based on the assumptionthat all the samples from different inputs are independent. These are all idealisticassumptions. The samples from which the SK value is estimated are not in the formof floating point values, because they are truncated to 4 bits before transferringto the GPUs. Furthermore, we might have a very bright celestial signal, e.g. theSun, that is measured by all of the antennas. In this case the signals measured atdifferent inputs are not totally independent. In this chapter I discuss the effect of4+4-bits truncation and the presence of a common mode signal on SK estimatesand on the variance of the estimator, respectively. I also introduce a correction forthe truncation bias through a python simulation.4.1 4+4 bits truncationSpectral kurtosis estimator is sensitive to the shape of the random variable dis-tribution. In our case, the random variable is the voltage and if there is no RFIcontamination, it follows a Gaussian distribution with mean zero. So accordingto the spectral kurtosis estimator properties discussed in the previous section, itsexpected value must be unity with a known variance. However, if for some reasonthe shape of the distribution changes, it could change the spectral kurtosis estimates.As I mentioned in the section 1.3.3, the real and imaginary parts of the voltageare separately truncated to 4-bits. Suppose that we have a Gaussian backgroundnoise with no RFI contamination. If the signal is truncated, then the shape of thesignal could change. This is because our samples are not floating point valuesanymore after the 4+4-bits truncation (4 bits real and 4 bits imaginary), but they areintegers in the interval [-7,7] (see figure 4.1). In this case, our SK estimates willdeviate from the expectation value and the RFI detection algorithm will mask thesignal. This is obviously not desired, as we don’t want to mask the celestial signal.I made a simulation in python to see the effect of 4+4-bits truncation on spectralkurtosis estimates:364.1. 4+4 bits truncation• A set of 푛푁 = 256 × 2048 complex random variable are drawn from acircularly symmetric Gaussian distribution with a known rms.• Real and imaginary parts of the random variable are separately truncated to4-bits. This means that they are rounded to their nearest integer and boundedbetween [-7,7]. All the values greater than 7 or smaller than -7 are replacedby 7 or -7, respectively. This is shown in the figure 4.1.• The spectral kurtosis with and without truncation are estimated.• This process is repeated 6000 times. To see the effect of truncation, themean and the variance of 6000 spectral kurtosis estimates with and withouttruncation are compared.Figure 4.1: Effect of 4-bits truncation on the shape of a Gaussian signal with anrms of ∼ 5.5.Looking at the figure 4.1, one can see that 4+4-bits truncation changes the shapeof the distribution significantly if the rms is high. Since kurtosis is a measure of theshape of the distribution, truncation introduces a bias in the estimator, because thespectral kurtosis estimates from truncated samples deviates from the expectationvalue. This is shown in the figure 4.2. The SK estimates generated from digitizedsamples tends to be smaller than those generated from non-digitized samples. Sotruncation introduces a negative bias in the SK estimates. This is expected, because374.1. 4+4 bits truncationat high rms, voltage samples beyond [-7,7] are replaced by -7, 7 which reduces the푆1 and 푆2 estimates, and 푆2/푆21 in the equation 2.27. For high rms values, theeffect of bounding of the samples to [-7,7] interval on 푆2/푆21 estimate is more thanrounding them to the nearest integer. For lower rms values (rms<2.2), there arenot many samples beyond [-7,7], but rounding of the samples to the nearest integerslightly reduces the 푆2/푆21 estimates. So truncation always leads to a negative biasin the SK estimates.Figure 4.2: Effect of 4+4-bits truncation on SK estimates. The SK estimates aregenerated from 6000 Gaussian dataset with rms ∼ 2.9. The SK estimates tend tobe smaller after truncation.Since the effect of truncation on the shape of a distribution is more significantfor high rms values (as there will be more sample with a value exceeding theinterval [-7,7]), we should expect a correlation between the rms and the bias in theSK estimate. Figure 4.3 shows the relationship between the truncation bias and therms of truncated samples in terms of nominal 휎 for CHIME. Note that the averagerms of the samples in the absence of bright sources in the sky is around 2.2. So theSK estimator is always biased, although the bias is less than 0.5휎 when there is nobright source in the sky.384.1. 4+4 bits truncationFigure 4.3: Bias of the SK estimator as a function of the rms of truncated samples.Each data point shows the deviation of the average SK value (over 6000 Gaussiandataset) from unity for a given rms.4.1.1 Effect of truncation bias on CHIME dataOn January 2020, we noticed an increase in the excision rate during the transit ofCygnus-A (hereafter Cyg-A) and Cassiopeia-A. This was obviously not desirable,as the RFI system should never excise the astrophysical signals. To understand thisfeature, spectral kurtosis values were dumped during Cyg-A transit on February2020. The goal was to study the spectral behaviour of the SK values during thetransit of bright sources, and to see how much the SK values deviates from theexpectations. Figure 4.4 shows the deviation of SK the values from unity duringthe Cyg-A transit. In the ideal case, we shouldn’t see the celestial sources in sucha plot, as they are Gaussian and SK estimator is not sensitive to them. However,we can see the trace of Cyg-A as a column which is darker than its neighbourhoodaround 9:29 AM. The lower part of the plot is darker during the transit, whichmeans that the negative SK fluctuations occur more at lower frequencies. In orderto see how much data are being flagged by the RFI system at lower frequencies,the average excised rate between 400MHz and 500 MHz is estimated. To find theaverage excised rate, an extreme threshold (푛, 푓 ) = (1.5, 0.4) was applied on 30msframes (see section 3.1.2). Then the number of samples exceeding our threshold isaveraged between 400 MHz and 500 MHz, and it is normalized by the total number394.1. 4+4 bits truncationFigure 4.4: Bias of SK values in terms of the nominal 휎 for CHIME. The SK valuesare more negatively biased during the transit of Cyg-a which is around ∼ 9:29 AM.Note that lower frequencies are more negatively biased than the higher frequencies.of samples in a 2-seconds frame. We can see this effect in the figure 4.5. Theincreases in the excision rate at the transit time of the Cyg-A is about %4. In otherwords, our kurtosis based RFI detection system mitigate %4 of useful data at lowerfrequencies during the transit of Cyg_A.On the other hand, we know that Cyg-A is brighter at lower frequencies and Ialso showed that 4+4-bits truncation introduces a negative bias in the SK estimator.This negative bias increases as the rms of the data increases. So, the observedcorrelation between the transit of Cyg-A and SK deviation from the expected value(which is higher at lower frequencies) is a hint that the increased excision rateduring the transit of bright celestial sources can be explained by the truncation bias.So correcting the estimator for the truncation bias should solve this issue.404.1. 4+4 bits truncation4.1.2 Correction for truncation biasThe truncation bias is defined as the deviation of the SK value from unity. Asdiscussed in section 4.1, the truncation bias depends on the rms of the samplesfrom which the SK value is estimated. Our goal in this section is to derive thetruncation bias as a function of rms. Suppose that we have a set of Gaussian datawith mean zero and known rms. If there is a table to relate the truncation bias withthe rms, the SK estimator can be corrected for the truncation bias for a given rms.To do so, I simulated 1000 SK estimates, where each estimate is generated from푁푛 = 2048 × 256 complex samples drawn from a circularly symmetric Gaussiandistribution with some rms. Then, I take the average of 1000 SK realizations,the rms of the truncated distribution is changed and the average SK value is re-estimated. The simulation steps are similar to what described in the section 4.1.The result is shown in the figure 4.6 (right).Once we have SK estimates as a function of rms, we can fit a polynomial to itto find the bias for a given rms. The bias for a given rms is:bias = SKfitted − 1 (4.1)and the SK estimate corrected for the truncation bias is:SKcorr = bias + SK, (4.2)Figure 4.5: The excised rate during Cyg-A transit . Each data point represents thefraction of excised samples in a 2 seconds frame averaged between 400 MHz and500 MHz.414.1. 4+4 bits truncationFigure 4.6: Left: Fitting residuals in terms of nominal 휎 as a function of rms fordifferent polynomials. Right: Spectral kurtosis estimated from truncated samplesas function of the rms of the truncated samples.where SK is the spectral kurtosis measured by the estimator given in the equation2.40. In other words, we continue to use equation 2.40 to estimate the SK value.This is a biased estimate. So we add one more step to correct the truncation bias:Having the instantaneous rms of the samples, we evaluate the bias for that rms usinga polynomial, and then it is added to the measured SK value. The result is the SKestimate corrected for the truncation bias.The major source of error in this process is the polynomial fit. Figure 4.6 (left)shows the residuals of the fitting in terms of the nominal sigma for 4 differentpolynomials. The residuals are the difference between the fitted SK value (usingpolynomial) and the simulated SK estimates. The lowest order polynomial whichgives a good fit within the range of ±0.5휎 is of degree 5. Polynomials of degree 8and more gives a much better fit with a lower residual, but fitting with a polynomialof degree 5 is just enough to test the system and to seewhether our idea for correctingthe truncation bias works as expected.4.1.3 Results of truncation bias correction on CHIME dataWe recorded the SK values together with the variances of 0.65 ms samples simul-taneously to correct the SK estimates for the truncation bias. Because the mean of424.1. 4+4 bits truncationthe random variable 푥 is zero, variance directly gives the rms:푉푎푟(푥) =〈푥2〉= rms(푥). (4.3)Using instantaneous rms values and the polynomial function I derived for thefitting (of degree 5), one can estimate the instantaneous truncation bias every 0.665ms. Adding the truncation bias to the measured SK value gives the unbiased SKestimate. The bias of the SK estimator before and after the correction for truncationbias is shown in the figure 4.7. Cyg-A transit is around 4:45 AM in these figures.I applied the threshold (푛, 푓 ) = (1.5, 0.4) to 30 ms frames. White spots are the30 ms frames exceeding the threshold. We can still see the trace of the Cyg-A onthe SK estimates after the correction, however, the bias is getting smaller and mostof the samples do not exceed the threshold after the correction. Moreover, the SKestimator is still sensitive to the other RFI events including LTE band, TV stationbands, repeating RFI around 470 MHz, and transient RFI.The effect of truncation is expected to be more clear at lower frequenciesbecause the sky is getting brighter at those frequencies during the transit of thegalactic plane. So, zooming into the lower frequencies can give us a better sense ofthe performance of the SK estimator with andwithout the truncation bias correction.The average excision rate for 2-sec frames between 400 MHz and 500 MHz beforeand after correction are shown in the figure 4.8. The bump around 4:45 in thefigure 4.8 (left figure) corresponds to the Cyg-A transit, showing that the excisedrate increases at lower frequencies without bias correction. The bump disappearsafter correction with a polynomial of degree 5. We can still see a %2 increase inthe excised rate at the transit time, however, this is a much better than %10 increasebefore correction.Although the correction with a polynomial of degree 5 improves the SK estima-tor, for better results we need higher order polynomials or a different fitting method.An alternative to the polynomial fitting is the cubic spline interpolation. It turnsout that this method gives a much better result. Figure 4.9 shows the SK deviationfrom the expected value in terms of nominal 휎 after correcting the SK values with acubic spline function during the transit of Cyg-A. The SK values are much closer tounity in this case, and the signature of Cyg-A in the plot is disappeared. Moreover,the SK values are much more uniform across the spectrum. The excision rate doesnot change with the cubic spline function, but in general it gives a much better fit tothe SK v.s rms plot. Currently, we are using the polynomial fitting in the pipeline.434.1. 4+4 bits truncationFigure 4.7: The bias of SK estimates in terms of nominal 휎 for CHIME without(Left) and with (Right) the correction for truncation bias.Cyg-A transit is around4:45 PT in this plot. Each data point is the average SKvalue over 30ms and thewhitespots are 30 ms frames that are masked by EOF algorithm with (푛, 푓 ) = (1.5, 0.4).The color scale is the same in both plots.Figure 4.8: Left: The excised fraction of the samples in an 2 second frame, averagedbetween 400 MHz and 500 MHz during Cyg-A transit (∼ 4:45 PT) without trunca-tion bias correction. Right: The same quantity as in the left figure, but corrected fortruncation bias with a polynomial of degree 5. The bump corresponding to Cyg-Adisappeared, i.e. the RFI system no longer excise the Cyg-A.444.1. 4+4 bits truncationFigure 4.9: Deviation of SKvalues fromunity after correcting the kurtosis values fortruncation biaswith a cubic spline function duringCyg-A transit around∼ 4 : 45푃푇 .The bias is closer to zero for all frequencies compared to the figure 4.7.454.2. Common mode signal4.2 Common mode signalWe derived the variance of the SK estimator for a compact array in the section 2.4assuming that the signals received by different receivers are independent. In generalthis is not the case, as the sky signal contributes to all the receivers. So there willbe common mode signal and this reduces the total number of independent samples,and increases the variance of the estimator. In other words, the sensitivity of thekurtosis-based RFI mitigation system which depends on the variance of the SKestimator is decreased in presence of a common mode signal.A simulation can reveal the effect of common mode signal on the variance ofSK estimator. Suppose that the receiver noise is described by the matrix N. A setof 256 × 2048 complex samples is selected randomly from a Gaussian distributionwith mean zero and variance 1. I put these samples in a 2 dimensional matrixN. This matrix has two axes: time axis which is represented by 256 samples; andreceiver axis representing 2048 inputs. Each row is independent from the otherone. This matrix represents the receiver noise of each input for 256 time samples.The component 푁푖 푗 represents the signal received by the 푖푡ℎ receiver at time 푗 . Soeach row is the time-stream for a single receiver (with a length of 256). Note thatthe noise of different receivers are independent from each other, and each of themfollows a circularly symmetric Gaussian distribution.푁 =푁1,1 푁1,2 ... 푁1,256푁2,1 푁2,2 ... 푁2,256... ... ... ...푁2048,1 푁2048,2 ... 푁2048,256 (4.4)Then, a common mode signal is added to the noise. It is described by thematrix C whose shape is same as matrix N. The only difference is that C follows acircularly symmetric Gaussian distribution along the time axis, but it is the sameover all of the receivers:퐶 =퐶1,1 퐶1,2 ... 퐶1,256퐶2,1 퐶2,2 ... 퐶2,256... ... ... ...퐶2048,1 퐶2048,2 ... 퐶2048,256 =퐶1 퐶2 ... 퐶256퐶1 퐶2 ... 퐶256... ... ... ...퐶1 퐶2 ... 퐶256Then the signal is:푆 = 휎푁푁 + 휎퐶퐶 (4.5)464.2. Common mode signalwhere 휎푁 is the noise amplitude and 휎퐶 is the common mode amplitude. Then,the SK value of the signal is estimated, and the variance of the estimator for a setof 1000 dataset is calculated. This process is repeated for different values of 휎퐶 .The fractional common mode amplitude shows how much of the amplitude of thesignal is dominated by the common mode:휎frac =휎퐶휎푡표푡푎푙=휎퐶√휎2푁 + 휎2퐶. (4.6)Figure 4.10 shows that the effect of the common mode signal on the variance ofthe SK estimator is negligible if the fractional common mode amplitude is lessthan 0.3. The signal is dominated by the common mode when 휎frac ≈ 1. In thiscase there are only 푛 = 256 independent time samples and the variance convergesto 4/푛 ≈ 0.015. Since the increase in the common mode amplitude reduces thenumber of independent samples, one can define the effective number of receiversusing the fact that for a compact array 푁 ≈ 4푛×푉 푎푟( ¯푆퐾 ) . When the fractionalcommon mode amplitude is unity the effective number of antennas becomes 1, asexpected. This is shown in the figure 4.11. Red horizontal line in both figurescorresponds to a signal that is totally dominated by a common mode.Note that the common mode signal does not alter the mean of the estimator,because mean of the SK estimator is an invariant quantity and does not depend onthe integration length. Therefore common mode is not counted as a bias. But itreduces the sensitivity of the SK estimator to the RFI as it reduces the number ofindependent samples.474.2. Common mode signalFigure 4.10: Variance of the SK estimator in presence of a common mode signal.The variance reaches to 10−2 if the signal is totally dominated by the commonmodesignal, as there are only 푛 = 256 independent time samples.Figure 4.11: Number of effective antenna as a function of fractional common modeamplitude.48Chapter 5ResultsIn this chapter, I summarize the results of the real-time RFI excision from May2019 to September 2020. In this period we have used all three excision algorithmsdescribed in the chapter 3. The final part of this chapter is a discussion aboutthe reliability of the Gaussianity test for the evaluation of the RFI system, theperformance of the system during the solar transit, and a brief discussion about theeffect of polarized RFI on the sensitivity of the RFI system.5.1 May to November 2019: Single stage algorithmWe started to use the kurtosis-based RFI detection system in May 2019. FromMayto November 2019, single stage algorithm with a detection threshold of 5휎 wasused. The estimator was not corrected for the truncation bias in this period.To compare the quality of the data before and after the RFI excision visually,I plotted the stacked autocorrelations normalized by the median across time foreach frequency. Stacked autocorrelations are the autocorrelations averaged over allinputs. Figure 5.1 shows the waterfall plot of autocorrelations before and after westart to use the kurtosis based RFI excision system. I zoomed in the range 400 MHzand 500 MHz, because this is the range that we see many repeating RFI lines. Onecan see that there are a few RFI lines in April 2019 (without RFI excision) that areabsent in May 2019 (after RFI excision). Most of these lines are not transient, forexample the repeating line around 468 MHz.To be more quantitative, I plotted the average excision rate as a function offrequency. Figure 5.2 shows the excision rate in each frequency bin averaged over1.5 hours. The shape of this plot is very similar to the figure 3.4. This shows thatthe excision algorithm is working as expected. Most of the RFI lines appearingbetween 400 and 500 MHz, can be recovered by removing less than 20% of thedata. On a 10-second cadence, this repeating line appears every 20 seconds. Butthe high cadence RFI excision system can clean the frequency without excising toomuch of the data. Figure 5.2 shows that this frequency can be recovered by excisingonly 8% of the data in a 10 second sample at this frequency.495.1. May to November 2019: Single stage algorithmFigure 5.1: Left: Normalized autocorrelations in April 10, 2019 without any RFIexcision. Right: Normalized autocorrelations in June 7, 2019 after deploying highcadence RFI excision system with single stage algorithm. Most of the repeatingRFI features disappeared.Figure 5.2: Excised fraction of 10 s samples averaged over 1.5 hours as a functionof frequency.505.1. May to November 2019: Single stage algorithm5.1.1 Moving blobIn July 2019, we detected a mysterious signal that was being removed by the RFIsystem. The signal was repeating every day with a regular pattern . Preliminaryevaluations showed that the signal appears around 480 to 490 MHz and it is lockedin sidereal time, i.e. it appears at the same time in each sidereal day, just like stars.Being locked in sidereal time was a hint that the signal source is not ground based.However, tracking the signal over a month showed that the blob is slowly driftingin sidereal time. Figure 5.3 shows how the blobmoves back in time. Each diagonalline corresponds to the excised rate at 483.98 MHz for 24 hours. We can see theincrease in the excised rate, two times per day, at specific right ascensions. Notethat CHIME is a transit radio telescope. So the right ascension of an object is equalto the local sidereal time when the object is in the CHIME beam. The questionwas whether this is a true RFI or the RFI system is removing a celestial source bymistake. There was no bright radio source at those right ascensions. Finally, withthe help from Scott Tilley, an amateur satellite observer, it turned out that the signalis coming from a Russian Meridian satellite and the RFI system did a good job indetecting the RFI.Figure 5.3: The signature of the Meridian satellite is visible in the plot of excisionrate of 10 s samples (in percent) as a function of time around frequency ∼484 MHz.Each diagonal line corresponds to 24 hours. The blob appears every day, but itslowly drifts in right ascension.515.2. December 2019 to February 2020: Excision on 30 ms frames with multiple thresholds5.2 December 2019 to February 2020: Excision on 30 msframes with multiple thresholdsThis section summarizes the overall performance of the RFI excision system fromDecember 2019 to February 2020. During this period the EOM algorithm with twothresholds (푛, 푓 ) = {(1.5, 0.4), (3, 0.13)} was used. This algorithm was initiallytested on SK/variance dataset in October 2019 (see section 3.2.2) .Figure 5.4 shows the excised fraction of 10 s samples (averaged over all fre-quencies) for 3 consecutive days. One can see that the excision rate averaged overall frequencies is around 16%, as expected from the analysis in chapter 3 for theEOF algorithm. We can also see 3 bumps which appear every day. Two of thesebumps are already identified as the meridian satellite (see section 5.1.1). Thesebumps are moving back in time every day. These three bumps have some similarcharacteristics: They are almost locked in sidereal day, and if one plots the excisedfraction as a function of frequency and time, they are all appearing as blobs between480 and 490 MHz (figure 5.8). However, unlike the previous blobs that are driftingslowly in RA, the new blob seems to be tightly fixed in RA, and it also seems to bemore broadband. Figure 5.8 shows that the new blob appears in 480-490 MHz and400-440 MHz at the time of Cyg-A transit. Therefore, it is obvious that the newblob has a different origin than the other two. In fact, the RFI system is excisingthe Cyg-A. This is not desirable and as already mentioned in chapter 4, it is due tothe 4+4-bits truncation. The effect of 4+4-bits truncation was not very clear withthe single stage algorithm because the detection threshold was at 5휎.To see how the excision work in different frequencies, I plot the average excisedfraction of 10 s samples at different frequencies over 1.5 hours. Figure 5.5 showsthat the LTE and TV station bands are completely masked by the EOF algorithm.However, the RFI lines around 450MHz are more excised than what the single stagealgorithm does. Note that removing more samples around 450 MHz by the EOFalgorithm does not necessarily mean that more RFI is removed. Higher excisionrate in those frequencies is mainly due to the fact that the EOF algorithm performsthe excision on 30 ms frames, while single stage algorithm masks the individualoutliers (0.65 ms).Another quantity of interest is the number of good frequency channels through-out the day. This is shown in the figure 5.6. Remember that good frequencychannels are those channels that pass the Gaussianity test (i.e. those with 퐺푇 <2).The maximum number of good frequency channels (∼ 700) occurs every nightbetween 2:00 and 4:00 AM. The minimum number occurs in the morning between8:00 AM and 12 PM. These numbers are consistent with the results of our test onsk/variance dataset.525.2. December 2019 to February 2020: Excision on 30 ms frames with multiple thresholdsFigure 5.4: Excised rate of 10 s samples averaged over all frequencies for 3 consec-utive days. The red box is the Cyg-A transit time and two green boxes correspond tothe meridian satellite. The excision rate increases every day during Cyg-A transit.This is due to 4+4-bit truncation.Figure 5.5: Average excised rate over 1.5 hours as a function of frequency for atypical day in February 2020. The excision algorithm completely masks the LTEand TV station bands.535.2. December 2019 to February 2020: Excision on 30 ms frames with multiple thresholdsFigure 5.6: Number of good frequency channels reported every 1.5 hours for 12consecutive days. We can see a repeatable pattern every day. Note that there arezero good frequency channels during solar transit, as the RFI system flags the Sun.This issue is discussed in section 5.4.2Figure 5.7: Change of the number of good frequency channels in 24 hours. Thenumber of good frequency channels is more at night and it goes down in themorning, when there are more RFI at site.545.3. March to September 2020: Two stage algorithmFigure 5.8: Excised fraction as a function of time and frequency. The green circlesand rectangle show the blob which is fixed in RA (which turned out to be Cyg-A).5.3 March to September 2020: Two stage algorithmAs of March 2020, we are using two stage excision algorithm: In the first stage, allindividual 0.65ms sampleswhose kurtosis exceed 5휎 limit aremasked. The secondstage excision is performed on 30 ms frames whose more than 13% of their samplesexceed 3휎. In June 2020, we corrected the SK estimator for truncation bias with afifth degree polynomial as described in chapter 4. The number of good frequencychannels as a function of time for 6 consecutive days in July 2020 is shown in thefigure 5.11. The maximum number of good frequencies is a little lower than thewinter run . This could possibly be due to the new RFI bands appearing around700 MHz and the change of excision algorithm and the threshold values. Figures5.9 and 5.10 show the effect of truncation bias correction on excision rate duringthe transit of Cyg-A. Since the truncation bias is more severe in lower frequenciesduring the Cyg-A transit, I averaged the excision rate over 400 and 430 MHz. Thebump in the figure 5.9 corresponds to the time of Cyg-A transit. This bump isdisappeared in June when the SK estimator was corrected for the truncation bias.555.3. March to September 2020: Two stage algorithmFigure 5.9: Excised rate of 10 second samples averaged between 400 and 430 MHzon May 19, 2020; without truncation bias correction. Red arrow shows the Cyg-Atransit time.Figure 5.10: Excised rate of 10 s samples averaged between 400 and 430 MHz onJune 16, 2020. The SK estimator is corrected for the truncation bias. The Thebump corresponding to the Cyg-A transit disappeared.565.4. DiscussionFigure 5.11: Number of good frequencies from July 9, 2020 to July 14,2020. Twostage algorithm was used in this period with 5휎 threshold on individual 0.65 mssamples and (3휎, 0.13) on 30 ms frames. Zeros correspond to the solar transit,when the RFI excision system is off.5.4 Discussion5.4.1 Reliability of Gaussianity testIn this analysis, the aim of gaussianity test was to compare the quality of the databefore and after the excision by identifying those frequency channels which areless contaminated by RFI. To see whether the gaussianity test selects only suchchannels, the standard deviation of the data over an hour or so can be used. The RFIshould appear as a spike in the standard deviation plot. So, if any RFI-contaminatedchannel pass the gaussianity test, we must see it as a spike in the plot of standarddeviation versus frequency.Figure 5.12 is the standard deviation of the frequency channels that have passedthe gaussianity test. Clearly, there are a few outliers: These frequencies have passedthe gaussianity test, but they have a large standard deviation and indeed they areRFI.Why the outliers pass the gaussianity test? Are they gaussian-like RFI? Notnecessarily. In the gaussianity test, we subtract the sky contribution by taking thedifference between neighbouring samples. But if the RFI is stable in time (over∼ 10seconds), the subtraction removes the RFI and the frequency channel might pass575.4. Discussionthe test. For example, figure 5.13 shows the histogram of one of the outliers beforeand after subtraction of neighbouring samples. After subtraction, the distributionlooks more Gaussian. This behaviour is not so much common, but it should betaken into account.Figure 5.12: Standard deviation of the data between 3:00 and 4:00 AM for thosechannels which pass the gaussianity test. Two of the outliers are shown with blackcircles.Figure 5.13: Left: Histogram of the autocorrelations one of the outliers at frequency419.14 MHz. Right: Histogram after subtracting the neighbouring samples.585.4. Discussion5.4.2 Solar transit issueThe RFI excision rate is extremely high during the solar transit. Ideally this shouldnot happen, because the Solar data are Gaussian and the RFI system should notremove the Gaussian signal. However, the two effects mentioned in the chapter 4affect the SK estimates and its variance. The truncation bias is very significant forthe sun, as it is the brightest object in the sky. The average rms of the solar data isgreater than 4. According to the simulations explained in the chapter 4, the bias ofthe SK estimates due to the truncation effects at rms∼ 4 is more than −40휎 . Thismeans that all the excision algorithms (single stage, EOM or two stage) will excisealmost all of the solar data. However, this issue persists even after the correction ofSK estimates for truncation bias. This is most probably due to the common modesignal: the common mode amplitude approaches to the unity in the case of solartransit. This reduces the number of independent samples from which the kurtosisestimates are generated. So the variance of the estimates will increase and the SKvalues exceed the detection thresholds. To avoid loosing solar data, the RFI systemmust be turned off during the solar transit.Figure 5.14: Deviation of SK values from unity in terms of nominal 휎 for CHIMEduring solar transit (around 13:00). Left: Without truncation bias correction SKvalues have a very high negative bias. Right: Since SK values sill highly deviatefrom expected value after truncation bias correction with a fifth degree polynomial,the Sun is still being flagged by the RFI system. This might be explained by thefact that the variance of the estimator is increased when the common mode signaldominates.595.4. Discussion5.4.3 RFI polarizationRunning theGaussianity test on all of the inputs reveals that in a few frequencies, onepolarization is more Gaussian than the other. This is possibly due to the polarizedRFI events. Since kurtosis is estimated by combining the time samples from all ofthe inputs, it can lead to a false estimate of RFI at some frequencies . The followingfigures show the result of gaussianity test over different inputs. It is clear that onepolarization is more non-Gaussian than the other one. One way to deal with thepolarized RFI is to estimate the kurtosis for different polarizations separately. Thismeans that instead of 2048 inputs, we will have 1024 inputs to combine. However,reducing number of inputs increases the variance of the estimator, and decreasesthe sensitivity of the excision algorithms to the RFI. Since it is important to havethe highest possible number of samples to keep the variance of the SK estimator aslow as possible, we ignore this effect for now.Figure 5.15: Result of the Gaussianity test for all feeds at a single frequency. TheNS and EW corresponds to two different polarizations. Note that the first 256 inputson each cylinder have a different polarization than the next 256.60Chapter 6Summary and future workIn this thesis I evaluated the performance of the high cadence kurtosis based RFIexcision system for CHIME. In general, RFI excision system is able to detect andflag many types of RFI, from repeating RFI lines between 400 MHz and 500 MHzto intermittent RFI from satellites and airplanes. A good example for this is theexcision of the Meridian satellite signal. Moreover, the RFI system detects andmasks the static RFI bands such as LTE and TV station bands. It is also possibleto recover a few frequency channels that were totally non-Gaussian by excisingless than %20 of millisecond samples in a 10 second sample. Before excision,these frequency channels were totally contaminated by the RFI. I also discussed theeffect of quantization on SK estimates and I corrected the SK estimator for 4+4 bitstruncation bias using polynomial fitting and cubic spline interpolation.Although the RFI detection system detects powerful non-Gaussian RFI events,it must be improved in many ways. First of all, we have to find new detectionthresholds and design new algorithms to increase the sensitivity of the RFI system;for example to excise 6 MHz TV channels more effectively. Moreover, truncationbias correction with a polynomial fitting must be replaced by the cubic splineinterpolation or by a polynomial of higher degree, as I showed that the cubicspline function gives a much better correction for the truncation bias, and fittingresiduals for a polynomial of degree 8 are much smaller than those for 5푡ℎ degreepolynomial. Another issue with the RFI algorithms that must be fixed is the solartransit issue. The number of independent samples from different inputs is reducedduring the solar transit. Therefore, the variance of the SK estimator increases andour detection algorithms flag the Sun. To avoid loosing the solar data, we turn offthe RFI system during the solar transit. Currently, there is no resolution for thisproblem, but another potential improvement of the system includes designing analgorithm which does not flag the Sun.61Bibliography[1] Taylor J., Denman N., Bandura K., Berger P., et al., 2019, JAI, 8, 1940004[2] Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009[3] Schmidt, B. P., Suntzeff, N. B., Phillips, M. M., et al. 1998, ApJ, 507, 46[4] Planck Collaboration, Ade P. A. R., Aghanim N., Arnaud M., et al., 2015,A&A 594, A13[5] Phillips, M. M., 1993, ApJ, 413, L105[6] Hamuy M., Phillips M. M., Suntzeff N. B. , et al. 1996, AJ, 112, 2398[7] Regnault, N., Conley, A., Guy, J., et al. 2009 A&A, 506, 999[8] Fukugita, M., Ichikawa, T., Gunn, J. E., et al. 1996, AJ, 111, 1748[9] Frieman, J. A., Bassett, B., Becker, A., et al. 2008, AJ, 135, 338[10] Sako, M., Bassett, B., Becker, A., et al. 2008, AJ, 135, 348[11] Sako, M., Bassett, B., Becker, A. C., et al. 2014, arXiv:1401.3317[12] Betoule, M., Marriner, J., Regnault, N., et al. 2013, A&A, 552, A124[13] Mosher, J., Guy, J., Kessler, R., et al. 2014, ApJ, 793, 16[14] Betoule. M., Kessler, R., Guy, J., et al. 2014, A&A, 568, A22[15] Planck Collaboration, Aghanim N., Akrami Y., Ashdown M., et al., 2020,A&A, 641, A6[16] Weinberg D. H., Mortonson M. J., Eisenstein D. J., et al., 2013, PhR, 530, 87[17] Landy, S. D. Szalay, A. S., 1993, ApJ, 412, 64[18] Eisenstein D. J., Zehavi I., Hogg D. W., et al., 2005, ApJ, 633, 560[19] Blake C., Davis T., Poole G. B., et al., 2011, MNRAS, 415, 289262Bibliography[20] Anderson L., Aubourg É., Bailey S., et al., 2014, MNRAS, 441, 24[21] Busca N. G., Delubac T., Rich J., et al., 2013, AA, 552, A96[22] Furlanetto S. R., Oh S. P., Briggs F. H., 2006, PhR, 433, 181[23] Peacock J. A., Schneider P., EfstathiouG., et al. 2006, arXiv:astro-ph/0610906[24] Bartelmann M., Maturi M., 2017, SchpJ, 12, 32440[25] de Jong J. T. A., Verdoes Kleĳn G. A., Kuĳken K. H., Valentĳn E. A.„Experimental Astronomy, 2013, 35, 25[26] Joudaki S., Mead A., Blake C., et al., 2017, MNRS, 471, 2, 1259–1279[27] Joudaki S., Mead A., Johnson A., et al., 2018, MNRS, 474, 4, 4894–4924[28] Martinet N., Schneider P., Hildebrandt, H., et al., 2018, MNRS, 474, 1,712–730[29] Troxel M. A., MacCrann N., Zuntz J., Eifler T. F., et al., 2018, Phys. Rev. D98, 043528[30] Aihara, H., Arimoto, N., Armstrong, R., et al., 2018, PASJ,70, SP1, S4[31] Hikage C., Oguri M., Hamana T., et al., 2019, PASJ, 71, 2, 43[32] Aihara, H., Armstrong, R., Bickerton, S., et al., 2018, PASJ, 70, S8[33] Pritchard J. R., Loeb A., 2012, RPPh, 75, 086901[34] V.D. Vrabie, P. Granjon, C. Servière, “Spectral kurtosis: from definitionto application”, IEEEEURASIP Workshop on Nonlinear Signal and ImageProcessing, Grado, Italy, June 8-11, 2003.[35] Nita G. M., Gary D. E., Liu Z., et al., 2007, PASP, 119, 805[36] Nita G. M., Gary D. E., 2010, PASP, 122, 595[37] Fridman P. A., Baan W. A., 2001, AA, 378, 32763
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- High cadence kurtosis based RFI excision for CHIME
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
High cadence kurtosis based RFI excision for CHIME Mirhosseini, Arash 2020
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | High cadence kurtosis based RFI excision for CHIME |
Creator |
Mirhosseini, Arash |
Publisher | University of British Columbia |
Date Issued | 2020 |
Description | This thesis describes the real-time Radio Frequency Interference detection system for the Canadian Hydrogen Intensity Mapping Experiment (CHIME). CHIME is a transit radio telescope located at Dominion Radio Astrophysical Observatory (DRAO) in Penticton, BC, and it is originally designed to map the large-scale structures in the redshift range 0.8 |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2020-10-28 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0394838 |
URI | http://hdl.handle.net/2429/76386 |
Degree |
Master of Science - MSc |
Program |
Astronomy |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2020-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2020_november_mirhosseini_arash.pdf [ 23.52MB ]
- Metadata
- JSON: 24-1.0394838.json
- JSON-LD: 24-1.0394838-ld.json
- RDF/XML (Pretty): 24-1.0394838-rdf.xml
- RDF/JSON: 24-1.0394838-rdf.json
- Turtle: 24-1.0394838-turtle.txt
- N-Triples: 24-1.0394838-rdf-ntriples.txt
- Original Record: 24-1.0394838-source.json
- Full Text
- 24-1.0394838-fulltext.txt
- Citation
- 24-1.0394838.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0394838/manifest