HEASARC Staff Scientist Position - Applications are now being accepted for a Staff Scientist with significant experience and interest in the technical aspects of astrophysics research, to work in the High Energy Astrophysics Science Archive Research Center (HEASARC) at NASA Goddard Space Flight Center (GSFC) in Greenbelt, MD. Refer to the AAS Job register for full details.

Search:

[Advanced Search]



RXTE Helpdesk/FAQ RXTE What's New HEASARC Site Map

RXTE
GOF		 Fourier Analysis of PCA Data and Searches for High-Frequency QPO
Recipes from the RXTE Cook Book 		RXTE
FAQ

Introduction

The analysis of X-ray time series is a large and involved topic with many subtleties. If you are unfamiliar with it, we recommend the following references:

  • Introductory: Bloomfield, P. 1976, Fourier Analysis of Time Series: an Introduction (John Wiley & Sons - New York).

  • General: Jenkins, G.M. and Watts, D.G. 1968, Spectral Analysis and its Applications (Holden Day-Oakland).

  • Searches for periodicities: Leahy, D.A. et al. 1983, Ap. J., 266, 160.

  • General: Press, W.H. et al. 1986, Numerical Recipes, (Cambridge Univ. Press).

  • Counting statistics corrections: Sutherland, P.G. Weisskopf, M.C., Kahn, S.M. 1978, Ap. J., 219, 1029.

  • Fourier techniques: Van der Klis, M. 1989, in Timing Neutron Stars, NATO ASI, Series C, ed. H. Ogelmann (Dordrecht: Reidel), p.27.

To exploit the full resolution of the PCA and to avoid the round-off errors associated with extraction it is better to work with the unbinned data, i.e. with data files that have not been put through saextrct or seextrct to produce a light curve. The xronos ftool powspec can handle these "raw" datafiles, but it can't clean or select the data beforehand. Consequently, this recipe is in two main parts. The first explains how to prepare your data for powspec; the second explains how to run powspec. A final section explains how to look for your QPOs in your power spectrum.

Note that this recipe is "tuned" for the study of QPOs in the range 10-2000 Hz, but can be quickly adapted for other frequency ranges of interest. It assumes that you have data in a Single-Bit mode (e.g. SB_125us_0_249_1s) or Event mode (e.g. E_125us_64M_0_1s, Good Xenon). It will also work for Binned mode, but only those with one channel (e.g. B_500us_1M_0_249_H).

Preparing your Data for Fourier Analysis

