Signal Processing techniques to optimize the detection of radio pulsar signals
Abstract
[ANGLÈS] Radio Pulsars are neutron stars that emits high polarized electromagnetic pulses with a very accurate and stable periodicity. Adding the fact that those pulses have a wideband nature, so they can be received almost everywhere, make Radio Pulsar signals a perfect candidate for navigation systems. The challenge, however, is that the Radio Pulsar signal is degraded and submerged in Additive White Gaussian Noise when it is propagated through the ISM, so that make them difficult to detect. The aim of the project is the design of the optimum receiver in order to estimate the time of arrival of the Radio Pulsar Signal and make possible the real-time navigation. In this contribution I introduce the detection theory applied to radio pulsars and show the different signal processing techniques to improve the detection with the minimum computational time possible. Although the previous work has centered in the Signal to Noise Ratio of the pulsar signals, our focus will be the improvement of the Energy to Noise Ratio. We will see as Epoch Folding and the increment of the Bandwidth of the receiver are the best solutions. Moreover, Integration in Time appears to be one of the most promising techniques to simplify the computation of the whole process. Also I demonstrate how the detection performance is not affected by any digital or analogue filter, so Low-Pass Filtering will not have any effect in the receiver but to eliminate spurious signals. Some experiments with simulated wideband signals are shown in order to prove the optimality of the signal processing techniques. Finally, summarizing all the results obtained in the Thesis, I propose two optimum receivers for pulsar-based navigation system applications. It is known that the interstellar medium has a frequency dependent transfer characteristic, so the higher frequencies of a signal arrive earlier than the lower frequencies. This will cause the pulse profile of the radio pulsar signal to appear dispersed in time. Up until now, the researchers have been performing de-dispersion techniques in order to obtain the original power pulsar profile spending more than the 80% of the actual processing time. The prove to avoid de-dispersion without decreasing the detection performance are presented. Furthermore, some simulations with data from the pulsar PSR B0329+54 recorded by the Westerbork Observatory are shown. PSR B0329+54 is one of the strongest pulsar signal visible in the northern hemisphere with one of the lowest Dispersion Measure, so it will be easier to perform the Signal Processing techniques and see the profile of the pulse. We will see how the results are not the ones expected, because the Radio pulsar signal happens to be very Narrowband. Finally, I am going to give a possible explanation of what is happening and in what part of the acquisition the wideband nature of the rotating star pulses are lost.
Full text
Signal and Information Processing (EEMCS, TU Delft) Signal Processing Techniques to optimize the Detection of Radio Pulsar Signals Daniel Hernando Portero Supervisors Specialisation Type of report Date Dr. Ir. Richard Heusdens Dr. Nikolay D. Gaubitch Signal Processing Master Of Science Thesis August 28, 2014
Signal Processing techniques to optimize the detection of Radio Pulsar Signals Master Of Science Thesis For the degree of Master of Science in Telecommunications Engineering at Polytechnic University of Catalonia under the Erasmus exchange programme at Delft University of Technology Faculty of Electrical Engineering, Mathematics and Computer Science (EEMCS) Delft University of Technology Escola Tecnica Superior d’Enginyeria de Telecomunicacions de Barcelona (ETSETB) Polytechnic University of Catalonia Daniel Hernando Portero August 22, 2014
Cover Image: Four antennas of the Atacama Large Millimeter/submillimeter Array (ALMA) gaze up at the star-filled night sky, in anticipation of the work that lies ahead. The Moon lights the scene on the right, while the band of the Milky Way stretches across the upper left. Polytechnic University of Catalonia Escola Tecnica Superior d’Enginyeria de Telecomunicacions de Barcelona Supervisor: Gregori Vazquez Delft University of Technology Intelligent SystemsDepartment Multimedia and Signal Processing (MSP) Group Signal and Information Processing (SIP) Lab Supervisors: Dr. Ir. Richard Heusdens Dr. Nikolay D. Gaubitch Copyright 2014 TU Delft, UPNA All rights reserved.
Abstract Key Words: ENR: Energy to Noise Ratio, GENR: Generalised Energy to Noise Ratio, TOA= Time of Arrival, Additive White Gaussian Noise, Antenna, De-Dispersion. Radio Pulsars are neutron stars that emits high polarized electromagnetic pulses with a very accurate and stable periodicity. Adding the fact that those pulses have a wideband nature, so they can be received almost everywhere, make Radio Pulsar signals a perfect candidate for navigation systems. The challenge, however, is that the Radio Pulsar signal is degraded and submerged in Additive White Gaussian Noise when it is propagated through the ISM, so that make them difficult to detect. The aim of the project is the design of the optimum receptor in order to estimate the time of arrival of the Radio Pulsar Signal and make possible the real-time navigation. In this contribution I introduce the detection theory applied to radio pulsars and show the different signal processing techniques to improve the detection with the minimum computational time possible. Although the previous work has centred in the Signal to Noise Ratio of the pulsar signals, our focus will be the improvement of the Energy to Noise Ratio. We will see as Epoch Folding and the increment of the Bandwidth of the receiver are the best solutions. Moreover, Integration in Time appears to be one of the most promising techniques to simplify the computation of the whole process. Also I demonstrate how the detection performance is not affected by any digital or analogue filter, so Low-Pass Filtering will not have any effect in the receptor but to eliminate spurious signals. Some experiments with simulated wideband signals are shown in order to prove the optimality of the signal processing techniques. Finally, summarizing all the results obtained in the Thesis, I propose two optimum receptors for pulsar-based navigation system applications. It is known that the interstellar medium has a frequency dependent transfer characteristic, so the higher frequencies of a signal arrive earlier than the lower frequencies. This will cause the pulse profile of the radio pulsar signal to appear dispersed in time. Up until now, the researchers have been performing de-dispersion techniques in order to obtain the original power pulsar profile spending more than the 80% of the actual processing time. The prove to avoid de-dispersion without decreasing the detection performance are presented. Furthermore, some simulations with data from the pulsar PSR B0329+54 recorded by the Westerbork Observatory are shown. PSR B0329+54 is one of the strongest pulsar signal visible in the northern hemisphere with one of the lowest Dispersion Measure, so it will be easier to perform the Signal Processing techniques and see the profile of the pulse. We will see how the results are not the ones expected, because the Radio pulsar signal happens to be very Narrowband. Finally, I am going to give a possible explanation of what is happening and in what part of the acquisition the wideband nature of the rotating star pulses are lost.
Acknowledgements Writing the Acknowledgement of my Master Thesis means that I am finishing my studies. I still do not realize the importance of that fact, that in a brief period of time I will have to face the first days of the rest of my life. I should feel scared, or doubtful with the things that await me in the future, but truth be told I am eager to start and create my own path. First of all I would like to thank Dr. Richard Heusdens for giving me the opportunity of working in that project with him. I have enjoyed all the moments of my research, even the hard ones. I have learned how to work alone and at the same time, being part of a group. I also want to thank him for let me being part of the weekly discussions about the project, where at the beginning I was a little bit lost but with the help of him and Dr. Nikolay Gaubitch I became one of them. I also want to express my gratitude to Dr. Nikolay Gaubitch for the good advices given during all my Thesis. Several people have contributed to make this semester an unforgettable one. I would like to start with the people from the Library Group. Although the first month we did not know each other, at the end we became a family. Of course I am talking about Roger, Paula, Benedetta, Cristina, Mikel, Pati and Riccardo. I hope we see each other soon. The second group of people are the Marcushof/Roland crew. It will be difficult to forget all the dinners, beers, conversations we have shared. At the end I do not want to forget the rest of the International Group that has contributed to the most amazing period of my life. In special, many thanks to Roger, my library mate and best friend in that Erasmus that has helped me to being constant in that long-distance race that is the Thesis. I will not end this Acknowledgements without saying thanks to my group of friends of Barcelona. Although we have not seen each other for 6 month, our friendship is so strong that the distance has not had any effect on us. I am referring to Xavi, Cesc, Julia, Carles, Gonzalo, Pablo, David, Pedro and Anna. Also I want to give special thanks to David, Pedro and Xavi for giving me advices whenever I needed and to Anna for always being there no matter what. Finally, I need to thank my parents, my sister, and the rest of my family, the most important people of my life. It is only thanks to them that I have managed to become the person I am today, for which I will be eternally grateful. Delft, The Netherlands Delft University of technology Daniel Hernando Portero August 28, 2013
Table Of Contents 1 Introduction 1 1.1 Introduction to the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Motivations..................................... 1 1.3 Thesisgoals..................................... 1 1.4 Thesiscontributions ................................ 2 1.5 Outline ....................................... 2 2 Radio Pulsar Signals 4 2.1 Pulsar description and emission properties . . . . . . . . . . . . . . . . . . . . 4 2.1.1 Pulsars ................................... 4 2.1.2 Radio Pulsar signal characteristics . . . . . . . . . . . . . . . . . . . . 6 2.2 Propagationeffects................................. 8 2.3 Pulsar-Based Navigation System . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3.1 Navigation system challenge . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4 PulsarsConclusions ................................ 10 3 Detection Theory applied to Radio Pulsar Signals 11 3.1 Basic Detection for deterministic signals . . . . . . . . . . . . . . . . . . . . . 11 3.2 Generalised Likelihood Ratio Test applied to Radio Pulsar Signals . . . . . . 14 3.3 Conclusions of the Detection Theory . . . . . . . . . . . . . . . . . . . . . . . 20 4 Radio Pulsar Signal Model 22 4.1 Filtered Analogue Signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4.2 A/Dconverter ................................... 24 4.3 Generalised Energy To Noise Ratio . . . . . . . . . . . . . . . . . . . . . . . . 26 5 Theoretical Signal Processing Techniques 28 5.1 EpochFolding ................................... 28 5.1.1 IntegrationinTime ............................ 30 5.2 Low-PassFiltering ................................. 32 5.3 Downsampling ................................... 36 5.4 Oversampling.................................... 41 5.5 Undersampling ................................... 48 5.6 Signal Processing Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 6 Radio Pulsar Signal PSR B0329+54 observations 53 6.1 PSRB0329+54features .............................. 53 6.2 Radio Pulsar Signal PSR B0329+54 data acquisition from WSRT . . . . . . . 54 6.2.1 Process to visualize the signal . . . . . . . . . . . . . . . . . . . . . . . 55 6.2.2 Representation of the radio Pulsar Signal . . . . . . . . . . . . . . . . 57
7 Signal Processing experiments 60 7.1 Simulateddata ................................... 60 7.1.1 EpochFolding ............................... 62 7.1.2 Low-PassFiltering............................. 64 7.1.3 WhiteningProcess............................. 66 7.1.4 Downsampling ............................... 68 7.1.5 Increase the Bandwidth of the signal . . . . . . . . . . . . . . . . . . . 69 7.2 Radio Pulsar signal B0329+54 from WSRT . . . . . . . . . . . . . . . . . . . 71 7.2.1 Epoch Folding + High Pass Filter . . . . . . . . . . . . . . . . . . . . 71 7.2.2 Downsampling + High-Pass Filter . . . . . . . . . . . . . . . . . . . . 75 7.2.3 IntegrationinTime ............................ 76 7.3 Experimental Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 8 Proposed Receptor 80 9 Summary and Future work 83 9.1 Summary ...................................... 83 9.2 FutureWork .................................... 85 Bibliography 86
Chapter 2 Radio Pulsar Signals 2.1 Pulsar description and emission properties 2.1.1 Pulsars Pulsars are highly magnetized, periodically rotating neutron stars that emit a beam of electromagnetic radiation. They were first discovered by Jocelyn Bell on November 28th, 1967. At first and for a short period of time scientists thought the electromagnetic emission was coming from an extra-terrestrial civilization, but this theory was soon rejected. Up to now over 1500 pulsars have been detected in our Galaxy [3], but it is expected that thousands more will be discovered during the next few years. In spite of more than four decades of intensive research there are still many open questions in pulsar astronomy, and thus it would be a fair statement that these neutron stars are understood only poorly. On the one hand, studies so far have allowed us to characterize the properties of the emitted signals after their travel through the Interstellar Medium (ISM), but on the other hand the complete description of the internal structure of a pulsar remains as a complex issue. At the best of our knowledge, the answer to questions such as how many pulsars are there in the Galaxy, what is their birth rate, how are isolated millisecond pulsars produced, how many pulsar planetary systems exist or many others are either unknown or simply there is a lack of general scientific agreement about them. These neutron stars have magnetic fields of the order of 108to 1015 G (Earths magnetic field magnitude at its surface ranges from 0.25 to 0.65 G), and as a result of Maxwells equations an electric field is induced. Charged particles are accelerated to the magnetic poles of the pulsar by this electric field, and as they are travelling through a magnetic field a beam of electromagnetic radiation of high magnitude is emitted alongside the magnetic axis. 4
Figure 1: Rotating pulsar model and its emission. Credit B. Saxton/NRAO/AUI. A representation of this phenomenon can be shown in Figure 1. It is clearly seen that the magnetic axis and the rotational axis are not necessarily the same. This misalignment causes the intensity of the electromagnetic radiation to vary in a periodic fashion when received from a fixed line of sight. Indeed, the beam is only seen from Earth as it sweeps past our line of sight once for every rotation of the neutron star, which leads to the pulsed nature of its appearance. In addition, both the angle between spin axis and magnetic axis and the frequency of the rotation is unique for each pulsar. So, this information becomes an exclusive signature. [17] Two parts of the spectrum of pulsar emissions are good candidates for navigation, mainly the radio spectrum and the high-energy spectrum, such as Xray and γ-ray. The choice of what kind of spectrum use can be made primarily using three criteria: quality of the received signal, equipment and pulsar availability. While high-energy photons, by definition, give a better SNR, emittance in the radio spectrum generally require much less from the receiver. Furthermore, pulsars emitting strong and usable signals in the radio spectrum are significantly more plentiful. Thus, while the SNR suffers in radio spectrum, the navigation system can be realized with less of a burden due to receiver equipment. [19] Furthermore, a navigation system based on radio pulsars could find use in vehicle navigation on Earth as well as rover tracking on other planets where high-energy pulsar signals are blocked by the atmosphere. Hence, the research of this thesis is targeted towards devising a navigation system based on radio pulsars. Such a system has useful possibilities that are not offered by the high-energy pulsar based system. 5
2.1.2 Radio Pulsar signal characteristics In this section it is going to discuss the main characteristics that make pulsars unique compared to other stellar entities. The first and more important one is the periodicity. Signals coming from pulsars are highly periodic. Each particular pulsar has its own particular periodicity which is different from the other ones. The most rapidly rotating neutron star currently known is PSR B1937+21 with a period of only 1.56 ms. In contrast, the longest period observed for any radio pulsar so far is 8.5 s for PSR J2144-3933. The pulse period of all pulsars slowly decays, supposedly until they come to a full stop. This decay is very slow, and quite stable, and has already been determined for most pulsars. It ranges from 10−13 s/s up to 10−19 s/s, which makes it a very slow decay. This gives the pulsars their characteristic frequency stability, but that decay can be used for navigational purposes as well, due to relativistic effects, as the pulse decay will appear faster or slower compared to the expected decay at earth. Nevertheless, the remarkable fact about the pulsar periodicity is that it is extremely precise. In some cases (millisecond pulsars), the regularity of the pulsation is as precise as an atomic clock. This stability allows millisecond pulsars to be used in establishing ephemeris time or building pulsar clocks. Due to this fact, pulsars are ideal for time-of-arrival (TOA) based navigation systems, as it will be explained in further sections. [18] The second main characteristic of the radio pulsar signals is their spectrum. Pulsar emissions are known to occupy a very wide band of the electromagnetic spectrum. However, based on the location of the frequency range of the emission it is common to classify them into X-ray pulsars (3 ∗1016 to 3 ∗1019 Hz) or radio pulsars (3 kHz to 3 THz), although some pulsars have been found to emit in visible light, gamma rays or even all of the frequency bands aforementioned, making a possible total emission spectrum from 3 kHz to 3 ∗1020 Hz. One interesting result of this issue is that no matter which frequency a receiver is tuned at it will still be able to receive the signal. That Wideband nature of the pulsar signals will be the first and most important assumption used in this Thesis to design the optimal receptor. Another interesting fact is that if you cut some part of the spectra (Pe: Filtering the signal), the integrated power profile of the pulsar will not change. However, the amplitude of that profile will be smaller. That will be an important fact to take into account in the Chapter 5. Note that for a particular pulsar the pulse shape varies as a function of the observing frequency, as stated in Figure 2. [17] 6
Figure 2. Multi-frequency pulse profile of two pulsars: (a) B1133+16 (b) J2145-0750. Credit Lorimer and Kramer, EPN database [5]. Pulsars emit the strongest signals at the lowest frequencies. At increasing frequencies, the signal levels will exhibit a decay, which differs from pulsar to pulsar. However, at lower frequencies the background noise temperature on earth is quite high, and even in space, interference caused by the planets and the sun are relatively strong. Selecting a high observing frequency is therefore beneficial from this point of view. Pulsars are one of the most polarized radio sources. They usually have linear polarizations, but in some cases they can be received with circular or elliptical polarizations. The Stokes Parameters describe the polarization of the signal, but for the pulsars application we are going to use the I parameter. I=E2 0=|Ex|2+|Ey|2 It is now clear that the I Stoke parameter is related to the total intensity or power of the electromagnetic radiation, i.e. is the actual pulse profile. Astronomical observations almost always record the whole four Stokes parameters so complete information about the state of polarization of the signal is achieved. There are a relation between the I Stoke Parameter and the total power of the electromagnetic radiation. That can be seen in the next equation where the total radiated power by the star is showed 7
PT=Z ZS |Eθ|2+|Eφ|2 ηdS(W) (2.1) Where θand φrefer to the spherical coordinates, ηis the characteristic impedance of the medium and S is a spherical surface emulating the the radio telescope antenna. To conclude, in order to obtain the shape of the pulse profile of a particular pulsar the radio telescope will acquire both Vx(t) and Vy(t) voltage signals. Then, they will be summed together following the expression |Vx|2+|Vy|2. [17] 2.2 Propagation effects A pulsar signal travels very large distances on its way to reaching our planet. Pulsars are located at several hundred or in other cases, several thousand light years away from Earth. The signals pass through the intergalactic space, which is scientifically known as the Interstellar Medium (ISM) and are affected by different effects: Dispersion, Scintillation, and Scattering. These effects are discussed and analysed in the following text. Furthermore, a brief review of the actual de-dispersion technique is presented in order to know its properties and analyse its influence in the detection. Scintillation is a process where inhomogeneities of the refractive index of the medium (caused by strong variations of electron densities) produce phase modulations on the propagating pulsar signal. That leads to a fluctuation of the intensity on a variety of bandwidths and time-scales. This effect is modelled as a thin screen of irregularities midway between the Earth and the pulsar [18]. It has been demonstrated that this effect is highly frequency dependent. Interference can occur only if the phases of the waves do not differ by more than about 1 radian. Then, as the phases are frequency dependent, there is a limitation in bandwidth of the interfering waves. This means that waves outside the scintillation bandwidth ∆f∞f4will not contribute [3]. The most powerful method to deal with the Scintillation is the average of different received periods. As state before, although the different pulses can arrive with a very different intensity, the integrated power profile happens to be very stable. Therefore, after some folds we can assume that the scintillation effect provoked by the ISM is gone. Scattering is basically a radiation effect related to multipath environments. In the thinscreen model introduced before this effect can be related directly to the variable path lengths. From the received point of view, the pulse shape will be broadened since not only the directpath component reaches it, but also several delayed versions of it that travel through different paths. This will cause the appearance of the characteristic exponential tails, with the consecutive reduction in the SNR. Note that this effect is also frequency dependent, with a much lesser impact when observing high frequencies. [17]. For a frequencies higher than 600 MHz the scattering disappear. Therefore, in order to avoid this effect we will avoid low-frequencies in the application of detect the radio pulsar signal. The interstellar medium has a frequency dependent transfer characteristic: higher frequency signals arrive earlier than lower frequency signals, even though the time of transmission was the same. This will cause the pulse profile in a broadband receiver system to appear smeared out in time, and will change the pulsar signal shape. The phenomenon is called Dispersion. The Dispersion depends on one term, and this is the Dispersion Measure (DM) 8
[17]. The dispersion measure is constant only for a certain measurement time and position. In other words, as the interstellar medium is not homogeneous, the dispersion measure will change depending on where and when the observer has taken the measurements. Therefore, every Pulsar has a different Dispersion Measure. Dispersion can be removed by the process of de-dispersion. There are two known methods to de-disperse the received signal, the first method de-disperses the signal in the time domain and is called incoherent de-dispersion, the second method employs frequency domain operations and is called coherent de-dispersion. The computational requirements of de-dispersing unit are very high [13]. Actually, all the acquisitions from Radio Pulsar Signals are de-dispersed in order to be processed with a better SNR. That techniques require more than the 80% of the computational time needed to receive, process and detect the signal. Later on and unlike a lot of researches think, we will see as the de-dispersion process can be avoided. That is due the fact that dispersion doesn’t change the energy of the Radio Pulsar Signal but only the shape. In fact, the effect of the Interstellar Medium (ISM) is described as a phase only filter, as by the Fourier Transform delay in time domain is equivalent to phase shift in the frequency domain. 2.3 Pulsar-Based Navigation System Pulsar based navigation is not a novel area of research, and might once offer the possibility to man of safely travelling distances much beyond Earth. Up until now the focus has been on X-ray based pulsar navigation, whereas recent studies focus on the possibility of using radio pulsars. The radio frequency range had been neglected in the past because the pulses were assumed to be too weak to detect with antennas of a reasonable size. Nowadays, however, due the really good performance of the Matched Filter as a detector [1] [16] and the faster evolution of the instrumentations on pulsar receivers, the goal of using the Radio Pulsar for real time navigation appear as a promising trend. Furthermore, not only can be used to localize a spacecraft but to localize targets on Earth. Therefore, pulsar-based navigation system can be a substitute of the actual navigation systems as GPS or Galileo. [7] provides an overview of the work that has been done on pulsar navigation and shows this new direction in pulsar-based navigation research. Since pulsar signals offer such a high stable periodicity, the idea is to use them as beacons for Time of Arrival (TOA) based navigation. There two kinds of navigation algorithms that use the extremely accurate periodicity of Pulsars. The first is the Doppler Shifted Navigation, which uses the Doppler effect in the estimated TOA’s in order to localize the target. The other technique is called Period Decay method and uses the period decay of the pulsars. These two techniques are explained in detail in [18]. Although some research has been made about building the actual pulsar navigation system still no practical implementation has been done. In following chapters we will propose a receptor that includes the optimal detector for the pulsar case and the signal processing techniques that improve the detection performance. 2.3.1 Navigation system challenge Although seemingly simple in principle, there are several hurdles that are needed to be overcome in realizing such a navigation system. The main challenges when using radio pulsars 9
for navigation are the following: 1. The extremely weak pulsar signal strength that is being used to navigate. The received signal is completely submerged in Additive White Gaussian Noise, so the Signal-Noise Ratios are very small on Earth. This is because the pulsar signals are emitted many light years away from Earth, leading to addition of noise and distortion due to the propagation channel as will be discussed later. 2. The requirement of receive, process and detect the radio pulsar signal in a few seconds in order to perform a real time navigation. Up to now the processing time to localize the signal it is higher than 10 minutes. Furthermore, it is needed really big antennas (+10 m diameter dish-antennas) to receive the radio Pulsar Signal with enough SNR. Nevertheless, due the improvement of the technology (computers and devices with higher computational costs), the antennas (possibility to reach higher Bandwidth) and the left of the receptor devices, someday it will be possible to achieve the goal of a pulsar-based navigation system. 2.4 Pulsars Conclusions In this chapter I have introduced a brief explanation about what is a pulsar and a classification of them depending on the signal they emit. We also have seen the main characteristics of the Radio Pulsar signals. That signals appear to be extremely periodic pulses that arrive to the Earth with a very weak intensity and submerged in Additive Gaussian White Noise. However, due their accurate periodicity they have been chosen as a perfect candidates for a real-time navigation system. Moreover, the Wideband feature of that pulses has been shown. Due that Wideband nature, we are able to receive them in a very big frequency range, from kHz to THz. That will be an important assumption in order to design the optimum signal processing techniques. Also it has been stated that Pulsar signals are highly polarized, hence an acquisition of two orthogonal polarizations is enough to obtain the power pulse profile. After that, the propagation effects has been explained. The ISM causes several frequency dependent unwanted effects on the wideband pulsar signal, such as dispersion, scattering and scintillation. However, some practical solutions have been given in order to avoid those effects. For example, choosing a adequate observation frequency. Moreover, it has been introduced two de-dispersion methods. Although it looks like they will be an important block of our receptor, we will see as this statement is not true. Finally, it has been explained the goal of implementing a radio real-time pulsar-based navigation system. The two main challenge of this aim has been introduced. That challenge are the low-intensity of the received pulsar signal and the requirement of process and detect the signals in few seconds. After that, two navigation algorithms has been introduced. That methods are the Doppler Shift and the Period Decay, and they both use the extremely periodicity of the radio pulsar signals as a key to localize a target. In the next section, a summary of the detection theory written by [1] and [22] will be explained in order to design a detector to estimate correctly the time of arrival of the pulsar. 10
Chapter 3 Detection Theory applied to Radio Pulsar Signals Detection theory deals with techniques to determine how good data obtained from a certain model corresponds to a given data set. An example of that can be the radars, where the presence of a target has to be detected. Another example could be to detect whether a 0 or 1 has been sent in a communication system. In this Thesis we will only deal with the detection of a pulsar signal over one period in presence of noise. Furthermore, we will assume that the pulsar signal is deterministic, so the detection will be easier to perform. Otherwise the process will be a detection with random processes. In this section I will summarize the detection theory applied to radio pulsar signals done by Richard Heusdens in [1] and [22]. As stated, Radio Pulsar emit a high polarised pulses extremely periodic. But, as these stars are located millions of km from the Earth they arrive with a very low intensity. Moreover, when we receive those signals only noise can be observed due the fact that they arrive submerged in Additive White Gaussian Noise with a very low Signal to Noise Ratio. In order to achieve our goal of using the Radio Pulsar signal for navigation applications we should assure that we estimate in a correct way the time of arrival of the pulses. First of all, in order to find the optimum detector for our applications a basic detection theory for deterministic signals is introduced. The fact that we know exactly how the radio pulsar signal is will be very important to choose a detector. Then, the detection theory for radio pulsars signal application will be explained. As we can guess, besides detect the pulsar, we will need to estimate the amplitude of the pulsar profile in order to implement a template and the Time Of Arrival to localize our target. 3.1 Basic Detection for deterministic signals The detection for deterministic signals is the simplest case because the prior we know about the signal play in our favour. The main idea behind the detection process is the statistical hypothesis testing. Given a data set and different hypothesis our aim will be determine which model fits the data best. Due the fact that we want to detect one signal (the one of the pulsar we want to use for localization), we will only consider in this Thesis two Hypothesis. The first hypothesis H0is the case when only random noise is received. In the second Hypothesis H1the deterministic signal is received in presence of the same random process. We will assume that the random process is an additive gaussian noise with 0 mean 11
and covariance matrix Φz. So, we can model our case as: H0:y(n) = z(n)n= 0,1,2, ..., N −1 H1:y(n) = s(n) + z(n)n= 0,1,2, ..., N −1 Being y(n) the discrete received signal, s(n) the discrete deterministic signal and z(n) the noise process modelled as N∼(0,Φz). The focus will be put in the probability density functions of the both hypothesis. That will be useful in order to choose one or the other hypothesis depending if the received belongs to the pdf of the first or second model. In the Hypothesis H0only noise is received. Then we can state that pdfois: pdf0(y) = 1 (2π)N 2|Φz|1 2exp(−1 2(y−µ)TΦ−1 z(y−µ)) = 1 (2π)N 2|Φz|1 2exp(−1 2yTΦ−1 zy) Being µthe mean of y when we only receive noise. As the H1will be the hypothesis when we receive the deterministic signal submerged in gaussian noise, the probability density function of these model will be also gaussian with a non-zero mean. Therefore, pdf1(y) = 1 (2π)N 2|Φz|1 2exp(−1 2(y−µ)TΦ−1 z(y−µ)) = 1 (2π)N 2|Φz|1 2exp(−1 2(y−s)TΦ−1 z(y−s)) In the next figure 3 we can observe the probability density functions of one both hypothesis. They are almost identical due the fact that they have the same random Gaussian process. The only difference is that the one of the hypothesis 1 is shifted s (being s the mean of the Hypothesis H1). Figure 3. Distributions of the Hypothesis H0and H1. Credit Dr. Richard Heusdens The detector will consist on choosing one of the distributions depending on the data set we receive, so a threshold will be needed. Looking at the figure 4 we can see how there are two possible mistakes we can make when we assess the detection with the threshold. The first one is called miss (II) and is produced when you choose for the hypothesis H0and, in fact, the deterministic signal is received. The other possible mistake is called false alarm 12
(I) and consist on deciding that you have received the deterministic signal (Hypothesis H1) although it is not true. These errors are unavoidable to some extent but may be traded off against each other by adjusting the detection threshold. It is not possible to reduce both errors at the same time once the probability density functions are set. As a consequence, a typical approach to design an optimal detector is to fix one error probability and minimize the other. In the next figure we can observe the Hypothesis testing errors and their trafe-off. Figure 4. Hypothesis testing errors and their trade-off adjusting the detection threshold. Credit Dr. Richard Heusdens To assess that probabilities we will assume a simplest model where we only receive one sample with amplitude s and that the noise process z is still a gaussian process and white with variance σ2 z1. So, z1∼N(0, σ2 z1). So the probability density functions will be: pdf0(y1) = 1 √2πσ2 z1 exp(−y2 2σ2 z1 ) pdf1(y1) = 1 √2πσ2 z1 exp(−(y−s)2 2σ2 z1 ) Now it is possible to compute the probability of false alarm as the probability that y1is bigger than the threshold and in fact we are not receiving the desired signal. That can be stated as: PF A =P(y1> γ |H0) = Z∞ γ 1 p2πσ2 z1 exp(−y2 2σ2 z1 )dy =Q(γ σz1 ) (3.1) Being Q(y) the Q-function or complementary cumulative distribution function related to the complementary (Gauss) error function by: Q(y) = 1 2erfc(y √2) 13
for the different shape transformations. The best results are obtained for the smoothed pulsar profile. All-pass filtering, however, has no effect on the performance, as expected. As the randomized signal is the one with a more independent variables (a really high and narrow peak in the correlation function), Dr. Heusdens has proved how the fact that the TOA is maximal for statistically independent random variables does not imply that this results in the best estimation. Moreover, as expected, the TOA estimation performance increases with the ENR. Figure 6. MLE of the TOA τo. Credit Richard Heusdens 3.3 Conclusions of the Detection Theory From that section we have obtained a lot of important conclusions in order to make feasible our aim of real time navigation using Radio Pulsar Signals. First, we have stated that the best approach for our case will be the performance of Neyman-Pearson detector. That is because the Neyman-Pearson gives the best detection performance once the probability of false alarm is fixed. However, as we need to estimate the Time Of Arrival of the signal, it is not possible to use that detector. So, it has been found a suboptimal detector, the GLRT. Assuming that the received signal is the radio pulsar pulse submerged in white noise, the GLRT will correlate N times (being N the length of on period) that received signal with the shifted replica or template of the radio pulsar Dτs. After doing that, we will estimate our 20
TOA as the shifted value τthat has given the highest correlation. Besides, we will decide that the signal has been detected if that correlation value is above a threshold γ. If the noise is colored, we have seen as the Detection Performance will depend on sΦ−1 zs. So, the implementation of the GLRT when the noise is colored will be a Whitening Process followed by N Correlators. Furthermore, it has been found the expressions for the detection performance and the estimation of the time of arrival, stating that we have to increase the ENR (if the noise process is white) or the GENR (if the noise process is colored) as much as possible. In the previous Thesis and Reports the aim has been to find the techniques that improve the SNR of the signal. Nevertheless, the SNR will not have any effect in the detection performance. Also it has been shown that for a ENR-GENR of 16 dB, the detection of the pulsar and the estimation of the time of arrival will be almost perfect. So, in the next chapters our aim will be to find which signal processing techniques improves the GENR of the process assuming that the Radio pulsar signals have a Wideband nature. The detection theory explained added to the signal processing techniques will lead to a theoretical optimum receptor for the detection of radio pulsar signals. That detector is proposed in the chapter 8. The last conclusion and more important of that section is the fact that the GLRT detector and MLE estimation of the TOA processes does not depend on any all-pass filtering operation. As the de-dispersion process can be modelled as a all-pass filtering, that means that performing the de-dispersion method will not have any effect in the detection performance neither the estimation of the TOA. That is a really important result since up to 80% of the computational cost dedicated to process the radio pulsar signal is spent in that process. 21
Chapter 4 Radio Pulsar Signal Model Figure 7. Block Diagram of the receiver and the Analogue to Digital Converter This section focus in the original properties of the Radio Pulsar signal and the Noise process. From now on I will assume that the received Radio Pulsar Signal is a deterministic wideband signal with unknown amplitude and Time Of Arrival submerged in Additive White Gaussian Noise. Therefore, we can model it as: yr(t) = asr(t−τ) + zr(t) Being s(t) the Radio Pulsar signal, z the noise process, a the unknown amplitude and τ the unknown Time Of arrival. In the Detection Theory section we have seen how to estimate the Amplitude and the time of arrival of the signal. In the next subsections we are going to assume that the amplitude and the time of arrival are known in order to compute in an easy way the ENR and GENR of the signal. Hence, the signal before passing through the analogue Low-Pass Filter will be: yr(t) = sr(t) + zr(t) 22
4.1 Filtered Analogue Signal Figure 8. Power Spectrum of the AWGN and its Autocorrelation function From now on, this Thesis will refer to the analogue filtered signal as ya(t) = sa(t) + za(t) Being yathe signal received and filtered with an analogue Low-Pass Filter, sathe deterministic wideband signal and zathe Additive White Gaussian Noise. As we assume that the Radio Pulsar signals are completely known, we model them as deterministic with an energy a=R∞ −∞ |sa(t)|2dt =1 2πfsR∞ −∞ |Sa(w)|2dw The second term of the equation is only true if we are working with angular frequency in rad/s. On the other hand, the noise is a stochastic White Gaussian signal N∼(0, σ2 za) with 0 mean and variance σ2 za. As state before, yais composed by the pulsar and noise signal, Band-limited to the Bandwidth of the analogue low-pass filter. In this first part of the Thesis we will assume that the antenna is ideal with a Bandwidth equal as the cut-off frequency of the Analogue Low-Pass Filter. Later on, we will observe how under that assumption the analogue Low-Pass Filter is useless from a detection performance point of view. As the noise is the principal problem why we can’t correctly detect the signal, we are going to focus on its properties. As stated, the received noise is a continuous time Additive White Gaussian process with 0 mean, variance σ2 zaand spectral density No. Therefore, the power spectrum of the noise is: Sza(w) = Nofor |w| ≤ 2πB 0 otherwise (4.1) where wis the angular frequency in rad/s. If we compute the inverse of the Fourier Transform of the Power Spectrum we will obtain the autocorrelation function of the noise: 23
Rza(t) = 1 2πZ∞ −∞ Sza(w)ejwtdw =1 2πZ2πB −2πB Noejwtdw =No πt sin(2πBt) = 2BNosinc(2πBt) (4.2) Therefore, the variance of the noise process is σ2 za=Rza(0) = 2BNo, with B the Bandwidth of the signal and Nothe Power Spectral density of the noise. Looking at the Power Spectrum and the Autocorrelation function of the noise processes after passing the signal through some processing techniques we will be able to know the variance of the noise and calculate the variations of the ENR. 4.2 A/D converter Just before starting to process the signal, it is passed through an the A/D Converter with sampling frequency fsin order to work in the digital domain. The reason of sampling the signal is to have more facilities to process it with a lower cost and a flexible digital hardware. At the end of the A/D converter we will have: y(n) = ya(nTs) = sa(nTs) + za(nTs) = s(n) + z(n) where Tsis the inverse of the sampling frequency fscalled sampling period. Hence, the discrete signals sand zare obtained by sampling saand zawith the sampling Nyquist rate fs= 2B. We use this sampling frequency in order to avoid aliasing and keep the white property of the noise process. Later on we will see how important is to have AWGN from a processing time point of view. As we have seen in the last section is used the un-normalised frequency was the variable for the representations of the signals in the frequency domain, being w= 2πf. To change between time and frequency domain it is used the Discrete time Fourier Transform: y(n) = 1 2πfsR2πfs 0Y(w)ejwnTsdw Y(w) = P∞ n=−∞ y(n)e−jwnTs First of all, to assess the ENR of the process ywe will focus on the deterministic signal s. As it is known, when a signal goes into an A/D converter, it loses information, and its visible spectra will be limited to the sampling frequency, 2πfs. Besides, in the frequency domain appears images of the signal every 2πkfs(with ka integer number, creating aliasing if the Nyquist Sampling Frequency, or a bigger rate is not used). As we stated before, the Nyquist sampling frequency is used in this section, so no aliasing will be obtained. Following that, the spectra of the discrete signal compared to the spectra of the previous analogue signal is: 24
S(w) = P∞ n=−∞ s(n)e−jwnTs=fsP∞ n=−∞ Sa(w+ 2πkfs) So we can find an expression of the energy of the deterministic sampled signal sin comparison with the energy of the analogue filtered signal saas: s=X n|s(n)|2 ∗ =1 2πfsZ2πfs 0|S(w)|2dw =fs 2πZ2πfs 0|∞ X k=−∞ Sa(w+ 2πkfs)|2dw ∗∗ =fs 2πZπfs −πfs|Sa(w)|2dw =fs 2πZ∞ −∞ |Sa(w)|2dw =fssa (4.3) * By Parseval Theorem that states that the energy of the signal in time is preserved in the frequency domain as well. ** Taking into account that Sa(w) do not have any contribution in w /∈[−πfs, πfs]. No aliasing. In this formula has been assumed that the bandwidth of the signal sadoes not exceed half of the total visible spectra 2πfs. As we can observe, the energy of the sampled signal increases with the sampling frequency in comparison to the analogue one. This result shows that increasing the fsthe energy of the process is improved, even when the bandwidth of the signal is not increased. On the other hand, to check the behaviour of the discrete time noise z, we are going to look to its Autocorrelation function: Rz(k) = z(n)∗z(n+k) = z(nTs)∗z((n+k)Ts) =Rza(kTs) = 2BNosinc(2πBkTs) = 2BNosinc(2Bπk 2B) = 2BNosinc(kπ) = 2BNoδ(k) (4.4) Being δ(k) the Kronecker delta. That happens because sampling the noise with the Nyquist sampling frequency preserve the Whiteness of the Noise. We can also observe that the variance of the discrete Noise will remain the same as before passing the signal through the A/D converter no matter the sampling frequency used. Therefore, σ2 z= Rz(0) = fsNoδ(0) = 2BN0=σ2 za. Finally, computing the Power spectrum of the Noise as the discrete-time Fourier Transform of the Autocorrelation Function, we have: 25
Sz(w) = fsNofor |w| ≤ πfs 0 otherwise (4.5) As we can observe, the shape of the Power spectrum of the noise remains equal, although the value of the discrete power spectral density is increased by the factor fs. Once obtained the values of the energy of the sampled signal sand the variance of the noise z, ENR is computed: ENR =s σ2 z =fssa σ2 za =fssa 2BNo =fssa fsNo =sa No (4.6) In the next subsection the term GENR is introduced. It is a Ratio to assess the ENR when the noise is not white. As we have seen in the Detection Theory, this Ratio assess the detection performance of the GLRT Detector. As we will see later, one way to improve the GENR is increasing the observation time of the signal to be detected, so sawill be increased. The problem is that we will focus on detecting the Pulsar signal over one period, so we will not be able to increase the energy of the signal integrating over a longer time span. 4.3 Generalised Energy To Noise Ratio To assess the behaviour and the improvement of the signals after performing some signal processing techniques, the scientific use different ratios. The most used is the SNR or Signal to Noise ratio, being the relation between the power of the desired signal and the variance of Noise. SNR =Ps σ2 z (4.7) Although up until now the researchers have been using this radio to assess the Radio pulsar signals behaviour, in this report we are going to use other two ratios, the Energy to Noise ratio and the Generalised Energy to Noise ratio. As it has been explained in the Detection Theory chapter, the detection performance only depend on the Generalised Energy to Noise ratio, and, in some cases, on the Energy to Noise Ratio. The term Generalised Energy to Noise Ratio has not been used before, but it describe the relation between the signal and the covariance matrix of the noise. Then, we are going to refer to Generalised Energy to Noise ratio of the signal y=s+zas: GENR =sTΦ−1 zs(4.8) 26
with Φz=E((z−z)(z−z)T) the covariance of the noise process and s the desired signal. We can see that if the noise is AWGN with 0 mean and covariance matrix Φz=σ2 zI, so Φ−1 z=1 σ2 zI, the GENR becomes the relation between the energy of the signal and the variance of the noise. Therefore GENR =sTΦ−1 zs=sTs σ2 z =s σ2 z =ENR (4.9) The GENR takes into account the color of the noise to assess the detection performance. From now on, Generalised Energy to Noise Ratio will be used to measure up whether the different signal processing algorithms increase or not the detection performance. 27
Chapter 5 Theoretical Signal Processing Techniques Figure 9. Block Diagram of the Signal processing techniques applied to the discrete-time signal The signal from the Radio pulsar is received submerged in AWGN and it is not possible to see or detect it without some processing. In this chapter we introduce some basic signal processing techniques to check if the ENR and GENR of the pulsar signal increase. The algorithms explained will be Epoch Folding, an average of the received signal at the exact period of the pulsar; Low-Pass Filtering, to eliminate the high frequencies of the signal; Downsampling, to change the sampling rate of the signals and therefore, the length of the data; decrease and increase the Bandwidth of the Antenna (and therefore, the cut-off frequency of the Analogue filter); and Oversampling/Undersampling, passing the analogue signal to the A/D converter with a sampling frequency higher/smaller than the Nyquist one fs. Downsampling can also be seen as a process that changes the bandwidth of the antenna and the sampling frequency of the A/D converter, an interesting feature that will help us to increase the ENR of the rotating star pulse without increasing the observing time. That result could be good to perform a real-time detection of the pulsar for navigation applications. 5.1 Epoch Folding Epoch Folding is a Signal Processing technique used to decrease the variance of the uncorrelated noise while keeping the energy/power of the desired signal. Epoch folding consist on choosing a range of periods, and average the data at those periods. The algorithm assumes that we know the periodicity, T, of the signal. The first step consist on breaking the received signal in intervals of time T. Then, sum all these clipped signals together and divide the resultant signal by the number of foldings, K. No matter if the signal is narrowband or 28
wideband, because we are not clipping any frequency spectrum of the signal provided the right period to perform the folding is used. So, the shape, amplitude, energy and power of the signal will remain the same no matter how many folds you do. In the case of Radio Pulsar signals, although they arrive with a very precise periodicity to the Earth, their amplitude can vary significantly over time due to the scintillation effects provoked by the ISM. Nevertheless, the averaged pulsar profile remains very stable, allowing us to perform the Epoch Folding without loosing any information of the signal not in time neither in frequency domain. As it has been explained, the noise will be additive, white and gaussian with 0 mean, variance σ2 zand uncorrelated. Is the last feature the important for the success of this technique, because averaging AWGN uncorrelated noise leads to a linear decrease of the noise variance with the number of folds. Furthermore, the noise is still white after passing through the averaging, so the GENR of the received data will increase. In the Figure 10 we can observe the process of epoch folding. Figure 10. Epoch Folding algorithm performed to a periodic signal submerged in uncorrelated noise. As we are going to perform this algorithm to a discrete signal, let me consider ras the discrete signal with KTfssamples, being K the number of folds and Tfsthe number of samples in one period. The next step is breaking the data in a sequence of discrete signals ykof length L, being L=fsT. As we have stated, yk=sk+zkwith zk∼N(0, σ2 z) an Additive White Gaussian Noise and skthe desired signal with length L. So, performing the Epoch Folding we have: x(n) = 1 K K−1 X k=0 yk(n) =1 K K−1 X k=0 sk(n) + 1 K K−1 X k=0 zk(n) ≈s(n) + z(n) (5.1) 29
From the last expression we can observe as the pre-whitening matrix will be: U=1 σzL−1. Then, taking into account that the filtered signal is p=Ls, the process Up can be written as Up =1 σzL−1Ls =1 σzs. That means that no matter what filter you have in the digital signal processing chain, it will not have any influence in the GENR of the filtered signal. As we have seen, that happens because the pre-whitening process Uis cancelling the effect of the filter once multiplied by the filtered signal. Finally, the energy of that process and as we stated before, the GENR will be: GENR =Up =X n|1 σz s(n)|2=1 σ2 zX n|s(n)|2 =1 σ2 z s=s σ2 z =sa No (5.11) To sum up, we have proved that filtering the discrete process ydoesn’t change its GENR. Therefore, the detection performance is not altered for this operation. So, assuming that the spectrum of the signal is almost flat, the relation between the loss of signal energy and noise variance will be the same. If we look at the de-dispersion section of the detection chapter, we can see as it has been explained that due to the dispersive nature of the interstellar plasma, lower-frequency radio waves travel through the medium slower than higher-frequency radio waves, which manifests itself as phase distortion. Then, as it has been proved that the filtering process does not affect the GENR and that the de-dispersion process can be assumed as a filter, the de-dispersion will not have any effect in the GENR of the radio pulsar signal. That fact agrees with the theory written by Richard Heusdens in [1] and proving that the de-dispersion does not affect the detection performance of the GLRT for any kind of noises. Another surprising thing of that section is the ability of the whitening matrix Uto recover the high frequencies cut to pby the Low-pass filter. Indeed, the matrix U=1 σzL−1 not only stretch the signal, but recover in a good way the shape of the high frequencies of s. So, we can state that Up =s σz. In the subsection 7.1.3 we will show the performance of the whitening matrix and how it recovers the signal spectra in the high frequencies. In the following chapter I am going to explain the process of changing the bandwidth and rate of the signal. The operations are called downsampling, that means decrease the bandwidth and the sampling frequency of the signal. 5.3 Downsampling Downsampling is signal processing operation that change the rate of the signal while keeping the relation Bandwidth-Sampling frequency stable. That means that, for example, downsampling a signal with a Bandwidth Band sampling frequency fs= 2Bby a factor of 36
Mwill be the same as Low-pass filtering the analogue signal with a bandwidth B/M and sampling it with a frequency fsM=fs/M = 2B/M. So if the signal before downsampling has been sampled with the Nyquist sampling frequency, the signal after the downsampling will also have a Nyquist sampling frequency. Hence, the property of whiteness of the noise after changing the rate will not be altered. As I will explain in the next section, this is very important to perform the detection without a lot of computational cost. In the next figure we can observe a representation of the Downsampling process. Figure 14. Downsampling process in the time domain applied to a sinusoid. It will be the same as decreasing the bandwidth of the receiver and the sampling frequency of the converter. If we perform the downsampling algorithm to the received signal y=s+z, we will obtain yd=r+W. Hence, ywill be the upsampled version of yd. The Downsampling process by a factor of M can be divided in two steps: 1) A low-pass filter with a cut-off frequency fc=B/M to eliminate the highest frequencies in order to avoid any aliasing in the next step. 2)Decimating yd(n) = y(nM) the signal by a factor of M As we have done the low-pass filtering with a factor of M before, now we only have to explain the effects of decimation. So, we will start applying the decimation block to the filtered signal yf(n) = p(n) + v(n). The filtered signal has a bandwidth of B/M, so it will not be aliasing after applying the decimation process. First, I am going to focus in the noise process W(n) = v(Mn) = z(Mn). Looking at the autocorrelation function of the downsampled noise: RW(k) = W(n)∗W(n+k) = v(nM)∗v((n+k)M) =Rv(kM) = fsNo Msinc(πkM M) =fsNo Msinc(πk) = fsNo Mδ(k) (5.12) 37
The variance of the process Wwill be σ2 w=RW(0) = fsNo/M =σ2 v=σ2 z/M, the same as the process vand Mtimes smaller than in the process z. So it is clear that the decimating algorithm doesn’t change the variance of the noise. Now, looking at the power spectrum of the noise we can see how the noise became white again after the decimation, although the power spectral density has decreased in a factor of M:SW(w) = PkRW(w)e−jwkTs= Nofs MPkδ(k)ejwkTs=Nofs Me0=Nofs M. Then, the power spectrum will be flat in all the bandwidth: SW(w) = fsNo Mfor |w| ≤ πfs M=πfsM 0 otherwise (5.13) Being fsM=fs/M the new sampling frequency and Bd=B/M the new Bandwidth of the Band-Limited Signal. As we can observe, the noise after the downsampling is still white, so the detection performance will depend on the ENR. In order to assess the energy of the pulsar signal after the downsampling, we have to take into account that downsampling by a factor of M a signal with a bandwidth Band a sampling frequency fscan be seen as the same process as decreasing the bandwidth of the receiver and the sampling frequency of the A/D converter by M. It is something trivial as the Decimation process can be seen as an A/D converter. In the next figure we can see a representation of a receptor with a processed signal r2with M times less Bandwidth than s. Therefore, it is easy to see how that discrete signal r2will be absolutely the same process as r, the downsampled signal. Figure 15. Scheme of the process equivalent to Downsampling the signals sand z. It consist on decreasing the bandwidth of the signal and the sampling frequency by a factor of M. We can see that in the Figures 15 and 9 that the energy of the downsampled pulsar signal rwill be the same as the energy of the new process r2(being r2the process with a analogue low pass filter cut-off frequency fcM= 2πB/M and sampling frequency fsM) (r=r2). In order to prove that, the power spectrum of the RW2(k) is computed: 38
RW2(k) = W2(n)∗W2(n+k) = W2(nTsd)∗W2((n+k)Tsd) =RW2a(kTsd) = 2BNo Msinc(2πBkTsd M) =2BNo Msinc(2BπkM 2BM ) =2BNo Msinc(kπ) = 2BNo Mδ(k) (5.14) Being Tsd the inverse of the sampling frequency fsd = 2B/M. So, the variance of W2will be σ2 W2=RW2(0) = fsNoδ(0)/M = 2BN0/M =σ2 W. Finally, computing the Power spectrum of the Noise as the discrete-time Fourier Transform of the Autocorrelation Function, we have: SW2(w) = fsNo Mfor |w| ≤ πfs M 0 otherwise (5.15) As we can see, W2is the same signal as the downsampled noise process W. That prove the thing that Downsampling by M is the equivalent process as reducing the Bandwidth of the receiver and the sampling frequency by a factor of M. So, knowing that the downsampled pulsar signal rwill be the same process as r2, we can compute the energy of r2. First, we start to see if what is the energy of ra2in comparison with the energy of sa: ra2=Z∞ −∞ |ra2(t)|2dt =Z∞ −∞ |Sra2(w)|2dw =Z2πB M −2πB M|Sra2(w)|2dw ∗∗∗ =1 MZ2πB −2πB |Ssa(w)|2dw =sa M (5.16) *** Assuming that the spectrum of the analogue pulsar signal is almost flat in a bandwidth of B, so clipping it by a factor of Mmeans having the total spectrum of the filtered analogue signal (ra) divided by M. Then, if we compute the energy of r2we have: 39
r2=X n|r2(n)|2 =1 2πfsMZ2πfsM 0|Sr2(w)|2dw =fsM 2πZ2πfs2 0|∞ X k=−∞ Sr2a(w+ 2πkfsM)|2dw =fsM 2πZπfsM −πfsM|Sra2(w)|2dw =fsM 2πZ∞ −∞ |Sra2(w)|2dw =fsMra2=fsra2 M =fssa M2=s M2=r (5.17) As I said before, r=r2, so r=p M=s M2. Due to the Wideband Nature of the pulsar signal, the energy will decrease in an order of M2when you downsample the signal. Therefore, it will decrease in an order of Mwhen you low-pass filter (as we have seen in the last section) and in an order of Mwhen you decimate. However, the variance of the noise has only decreased in an order of Mduring the downsampling. That is due to the fact that the variance is power and it takes into account the energy per sample. If we compute the ENR of the downsampled process yd=r+Wwe will obtain: ENR =r σ2 W = p M σ2 v = s M2 σ2 z M =fssa MfsNo =sa MNo (5.18) Then, after downsampling the pulsar signal the ENR becomes smaller. To assess the detection performance we are going to look at the GENR. GENR =rTΦ−1 Wr=rTr σ2 W =r σ2 W =sa MNo =ENR (5.19) As we were expecting, the value of GENR is the same as the ENR because the noise process Wis white. So, we can state that the downsampling and/or decreasing the observing bandwidth of the received signal makes the GENR and ENR decrease. Hence, deteriorate the detection performance. If we realize that increasing the observing bandwidth (increase the cut-off frequency of the analogue low pass filter and the sampling frequency) is the inverse operation of the downsampling, we can state that increase the observing Bandwidth improves the GENR and the detection performance. So, from now on we have two different ways to improve the detection performance: 1) Epoch Folding 40
2) Increase the observing bandwidth of the receiver and the sampling frequency of the A/D converter Some experiments about the improvement of the GENR due the increment of the Bandwidth of the signal will be shown, as well as the performance of the other explained algorithms. As we can’t perform a Epoch Folding with a really big number of folds because we need to make feasible the real-time navigation, the main solution to improve the detection performance will be increase the bandwidth of the antenna. Now that solution doesn’t allow us to increase the GENR too much because of the technological limits. Nevertheless, every year the bandwidth of the receiver is increasing with technological improvements, so the detection performance will be increasing along time until it will be feasible to do a real-time navigation. To finish this signal processing background I am going to introduce the oversampling and undersampling techniques. That consist in increasing or decreasing the sampling frequency of the A/D converter without changing the bandwidth of the antenna. So, as we can expect, the noise will be colored although it will keep constant the gaussian and 0 mean features. We will see how undersampling doesn’t vary the detection performance under some assumptions, and how depending on the way you perform oversampling, you can improve the detection performance. 5.4 Oversampling Figure 16. Scheme of the process of Oversampling. It consist on keeping the bandwidth of the receiver while increasing sampling frequency by a factor of M over the Nyquist one. Oversampling is a signal processing technique that consist in increasing the sampling frequency while keeping the Bandwidth of the signal. That means sampling the filtered analogue signal sao, with Bandwidth B, with a sampling frequency higher than 2B. As we can see in the figure 16, oversampling will increase the sampling frequency above the Nyquist one. One of the reasons to perform oversampling is that the energy of sowill increase in comparison with swhereas the variance of the oversampled noise zowill not change. In that 41
section we will assume that the Bandwidth of the Antenna is BM. Later on, we will see how important is that in a detection performance point of view. First of all, to assess the energy of the signal and the variance of the noise, we shall notice that the energy of the oversampled analogue signal yaowill be the same as the analogue signal of the first chapters ya. That is because we are band-limiting it with the same bandwidth B. Hence, sao=sa. The variance of the noise zaowill also be the same as the variance of za. So, focusing in the oversampled process yo=so+zoafter the sampling, we will have a signal with a visible spectrum of fso=fsM= 2BM band-limited with a bandwidth B. Before calculating the energy of the oversampled pulsar signal and the variance of the oversampled noise, let’s compute how the autocorrelation function and power spectrum of zrand analogue filtered signal zaoare: Szao(w) = Sza(w) = Nofor |w| ≤ 2πB 0 otherwise (5.20) Rzao(t) = Rza(t) = 1 2πZ∞ −∞ Sza(w)ejwtdw =1 2πZ2πB −2πB Noejwtdw =No πt sin(2πBt)=2BNosinc(2πBt) (5.21) So, the variance of the noise process zaois σ2 zao=Rzao(0) = 2BNo=σ2 z. In order to compute the autocorrelation function and the Power Spectrum of the analogue unprocessed noise zr, we have to take into account that is the same process as zao, but with M times more Bandwidth. Szr(w) = Nofor |w| ≤ 2MπB 0 otherwise (5.22) Rzr(t) = 1 2πZ∞ −∞ Szr(w)ejwtdw =1 2πZ2πMB −2πMB Noejwtdw =No πt sin(2πMBt)=2MBNosinc(2πMBt) (5.23) As we can see, the variance of that noise process is σ2 zr=Rzr(0) = 2MBNo=Mσ2 z. Hence, M times higher than the variance of the Low-Pass Filtered analogue process. That is due the effect of the Low-Pass Filter that is cutting the Bandwidth of the received analogue signal by a factor of M . In the figure 17 we can observe the shape of the Power Spectrum and the autocorrelation function of the noise zr. 42
Figure 17. Representation of the Power Spectrum and the autocorrelation function of the analogue noise process zr. Now, if we sample the signal yaowith a sampling frequency fso=fsM= 1/Tso, the resultant process will be yo(n) = yao(nTso) = sao(nTso)+zao(nTso) = so(n)+zo(n). In order to assess the energy of the oversampled pulsar signal, its spectra is shown: So(w) = ∞ X n=−∞ so(n)e−jwnTso=fso ∞ X k=−∞ Sao(w+ 2πkfso) (5.24) So we can find an expression of the energy of the deterministic signal soin comparison with the energy of the analogue filtered signal saoand sr: so=X n|so(n)|2 =1 2πfsoZ2πfso 0|So(w)|2dw =fso 2πZ2πfso 0|∞ X k=−∞ Sao(w+ 2πkfs)|2dw =fso 2πZπfso −πfso|Sao(w)|2dw ∗ =fso 2πZπfs −πfs|Sao(w)|2dw =fso 2πZ∞ −∞ |Sao(w)|2dw =fsosao=Mfssa=Ms (5.25) * Assuming that the spectrum of the discrete-time oversampled signal is has a bandwidth of πfsbecause of the analogue low-pass filter. ** Because the energy of the analogue process sris M times higher than the energy of saoand sadue the effect of the Low-Pass Filter and assuming that the receiving pulsar signal has a Wideband nature. Also fso=Mfs 43
Figure 18. Power spectrum of the Oversampling Noise. Hence, as the sampling frequency fsois bigger than the Nyquist one, the energy of the oversampled pulsar signal is also higher by a factor of M. To check the behaviour of the discrete time noise zo(n), we are going to look to its Autocorrelation function: Rzo(k) = zo(n)∗zo(n+k) = zo(nTso)∗zo((n+k)Tso) =Rzao(kTso)=2BNosinc(2πBkTso) = 2BNosinc(2Bπk 2BM ) = 2BNosinc(kπ M) (5.26) Being Tso= 1/2BM. That happens because sampling the noise with a sampling frequency higher than the Nyquist one eliminates the whiteness property of the noise. Nevertheless, the variance of the sampled noise will remain the same as before passing the signal through the A/D converter no matter the sampling frequency used. So, the oversampled noise zowill be the same as the variance of the normal process z. σ2 zo=Rzo(0) = fsNo= 2BN0=σ2 z(5.27) Finally, computing the Power spectrum of the oversampled noise as the discrete-time Fourier Transform of the Autocorrelation Function, we have: Szo(w) = fsoNo=MfsNofor |w| ≤ πfs 0 otherwise (5.28) 44
Where we can observe as the noise is not white anymore. Moreover, we can see that in the Figure 19 where it is shown the power spectrum of the colored oversampled noise. Now that we have the energy of the pulsar signal and the variance of the noise, we can compute the ENR: ENR =so σ2 zo =fsosao fsNo =Mfssao fsNo =Msao No =Msa No (5.29) But as I explained before that does not mean that the detection performance is improved. To assure that, we have to compute the GENR. As the noise is not colored, we will do the same procedure as in the Low-Pass Filtering section to calculate it. In the Low-Pass Filtering section we have seen as the digital filters had no effect in the GENR of the process. So, in order to being able to compute the covariance matrix of the oversampled noise, needed to assess the GENR, some mathematical stuff will be shown. Knowing that sampling and Low-Pass filter are lineal operations and looking at the equations stated in that section, we can write the oversampled signal as: yo(n) = so(n) + zo(n) = sao(nTso) + zao(nTso) = Lsr(nTso) + Lzr(nTso) Being L the filtering matrix taken from the filter h(nTso), and h(nTso) the sampled version of the analogue low-pass filter h(t) with cut-off frequency of B Hz. Do not confuse the processes sr(nTso) or zr(nTso) with sr(t) and zr(t), as the first ones are the sampled version of the second ones. That means that sr(nTso) and zr(nTso) are discrete signals. Hence, as proved before, the energy of sr(nTso) will be fsotimes the energy of r(t). Another important feature of sr(nTso) or zr(nTso) is that they have a Bandwidth MB, so the sampling frequency fsowill be, in fact, the Nyquist sampling frequency of theses processes. As the reader already knows from the previous chapters, that mean that the noise process zr(nTso) is white with the same variance as the analogue process zr(t). That is due the fact that sampling the noise does not vary the variance. So, the variance of that process will be σ2 zr(nTso)=σ2 zr(t)= 2MBNo. That formulation will be very useful when computing the GENR of the oversampled signal. Then, taking into account the results we had in the section 5.2, the covariance matrix of the oversampled noise and pre-whitening matrix U of the noise process zowill be: Φzo=E(zozT o) = E(zao(nTso)zao(nTso)T) = E(Lzr(nTso)zr(nTso)TLT) =LΦzr(nTso)LT∗ =σ2 zr(nTso)LLT(5.30) * Knowing that the process zr(nTso) is white gaussian noise with variance σ2 zr(nTso)= 2BMNo. 45
real pulsar signal will not have a completely flat spectrum, so the aliasing will not be added linearly to the spectrum and the performance will go down. . Finally, the Integration in Time has been introduced. That technique will reduce the detection performance. However, it also will highly decrease the number of samples of our data. Therefore, Integration in Time may be an interesting technique in order to decrease the computational complexity of the whole receptor. 52
Chapter 6 Radio Pulsar Signal PSR B0329+54 observations 6.1 PSR B0329+54 features In order to prove the theoretical background stated in the last chapter, the results of some experiments with simulated and real data are shown. But, first of all, I will explain the features of the Radio Pulsar signal PSR B0329+54 and the way it has been recorded. Table 2. Parameters of the Radio Pulsar B0329+54 PSR B0329+54 is a neutron star situated approximately 2,643 light-years away from the Earth in the constellation of Camelopardalis and it was created 6.74 millions of years ago. In 1979, two extrasolar planets were announced to be orbiting the pulsar (being classified as pulsar planets). Later observations however ruled out this idea. These radio pulsar emits one of the strongest polarized pulses received in the north hemisphere with a periodicity of 0.71451866398 s. Furthermore, the fact that it has a really low dispersion and an almost negligible spin down, makes it a good candidate to perform experiments. The Dispersion Measure of that periodic signal is 26.776 cm−3pc, a low enough value that allow us to avoid performing the de-dispersion process. Due the intensity of that Radio Pulsar we can make out its pulsar profile after some foldings. The pulsar has an average flux density at the observation frequency fo= 1400 MHz of 203 mJy and an average flux density of 1650 mJy in a fo= 400 MHz. As we can observe in the next figure, the star has three nested cones of emission and a central core emission. Also we can see how the pulsar is visible and almost identical in all the observed frequencies, from 117 MHz to 1170 MHz. 53
Figure 21. Profiles of the Radio Pulsar B0329+54 with different observed frequencies. This neutron star can be classified as a normal pulsar as it is not a millisecond pulsar. However, the period of that pulse is enough to being able to perform hundreds of folds without losing many time. Although this pulsar is not strongly affected by dispersion, it is known that it scintillate a lot. Therefore, the amplitude of the pulsar signal will vary over time. That is not something we have to worry about due the fact that the integrated pulsar profile (after folding) is quite stable. In the EPN database we will be able to find integrated radio pulsar profiles of the B0329+54 observed in a different frequencies and recorded with a different Bandwidth. 6.2 Radio Pulsar Signal PSR B0329+54 data acquisition from WSRT On February 2nd 2012, the Westerbork Synsthesis Radio Telescopes observatory recorded data from the radio pulsar PSR B0329+54 for a total observation time of 140 s. The Westerbork Synthesis Radio Telescope (WSRT) is an aperture synthesis interferometer near camp Westerbork, north of the village of Westerbork, Midden-Drenthe, in the northeastern Netherlands. It consists of 14 dish-shaped antennas. The operator in the control room has a good view of the dishes in the array. By means of a variety of computers it is possible for the operator to control the telescopes, receivers, and everything in the observing system. In the control room are instruments, which convert the signals to digital information to be read and processed by a computer. The software that has been specially developed for this purpose is so clever that it makes the 14 dishes look like one large dish. The acquisition was a demo observation planned in order to obtain test data with a high bandwidth. The signal was recorded at the observing frequency of 1330 MHz with a Bandiwdth of 20 MHz. After that, the signal was sampled with a sampling frequency of 40 MHz 54
in order to use the Nyquist sampling frequency and preserve the whiteness property of the received noise. Taking into account the period of the pulsar, 196 complete periods can be extracted from the data. This acquisition was stored into 14 files in .dada format with a total weigh of 10.4 GB. Each file consists of 4096 bytes of header and then 800000000 bytes of X and Y polarization real voltage signal samples interleaved. The format used for each sample is signed integer with little endian byte ordering: the digital dynamic range goes from -127 to 127, but no information about the amplitude of the voltage signal can be extracted from there. The header contains all the information about the acquisition besides the number of the file being open. Each file has 10 seconds of data and the acquisition is consecutive, meaning that from the ending of one file to the beginning of the next no data is lost. The signal is formed as samples of real voltage signals in X and Y polarization interleaved (XYXYXYXYXY). Those signals has been collected by the WSRT with 14 Telescopes. As it has been stated before, some beamforming techniques has been applied in the control room in order to have only one recorded output signal. The purpose of that is to take profit of the features of the Telescopes of WSRT and to increase the SNR and ENR of the signal with the beamforming techniques. Apart from that, the sampled voltage signal is not processed in any other way. That means that it is an almost pure, raw signal acquisition. In the next sections we will assess the GENR of the recorded pulsar signal with the help of the template taken from the EPN database. That template will allow us to compute the variance of the noise subtracting the pulsar signal from the acquisition. But before that, some adjustments to the template have to be done in order to have the same amplitude and time of arrival than the received pulsar signal. 6.2.1 Process to visualize the signal As I have stated in the last paragraph, we have 14 files .dada with 140 s of recording of the pulsar B0329+54. Taking into account to the period, we can extract 196 complete periods of that data. However, as we are going to perform Epoch Folding with all those periods, is important to check if the spin down of the pulsar signal will vary our period over time. Before it has been stated that the spin down of the radio pulsar is very small (in the order of 10−15 ss−1), but not for that reason we have to forget it. In the case we perform Epoch Folding over a long time without changing the initial period, we will have a misaligned of pulsar signals added together in the wrong position. That will lead to a bad performance of Epoch Folding and a change of the pulsar profile. As we do not want that to happen due the fact that we need a clean pulsar profile to perform the Matched Filter, we are going to check the true period of the pulsar signal. [17] explains the procedure to calculate the correct period of the Radio Pulsar Signal depending on how many folds you perform. During the first folds the Ptrue will be almost the same as the original period. But, after 50 folds the folded signal start drifting to the right, the profile starts to broaden and the pulsar signal profile starts to blur. So, as we can see in [17], the correct period for 196 foldings should be Ptrue = 0.7145579 s and the number of samples per period will be 28582316. That change in the period will be a problem in order to read the data. The initial data has been recorded in order to have a period of 0.71451866398 s. So, in every file there are 14 exact periods of 28580746 samples. Now, if we change the length of the periods, we will have to split the reading of the 14th period of every file in two steps. Finally, that will provoke to have 195 whole periods in all the recorded data, having to left without reading a little portion of the 196th period. Although it look 55
like a problem, at the end we will see as the final folded pulsar profile will be very accurate. First, the two polarized voltage shas been integrated in order to have the pulsar power profile. To do that, both voltages signal will be combined in this way: Pprofile =|Vx|2+|Vx|2 That profile will be a signal with a range of values between [0,16128]. As it is very complex from a computational point of view to work with a signal with 195 periods, we will perform the next techniques to one period. So, in case we want to perform the Epoch Folding, that will be the moment. That is to say, if we want to perform Epoch Folding, we will perform it in this part of the process to adequate the signal. So, from now on, our signal will be modelled as: y(n) = s(n) + z(n) Being s the pulsar signal and z Additive White Gaussian noise with unknown mean over one period. Although is not something critic, it will be clearer to work with the exact mean of the radio pulsar signal, so we have to find the mean of the process y. To do that, we have to take into account that the received noise process have 0 mean, so the mean of the total signal will be: E(y) = E(s+z) = E(s) + E(z) = E(s) Hence, we need to know the mean of the radio pulsar signal in order to compute the mean of the whole process. To do that, we need to use the template of the pulsar signal B0329+54 taken from the EPN database. That template is a signal with 1024 samples with an amplitude Aand a time of arrival τ. In order to know the exact amplitude, and therefore, the mean of the received signal, first we have to estimate the values of the amplitude and time of arrival of the template. Furthermore, as our signal yhave 28582316 samples, first we have to upsample the template. The next step was estimating the time of arrival of the radio pulsar signal. As it has been stated in the detection theory, the best way to do it is performing N matched filters of the templates and the radio pulsar signal and see in what position the value is maximum. After doing that, it was clear that the signal was shifted 335000 samples, so I shifted the template that number of samples to being able to compute the better ENR possible. The last step to being able to have the real mean of the signal was estimating the amplitude of the radio pulsar received signal. As stated in the detection theory, the value of the MLE of the amplitude will be ˆa=<y,p> <p,p> With that processes, we can state that the resulting template pis the best estimation of the received pulsar. This process has to be repeated if we use a folded signal or if we change the number of foldings. The time of arrival will not vary with the folding, though, so we do not have to compute again the matched filter. So, if instead of working with the unprocessed signal we start working with a K folded signal, we have to estimate the amplitude of the template again. 56
From now on, we can visualize and work with the signal yand the template p. In the next section I will show the visualization of the radio pulsar signal, as well as compute the amplitude, mean, energy of the pulsar signal, the variance and spectra of the noise and the ENR of the whole process. 6.2.2 Representation of the radio Pulsar Signal In this section I am going to show the plots and the features of the radio pulsar signal B0329+54 described before. First of all, we are going to assume that our unprocessed signal is y=s+z, being s the pulsar signal and z the Additive White Gaussian Noise. Also we have the Template p, with the same time of arrival, shape and amplitude as the received pulsar signal. The process described in the last sections have been performed in this template in order to have a replica of the received pulsar signal. So, we can state that p=s. Therefore, the GENR and ENR of the process can be assessed. In fact, the value of the GENR will be the same as the ENR due the fact that the noise is white. We will see that in the next figures, where the Power Spectrum of the noise process will be shown. So, finally we can see as the GENR: GENR =ENR = 39.9dB In the next figure we can observe the plots of the noisy signal and the Template in the time domain. Figure 22. Representation of the unprocessed noisy pulsar signal (left) and the radio pulsar signal Template (right). From that unprocessed noisy pulsar signal we can see as the radio pulsar signal is completely submerged in noise. Moreover, there are 4 interference with a big amplitude that has not interest. As we will see in the results of the experiments in the Chapter 7, those interference will disappear with the performance of the Epoch Folding technique. Then, if we look at the Template we can see as the estimating amplitude is around 30 while the maximum amplitude of the noise is about 2000. That shows how small are the unprocessed pulsar signal if we compare it with the noise. If we keep looking to the template, we can 57
see the power pulsar profile, and how it has the same features that have been described in the last sections. Now, after looking the signal in the time domain, the power spectrum of the noise process zand the noisy signal yis computed in order to see the frequency domain characteristics. Due to the fact that with Matlab is difficult to show the real spectrum of the signal, I am going to compute the power spectrum of the real data. To compute it, the periodogram approximation is used: 1) Compute the FFT of the time-domain signal and take the absolute value 2) Divide it by the square root of the number of samples 3) Apply an average filter 4) Do the square of the filtered signal Finally, we have 1 L|Sxx(f)|2, being Sxx the average FFT of the signal s. This is a good approximation of the power spectrum, and therefore, of the spectrum of the signal. First of all, I will show the noise spectrum in order to see the whiteness property and the value of the power spectral density. Figure 23. Representation of the Power Spectrum of the unprocessed noise. Looking at the Figure 23 we can take some conclusions. The first one is that, as we were expecting, the noise power spectrum is completely flat, so the unprocessed noise process is white. However, we can see a small peak in the center of the spectrum due the little difference between the Template and the Real Pulsar Signal. That is because the noise have been computed as y−p=y−s=z, so any small difference between the real pulsar signal and the template will be reflect in the noise power spectrum. Anyway, this little interference is negligible in the computation of the variance of the noise. Another conclusion we can observe is the fact that the amplitude of the power spectrum is equal to the variance of the noise. That proves the fact that the power spectrum of the white sampled noise is fsNo, exactly the same value as the variance if the noise is white. Now, the power spectrum of the whole 58
signal is shown in order to see if the signal is wideband as we are expecting. Computing it with the same process as the used for the noise power spectrum, we have: Figure 24. Representation of the Power Spectrum of the Radio Pulsar noisy signal. Looking carefully at the figure 24, we can see as the bandwidth of the pulsar signal is 200 Hz and with an amplitude higher than the noise spectrum . So, the signal received from the WSRT looks very narrowband. That result is the opposite that the one we were expecting due the fact that the radio pulsar signals has a wideband nature. Moreover, looking at those kind of Bandwidth (the order of MHz) we should see an almost completely flat spectrum. From this figure we cannot still conclude that the radio pulsar signal recorded by the WSRT has a narrowband nature because we are looking at the power spectrum of the signal and not the real spectrum. However, I will prove in the Chapter 7 the fact that the pulsar signal from the WSRT has an almost narrowband nature. 59
Chapter 7 Signal Processing experiments In this section, experiments proving the performance of the algorithms explained in the Chapter 5 with Wideband Signals submerged in Additive White Gaussian Noise will be shown. To do that, I am going to divide it in two subsection. First, I will perform simulations with fake data created by Matlab. In that part I am going to prove that the GENR of the Wideband Signal increases when you perform Epoch Folding and when you improve the bandwidth of the receiver. As stated in Detection Theory, increasing the GENR will mean increasing the detection performance and the estimation of the Time Of Arrival of the Wideband Signal. Furthermore, I will show that Low-Pass filtering doesn’t change the detection performance and that it decreases with Downsampling. After that, I will perform the same algorithms to a Real Radio Pulsar Signal PSR B0329+54 provided by the Westerbork Synthesis Telescope. We will see as the results with that real data will not be the ones we are expecting. I am going to provide an explanation of those results. 7.1 Simulated data In that subsection we are going explained how I have created a Wideband Fake Radio Pulsar Signal submerged in white noise to apply the 5 algorithms explained in that report. After that, we are going to assess the ENR and GENR to the signal after every processing step to prove the theory stated before. The data has been implemented with Matlab. Furthermore, the other signal processing techniques also has been performed with that Software. First of all, it has been created a simulated data with features similar to a Radio Pulsar signal swith length of 2000 samples and a bandwidth of 10 kHz. The simulated data that reach the properties of the star pulses (Wideband, Time-Limited and accurate periodicity) is the Additive White Gaussian Noise with a duration of few samples. In this case, a AWGN signal of length 100 samples over a 2000 samples background has been implemented with an energy of s= 1.01 ∗106. The fact that this fake pulsar signal is an additive white gaussian noise makes the spectrum Wideband. Moreover, the spectrum is almost flat, important property to show the linear improvement of the GENR of the signal after applying the signal processing techniques. In the next figure we can observe the fake rotating star sampled pulse in the time domain. 60
Figure 25. Representation of the simulated radio pulsar signal in time domain with a total length of 2000 samples and a sampling frequency of 20 kHz. As I said, this signal has a wideband nature. Due to the fact that with Matlab is difficult to show the real spectrum of the signal, I am going to compute the power spectrum of the simulated data. So, we will have 1 L|Sxx(f)|2, being Sxx the average FFT of the signal s. In the next figure we can see the Power Spectrum of the simulated data. The red line shows the averaged Power Spectrum while the blue one shows the same signal but without performing the average filter. The total bandwidth of this signal is 10 kHz in order to have the Nyquist sampling frequency of 20 kHz. Figure 26. Representation of the Power Spectrum of the simulated radio pulsar signal with a total bandwidth of 10 kHz. 61
Looking at the Figure 32 we can state that the whitening matrix is recovering the high frequencies of the signal s, eliminated by the filter. Also we can observe that the amplitude has decreased in σ2 z. That is due the whitening process have a normalization by a factor σ2 z, the variance of the noise before applying the filter. 7.1.4 Downsampling In that section we are going to assess the GENR when you downsample a signal. To do it easier, I am going to perform the downsampling to the simulated noisy signal of L= 2000 with the same factors Mas the ones in the Low-Pass Filtering experiment. In the chapter 5.3 I have concluded that downsampling decreases the ENR and the GENR by a factor of M. Hence, the detection performance is decreased. So, if we take into account that yd=r+W, being ydthe downsampled noisy signal, r the fake downsampled pulsar signal and Wthe downsampled noise process, the ENR of the downsampling signal will be ENRd=r/σ2 W=s M2/σ2 z M=ENR/M. Moreover, as the noise process zis Additive White Gaussian Noise, the downsampling will not change the white feature on the process W. Therefore, the value of the GENR will be the same as the value of ENR. Figure 33. Representation of the Generalised Energy To Noise Ratio for a different values of the downsampling factor M in the linear domain As we were expecting, the detection performance decrease when M increases. Is not decreasing in a linear way because the spectrum of the simulated pulsar signal is not completely flat. Moreover, as we are looking at signals with a length up to 200 samples, the accuracy of the ratios in not perfect. However, if you look carefully, the value of the GENR is around 10 when M= 10 and around 120 when M= 1, so the decrement is almost lineal. 68
7.1.5 Increase the Bandwidth of the signal In the last section it has been concluded that decreasing the Bandwidth of the receiver/antenna decreases the ENR, GENR and the detection performance. So, we can demonstrate that increasing the Bandwidth of the antenna improves the detection performance knowing that the opposite technique, downsampling, decrease the GENR. Anyway, in this section I will show the results of how increasing the Bandwidth of the signal, and therefore, the samples per period, improve the ratios used in this Thesis. Before showing the result of the experiments done with simulated data, first I will sum up the theory of how increasing the Bandwidth of the receiver/antenna improves the GENR. If we take a look at the next figure, we can see the schemes of both cases. The one with a Bandwidth Band the one with a Bandwidth MB. The sampling frequency will always be the Nyquist one. Hence, two times the Bandwidth of the signal. Figure 34. Representation of the receiver scheme of the signal with a Bandwidth BHz (left) and the one with the signal with a Bandwidth MB Hz (right). As we can observe, the analogue Low-Pass filter will not change anything as its cut-off frequency is the same as the Bandwidth of the antenna. However, it is shown in order to eliminate the spurious. Looking at the left figure, we can see as is the same scheme as the unprocessed signal shown in the chapter 4. Therefore, we know that the ENR and the GENR of that signal will be: GENR =ENR =sa No Now, looking at the second Scheme, we can see as the energy of the signal saMwill be M times higher than the energy of sadue it has M times more Bandwidth. In addition, as the sampling frequency is also M times bigger in the second scheme, the energy also will increase in a factor M. So, finally we can state that: sM=fsosaM=MfssaM=M2fssa The noise process zMthough, only will increase its variance in a factor of M if we compare to the noise process z. That is due the fact that the A/D converter doesn’t increase the variance of the noise no matter the sampling frequency used. Then, it will be increased by Mdue the bigger Bandwidth of zM. So, the variance of the noise process zMwill be: 69
σ2 zM=Mσ2 z=MfsNo Finally, and taking into account that the process zMis white with autocorrelation function RzM(k)=2MBNoδ(k), the GENR will be: GENR =ENR =sM σz2 M =M2fssa MfsNo=Msa No So, it can be seen as the GENR is increased by a factor of M. Now, to prove that I will show the results of the experiments done. First of all we will consider the signal without the Bandwidth increased. That will be a simulated signal with Bandiwdth 1 kHz, a sampling frequency of 2 kHz and length 200 samples. Then, it has been created the signals with the same properties but with a increased Bandwidth. The new signals will have a Bandwidth of B=MkHz with a sampling frequency of 2MkHz and length 200Msamples. To do that experiment, it has been chosen different values of M between 1 and 10, being the signal with M= 10 the one I have shown in the section 7.1.1. So, after computing the GENR of the different signals, we can see the evolution of it with the increment of M. Figure 35. Representation of the Generalised Energy To Noise Ratio for a different values of M. Being M the factor of increment of Bandwidth Looking at the figure 35 we can see how increasing the Bandwidth of the receiver improves in a linear way the GENR. Hence, the detection performance will also be increased. The improvement is not completely linear due the fact I have done the simulations with signals of few samples, so the accuracy of the process is not perfect. Anyway, it can be easily seen as the ratio is improved. 70
7.2 Radio Pulsar signal B0329+54 from WSRT In that section I will show the results of the experiments with the Real Radio Pulsar Signal data from the Westerbork Synthesis Radio Telescope. The problem is that this pulsar looks behave like a narrowband signal, defying its wideband nature. However, only observing the power Spectrum of the real data we cannot conclude that the pulsar signal data is narrowband due the fact that the Power spectrum we compute is not the same as the real spectrum of the signal. So, the experiments done in this subsection will be shown in order to prove the ”narrowband” nature of this recorded pulsar signal. The experiments done to the real data has been: 1) Epoch Folding, in order to see how this technique increase the SNR and GENR no matter the spectrum of the signal. Also to being able to see the blurred pulsar profile and how it dissapear after performing a High-Pass Filter. 2) Downsample the signal to see how the SNR of the signal is improved and the GENR is not changed. In fact, a downsampling by a factor of 27885 has been performed in order to see how the noise disappear and only the signal remains. Finally, a high pass filter will prove the almost narrowband nature of the pulsar signal taken from the WSRT. 3) Integration in Time, in order to see the Pulsar Profile for every frequency channel. We will see as the shape of the pulsar can be seen for a high frequencies although really attenuated. So, the behaviour of the pulsar is closer to a Narrowband signal. 7.2.1 Epoch Folding + High Pass Filter First of all, I performed the Epoch Folding with different values of K in order to assess the GENR. The periods have been added every 28582316 samples. This number of samples is used in order to perform the Epoch Folding with the new period computed in the section 6.2.2. In the next figure we can observe the evolution of the GENR of the signal depending the number of folds. The assessment of the ratios has been done every 14 folds, finishing with K=195. The results we are expecting is an improvement of 10log(195) = 22.90dB in both ratios. 71
Figure 36. Representation of the the GENR of the pulsar signals after performing Epoch Folding with K folds. The ENR has not been computed because as we have seen, the noise process zis completely white, so the GENR will be equal as the ENR. From that figure we can conclude that Epoch Folding increase the ratios in a linear way. If we check at the GENR after 195 folds we can observe as: GENRk=195 = 62.95dB =GENR + 23.05dB So we can prove as the increment is almost linear as we were expecting. In the next figure is shown the folded pulsar signal after 195 folds in the time and frequency domain. We can see how the pulsar is visible although very blurred by the noise. Furthermore, in the frequency domain we can still see how the pulsar is narrowband and its power spectrum is much higher than the noise one. 72
Figure 37. Representation of the Folded Radio pulsar signal after 195 folds in time domain (left) and in the frequency domain (Power Spectrum)(right) We can observe as the signal is still narrowband with a Bandwidth of 200 Hz after the folding. From now on, to make easy and intuitive the experiments we are going to work with the Folded signal. After showing the effects of Epoch Folding , a high-pass filter with a cut-off frequency of 80 kHz has been performed in order to show how the pulsar signal of the Figure 37 disappear. To start with the experiment, first it has been created a High-Pass Filter h(n) with the Software Matlab. The filter have has an amplitude of 1 in the frequencies between 80 kHz and 20 MHz and almost 0 from 0 Hz to 800 kHz. The next figure shows the transfer function of this filter. Figure 38. Representation of the frequency response of the High-Pass Filter created with Matlab. 73
The following step was filtering the folded signal y195 with that filter. So, in order to see if the high pass-filter has had any effect on the pulsar signal, the representation of the filtered signal is shown. Figure 39. Representation of the folded radio pulsar signal after high-pass filter. From that plot we can state that the pulsar signal visible in the Figure 39 is almost eliminated. We also can see that the noise level has decreased. That can be explained if we take into account that the filter h(n) is eliminating part of the noise spectra, so the variance of the noise (power/sample) is also reduced. Another simulation is performed to those signals. As we have stated, the template pis the same process as the Folded Radio pulsar signal s195. Also, we know that if we correlate that Template with the whole noisy signal the result should be the autocorrelation function of the radio pulsar signal. Hence, a really narrow peak in the middle of the radio pulsar B0329+54 period. Now, if we perform this correlation between the High-Pass Filtered signal and the Template and we do not see anything, that will mean that the radio pulsar signal has been eliminated with the filtering. In the next figure we can observe that fact. 74
Figure 40. Representation of the correlation between the folded signal and the Template (left) and correlation between the folded signal after high pass filter it and the Template(right) The correlation between the high-pass filtered signal and the template is almost 0. That means that the signal does not have nothing in common with the template. 7.2.2 Downsampling + High-Pass Filter In this section it is showed the performance of Downsampling to the folded signal y195 by a factor of 27885. That is done in order to obtain a resultant signal with 1025 samples. Also it has been done to compute the extreme case when the signal is low-pass filtered with a cut-off frequency of 717 Hz. In the next Figure we can see the Downsampled pulsar signal in the time domain. Figure 41. Representation of the folded radio pulsar signal after downsampled it by 27885. 75
The Figure shows the pulsar power profile with the same level as the template before Downsampling. That means that the Low-Pass filter is not cutting any spectra of the pulsar signal. With that result we can take two conclusions. First, the low-pass filter is not decreasing the energy of the pulsar signal. The second is that the estimation of the amplitude of the template is really accurate. Now, if we compute the SNR of that downsampled signal we can see as SNR = 24.60dB. That means that the SNR has increased in 36.2 dB, where it should be 0 dB if the signal had a wideband nature. Finally, the a high-pass filter has been applied in order to see how the Radio Pulsar Signal is almost eliminated. First, we have to take intro account that the downsampled signal will have a Bandwidth of 717 Hz and a sampling frequency of 1434 Hz. So, the cut-off frequency of the high pass filter is fc= 250Hz. In the next figure we can see the downsampled process after high pass filter it. As it can be observed, the signal has been almost eliminated by the filter, another prove that the data from WSRT has been stored eliminating its wideband feature. Figure 42. Representation of the downsampled radio pulsar signal after high-pass filter it. 7.2.3 Integration in Time The final experiment has been the Integration in Time separating the signal in 33 frequency channels. Doing that we will be able to see the Integrated Pulsar Profile of the signal for a different frequencies. The channel 1 shows the lowest frequencies and the channel 33 shows the highest ones. In the next figure we can see the Power Pulsar Profile without noise in the first channel. That is due the fact that Integration on time highly increases the SNR of the signal. 76
Figure 43. Representation of the Power Pulsar Profile in the first Frequency Channel (of 33). The most important thing of that figure in order to check if the received signal is narrowband is the amplitude. If the data behave as a real Pulsar, the other frequency channels should show the same Pulsar Power Profile with almost the same amplitude. This is due the fact that Radio Pulsar Signals are really Wideband, so the intensity of the pulsar should not vary in a observation Bandwidth of 20 MHz. In the next figure we can see the Pulsar Profile for the 4th frequency channel. Figure 44. Representation of the Power Pulsar Profile in the fourth Frequency Channel (of 33). We can see as the Pulsar Profile can still be observed but with a lot of noise. Furthermore, 77
of all-pass filtering. Then, as de-dispersion can be modelled as an all-pass filtering process because it only changes the phase of the received signal, it can be proved that we can avoid the de-dispersion without decreasing the detection or TOA performance. After that, a review of the different signal processing algorithms has been shown in order to assess if the ENR/GENR of the signals are improved with their application. Some interesting conclusions can be extracted. The first one and more important is that the detection performance will not be affected by the application of any kind of analogue or digital filtering. Therefore and assuming the Radio Pulsar Signal as Wideband, nor low-pass filtering neither the de-dispersion method will improve the detection or TOA performance. A second conclusion is that to improve the GENR/ENR it is necessary to increase the Bandwidth of the receiver. Furthermore, it has been found that the limitation of Bandwidth in order to increase the ENR/GENR is the Antenna since the filtering does not affect the detection performance. That can be understood due the fact that it is not possible to reconstruct any signal that the antenna has not received. Then, in order to improve as much as possible the GENR and therefore, the detection and TOA performance, the aim will be to increase as much as possible the Bandwidth of the Antenna. Other signal processing techniques as Downsampling, Oversampling and Undersampling has been rejected. Nevertheless, as known from other Thesis, the Epoch Folding improves the GENR and ENR linearly with the number of folds. Some experiments with simulated signals have been performed in order to prove the Signal Processing theoretical techniques stated in the Thesis. Another Signal Processing technique has been introduced in order to assess if its application improves the detection performance. That method has been the Integration in Time formed by an average and a Downsampling process. That algorithm increases the SNR of the signal but we concluded that it does not improve the detection performance. However, it can be useful due the fact that it reduces the length of the signal. That may decreases the computational complexity of the whole receptor. In this Thesis is not proved the efficiency of that technique, so it will not be included in the first proposed detector but in the second. Later on, a description of the WSRT observatory is presented, along with a PSR B0329+54 characterization. After the performance of some experiments with that data I have concluded that the signal has been recorded in a way that the Wideband Nature of the Radio Pulsar Signal has disappeared. After looking at the way they usually record this signals, it may be that the loss of that high frequencies has been due the fact that they record the signals with 14 antennas. So, they compute different Beamforming techniques to convert the 14 received signals into 1. Therefore, we can’t use that data for prove the theoretical things explained in that Thesis since the main assumption has been the Wideband nature of the Radio Pulsar signals. Finally two receptors are proposed to detect the radio pulsar signals. That analogue part of that detector is formed by an antenna with the Bandwidth as large as possible, a Low-Pass Filter with a cut-off frequency equal as the Bandwidth of the antenna and a A/D converter with an sampling frequency twice the cut-off frequency of the filters (antenna). Then, in the digital domain two different approaches are presented depending on the efficiency of the Integration in Time technique. The digital blocks of the first receptor will be formed by an Epoch Folding Block and a GLRT detector. The implementation of that detector will be N different correlators in parallel with different shifted templates. The digital part of the second receptor will be exactly the same as the first one but including the Integrator between the A/D Converter and the Epoch Folding in order to reduce as much as possible the computational complexity of Epoch Folding and the detector. In the Figures 59 and 60 84
we can observe the proposed detectors. As expected, the de-dispersion block is not included since it will not improve nor the TOA neither the detection performance. 9.2 Future Work Since the detection and signal processing theory to detect and make feasible the real time navigation with pulsar signal has been presented in that Thesis, some practical experiments are needed in order to prove them with real Radio Pulsar Signals. Therefore, the future work to do for the next researchers will be: 1) Prove with real data that Filtering and de-dispersion does not have any effect in the detection performance using the GLRT detector. 2) Prove with real data that increasing the Bandwidth of the receiver leads to a better detection and TOA performance due the fact that the GENR/ENR is increased. 3) Perform some experiments to assess the performance of the Integration in Time technique in order to see if it decreases the time needed to receive, process and detect the radio pulsar signal. 4) Creation of a dispersed pulsar profiles database: astronomers have always corrected dispersion in pulsar recordings, so all the pulsar databases show templates for the de-dispersed versions of them. [5] 5) Try to record and process (with the proposed receptors) a radio pulsar signal with an 2-3 m diameter dish-shape antenna in order to assess the actual processing time required to correctly detect Radio Pulsar Signals. 85
Bibliography [1] Dr. R. Heusdens, Statistical Decision Theory, technical Report, Delft University of Technology, 2014. [2] Jeongmin Lee, Effects of Low-Pass Filtering on inversion of airbone gravity gradient data Master Thesis, Colorado School of Mines, 2006. [3] D. Lorimer and M. Kramer, Handbook of Pulsar Astronomy. Cambridge University Press, Cambridge, 2005. [4] P. P. Vaidyanathan, Multirate Systems and Filter Banks. Prentice Hall Signal Processing Series. Prentice Hall, Inc., Englewood Cliffs, New Jersey, 1993. [5] European Pulsar Network (EPN) database. http://www.mpifrbonn.mpg.de/div/pulsar/data/browser.html, Retrieved October 13, 2005. [6] A. Lyne and B. Ricket, Measurements of pulse shape and spectra of pulsating radio sources, Nature, no. 218, 1968. [7] P. J. Buist, S. Engelen, A. Noroozi, P. Sundaramoorthy, A. A. Verhagen, C. Verhoeven, Overview of pulsar navigation: Past, present and future trends, avigation, vol. 58, no. 2, pp. 153164, 2011. [8] P. A. G. Scheuer, Amplitude variations in pulsed radio sources, Nature, no. 218, 1968. [9] O. Lhmer, M. Kramer, D. Mitra ,D. R. Lorimer , A. G. Lyne, Anomalous scattering of highly dispersed pulsars, The Astrophysical Journal, pp. 157161, 2001. [10] A. Hewish, S. Bell, J. Pilkington, P. Scott, R. Collins, Observation of a rapidly pulsating radio source, Naturel, pp. 709713, 1968. [11] P. A. G. Scheuer, Amplitude variations in pulsed radio sources, Nature, no. 218, 1968. [12] J. Sala, A. Urruela, X. Villares, R. Estalella, J.M. Paredes, HFeasibility study for a spacecraft navigation system relying on pulsar timing information, Tech. Rep. Ariadna 86
Study 03/4202, Universitat Politecnica de Catalunya and Universitat de Barcelona, June 2004. [13] V. K. Chaudhri, Fundamentals, Specifications, Architecture and Hardware Towards a Navigation System Based on Radio Pulsars, Master thesis, Delft University of Technology, 2011. [14] V. K. Chaudhri, A Framework for Designing and Testing the Digital Signal Processing unit of a Pulsar Based Navigation System, Master thesis, Delft University of Technology, 2012. [15] ASTRON, Guide to observations with the westerbork synthesis radio telescope. http: //www.astron.nl/radio-observatory/astronomers/wsrt-guide-observations/ wsrt-guideobservations, 2010. [16] R. Heusdens, S. Engelen, P.J. Buist, A. Noroozi, P. Sundaramoorthy, C. Verhoeven, M. Bentum, E. Gill, Match filtering approach for signal acquisition in radio-pulsar navigation, 63rd International Astronautical Congress, Naples, Italy, 2012. [17] Yuri M. Mange, Detection of radio frequency pulsar signals using a matched filtering approach, Master thesis, Delft University of Technology, 2013. [18] S. Engelen, Deep space navigation system using radio pulsars, FRONT-END Master thesis, Delft University of Technology, 2009. [19] A. A. Kestil, Deep space navigation system using radio pulsars, BACK-END Master thesis, Delft University of Technology, 2009. [20] Harald Martens, Martin Hoy, Barry M. Wise, Rasmus Bro and Per B. Brockhoff Pre-whitening of data by covariance-weighed pre-processing Journal of Chemometrics, Published online in Wiley InteScience (www.interscience.wiley.com), 2002. [21] Marianna V. Ivashina, Oleg Iupikov, Rob Maaskant, Wim A. van Cappellen, and Tom Oosterloo An Optimal Beamforming Strategy for Wide-Field Surveys With Phased-Array-Fed Reflector Antennas IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 59, NO. 6, JUNE 2011 [22] Alan V. Oppenheim and George C. Verghese Detection Theory MIT OpenCourseWare. Introduction to Communication, Control, and Signal Processing Spring 2010 87