Timing tools

Foucalc

[In work! You might use it as a reference in the meantime, but some things are not yet covered. Also let me know if things are not clear! - VG]

The basic tool for X-ray timing analysis with the Remeis scripts is foucalc , which calculates power-, cross power-, coherence- and timelag-spectra . It requires equally spaced, gapless lightcurves as well as the length of the lightcurve segment on which the Fourier transformation is to be performed, given in number of bins. Note that because the underlying algorithm is FFT, the segment length should be a power of two to increase the performance.

When only one lightcurve is given, foucalc calculates the power density and related quantities only. When more than one lightcurve is given, the cross power-, coherence- and timelag-spectra are calculated for all pairs of lightcurves. The lightcurves have to be in the format:

lc = struct { time=[t_1, t_2, ...], rate1=[r1_1, r1_2, ...], rate2=..., rate3= };

The time has to be given in seconds (so be cautious when using fits_read_lc, which usually returns the time in MJD).

If you use RXTE/PCA lightcurves, don't forget to set the numinst keyword, which defines the number of instruments (=PCUs) on, since otherwise your noise correction will be wrong. If you are not using RXTE/PCA set the deadtime keyword to whatever is appropriate for the instrument you are using (e.g. 0).

Given a single PSD value before any further calculations, its uncertainty will be of the same order of magnitude as the PSD itself. Averaging therefore crucial. This is build in into foucalc: given a lightcurve of say 1288192 bins and a segment length of 16384, foucalc will calculate a PSDs on every segment and average every value over all 1288192/15385 = 78 segments.

PSDs

Calculating the PSD

Given two PCA lightcurves,

50110-01-45-00_437_AND_good_14off_excl_0-10.lc

and

50110-01-45-00_437_AND_good_14off_excl_11-13.lc

(which have been checked for gaps and found to be gapless ) we want to calculate and bin the PSD, lags and coherence.

%% read in lightcurves
variable lc1 = fits_read_table("50110-01-45-00_438_AND_good_14off_excl_0-10.lc");
variable lc2 = fits_read_table("50110-01-45-00_438_AND_good_14off_excl_11-13.lc");

%%define variables to store average rate of lightcurve & rms
variable rms,avg;

%%calculate timing quantities
variable result = foucalc(struct{time=lc1.time,rate1=lc1.rate,rate2=lc2.rate},16384;
                          numinst=3,RMS=&rms,avgrate=&avg,verbose);

This produces the following output - thank to the verbose keyword:

light curves: 2516 s with 1288192 bins of 0.00195312 s  (0 s gaps)
  segmentation: 78 segments of 16384 bins
=> frequencies: 0.03125-256 Hz.
Power spectra will be calculated in Miyamoto-normalization.
Calculating power spectrum for rate1.  RMS(0.03125-256 Hz) = 27.1%
Calculating power spectrum for rate2.  RMS(0.03125-256 Hz) = 25.7%
Calculating cross power spectrum (=> coherence and timelags) for rate1 and rate2.

This is right now the only place at which foucalc will complain if the lightcurve contains gap. It will calculate the quantities even if the lightcurve contains gaps - only you can't trust your results. Therefore: use verbose, check that you've done everything right, at least while building up your scripts.

Plot the PSD for the first lightcurve with:

xlog;ylog;plot_with_err(result.freq,result.signormpsd1,result.errnormpsd1);

Rebinning the PSD

The uncertainties on the individual PSDs values are still large, especially at high frequencies. To further improve the SNR, the PSD is rebinned to a new frequency grid and the final PSDs is then the average over the individual PSDs at the frequencies corresponding to the frequencies included into the new frequency bin. Mostly, a logarithmic rebinning with df/f = const. is used, a typical value would be df/f = 0.15 .

Quantities calculated with foucalc are rebinned with the rebin_fouquan tool - the quantity and it's error have to be rebinned in two individual steps:

variable reb_psd = rebin_fouquan(result.freq,result.signormpsd1,
                                 result.numavgall;logfreq=0.15);
variable reb_err = rebin_fouquan(result.freq,result.errnormpsd1,
                                 result.numavgall;logfreq=0.15);

reb_err does not yet contain the real errors of the psds, but needs to be weighted with the square root of the number of values the averaging was performed over. The final PSD is therefore:

variable psd1 = struct{freq_lo = reb_psd.freq_lo,freq_hi=reb_psd.freq_hi,
  value = reb_psd.value, error= reb_err.value/sqrt(reb_err.n)};

Once again plot you PSD - either as PSD vs. frequency or PSD x frequency vs. frequency:

%%PSD vs. frequency
xlog;ylog;plot_with_err((psd1.freq_lo+psd1.freq_hi)/2,psd1.value,psd1.error);

%%PSDxfreuency vs. frequency
xlog;ylog;plot_with_err((psd1.freq_lo+psd1.freq_hi)/2, psd1.value*(psd1.freq_lo+psd1.freq_hi)/2,
psd1.error*(psd1.freq_lo+psd1.freq_hi)/2);

Lag-energy spectrum for archival XMM-Newton data

[--- under construction --- written by OK, please ask if something is unclear!]

The following is an example how to compute PSDs, a lag-frequency spectrum, and a lag-energy spectrum on archival XMM-Newton data of GX 339-4 .

We start by extracting the XMM-Newton PN lightcurve in two energy bands (500-900eV and 2000-3000eV)

