Loading...

Summer Research Fellowship Programme of India's Science Academies

Application of wavelets for the detection of the Global 21cm signal from the Epoch of Reionisation

Akshatha K V

R V College of Engineering, Mysore Road, Bengaluru 560059

Guided by:

Prof. Udaya Shankar N

Astronomy & Astrophysics, Raman Research Institute, Bengaluru- 560080

Abstract

According to the standard model of Cosmology, the expansion of the universe after the big bang resulted in its cooling and the drop in temperature allowed combination of neutrons and protons leading to the formation of atomic nuclei. The initial inhomogenity of matter density resulted in the formation of the first stars in the universe, after 400 million years from the big bang. The radiation and heat from the first stars started ionising the matter, marking the beginning of the Epoch of Reionisation (EoR). The most promising approach so far for observing this early epoch of the universe is the 21cm hyperfine transition line of neutral hydrogen, red-shifted in the observer’s frame to the frequencies in the range of 40-200 MHz. Recently an absorption profile centred at 78 MHz in the sky averaged spectrum was reported by researchers at ASU (Arizona State University, in their Experiment to Detect the Global Epoch of Reionisation Signature, called EDGES Bowman, et al, 2018 ). This absorption spectrum with a width of ~18 MHz and an amplitude of 530 mK is buried in foreground of the order of 2000K. We are exploring Wavelet transform technique to find the domain which gives the best signal to noise ratio of this absorption profile to its foreground. Initially, we studied the behavior of the signal to the Fourier transform (FT) technique. Then we explored the Continuous Wavelet Transform (CWT) to find a suitable technique for feature extraction in spectral analysis and signal reconstruction. The Morlet, Mexican Hat (Ricker) and Symlet wavelets were explored in CWT. We also explored discrete wavelets like Daubechies and Haar. It was found that the absorption profile is very broad, making the available baseline in the foreground estimation too small for applying linear transformation techniques. In view of this, we explored the linear transformation technique by compressing the absorption profile in the frequency domain. We studied the efficacies of using FT, CWT and DWT techniques for extracting the Global 21cm feature in the presence of the foreground signal of sky averaged radio spectrum. We also implemented DWT and reconstructed the signal on GPU (GeForce GTX 1050) using PyCuda libraries.

Keywords or phrases: Absorption profile, Wavelet, Continuous Wavelet Transform, Discrete Wavelet Transform, Graphics Processing unit

Abbreviations

EoREpoch of Reionisation
CMBRCosmic Microwave Background Radiation
FFTFast Fourier Transform
CWTContinuous Wavelet Transform
DWTDiscrete Wavelet Transform
EDGESExperiment to Detect the Global EoR Signature
GPGPUGeneral-purpose Graphics Processing Unit
CUDACompute Unified Device Architecture
SNRSignal to Noise Ratio
RFIRadio Frequency Interference

INTRODUCTION

Background/Rationale

The Universe is believed to have begun in an extremely hot and dense state, followed by a sudden expansion after about 10-32sec, and subsequent slow cooling and expansion. This aided for a gradual drop in temperature and density of the Universe, facilitating for the combination of neutrons and protons to form the atomic nuclei. After a few hundred thousand years, when the Universe was further cooled, these atomic nuclei captured the electrons to form singly ionized Helium, neutral Helium and neutral Hydrogen (​Ganjam, 2018​). This phase of the universe, widely known as Epoch of Recombination contributed to the Cosmic Microwave Background Radiation. When the first stars were formed, the heat and radiation ionized the neutral matter leading to the formation of bubbles of plasma surrounding the source (MIT Haystack Observatory: EDGES). As the bubbles grew bigger, they started to overlap, leading to the formation of much larger bubbles, accelerating the ionisation process. This era of Universe is termed as Epoch of Reionisation. The most powerful probe to study the EoR signal so far is the redshifted Global 21-cm spectral line resulting from the spin-flip hyperfine transition of the hydrogen atom, red-shifted in the observer’s frame to the frequencies in the range of 40-200 MHz (​Bowman, et al, 2018​).