For the sake of clarity, the recipe below is for the specific example of four data files in the Event configuration E_125us_64M_0_1s. If you have Single-Bit data or if you have an Event configuration that doesn't have channel information, e.g. E_1us_1M_0_1s, then simply skip the third step and proceed to the powspec section.

  1. Identify continuous sections of clean data: The best way to do this is to make a light curve from your Standard-1 or Standard-2 data with GTI derived from the corresponding filter file.

    1. Using fplot (with offset=yes), examine the light curve and note the relative time ranges encompassing the "good" parts of the light curve.

    2. Since you'll need to know which of your four Event files these time ranges fall in, run timetrans to generate the corresponding absolute times which will appear, in the example below, in the ASCII file good.tint: 
      
   kaah[4]% timetrans

   Running TIMETRANS v3.5.2
   #######################################

   Input file name or @file-of-filenames:[]standard2.lc
   Output filename:[outfile_11847]good.tint
   Input relative time range:[0-250.0, 350.0-500.0, 750.0-1000.0]
   0-3020,4400-6210,7900-8300

  2. Create individual files containing only these good data:

    1. On each of your four Event files run fkeyprint to derive the values of the header keywords TSTART and TSTOP.

    2. Compare these with the absolute time ranges to match the relative time ranges to the individual files.

    3. For each relative time range, run grosstimefilt on the appropriate Event data file to create new Event files containing only the "good" data. We'll call them, say, good_01.se etc.

  3. If desired or applicable, select energy channels: Some Event configurations (but no Single Bit configurations) contain channel information - the configuration E_125us_64M_0_1s, for example, has 64 energy channels. Note that unlike the case when using the extractors, the channel selection on data files must be done with respect to the binning scheme of the data themselves and not with respect to the "absolute" 256-channel scale. Here, we'll be using xspec to determine the channel range, sebitmask to make the bitmask and fselect to create the channel-selected data files:

    1. Assuming that you don't know a priori the energy-to-channel conversion for the configuration, first make a spectrum with seextrct from a sample of the data.

    2. Generate the corresponding response matrix with pcarsp.

    3. Read the spectrum into xspec and examine it. If you want, say, to include only the range 5-20 keV in your Fourier analysis, type: 
      
   XSPEC> ignore 0.0-5.0 20.0-**
   XSPEC> show
      The output of the "show" command will tell you which binned energy channels correspond to the range 5-20 keV. In this case it's: 
      
   Noticed channels    10 to    28

    4. Create the corresponding bitmask with sebitmask. For Event configurations, the channel information is contained in the C-token of the DDL string occupying the TEVTB2 header keyword. In this example, we'll create the bitmask e_5_10_keV.bit, using the channel boundaries derived above, for the E_125us_64M_0_1s file good_01.se. User input is emboldened: 
      
   rufus [74] [day] herx1: sebitmask

   Running SEBITMASK version 3.6
   ==============================================
   Filename containing SE data to be processed:[] good_01.se
   Enter filtering expression (or @filename):[] C[] >= 10 && C[] <= 28
   Enter the column NAME to be processed (Event):[Event]
   Output filename to contain boolean expression
    (acts as input for FSELECT):[bitfile] e_5_20_keV.bit
      The key input here is C[] >= 10 && C[] <= 28 which tells sebitmask to create a bitmask corresponding to values of the C-token between 10 and 28 inclusive. [By the way, you should ignore the rather fearsome warnings!] To view the contents of the newly created bitmask: 
      
   rufus [74] [day] herx1: more e_5_20_keV.bit
   Event >= bxxxx001010xxxxxx &&
   Event <= bxxxx011100xxxxxx &&
   Event == b1xxxxxxxxxxxxxxx
      The last line removes the non-event markers.

    5. Apply the bitmask with fselect to each of the "good" Event files created in step-2: 
      
   rufus [82] [day] herx1: fselect
   Name of FITS file and [ext#][] good_01.se
   Name of output FITS file[] good_01_5_20_keV.se
   Selection Expression[] @e_5_20_keV.bit
      Here, we've used e_5_20_keV.bit to create good_01_5_20_keV.se from good_01.se.
Important warning: the filtered versions of datafiles produced following this recipe cannot be properly processed by saextrct and seextrct (or any other software that uses the GTIs or channel information contained in the header) unless the user inputs the same GTI information used in grosstimefilt into the timeint parameter and the channel filtering information be given using either the chmin and chmax or chint parameters.

Running Powspec

After creating cleaned and selected data files, we're now ready to put them into powspec, the xronos task that generates a power density spectrum. At this point, it's worth repeating that this sort of temporal analysis should not be thought of as a simple recipe. Rather, it's an analysis technique that requires the user to make several choices based on the scientific priorities of the investigation and on the properties of the source.

This example gives the most convenient and common usage of powspec, but note that the tool has a long and comprehensive help file ("fhelp powspec") which the advanced user will wish to study. First type "powspec", then:

  1. Ser. 1 filename +options (or @file of filenames +options)[file1]: Here, give the name of the one of the data files you've prepared, e.g. good_01.se. But note:

    • You can also enter more than one datafile here - either by name or via an ASCII list prefaced with an @-sign). However, there are two good reasons to be cautious. The first is speed: 12 minutes of data at 122-microsecond resolution - i.e. 6 million points - will take powspec an hour to work through on a Sparc10. The second is scientific: calculating a power spectrum over data gaps will introduce spurious signals which may or may not be important.

    • If you are reading Single-Bit or Binned data directly into powspec, then the string VY2 should follow each filename (the same is true if the filenames are in an ASCII list). This tells xronos which column in the FITS file to use, e.g. 
      
   good_singlebit.sa VY2
      In the case of Event files, VY2 isn't needed.

  2. Name of the window file ('-' for default window)[]: Type in a hyphen for the default window (experienced users can supply additional time filtering here - see 'fhelp xronwin' for help on how to create a window file).

    [If you get a message saying that the defaults_win.wi file is not found, this means you don't have the xronos XRDEFAULTS environment variable set correctly. You'll need to set it by issuing a command like: 

    
   setenv XRDEFAULTS '/ftools/SUN/develop/xronos/defaults/'
    With the path correctly set for your ftools installation, run powspec again. Be careful: the last / is important.]

  3. Newbin Time or negative rebinning: To obtain a power spectrum extending up to 2000 Hz - probably a good idea if you're hunting for kHz QPO - your bins should be about 1/4096 seconds, i.e. 250 microseconds (The Nyquist frequency = 1/2*newbin.) Although powspec can accept arbitrary newbin sizes, it's better to specify an integral factor by which the original resolution is to be multiplied. For example, in the case of the Event configuration E_62us_64M_0_1s_L1R1, the time resolution is 1/2**14 seconds. So, to obtain bins of 250 microseconds, the original resolution should be multiplied by a factor of four. To do this in powspec, give the newbin size as -4. Note that powspec prints the time resolution - but often truncated - on the screen after reading the input file. Look for a line like: 
    
   Bin Time (s) ......   0.6104E-04

  4. Number of Newbins/Interval: Powspec divides data into "intervals" over which it calculates individual power spectra. The final power spectrum is the average of these. Choosing the interval size is scientific question: the longer the interval, the larger the frequency range of the power spectrum. But the shorter the interval, the lower the noise. If, for example, you're not interested in frequencies below about 1 Hz, then your intervals needn't be longer than 1 s, that is, 4096 newbins of 250 microseconds.

  5. Number of Intervals/Frame: As you probably want to create the summed power spectrum for the whole dataset, make sure that the number you enter is equal to or bigger than the number given for "default intervals per frame" supplied by the program two lines above this prompt.

    Make a note of the number you enter here, known to the code as NINTFM. It is the number of 'processing intervals' the program will later work through before giving a result.

  6. Rebin results?: You should definitely do this, or your power spectrum will be very noisy at high frequencies. For kHz QPO searches, geometric rebinning of -1.005 is good. If you're more interested in the 10-100 Hz range, use -1.01 or even -1.02. Experimentation will show what level of binning is ideal for your purposes.

  7. After this you'll go through a series of prompts for: output filenames (powspec will add .fps to whatever you specify); whether you want to plot; and the name of the plotting device. (If you say "no" to the "Want to Plot" question, the data will still be saved in the file for your later use.

