COMPUTER PROJECT
IN
COHERENT DE-DISPERSION OF PULSAR DATA: PYTHON IMPLEMENTATION
AT
Birla Institute of Technology and Sciences, PilaniGoa Campus
1
Report Submitted by:
Adithya Reddy – 2005A3PS055G
Report Done in the Partial Fulfillment of Computer Project BITS GC 335 at
BITS – Pilani, Goa Campus
Date of Submission: Submission: April 30th, 2009
2
CERTIFICATE
This is to certify that project work on “Coherent De-dipsersion of Pulsar Data: completed successfully successfully working Python Implementation” by Adithya Reddy was completed with Dr. Ashish Asgekar during the 8th Semester(January to May), 2009 at BITS Pilani Goa Campus .
Project Guide: Dr. Ashish Asgekar
Date: 30th April, 2009
3
ACKNOWLEDGEMENTS
I thank Dr. Ashish Asgekar, giving me an opportunity to work on this project. Work on this project would not have been possible without his constant support, guidance and motivation. I also would like to thank other Bitsians: Nishant who helped me out beginning with python and Amit for our discussions about the project.
4
ABSTRACT
Pulsars are highly magnetized, rotating neutron stars that emit a beam of electromagnetic radiation. The radiation can only be observed observed when the beam of emission emission is pointing towards the Earth. Due to the dispersive nature of the interstellar plasma, lower-frequency radio waves travel through the medium slower than higher-frequency radio waves. Before further analysis of pulsar data there is a need to correct or De-disperse the band limited signal. Coherent De-dispersion being one of the efficient methods for De-dispersion, this project aims at a software implementation for the same. A python code using scientific computing tools in this platform is written and tested on pulsar data. This report explains the usage of code and an overview of Software Implementation. It also compares plots of original pulsar data and the data after Coherent De-dispersion. Keywords: Pulsars, Python, Phase-Correction filter, Folding, De-dispersion
5
CONTENTS:
I. Introduction II. Important Expressions and Definitions: a. Dispersion Measure b Dispersion Delay c. Base-band Voltage III. Coherent De-Dispersion a. Basis of Method b. Overlap save algorithm used IV. Block Diagram of System V. Software Resources needed VI. Stages in Software Implementation VII. Usage of code VIII. Plots obtained IX. Current work X. Conclusion XI. Appendix A – Installing FFTW bindings XII. Appendix B – Full Code XIII. References
7 8
9 9 10 11 12 13 14 15 19 20 21 22 29
6
I. Introduction Pulsars are highly magnetized, rotating neutron starts that emit a beam of electromagnetic radiation. The radiation can only be observed when the beam of emission is pointing towards the Earth. This is called the lighthouse effect and gives rise to the pulsed nature that gives pulsars their name. The radiation from pulsars passes through the interstellar medium (ISM) before reaching Earth. Due to the dispersive nature of the interstellar plasma, lower-frequency radio waves travel through the medium slower than higher-frequency radio waves. The dispersion within the bandwidth causes the pulse to be broadened once the signal of entire band is put together without special measures used to de-disperse the signal. Hence, before further analysis of pulsar data there is a need to correct or De-disperse the band limited pulsar signal. Coherent De-dispersion being one of the efficient methods for De-dispersion, this project aims at a software implementation for the same. It involves phase correction of the received data. A python code using scientific computing tools in this platform is written and tested on pulsar data. This report explains the usage of code. FFTW algorithm is used for implementing all Fourier transforms involved. Bindings to call the FFTW (C library) f rom python are created and used in the project.
Illustration 1: Schematic view of a pulsar. The sphere in the middle represents the neutron star; the curves indicate the magnetic field lines and the protruding cones represent the emission beams. Source: Wikipedia
7
II. Important Expressions and Definitions: a. Dispersion Measure:
The delay in the arrival of pulses at a range of frequencies is directly measurable as the dispersion measure (DM) of the pulsar.
Where D is the distance from the pulsar to the observer and ne is the electron density of the ISM
b. Dispersion Delay:
The propagation delay introduced by dispersion between the band at frequency v and at an infinitely high frequency band is called Dispersion delay.
Equation 1
DM=Dispersion Measure value in cm -3pc v = Frequency value in MHz
c. Base-Band Voltage:
The pre-detection baseband voltage voltage signal of a pulsar is given by
Equation 2 Øp = Phase acquired by probation vb = Baseband Frequency
8
III. Coherent De-Dispersion: Basis of Method Coherent De-dispersion was deviced by Hankins. As explained detailedly in Reference [1] Ø P is split and written as sum of three factors. ØP=ØP0 +2πvb ĩ 1 + ØD
Equation 3
ØP0 ØD vb ĩ 1 2πvb ĩ 1
= Constant Phase factor = Differential delay between the frequencies with in the band = Base band frequency = Propagation Delay at Frequency v, at Lower edge of RF band = Linear Phase gradient through the band
The expression for Ø D is given by
D1
= Dispersion Delay at frequency v 1 i.e. at lower edge of RF band
Using Equation 2 and 3 in frequency domain
Equation 4
F (V(t)) =
Multiplying above spectrum by a factor f actor e –iØD, We obtain the corrected spectrum
Equation 5
F (Vc(t)) =
To convert it back to time domain we perform an IFFT on the corrected spectrum. Hence, Coherent De-dispersion mainly involves 3 steps: 1. Converting the block from pulsar signal to frequency domain by performing FFT with overlap from adjacent blocks 2. Multiplying with the factor e –iØD 3. Getting back the corrected signal in time domain by performing IFFT. Taking care of edge effects and hence rejecting the polluted points and storing the middle values in the output block.
9
Overlap Save method:
During overlap-save, the input data (top) are divided into over-lapping segments of N points. Each N-point segment is separately transformed into the Fourier domain, multiplied by the Phase correction filter transformed back into the time domain (middle). Polluted points from the wrap-around region (shaded grey) are discarded before copying the result into the output dataset (bottom). Overlap percent chosen equals 50% of the considered Input block size on either side.
Above Figure shows the Overlap save Algorithm used to account for the edge effects
For more Information on the algorithm in reference r eference with Coherent De-dispersion refer to Willem van Straten [4]
10
IV. Block diagram of the System: a. Bloc Block k diag diagra ram: m: Based on above discussion the Phase correction system to be implemented in Python has following blocks involved:
Pulsar Data
Read data in blocks
Execute FFT with Overlap
Apply Phase Correctio n Filter
Execute IFFT and reject Polluted Points
For detailed description of block diagram refer to Ithaki and Gupta, 1995 [3].
11
Folding to decrease time Resolution
De-dispersed data
V. Software Resources needed: Python:
Python is a general-purpose high level programming language. Its design philosophy emphasizes code readability. Python's core syntax and semantics are minimalistic, while the standard library is large and comprehensive. comprehensive.
Packages that needed to installed before executing the code: NumPy is an extension to the Python programming language, adding support for large, multidimensional arrays and matrices, along with a large library of high-level mathematical functions to operate on these arrays
f or the Python programming language and its NumPy numerical Matplotlib is a plotting library for mathematics extension. It provides a "Pylab" API designed to closely resemble that of MATLAB, thereby making it easy to learn for experienced MATLAB users. FFTW bindings for python:
As implementation of Coherent De-dispersion involves multiple FFT and IFFT calculations, using the inbuilt FFT in python is slower, Python being a high level language. Hence the bindings are generated for the faster C library i.e. FFTW3 used for execution of FFT. The installation procedure for the bindings is in Appendix A Once the installation is complete we can import the binding into Python as >>import fftw3 We need to first f irst create a forward and backward FFT plan depending upon requirement. For FFT >>fft1 = fftw3.Plan(inputa,outpu fftw3.Plan(inputa,outputa, ta, direction='forward', flags=['measure']) For IFFT: >>ifft1 = fftw3.Plan(outputa, inputa, direction='backward', flags=['measure']) Where inputa and outputa are 1-Dimensional Numpy arrays For computing the FFT at a point of time you wish >>fft1.execute() Similarly for computing IFFT >>ifft1.execute()
12
VI. Stages in Software Implementation: Implementation: The steps involved in writing the code for the Coherent De-dispersion Algorithm: 1. Reading the pulsar Header and unpacking the information regarding the data format. 2. Forcing some values in header from Config.das file: 3. Some of the values read into the code: code: a. Bandwidth b. Sampling frequency c. RF channels d. Dispersion Measure 3. Reading the 2-bit data for 2 channels alternatively present (i.e. channel1_2bits, channel2_2bits, channel2_2bits, channel1_2bits channel1_2bits .. ), for each channel 4096 samples are read at once. The 2 bit values are unpacked according to the thresholding look-up table 4. Depending upon the choice 4096 samples of specific channel converted to frequency domain. The FFT is implemented in python with FFTW bindings to be faster. 6. The % of data to be considered for De-dispersion is given and an input. Also prefix block value where the pulsar observation starts is given as input to code. 7. 50 % Overlap and save on either sides is used while computing the spectrum as explained in Section VI. 8. Phase correction is done by multiplying by a factor e-iØd also called as Chirp function in [4] 9. Corrected or De-dispersed spectrum is Inverse Fourier transformed to get back to time domain. 10. The time data is plotted after decreasing the resolution by a factor (say 16: Multiple of 2 and <=512 as limited by the code).
13
VII. Usage of code: Code can be executed by following command >>python filename Or >>python code.py Code requires three three inputs for execution: execution: 1. The percentage of total signal present in the pulsar data file to be considered for De-dispersion. 2. The channel number whose data you would ask to de-disperse (1 for 1 st channel and 2 or any other character for 2 nd channel) 3. Number of blocks (each 2048 bytes long) after which pulsar observation should start. For channel 1 pulsar peaks can be observed for a value of 750. 4. Factor by which time resolution of the output is reduced.
Uncompress the code.tar.gz file to a folder named 'code': It would contain folder 1. “FFTW”: Containing the FFTW bindings to be installed. Refer to Appendix A for installation details. and files: 1. Pulsar data file with markers removed: 'data_fixed' 2. Configuration file: 'config.das' 3. Code file named 'code.py' 4. Upon execution will create a file “disp.out” containing de-dispersed de-dispersed data.
14
VIII. Plots Obtained:
Illustration 2: Power Spectrum of channel 1 at several frequency with 50 % of data taken to compute spectra [Average Power-Y axis vs Channel frequcny X axis]
Average Power
Illustration 3: Power Spectrum Spectrum of channel channel 2 at different frquencies frquencies with 50 % of data taken to compute spectra[ Average Power-Y axis vs Channel frequcny X axis ]
15
Illustration 4: Pulsar Peak observed after decreasing the Time Resolution by 16 and plotting 50% of total data
16
Illustration 5: Denser Picture of Peak observed at 750 th block in channel 1
Illustration 6: Same position in Illustration 5 Observed after accommodation for OVERLAP SAVE
17
Illustration 7: Original Pulsar data plotted at point of peak observation
Illustration 8: Pulsar peak observed with folding/Time resolution factor 64 and Overlap save
18
Illustration 9: After final analysis 32K point FFT has been taken to obtain the beautiful pulsar peak at ms resolution
19
IX. Current Work:
Interference rejection step is being accommodated before de-dispersion. Work in [2] which deals with median analysis is being used for the interference rejection. The channels which have heavy interference i.e. above a particular threshold are not considered for coherent De-dispersion process. The code median_analysis.py median_analysis.py from [2] will be used for the same.
20
IX. Conclusion:
Coherent De-dispersion had been implemented in Python with using Numpy and Matplotlib packages. FFTW is used for fast implementation of Fourier transform. The code caters for 2 channel pulsar data plotting both power spectrum and normal spectrum. Efforts are being made to include Interference rejection with help of median analysis in the code. Overlap save method accounts for the edge effects that arise due to finite point FFT. The De-dispersed and normal time series when folded by different factors have been observed and compared. Plots showing the effect of Overlap save method have also been accommodated accommodated in the report.
21
X. APPENDIX A: Steps to install FFTW3 bindings: 1. Check if FFTW3 package and Numpy are installed on your computer
2. Unzip the code.zip to form a folder “code” 3. Make code as current directory 4. Execute the following commands in termainal $cd FFTW $su root $python setup.py
5. After executing the above commands, go to Python work environment and execute the following command to test the successfull installation >>import fftw3 6. FFTW3 importing to Python would make sure that bindings are successfully uploaded
22
XI APPENDIX 2: Full Code #Code Function: Observation of De-dispersion on pulsar data #Author: Adithya Reddy G #Supervisor: Dr. Ashish Asgekar #Comments: for ready observation of de-disperse pulsar peak first input aver=1 #2nd input chan=1 #3rd input prefix=755 (number of blocks(each 2048 bytes long the pulsar starts coming)) #4th input factor=16 (factor by by which time resolution resolution should be decrease) decrease)
import numpy from numpy import * from pylab import * import os import fftw3 ############################# ####################### ######
Reading File - Header Part
fp=open('0809_mod_fixed','r') fp=open('0809_mod_fixed ','r') #if NULL file does not exist g=fp.read(2048) #Unpacking Header dh0= numpy.frombuffer(g[0:8], dtype=numpy.float64) fh2=numpy.frombuffer(g[8:12], dtype=numpy.float32) fh2=numpy.frombuffer(g[8:12], fh3=numpy.frombuffer(g[12:16], fh3=numpy.frombuffer(g[ 12:16], dtype=numpy.float32) ih4=numpy.frombuffer(g[16:20], dtype=numpy.int) ih4=numpy.frombuffer(g[16:20], ih5=numpy.frombuffer(g[20:24], ih5=numpy.frombuffer(g[ 20:24], dtype=numpy.int) ih6=numpy.frombuffer(g[24:28], ih6=numpy.frombuffer(g[ 24:28], dtype=numpy.int) ih7=numpy.frombuffer(g[28:32], ih7=numpy.frombuffer(g[ 28:32], dtype=numpy.int) ih8=numpy.frombuffer(g[32:36], ih8=numpy.frombuffer(g[ 32:36], dtype=numpy.int) ih9=numpy.frombuffer(g[36:40], ih9=numpy.frombuffer(g[ 36:40], dtype=numpy.int) ih10=numpy.frombuffer(g[40:44], ih10=numpy.frombuffer(g [40:44], dtype=numpy.int) mark_gap=ih4 if mark_gap != 4096: print " The header has wrong entry for marker_interval. May be header is absent\n"
bw = dh0 invert = 1 if bw < 0.0: bw = -bw; invert = -1; f_clk = bw*(double)(2./1000.) f_0 = fh2
#/* f_clk in MHz */
dm = fh3
23
dd = ih5 mm = ih6 yy = ih7 hr = ih8 mt = ih9 sec = ih10 #memcpy(psrname, ph+44, 32); psrname=g[44:76] print "PSR name: ", psrname data_file=g[76:108] print "Data file: ",data_file site_name=g[108:140]; print "Site Name: ",site_name
print " ======================= ===================================\n" ============\n"
############ Opening "DASCONFIG" file for forcing the values which might be wrong in the header print "Forcing with Config.das" hand_cfg = open('config_Nov_23_2002 open('config_Nov_23_2002.das',"r") .das',"r")
f_0=float(hand_cfg.readline()) invert=int(hand_cfg.readline()) f_clk=float(hand_cfg.readline()) div_panel=int(hand_cfg.readline()) site_name=str(hand_cfg.readline()) mark_gap=int(hand_cfg.readline()) dm=float(hand_cfg.readline()) hand_cfg.close()
bw = (1000.0*f_clk)/(2*div_ (1000.0*f_clk)/(2*div_panel+2) panel+2)
print print print print print print print print
"bandwidth of obs.n :",bw," kHz" "Centred at :",f_0," MHz" "header DM :",dm," pc/cc" "marker interval :",mark_gap," words" "date :",dd,mm,yy "Field name :",psrname "Original data file name :",data_file "Telescope site_code :",site_name
print " ======================= ===================================\n" ============\n"
print "Reading Header extension"
24
ih35=numpy.frombuffer(g[140:144], ih35=numpy.frombuffer(g [140:144], dtype=numpy.int)
n_rfch = ih35; if (n_rfch > 0) & (n_rfch < 3): nbits=numpy.frombuffer(g[144:148], dtype=numpy.int) nbits=numpy.frombuffer(g[144:148], f_1 = numpy.frombuffer(g[14 numpy.frombuffer(g[148:152], 8:152], dtype=float32) bw1 = numpy.frombuffer(g[15 numpy.frombuffer(g[152:160], 2:160], dtype=double) #memcpy(observer, ph+160, 32); /* ih[40] */ observer=g[160:192] #printf(" Observer %s \n",observer); #printf(" ===================== ===================================\n"); ==============\n"); hund_sec = numpy.frombuffer(g[192: numpy.frombuffer(g[192:196], 196], dtype=int) # ih48; # /* next usable: usable: ih[49] or fh[49] or dh[25-] or ph+196 */ print "2nd channel details" print "no of RF_channels : No of bits : 2nd banddwidth : observer" print n_rfch, " ",nbits," ", bw, " KHz "," ",observer else: nbits = 2 n_rfch = 1 bw1 = bw f_1 = f_0 hund_sec = 0 print "Observer Unknown (due to old header)\n"
######################################
print "Creating split codes \n"
#Initialising Variables for carrying out FFT #accum=zeros(4*2048) accum = append(numpy.zeros((409 append(numpy.zeros((4096*2,1), 6*2,1), dtype=complex),0) outputa=append(numpy.zeros((4096*2,1), outputa=append(numpy.ze ros((4096*2,1), dtype=complex),0) modspect=append(numpy.zeros((4096*2,1), modspect=append(numpy.z eros((4096*2,1), dtype=complex),0) e1 = numpy.zeros((4096*2,1), dtype=complex) #For holding overlaped FFT out1=accum f1=accum overlap=0.5 #for 50% overlap on either sides of data
# Creating FFT plan to use FFTW3 library and call FFT, IFFT functions fft1 = fftw3.Plan(f1,out1, direction='forward', flags=['measure']) flags=['measure']) #fft2 = fftw3.Plan(f1,out2, direction='forward', direction='forward', flags=['measure']) #fft3 = fftw3.Plan(junk,out3, direction='forward', flags=['measure'])
25
ifft1 = fftw3.Plan(modspect,out fftw3.Plan(modspect,outputa,direction='backward puta,direction='backward', ', flags=['measure'])
#byte_split method borrowed from Asgekar Sir's code 2_compute_spectra
nfft=4096 byte_split=zeros([256,4]) nfft_4 = nfft/4 limit =nfft_4*2 # For reading two channels of data each 2048 bytes (1024 short ints) at once offset = 0 j=0 while j < 256: if nbits > 1: n = 1 o = 0 while o <= 3:
#//
deal with only one byte
k = j; temp = ((k/n) & (0x03)) byte_split[j][o] = (2*temp - 3) n = n * 4; o=o+1 j=j+1 #print byte_split[0:255] #for debugging print "Unpacking data with split codes \n" #byte_split method borrowed from Asgekar Sir's code 2_compute_spectra
print "How much % of data do u wish to fold the spectrum for ? \n" aver=input() print "Press 1 to correct the 1st channel and anything else for the 2nd chanel correction \n" chan=input() print "enter prefixed value with which pulsar observation starts prefix=input()
\n"
print "enter a factor by which output time resolution should be decreased *SHOULD be multiple multiple of 2 : and <=512 <=512 \n" factor=input() nblocks=int(aver*int(os.stat('0809_mod_fixed'). nblocks=int(aver*int(os .stat('0809_mod_fixed').st_size)/(100*2048)) st_size)/(100*2048))
26
bw= bw/1000.0 # To write the value in MHz nu1= f_0-(bw/2) #nu1=f_0 Frequency at lower edge of RF band
#dm=50.877 rblock=2048.0 #The number of bytes I deal with in 1 loop traversal dd=(4150*dm*pow(nu1,-2))
fi=zeros(2*rblock+1) vi=zeros(2*rblock+1) vb=zeros(2*rblock+1) ci= append(numpy.zeros((rbl append(numpy.zeros((rblock*2,1), ock*2,1), dtype=complex),0) e= numpy.zeros((4096*4,1), dtype=complex)
i=1 while i<(2*rblock+1): vb[i]=((bw/(rblock*2*2))*i) fi[i]=((vb[i])/(nu1+(vb[i])))*2*pi*dd*vb[i]*1000 fi[i]=((vb[i])/(nu1+(vb [i])))*2*pi*dd*vb[i]*1000000 000 cos it would however cancel ci[i]=complex(cos(fi[i]),sin(fi[i])) i=i+1
# no minus sign
#Function to look at the output absolute voltage sum by decreasing the time resolution by input "Factor" ac=0 def res(n,b): for df in range(0,4096/factor): b=append(b,(sum(pow(n[df*factor:factor*(df+1)],2)))) return b
n_1=prefix while n_1<(nblocks+prefix): leave=0 fp.seek((n_1*2048)+2048+leave) p=fp.read(4096) a=[0,0,0,0] b=[0,0,0,0] d=[]
e[:] = numpy.zeros((4096*4,1), dtype=complex)[:]
27
z=numpy.frombuffer(p, dtype=numpy.uint16)
m = 0; # l = offset while
// Initialize the Dest. index for d array
l < limit:
#// assumed 4k words
temp = z[l]/256; n = z[l] - temp*256; #// n = 0; // only to test e[m] = byte_split[n][0]; # // hopefully not reversed m=m+1 e[m] = byte_split[n][1]; m=m+1 e[m] = byte_split[n][2]; m=m+1 e[m] = byte_split[n][3]; m=m+1 n = temp; #// n = 255; e[m] = byte_split[n][0]; m=m+1 e[m] = byte_split[n][1]; m=m+1 e[m] = byte_split[n][2]; m=m+1 e[m] = byte_split[n][3]; m=m+1 #print m #print limit l=l+1
// only to test
#print e[1:100] #print d[1:100] cn=0
print e if chan==1: cn=0 while cn<4096: e1[cn]=e[cn*2] cn=cn+1 else: cn=0 while cn<4096: e1[cn]=e[cn*2+1] cn=cn+1
f1[:]=append(e1,0)[:]
28
fft1.execute()
modspect[0]=out1[0] i=1 while i<(2*rblock+1): modspect[i]= out1[i]*ci[i] modspect[4*rblock+1-i]=conj(modspect[i]) i=i+1
ifft1.execute() ac=res(outputa[2048:6144]/8192,ac) n_1=n_1+1
print "Calculating and Plotting Spectrum u choose with 2048 points"
plot(ac) show()
29
XII. References:
[1] Detection of Radio Emission from pulsars, Dipankar Bhattacharya, Raman Research Institute, India
[2] Identification and Removal of Interference from Time-Frequency Data , S. Vaishnavi, Unpublished report, BITS Pilani Goa Campus
[3] De-dispersion System for GMRT, Ithaki, S., Gupta, Y., & Subrahmanya, C. R. ,JRNL. ASTROPHYSICS & ASTRON. V. 16, NO. SUPPL, P. 243, 1995 [4] High-Precision Timing and Polarimetry of PSR J04374715, Willem van Straten, Swinburne University of Technology [5] Python Tutorial: http://docs.python.org/tutorial/ http://docs.python.org/tutorial/
30