bowman.png
    a, Measured spectrum for the reference dataset after filtering for data quality and radio-frequency interference.The spectrum is dominated by Galactic synchrotron emission.b, c, Residuals after fitting and removing only the foreground model (b) or the foreground and 21-cm models (c). d, Recovered model profile of the 21-cm absorption, with a signal-to-noise ratio of 37, amplitude of 0.53 K, centre frequency of 78.1 MHz and width of 18.7 MHz. e, Sum of the 21-cm model (d) and its residuals (c)​Bowman, et al, 2018​.

    The 21-cm signal from EoR can be studied as fluctuations of temperature in the global signal varying as a function of frequency and can be seen as the absorption profile in the sky-averaged radio spectrum. Recently, an absorption profile centered at 78 MHz in the sky-averaged spectrum was reported by researchers at ASU (Arizona State University), in their Experiment to Detect the Global Epoch of Reionisation Signature, called EDGES (​Bowman, et al, 2018​). This absorption spectrum has a width of ~18 MHz and an amplitude of 530 mK and is buried in the foreground of the order of 2000K. The reported profile was consistent with the expectations of the 21-centimeter signal induced by early stars, except that the amplitude of the profile was twice the predicted amplitude. This large variation between the observed and predicted profiles suggests that either the primordial gas was much colder than expected or the background radiation temperature was hotter than expected (​Bowman, et al, 2018​).

    In this summer project, an attempt is made to analyze the spectral data using Wavelet transform technique to find a contrast between the temperature dip, and the foreground. As a second phase of the project, the wavelet transform technique was implemented on the GPU using PyCuda to explore and appreciate the properties and advantages of GPU.

    Problem Statement

    The EoR signal, seen as the absorption profile is of the order of mK, while the foreground of Global 21cm sky-averaged radio spectrum is of the order of 2000k. This fluctuation of mK is hidden in the foreground, making it challenging to detect. This signal is predicted to be extremely weak, spread over a wide frequency range from 40-200MHz and is buried in the foreground having cumulative effect of Galactic and Extragalactic foregrounds as well as Radio Frequency Interference ​(​Singh, et al, 2017​)​. In order to extract this weak signal from the enormous foreground, it is paramount to find the contrast between the temperature dip and the foreground itself. In this project, an attempt is made to find a domain which gives the highest signal to noise ratio. Since wavelet transform technique involves localisation in the time-frequency domain, this technique is explored to analyse the possible ways of extracting the Bowman dip from the foreground data.

    data_2.png
    Foreground
      Signal_1.png
      Signal
        Comparision of magnitude of Foreground and EoR signal

        For this project, the foreground, as measured by the instrument, data ranging from 63MHz to 103MHz was considered.

        Objectives of the Research

        • To study the Fourier and Wavelet transform of the EoR signal and the foreground and to verify if they can be used as a tool for extracting the signal from the foreground
        • To attempt to find the domain of the highest signal-to-noise ratio for foreground and EoR signals
        • To understand and analyse both continuous and discrete wavelet transform technique for the case of foreground data and the temperature dip.
        • In addition, to understand and appreciate the scientific computational capabilities of the GPU and study it for the case of spectral data.

        In addition to the above objectives, a study was also carried out to understand and appreciate the scientific computational capabilities of the GPU and study it for the case of spectral data analysis.

        LITERATURE REVIEW

        Mathematically, wavelet transform of a signal s(t) with respect to a mother wavelet g(t) is,

        S(τ,a)=12πaS(ω)G(aω)ejωτdω S(\tau,a)\;=\frac1{2\mathrm\pi}\sqrt a\;\int S(\omega)G^\ast(a\omega)e^{j\omega\tau}\;d\omega

        where S(ω) is a Fourier transform of the signal, a>0 is a dilation parameter, τ is a dilation parameter and G*(ω) is a complex conjugate of the Fourier transform of g(t) (​Suvichakorn, et al, 2008​).

        In the discrete wavelet transform, the scale parameter is always discretized to integer powers of 2, ​2j​ , j=1,2,3,..., so that the number of voices per octave is always 1. The DWT provides a sparse representation for many natural signals. In other words, the important features of many natural signals are captured by a subset of DWT coefficients that is typically much smaller than the original signal. This compresses the signal. ​ Thus, DWT always gives the same number of coefficients as the original signal, but many of the coefficients may be close to zero in value. As a result, it is often practical to throw away those coefficients and still maintain a high-quality signal approximation ​Vijendra, 2016​. This forms the basis of signal compression.

        In feature extraction of Heart rhythms, wavelet transform was used for noise removal and peak detection. The major source of noise, in this case, were muscle contraction and movement of the patient while recording ECG. The Heart Rate Variability, HRV, signals were obtained after pre-processing, which were then anaysed in time and frequency domains. Wavelets have been utilized and performed well for signal processing of ECG signals and then the HRV signal is generated and analyzed on the basis of features extracted (​Gautam, et al, 2017​).

        Wavelet transform technique has also been used widely in feature extraction in signal processing. Sea surface elevation and ocean wave spectrum, both significant in marine science research and ocean engineering, have been calculated using the wavelet transform algorithm. The wave scalogram obtained from the data collected by ultrasonic wave gauges on the Cigu pile station from the Taiwan Strait was used to analyse instant nonlinear features at short-time duration (​Wu, et al, 2010​).​

        DWT has been effectively used in the feature extraction of ultrasonic flaw data. The signal was decomposed into detail and approximation coefficients using successive high pass and low pass filters respectively. The approximation coefficients obtained from level-1 of DWT operation was subjected to level-2 DWT operation giving sets of approximate and detail coefficients of level-1 coefficients. The authors considered four levels of decomposition for various feature extraction such as Mean, Variance, Mean of Energy, Maximum and Minimum amplitude among others (​Shodhganga​).

        DWT is a non-redundant transform. It was developed so there would be a one to one correspondence between the information in the signal domain and the transform domain. This tight correspondence in the DWT allows signal reconstruction (​Jorgensen, et al, 2009​).

        As a second phase of the project, discrete wavelet transform was implemented on a GPU. A graphics processing unit (GPU) is similar to a computer's CPU but is designed specifically for performing the complex mathematical and geometric calculations that are necessary for graphics rendering. However, it is interesting to note that GPUs are widely used in Scientific Computations. Some of the fastest GPUs have more transistors than the average CPU. GPUs are highly parallel and architecture sensitive (​Andreas Kl ̈ockner, 2009​) and are capable of handling as many threads as 10 times the number of cores (​Vincent Hindriksen, 2017​). For a computational activity in GPU, Host to Device and Device to Host memory transfer is essential. For this, GPU receives strings of data and instructions through PCI-Express from CPU.

        METHODOLOGY

        Fourier Transforms

        As an initial exercise, the foreground and the signal were analysed in Fourier domain. A real-valued function or the analytic signal gives a Hermitian symmetric signal or complex-valued function in Fourier domain. Since the foreground and the signal are both analytic signals, a complex-valued function is obtained in the Fourier domain. The amplitude is thus calculated as

        Amplitude=(Realpart)2+(Imaginarypart)2 Amplitude=\;\sqrt{(Real\;part)^2\;+\;(Imaginary\;part)^2}

        and is plotted as a function of frequency.

        fft_fore_1.png
        FFT of Foreground
          fft_dip_1.png
          FFT of Dip
            Foreground and Signal in fourier Domain

            It can be clearly seen that the FFT of foreground and Dip vary on a wide dynamic range. Therefore, the FFT results were plotted on a log scale. Since FFT of both foreground and the dip is Hermitian symmetric, the results were shifted to the center of the band.

            fft_fore_shift_1.png
            Shifted FFT of foreground on log scale
              fft_dip_shift_1.png
              Shifted FFT of Signal on log scale
                Shifted FFT of foreground and Signal on log scale

                A similar approach was followed for the sum of foreground and signal. FFT was computed, shifted to the center and was plotted on a log scale.

                fft_sum_1.png
                FFT of Sum of foreground and signal
                  fft_sum_log_2.png
                  Shifted FFT of sum of foreground and signal on log scale
                    FFT of sum of foreground and signal on linear and log scale

                    The purpose is to find the domain with the best signal to noise ratio. For example, FFT of sine wave is a delta peak in frequency domain. Now, consider a noisy sine wave. Extracting the information in the presence of noise in the time domain is difficult. Taking the Fourier transform of the noisy sine gives two peaks with a small Gaussian noise. The sine wave is distinctively seen in Fourier domain.

                    pure_sin.png
                    Pure Sine wave
                      fft_sin_1.png
                      FFT of Pure Sine wave
                        Sine wave and its Fourier transform

                        Gaussian white noise (random noise at all frequencies) with mean=0 and standard deviation=1 was added to the above sine wave.

                        noisy_sin.png
                        Noisy Sine wave
                          fft_noisy_sin_2.png
                          FFT of Noisy Sine wave
                            Noisy sine wave and its Fourier transform.

                            In the case of a Fourier transform of a noisy sine wave, the basis functions are sine and cosine functions, orthogonal to the sine wave. Therefore when FFT was computed, characteristic delta functions were obtained. Noisy sine in time domain is indistinguishable, due to low signal to noise ratio. Transforming it to frequency domain gives characteristic delta peaks with amplitude of ~128 (for 256 point FFT) and noise of amplitude 1, making SNR very high. From Parsival’s theorem, power in time and frequency domain is identical. Thus a domain which gives highest signal to noise ratio is preferred.

                            Continuous Wavelet Transforms

                            In order to study the foreground and the temperature dip using wavelet transform technique, a suitable mother wavelet is to be selected. Results of two important wavelets, the Morlet and the Mexican hat (or Ricker) wavelet were compared. Morlet Wavelet has been used in signals containing Lorentzian function, described by full width half maximum and Voigt function (obtained by convolution of two broadening mechanisms, one of which would produce a gaussian profile) (​Suvichakorn, et al, 2008​).

                            Comparision of Morlet and Ricker Wavelets

                            The wavelet transform was computed using SciPy library in Python. Wavelet transform of the temperature data, Bowman dip, Sum of temperature data and dip, difference of CWT coefficients of sum of dip and foreground & the foreground was computed. Each of these results was compared for both the wavelets.

                            cwt_temp.png
                            CWT of Foreground using Morlet
                              cwt_data.png
                              CWT of Foreground using Ricker
                                cwt_dip_1.png
                                CWT of Signal using Morlet
                                  cwt_dip.png
                                  CWT of Signal using Ricker
                                    cwt_sum.png
                                    CWT of sum using Morlet
                                      cwt_sum_1.png
                                      CWT of sum using Ricker
                                        diff_coeff.png
                                        Difference in coefficients foreground+signal & foreground for the case of Morlet wavelet
                                          diff_coeff_1.png
                                          Difference in coefficients foreground+signal & foreground for the case of Ricker wavelet
                                            int_trans_s1.png
                                            Scale 1 Intensity curve of Signal using Morlet wavelet
                                              int_trans_s1_1.png
                                              Scale 1 Intensity curve of Signal using Ricker wavelet
                                                int_trans_log_s1.png
                                                Scale 1 Intensity curve of Signal on log scale using Morlet wavelet
                                                  int_trans_s1_log.png
                                                  Scale 1 Intensity curve of Signal on log scale using Ricker wavelet
                                                    Comparision of Morlet and Ricker wavelets

                                                    In the scalogram, Fig 8a-8h​​, values on X-axis correspond to the 592 data points (ranging from 63-103MHz), while that on Y-axis corresponds to scales ranging from 1 to 4. For the CWT of foreground using Morlet and Ricker wavelets, there is visible intensity on the scalogram at higher and lower frequencies, while there are no significant values in the mid- frequency range.

                                                    From the above comparison table, it can be observed that, the wavelet transform coefficient for Ricker wavelet in the translation range of 0 to 500 ranges between -0.00527 and 1.51e-09 with the minimum of 1.51e-09 at (591,2) and maximum of 0.00234 at (464, 1.16) , while for the Morlet wavelet, it ranges between -0.000126 and -4e-06 with the minimum of -2.45e-07 at (206, 2.16) and maximum of 0.00162 ,at (457, 1.6). It can be seen from the corresponding intensity plots of scale 1 that the Ricker wavelet gives better results for the given data. Thus, for the rest of the analysis, Ricker wavelet was considered.

                                                    Close observation of Fig 2a & 2b​​ and ​Fig 8l​​ shows that

                                                    • the temperature dip and the foreground are in the same frequency range (dip is not sharp)
                                                    • the foreground is not smooth and has few RFIs
                                                    • there is an edge effect on either side of the CWT coefficients.

                                                    In order to overcome the above problems, either the foreground data set is to be extrapolated over a wide range of frequency or the dip is to be compressed. In the extrapolation method, if the polynomial is of a lower order, there would be a larger error between actual foreground and the polynomial. However, if the polynomial is of the higher order, it would have higher degrees of freedom and the polynomial tends to deviate from the foreground. Thus, the idea of extrapolating the foreground was dropped. The latter method did not appreciably reduce the edge effect. However, higher order (higher scale) CWT gave lesser edge effects (​Fig 10​)​.

                                                    It is well understood that, if a particular frequency data in a signal is to be found, the wavelet should look like the signal we are looking for, and that the basis of the transform must be orthogonal. For everything else but the identical signal, the output will be zero (​Lathi, 2009​). This was verified by adding a Morlet wavelet to an exponential function and finding morlet wavelet transform of the data over a large scale ranging from 1 to 256. A clear intensity peak was observed in the scalogram of the CWT coefficients.

                                                    double_mor_exp.png
                                                    Morlet wavelet added to an exponential signal
                                                      cwt_double_mor_noisyexp.png
                                                      Scalogram of the CWT coefficients
                                                        Feature extraction using Morlet wavelet

                                                        Wavelet transform technique is based on the match filtering concept. The idea of match filtering worked well when Morlet function was to be found from an exponential signal, where the feature to be extracted was of the order of the signal itself. The scale variable provides analysis in the frequency domain: a compressed version of the wavelet function corresponds to the high-frequency components of the original signal while the stretched version corresponds to the low-frequency components; and the translation variable provides analysis in the time domain (​R. E. J. Yohanes, et al, 2012​).

                                                        The scalogram gives a 3-dimensional intensity plot of Wavelet transform coefficients. To study these coefficients, intensity vs translation of CWT of Foreground and the Signal was plotted corresponding to each scale ranging from 1 to 6.

                                                        scaleall.png
                                                          CWT of Foreground and the Signal using Ricker wavelet over the scale 1 to 6 (as described in step 2)

                                                          In order to find the contrast between the signal and the foreground, the following procedure was followed:

                                                          • Compute CWT of the Signal and the foreground
                                                          • Compute and plot the coefficients on a logarithmic scale (Consider absolute values, as logarithm can be computed only for positive values)
                                                          • Compute FFT of data obtained in 2
                                                          • Plot the results on a log scale to cover the Dynamic range
                                                          fft_scallall.png
                                                            FFT of CWT coefficients (as described in step 3)
                                                            cwt_log_fft_log_exp.png
                                                              FFT plotted on a log scale, maximised near the peak (as describes in step 4)

                                                              It can be clearly seen that the dip shoots up at about 81.8 MHz. Although this method gives an excellent contrast between the signal and the foreground, in real-time data, the Dip is uncertain and is not known beforehand. Thus, this cannot be a generalised solution.

                                                              Discrete Wavelet Transforms

                                                              The important feature to be noted in DWT coefficients is that it always gives the same number of coefficients as the original signal, unlike CWT, which gives M-by-N matrix of N samples for an N-length signal and M number of scale. Also, in DWT, many of the coefficients may be close to zero in value, making it convenient to throw away those coefficients and still maintain a high-quality signal approximation. The CWT is a highly redundant transform. There is significant overlap between wavelets at each scale and between scales. It also requires much larger computational resources than that in DWT Wavelet Toolbox, MathWorks.​

                                                              Reconstruction of Signal from the foreground data

                                                              In order to complete the feature extraction, signal reconstruction is necessary. Inverse CWT can be implemented, but it is often the case that the reconstruction is not perfect. Because the CWT is a redundant transform, there is no unique way to define its inverse Wavelet Toolbox, MathWorks. The inverse CWT requires an infinite number of coefficients. Due to the fact that CWT takes in non-base 2 values also, there can be many values in an octave, unlike in DWT, which has only one value in an octave. This makes CWT an overlapping transform with 2D coefficient matrix while DWT gives one-dimensional coefficients corresponding to each data point.

                                                              The DWT provides a perfect reconstruction of the signal upon inversion. This implies that DWT coefficients of a signal can be used to synthesize an exact reproduction of the signal within numerical precision. The outputs of the DWT algorithm are detailed coefficients, corresponding to high-frequency detail signals, and approximation coefficients, corresponding to the coarse approximation of the original signal in the time domain (​Lee, et al, 1994​). They represent the degree of correlation between the analyzed signal and the wavelet function at different instances of time, making it suitable for feature extraction in spectral analysis.

                                                              The procedure followed for reconstruction of Signal using DWT is:

                                                              • Compute DWT of Foreground
                                                              • Compute DWT of foreground+dip
                                                              • Find the difference between DWT coefficient of (2) and (1), both detailed and approximate
                                                              • Compute IDWT of Difference
                                                              dwt_fore.png
                                                              DWT of Foreground using Daubechies2 wavelet
                                                                dwt_sum_1.png
                                                                DWT of Foreground+Signalusing Daubechies2 wavelet
                                                                  difference.png
                                                                  Difference in Coefficients from (b) and (a)
                                                                    recons.png
                                                                    Reconstructed signal
                                                                      error.png
                                                                      Error in the reconstructed signal
                                                                        db2.png
                                                                        Daubechies 2 wavelet 
                                                                          Signal reconstruction using db2 wavelet in DWT

                                                                          The above results are obtained from asymmetric, orthogonal Daubechies 2 wavelet. Considering both approximate and detailed coefficients, the dip can be satisfactorily reconstructed. The error in the reconstructed signal is observed to be in the order of 10-13, which is an acceptable margin of error in signal reconstruction. Discrete wavelets like Haar, Daubechies 2 and Symlet 2 wavelets were also used in the signal reconstruction. All of them gave similar results.

                                                                          The DWT analysis windows are fixed in both the time and scale directions, so if the coefficients are plotted, a grid of boxes that start out large at one end of the scale axis and end as smaller ones will be obtained (​Jorgensen, et al, 2009​). To verify the details of the DWT, second level decomposition is applied to the obtained approximate coefficients. The results were similar, upholding the fact that the approximate coefficients give similar plot as that of the input signal while detailed coefficients are much smaller.

                                                                          dwt_dip_log_level12.png
                                                                            Level 1 and Level 2 DWT of Signal on a log scale

                                                                            The difference in wavelet transform coefficients of the sum of foreground and signal & the foreground gives the signature of the Dip.

                                                                            cwt_scale6.png
                                                                              Signature of the Signal in the wavelet domain

                                                                              This signature is clearly due to Dip, as the effect of the foreground has been removed. But when there is variation in the foreground, this signature may not remain the same. Thus finding a basis function which makes the signal stand out in the foreground is essential.

                                                                              Computing with GPU

                                                                              In a second phase of the project, DWT was implemented on a GPU available in the lab, Nvidia GeForce GTX 1050. A GPU uses special programming tool to help it analyze and use data. In this project, PyCuda, which gives Pythonic access to Nvidia’s CUDA, parallel computation API was used. The pycuda module was used to fetch the details of GPU being used. It was also a verification if the GPU is active in computation. The various device properties that were verified were:

                                                                              Device Properties for GeForce GTX 1050
                                                                              CUDA Driver Version / Runtime Version9.1 / 9.1 C
                                                                              CUDA Cores128
                                                                              GPU Max Clock rate1569 MHz (1.57 GHz)
                                                                              Memory Bus Width128-bit
                                                                              Total number of registers available per block65536
                                                                              Warp size32
                                                                              Maximum number of threads per multiprocessor2048
                                                                              Maximum number of threads per block1024
                                                                              Max dimension size of a thread block (x,y,z)​(1024, 1024, 64) ​
                                                                              Multiprocessors 5

                                                                              To implement DWT on GPU, PyCuda's PyGASP package was used. It provides with the dwtCuda function for the computation of DWT using GPU. DWT was implemented on both CPU and GPU and the results were compared.

                                                                              cpu_foreground.png
                                                                              Level-1 DWT of  foreground using db2 wavelet computed in CPU
                                                                                gpu_foreground.png
                                                                                Level-1 DWT of  foreground using db2 wavelet computed in GPU
                                                                                  Comparision of CPU and GPU results

                                                                                  The error in the DWT coefficients was calculated. It can be seen that the GPU and CPU give similar results.

                                                                                  error_1.png
                                                                                    Difference in DWT coefficients as calculated in CPU and GPU

                                                                                    When the time duration of the operation was seen, it was noted that

                                                                                    CPU Time= 117 μsec

                                                                                    GPU time= 1901 μsec

                                                                                    Enabling nvprof to see the detailed time scale of computation for the DWT of foreground using db2 wavelet, GPU activity of about 51μsec was observed.

                                                                                    nvstat.png
                                                                                      Profiling of the GPU activity in the computation of DWT

                                                                                      This could be because of the fact that the GPU is configured for the floating point operations, while the data given is a double precision array. However, Signal reconstruction was done in GPU using DWT. The procedure followed was the same as in the case of CPU explained in ​Sec. 3.3.1​.

                                                                                      reconstruct.png
                                                                                        Reconstructed Signal using DWT module in PyGASP

                                                                                        CONCLUSION AND RECOMMENDATIONS

                                                                                        When the signal is a sinusoid buried in the noise, analysis of the Fourier transform technique very clearly indicated that a sinusoid shows up with a higher signal to noise ratio in the frequency domain. However, if the signal is of the nature of the predicted EoR signal in the frequency domain, Fourier transform to the other domain in no way improves the SNR. Hence a simple Fourier analysis is not helpful in detecting EoR signal buried in the foreground. We wanted to investigate and find out whether there are wavelets which are more suitable for separating out EoR signal from the foreground. Thus we started with CWT and found certain disadvantages (as described in ​​Sec. 3.3.1​​). So we moved on to DWT and explored its properties. We found that by applying the procedure (as described in ​Sec. 3.2.1​​), the contrast of the EoR signal could be significantly improved. However, it has not been verified if this method is independent of the foreground and noise. As the foreground changes with the direction of observation, this method cannot be generalized. Therefore, in an overall sense, one can continue with different foregrounds and verify if it works with all foregrounds. It is to be noted that, if a particular frequency data is in the signal and it is to be found, the wavelet should be similar to the signal we are looking for. Therefore, in order to extract the signal, a custom wavelet similar to the signal and of the comparable magnitude forming orthogonal basis should be developed. Alternatively, we can also take up the job of constructing a set of wavelets which would remove all the foreground but remain orthogonal to the EoR signal.

                                                                                        The results from the GPU were as seen in the CPU, except for the time duration of computation. This time duration of computation can be reduced considerably with various optimization techniques and by upgrading to higher end GPU. However, successful signal reconstruction suggests that GPUs can be used as a potential computational tool in spectral analysis.

                                                                                        ACKNOWLEDGEMENTS

                                                                                        First and foremost I express my deep sense of gratitude to the Indian Academy of Sciences (IASC-INSA-NASI) for giving me an opportunity to carry out this project in the Summer Research Fellowship Programme-2018.

                                                                                        I owe my heartiest gratitude to respected Prof. Udaya Shankar N, Raman Research Institute, Bangalore for the continuous support, motivation, and immense knowledge. I could not have imagined having a better advisor and mentor for the project.

                                                                                        My sincere thanks to Prof. R Subrahmanyan for his valuable guidance and encouragement throughout the project.

                                                                                        My sincere thanks to Srivani K S, Girish B S and Somshekhar R for their valuable guidance, consideration and always inspiring and supporting me during the internship.

                                                                                        I thank Shruti Bhatporia, Mugundhan V, Saurabh Singh, Jishnu Nambissan T and Sindhu Gaddam for their continuous guidance, encouragement, moral and academic support throughout my project.

                                                                                        It is my privilege to express regards to Dr. Bhuvaneshwar Babu and Dr.Karthik Shastry for providing the letter of recommendation to the Indian Academy of Sciences.

                                                                                        I also thank my institute R V College of Engineering, Bengaluru for allowing me to pursue this internship.

                                                                                        I am glad I used Authorcafe for writing my report and I would definitely recommend this to my peers, provided some minor bugs are fixed.

                                                                                        I express my deepest sense of gratitude towards my beloved parents, fellow interns and friends for their constant support and encouragement.

                                                                                        APPENDICES

                                                                                        The various libraries used in this project are:

                                                                                        • Numpy : A fundamental package for scientific computing with Python, mainy used for numerical opreations of arrays. (http://www.numpy.org/)
                                                                                        • Matplotlib: A Python 2D plotting library used extensively in this project for visualising the results. (https://matplotlib.org/)
                                                                                        • FFT package in Numpy: Discrete Fourier Transform package of numpy used for computation of FFT.(https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.fft.html)
                                                                                        • SciPy: A scientific computation tool in Python. The package scipy.signal was used extensively for the computation of CWT, using available wavelets. (https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.fft.html)
                                                                                        • PyWavelet: Package used for computation of DWT of data using a wide range of available wavelets. (https://pywavelets.readthedocs.io/en/latest/)
                                                                                        • PyGASP: Package used for computation of DWT and IDWT on GPU. This package contains several transformations used for signal processing, along with associated tools (https://pypi.org/project/PyGASP/).

                                                                                        References

                                                                                        • Bowman, Judd D., et al. "An absorption profile centred at 78 megahertz in the sky-averaged spectrum." Nature 555.7694 (2018): 67.

                                                                                        • Ganjam, Samskruthi. (2018)., "Satellite Mission for the study of Epoch of Reionisation", submitted as B.Tech thesis to Ramaiah Institute of Technology, Bangalore

                                                                                        • https://www.haystack.mit.edu/ast/science/epoch/

                                                                                        • Singh, Saurabh, et al. "SARAS 2: a spectral radiometer for probing cosmic dawn and the epoch of reionization through detection of the global 21 cm signal." arXiv preprint arXiv:1710.01101 (2017).

                                                                                        • Suvichakorn, Aimamorn, and Jean-Pierre Antoine. "Analyzing NMR spectra with the Morlet wavelet." Signal Processing Conference, 2008 16th European. IEEE, 2008.

                                                                                        • Vijendra, V., and Meghana Kulkarni. "ECG signal filtering using DWT haar wavelets coefficient techniques." Emerging Trends in Engineering, Technology and Science (ICETETS), International Conference on. IEEE, 2016.

                                                                                        • Gautam, Desh Deepak, V. K. Giri, and K. G. Upadhyay. "Feature extraction of HRV signal using wavelet transform." Convergence in Technology (I2CT), 2017 2nd International Conference for. IEEE, 2017.

                                                                                        • Wu, L. C., C. C. Kao, and D. J. Doong. "Application of the wavelet transform and inverse wavelet transform to analyze the ocean wave signals from data buoys." OCEANS 2010 IEEE-Sydney. IEEE, 2010.

                                                                                        • shodhganga.inflibnet.ac.in/bitstream/10603/26360/12/12_chapter 7.pdf

                                                                                        • Jorgensen, Palle ET, and Myung-Sin Song. "Comparison of discrete and continuous wavelet transforms." Encyclopedia of Complexity and Systems Science (2009): 1-30.

                                                                                        • Andreas Kl ̈ockner., "GPU Metaprogramming using PyCUDA: Methods & Applications", Division of Applied Mathematics Brown University Nvidia GTC · October 2, 2009

                                                                                        • https://streamhpc.com/blog/2017-01-24/many-threads-can-run-gpu/

                                                                                        • Lathi, B. P. Principles of Signal Processing and Linear Systems. ser. 9780198062288, Oxford International, 2009.

                                                                                        • R. E. J. Yohanes, W. Ser and G. b. Huang, "Discrete Wavelet Transform coefficients for emotion recognition from EEG signals," ​ 2012 Annual International Conference of the IEEE Engineering in Medicine and Biology Society, San Diego, CA, 2012, pp. 2251-2254. doi: 10.1109/EMBC.2012.6346410

                                                                                        • https://in.mathworks.com/help/wavelet/gs/continuous-and-discrete-wavelet-transforms.html

                                                                                        • Lee, Daniel TL, and Akio Yamamoto. "Wavelet analysis: theory and applications." Hewlett Packard journal 45 (1994): 44-44.

                                                                                        Source

                                                                                        • Fig 1: Bowman, Judd D., et al. "An absorption profile centred at 78 megahertz in the sky-averaged spectrum." Nature 555.7694 (2018): 67.
                                                                                        More
                                                                                        Written, reviewed, revised, proofed and published with