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:
- 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.
- 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.]
- 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
- 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.
- 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.
- 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.
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:
Limit your plot to the frequency range of interest,
using 'r x {lowfreq} {highfreq}
Fit a model to the 'background' or continuum power
spectrum: 'model pow' or 'model cons linr' might be good first
choices.
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.
When you're happy with the fit you can determine the
formal error on the fit using the UNCER command.
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