Powspec will now trundle along, performing NINTFM fast Fourier transforms, writing to the screen the number and start time of the data segment transformed, plus the data maximum, minimum, variance etc.

If you answered "yes" to the plot question, the power spectrum will appear on your screen when the number crunching completes, and will leave you with a PLT> prompt. It only writes the data to a file once you exit the plot by typing "q" - in other words, don't control-Z out of the plot interface, or you'll lose everything.

If you're analyzing LMXB data and searching for high-frequency QPO, you'll look with interest at the region of the power spectrum between (say) 300 and 1200 Hz. You may also see a sloping low frequency red noise component, or QPO in the 1-60 Hz range.

Plotting your Powspec Output to Find QPOs

Once the file is made and you've exited powspec, you can plot the data in the output.fps file by typing "fplot output.fps" and giving the following columns for the axis names: 


   X Axis Parameter[error]: FREQUENCY[XAX_E]
   Y Axis Parameter[error]: POWER[ERROR]
You'll probably want to look at the plot in log-log space: 

   PLT> log x
   PLT> log y
   PLT> p
See the plotting recipe for further details.

If you discover one or more QPO peaks, you'll want to do some fitting to determine some scientifically interesting parameters. The simplest way to do this is to plot the data using FPLOT and the do your fits within the PLT interface. There is extensive help on how to do this in the plotting recipe and also within the PLT interface itself (just type "help" at the PLT> prompt) but, stated briefly, the steps you'll do are as follows:

  1. Limit your plot to the frequency range of interest, using 'r x {lowfreq} {highfreq}

  2. Fit a model to the 'background' or continuum power spectrum: 'model pow' or 'model cons linr' might be good first choices.

  3. Add a Lorentzian to this model, e.g. 'model pow lore'. The three parameters of the Lorentzian are LC, the line center (centroid frequency of the QPO); LW, the line width; and LN, the line normalization.

  4. When you're happy with the fit you can determine the formal error on the fit using the UNCER command.

  5. The centroid frequency and line width can be used in your IAU Circular 'as is'. To get the rms amplitude of the QPO, work out the integral of the Lorenzian: 
    
   I = pi*LN*LW/2
    Finally, divide this by the mean count rate of the source, take the square root, and multiply by 100 to get this as a percentage, i.e. 
    
   rms amplitude = 100*sqrt( I/mean )

Need more info? Need help? In the first instance, feel free to send any questions you have about using powspec or other xronos ftools for RXTE analysis to xtehelp@athena.gsfc.nasa.gov, but note that these tools are provided by (and maintained by) the HEASARC under the guidance of Lorella Angelini, and if it's something we can't answer, we'll refer you to them/her.

If you have a question about RXTE, please send email to one of our help desks.

This page is maintained by the RXTE GOF and was last modified on Wednesday, 24-Aug-2022 11:10:30 EDT.