#!/bin/tcsh
source $SOFTDIR/sas_init.csh # initialize the SAS software
set xmmscripts=${XMMTOOLS} # make sure the proper xmmscripts are loaded
set obsid='0204730201'
set datadir=/userdata/data/koenig/GX339-4/XMM-2004-03-16/odf
set dt=0.1

${xmmscripts}/xmmprepare --datadir=$datadir --prepdir=$obsid --pn --noapplyflaregti --timing --notimelog # prepare the data

${xmmscripts}/xmmextract --prepdir=$obsid --pn --timing --full --energyband="500,900" --noflarescreen --defaultpattern --dt=$dt --clobber
${xmmscripts}/xmmextract --prepdir=$obsid --pn --timing --full --energyband="2000,3000" --noflarescreen --defaultpattern --dt=$dt --clobber

Now we can start the analysis. Let's load the data and correct them for the telemetry drop-outs (encoded in the Good Time Intervals):

%% Load the lightcurve and GTI file
a=fits_read_table("full_sd_0.1s_0500-0900.lc");
b=fits_read_table("full_sd_0.1s_2000-3000.lc");
gtia=fits_read_table("full_sd_0.1s_0500-0900.lc[GTI0006]");
gtib=fits_read_table("full_sd_0.1s_2000-3000.lc[GTI0006]");

%% Get the start time such that LC starts at 0s
tstart=fits_read_key("full_sd_0.1s_0500-0900.lc[1]", "TSTART");

%% Subtract TSTART and do GTI correction
(filta,filtb)     = filter_gti(a.time, gtia), filter_gti(b.time,gtib);
(timea,timeb)     = a.time[filta]-tstart,b.time[filtb]-tstart;
(countsa,countsb) = a.counts[filta],b.counts[filtb];

If we plot the (re-binned) lightcurve, there is not much to see.

File:Gx399m4 xmm lc.pdf

Fortunately, we can do a Fourier analysis in order to reveal the variability and correlations of these lightcurves. As described above, this can be done with the foucalc routine. I choose the Miyamoto normalization (see ) as in the publication. This normalization yields the PSD in units of fractional root mean square (rms/R_signal)^2.

%% Choose segment length
dimseg=8192; %% (bins), should be a power of 2 for FFT

%% Compute PSDs, CPD, and time lags with foucalc
r=foucalc( struct{time=timea, rate1=countsa/dt, rate2=countsb/dt}, dimseg ; normtype="Miyamoto");
xlog;ylog;plot_with_err(r.freq,r.signormpsd1*r.freq,r.errnormpsd1*r.freq);

%% There is a warning that there are gaps in the LC if verbose is set in foucalc. We ignore this...

File:Gx399m4 xmm psd unbinned.pdf

This really needs some re-binning to increase the S/N-ratio. All foucalc quantities can be re-binned with the rebin_fouquan routine (note that the PSD has already been binned when we set the length of the segments.).

%% Choose Delta f/f, the logarithmic frequency resolution
dff=0.15
reb_psd1 = rebin_fouquan(r.freq, r.signormpsd1, r.numavgall ; logfreq=dff);
reb_err1 = rebin_fouquan(r.freq, r.errnormpsd1, r.numavgall ; logfreq=dff);
reb_psd2 = rebin_fouquan(r.freq, r.signormpsd2, r.numavgall ; logfreq=dff);
reb_err2 = rebin_fouquan(r.freq, r.errnormpsd2, r.numavgall ; logfreq=dff);

%% Need to weight the PSD error with the square root of the number of values the averaging was performed over
hplot_with_err( reb_psd1.freq_lo, reb_psd1.freq_hi, reb_psd1.value*reb_psd1.freq_lo, reb_err1.value/sqrt(reb_err1.n)*reb_psd1.freq_lo );

File:Gx339m4 xmm psd binned.pdf

What can we deduce from this PSD? We observe that the minimum sampled frequency is 1/dimseg - if you choose a long segment length, small quantities can be sampled but then also the S/N-ratio is smaller. On the other hand, choosing a small segment length can lead to an effect called red-noise leakage, which is due to the finite extend of the data (also called windowing, see p.8ff or p.14f.). The maximum sampled frequency is 1/(2*dt), where dt is the timing resolution of the lightcurve. This is called the Nyquist frequency.

foucalc also computes the time lag between the two lightcurves. Let's first rebin and then plot them together with the coherence:

%% Rebin time lag results
variable reb_lag12 = rebin_fouquan(r.freq, r.lag12, r.numavgall ; logfreq=dff);
variable reb_errlag12 = rebin_fouquan(r.freq, r.errlag12, r.numavgall ; logfreq=dff);
hplot_with_err(reb_lag12.freq_lo, reb_lag12.freq_hi, reb_lag12.value, reb_errlag12.value/sqrt(reb_errlag12.n) );

%% Rebin coherence function
variable reb_cof12 = rebin_fouquan(r.freq, r.cof12, r.numavgall ; logfreq=dff);
variable reb_errcof12 = rebin_fouquan(r.freq, r.errcof12, r.numavgall ; logfreq=dff);
hplot_with_err( reb_cof12.freq_lo, reb_cof12.freq_hi, reb_cof12.value, reb_errcof12.value/sqrt(reb_errcof12.n) );

File:Gx339m4 xmm tlag.pdf File:Gx339m4 xmm coherence.pdf

Explain large error bars, explain decay, explain coherence>1, explain how to build lag-energy spectrum