<jupyter_start><jupyter_text>Using Fourier Analysis to Analyze Quasi-Periodic OscillationsBy <NAME>  Problem 1: damped harmonic oscillator example Generating a light curve<jupyter_code>a = Table()
a.meta['dt'] = 0.01  # time step, in seconds
a.meta['duration'] = 200 # length of time, in seconds
a.meta['omega'] = 2*np.pi  # angular frequency, in radians
a.meta['phi'] = 0.0  # offset angle, in radians<jupyter_output><empty_output><jupyter_text>1a. Compute the time steps and a cosine harmonic with the above-defined properties.<jupyter_code>time_steps = np.arange(0, a.meta['duration'] + a.meta['dt'], a.meta['dt'])
amplitudes = np.cos(a.meta['omega'] * time_steps + a.meta['phi'])<jupyter_output><empty_output><jupyter_text>1b. Compute four exponentially damped versions of the harmonic oscillation. $$D(t_i) = e^{-\zeta t_i}H(t_i)$$Pick your own four $\zeta$s! I recommend values between 0.01 and 1.<jupyter_code>def dampen(time_steps, amplitudes, zeta):
    damp_factor = np.exp(-1.0 * zeta * time_steps)
    return amplitudes * damp_factor

zeta_dict = {}
for zeta in [0.01, 0.1, 0.4, 1.0]:
    zeta_dict[zeta] = dampen(time_steps, amplitudes, zeta)<jupyter_output><empty_output><jupyter_text>1c. Plot them all on top of each other.<jupyter_code>plt.figure()
for zeta, values in zeta_dict.items():
    plt.plot(time_steps, values, label=r'$\zeta =$ %.2f' % zeta)
plt.legend()
plt.xlabel('Time (sec)')
plt.ylabel('Amplitude')
plt.show()<jupyter_output><empty_output><jupyter_text>1d. Take the power spectrum of the harmonic and 4 damped harmonic time series.The power $P$ at each frequency $\nu_i$, for the Fourier transform $F$, is $$P(\nu_i)=|F(\nu_i)|^2$$<jupyter_code>power_dict = {}
for zeta, amp in zeta_dict.items():
    power_dict[zeta] = np.abs(fftpack.fft(amp))**2<jupyter_output><empty_output><jupyter_text>Computing the Fourier frequencies and the index of the Nyquist frequency.<jupyter_code>freq = fftpack.fftfreq(len(time_steps), d=a.meta['dt'])
nyq_ind = int(len(time_steps)/2.)  # the index of the last positive Fourier frequency<jupyter_output><empty_output><jupyter_text>1e. Plot them! Notice the trend between the width of the peak in the power spectrum, and the strength of the damping factor.For bonus points, put in a key/legend with the corresponding $\zeta$ value for each curve.<jupyter_code>plt.figure()
for zeta, values in power_dict.items():
    plt.plot(freq, values / np.max(values), label=r'$\zeta =$ %.2f' % zeta)
plt.legend()
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
plt.xlim(0.9, 1.1)
plt.show()<jupyter_output><empty_output><jupyter_text>Problem 2: Analyzing *NICER* data of the black hole X-ray binary MAXI J1535-571  Import it with astropy tables from the file "J1535_evt.fits".<jupyter_code>! ls DSFP_Day4_data
from astropy.io import fits
filename = '/Users/rmorgan/Downloads/DSFP_Day4_data/J1535_evt.fits'
hdu = fits.open(filename)
data = hdu[1].data
hdu.close()
data[0:5]<jupyter_output><empty_output><jupyter_text>The data have come to us as an 'event list', meaning that it's a list of the time at which a photon was detected (in seconds, in spacecraft clock time)  and what the energy of the photon was (a channel number; channel/100=photon energy in keV).   2a. Turn this rag-tag list of photons into an evenly-spaced light curve 2a.i. First, clean it up a little by only keeping photons with energies greater than 1 keV and less than 12 keV, using a mask.<jupyter_code>#Note that the units in the table are in 10*ev
tab = Table(data)
tab = tab[(tab['ENERGY'] > 100) & (tab['ENERGY'] < 1200)]<jupyter_output><empty_output><jupyter_text>2a.ii. Then, make an evenly-spaced light curve array with `np.histogram`. Pick a light curve time resolution of `dt=1/8`seconds to start with. To put your light curve in units of count rate, divide the histogram (counts per bin) by `dt` (seconds per bin); to avoid typecasting errors, do this by multiplying by int(1/dt).<jupyter_code>bins = np.arange(np.min(tab['TIME']), np.max(tab['TIME']) + 1./8, 1./8) 
values, bins = np.histogram(tab['TIME'], bins=bins)
counts = values * 8<jupyter_output><empty_output><jupyter_text>(yes, it takes a second or two; you're using half a million time bins in your light curve!)<jupyter_code>plt.figure()
plt.step(bins[0:-1], counts)
plt.yscale('log')
plt.show()<jupyter_output><empty_output><jupyter_text>2b. Let's try taking the power spectrum of it.$$P(\nu_i)=|F(\nu_i)|^2$$where $F$ is the FFT of your light curve.<jupyter_code>power_spectrum = np.abs(fftpack.fft(counts))**2
freq = fftpack.fftfreq(len(power_spectrum), d=1./8)
nyq_ind = int(len(power_spectrum)/2.)  # the index of the last positive Fourier frequency<jupyter_output><empty_output><jupyter_text>Plot it!<jupyter_code>plt.figure()
plt.plot(freq, power_spectrum)
#plt.xlim(-0.005, 0.005)
plt.yscale('log')
plt.show()<jupyter_output><empty_output><jupyter_text>It's ugly! But more importantly, you can't get useful science out of it.  What's going on?1. There are gaps in the light curve due to the orbit of the spacecraft (and occasionally stuff gets in the way). This has the effect of inserting top-hat windows into our function, which give the lumpy bumps at ~0.25 Hz. So, we need to break the light curve up into shorter segments that won't have weird drop-outs.2. There is a giant DC component at $\nu=0$. This is not astrophysical in origin, but from the mean of the light curve.3. Power spectra are often plotted on log-log scales, but the power gets really noisy and 'scatter-y' at higher frequencies. 4. The astute observer will notice that we can only go up to a Nyquist frequency of 4 Hz. There's interesting astrophysical signals above 4 Hz, but if we did smaller `dt` with keeping the very long segment length, we'd have >1 million time bins, which can be asking a lot of a laptop processor.   2c. Segments!2c.i. Turn your light curve code from 2a.ii. into a function:<jupyter_code>def make_lc(events, start, end, dt):
    bins = np.arange(start, end + dt, dt) 
    values, bins = np.histogram(events, bins=bins)
    counts = values * 8
    return counts, bins<jupyter_output><empty_output><jupyter_text>2c.ii. Sometimes, the detector is on and recording photons, but it's pointed too close to the Earth, or a structure on the spacecraft is occulting part of the view, or the instrument is moving through a zone of high particle background, or other things. The times when these things happen are recorded, and in data reduction you make a list of Good Time Intervals, or GTIs, which is when you can use good science data. I made a list of GTIs for this data file that are longer than 4 seconds long, which you can read in from "J1535_gti.fits".<jupyter_code>filename = '/Users/rmorgan/Downloads/DSFP_Day4_data/J1535_gti.fits'
hdu = fits.open(filename)
gti_data = hdu[1].data
hdu.close()
gti_data[0:5]<jupyter_output><empty_output><jupyter_text>2c.iii. Not only do we want to only use data in the GTIs, but we want to split the light curve up into multiple equal-length segments, take the power spectrum of each segment, and average them together. By using shorter time segments, we can use finer `dt` on the light curves without having so many bins for the computation that our computer grinds to a halt. There is the added bonus that the noise amplitudes will tend to cancel each other out, and the signal amplitudes will add, and we get better signal-to-noise!As you learned in Jess's lecture yesterday, setting the length of the segment determines the lowest frequency you can probe, but for stellar-mass compact objects where we're usually interested in variability above ~0.1 Hz, this is an acceptable trade-off.<jupyter_code>time = np.asarray(tab['TIME'])  ## Doing this so that we can re-run
seg_length = 32. # seconds
dt = 1./128.# seconds
n_bins = int(seg_length/dt) # Number of time bins in a segment of light curve
psd_avg = np.zeros(n_bins)  # initiating, to keep running sum (then avearge at end)
n_seg = 0
for (start_gti, stop_gti) in zip(gti_data['START'], gti_data['STOP']):
    start_time = start_gti
    end_time = start_time + seg_length
    while end_time <= stop_gti:
        ## Make a mask of events in this segment
        events = time[np.where((time > start_time) & (time < end_time))]
        
        ## Keep the stuff not in this segment for next time

        ## Make the light curve
        counts_per_sec, _ = make_lc(events, start_time, end_time, dt)

        ## Turn that into a power spectrum
        psd = np.abs(fftpack.fft(counts_per_sec))**2

        ## Keep a running sum (to average at end)
        psd_avg += psd

        ## Print out progress
        if n_seg % 5 == 0:
            print(n_seg)
        ## Incrementing for next loop
        n_seg += 1
        start_time += seg_length
        end_time += seg_length
        
## Divide summed powers by n_seg to get the average
psd = psd_avg / n_seg<jupyter_output>0
5
10
15
20
25
30
35
40
45
50
55
60
65
70
75<jupyter_text>Plot it! Use similar code I gave you above to make the array of Fourier frequencies and get the index of the Nyquist frequency.<jupyter_code>freq = fftpack.fftfreq(len(psd), d=dt)
plt.figure()
plt.plot(freq, psd)
plt.xlim(-1, 20)
plt.yscale('log')
plt.show()<jupyter_output><empty_output><jupyter_text>2d. Mean-subtractedSo, you can see something just to the left 10 Hz much clearer, but there's this giant whopping signal at the lowest frequency bin. This is what I've heard called the 'DC component', which arises from the mean count rate of the light curve. To get rid of it, subtract the mean from your light curve segment before taking the Fourier transform. Otherwise, keep the same code as above for 2c.iii. You may want to put some of this in a function for future use in this notebook.  2e. Error on average powerThe average power at a particular frequency has a chi-squared distribution with two degrees of freedom about the true underlying power spectrum. So, error is the power value divided by the square root of the number of segments. A big reason why we love power spectra(/periodograms) is because this is so straight forward! One way to intuitively check if your errors are way-overestimated or way-underestimated is whether the size of the error bar looks commeasurate with the amount of bin-to-bin scatter of power at neighbouring frequencies. Plotting, this time with `ax.errorbar` instead of `ax.plot`. The thing at ~8 Hz is a low-frequency QPO, and the hump at-and-below 1 Hz is broadband noise (which we'll discuss in detail this afternoon)!! Now that you've got the basic analysis step complete, we'll focus on plotting the data in a meaningful way so you can easily extract information about the QPO and noise.  2f. Re-binningWe often plot power spectra on log-log scaled axes (so, log on both the X and Y), and you'll notice that there's a big noisy part above 10 Hz. It is common practice to bin up the power spectrum `geometrically` (which is like making it equal-spaced in when log-plotted). For this example, I'll use a re-binning factor of 1.03. If new bin 1 has the width of one old bin, new bin 2 will be some 1.03 bins wider. New bin 3 will be 1.03 times wider than *that* (the width of new bin 2), etc. For the first couple bins, this will round to one old bin (since you can only have an integer number of bins), but eventually a new bin will be two old bins, then more and more as you move higher in frequency. If the idea isn't quite sticking, try drawing out a representation of old bins and how the new bins get progressively larger by the rebinning factor.For a given new bin `x` that spans indices `a` to `b` in the old bin array: $$\nu_{x} = \frac{1}{b-a}\sum_{i=a}^{b}\nu_{i}$$$$P_{x} = \frac{1}{b-a}\sum_{i=a}^{b}P_{i}$$$$\delta P_{x} = \frac{1}{b-a}\sqrt{\sum_{i=a}^{b}(\delta P_{i})^{2}}$$<jupyter_code>def rebin(freq, power, err_power, rebin_factor=1.05):
    """
    Re-bin the power spectrum in frequency space by some re-binning factor
    (rebin_factor > 1). This is sometimes called 'geometric re-binning' or 
    'logarithmic re-binning', as opposed to linear re-binning 
    (e.g., grouping by 2)

    Parameters
    ----------
    freq : np.array of floats
        1-D array of the Fourier frequencies.

    power : np.array of floats
        1-D array of the power at each Fourier frequency, with any/arbitrary
        normalization.

    err_power : np.array of floats
        1-D array of the error on the power at each Fourier frequency, with the
        same normalization as the power.

    rebin_factor : float
        The factor by which the data are geometrically re-binned.

    Returns
    -------
    rb_freq : np.array of floats
        1-D array of the re-binned Fourier frequencies.

    rb_power : np.array of floats
        1-D array of the power at the re-binned Fourier frequencies, with the
        same normalization as the input power array.

    rb_err : np.array of floats
        1-D array of the error on the power at the re-binned Fourier
        frequencies, with the same normalization as the input error on power.
    """
    assert rebin_factor >= 1.0

    rb_power = np.asarray([])  # Array of re-binned power
    rb_freq = np.asarray([])   # Array of re-binned frequencies
    rb_err = np.asarray([])    # Array of error in re-binned power
    real_index = 1.0           # The unrounded next index in power
    int_index = 1              # The int of real_index, added to current_m every iteration
    current_m = 1              # Current index in power
    prev_m = 0                 # Previous index m

    ## Loop through the length of the array power, new bin by new bin, to
    ## compute the average power and frequency of that new geometric bin.
    while current_m < len(power):



    return rb_freq, rb_power, rb_err<jupyter_output><empty_output><jupyter_text>Apply this to the data (using JUST the frequency, power, and error at positive Fourier frequencies). Start with a rebin factor of 1.03. Play around with a few different values of rebin_factor to see how it changes the plotted power spectrum. 1 should give back exactly what you put in, and 1.1 tends to bin things up quite a lot. Congratulations! You can make great-looking power spectra! Now, go back to part 2d. and try 4 or 5 different combinations of `dt` and `seg_length`. What happens when you pick too big of a `dt` to see the QPO frequency? What if your `seg_length` is really short?One of the most important things to notice is that for a real astrophysical signal, the QPO (and low-frequency noise) are present for a variety of different `dt` and `seg_length` parameters.   2g. NormalizationThe final thing standing between us and a publication-ready power spectrum plot is the normalization of the power along the y-axis. The normalization that's commonly used is fractional rms-squared normalization, sometimes just called the rms normalization. For a power spectrum created from counts/second unit light curves, the equation is:$$P_{frac.rms2} = P \times \frac{2*dt}{N * mean^2}$$`P` is the power we already have,`dt` is the time step of the light curve,`N` is `n_bins`, the number of bins in one segment, and`mean` is the mean count rate (in counts/s) of the light curve (so, you will need to go back to 2d. and re-run keeping a running sum-then-average of the mean count rate).Don't forget to apply the normalization to the error, and re-bin after!  2h. Poisson noise levelNotice that the Poisson noise is a power law with slope 0 at high frequencies. With this normalization, we can predict the power of the Poisson noise level from the mean counts/s rate of the light curve! $$P_{noise} = 2/mean$$Compute this noise level, and plot it with the power spectrum. Your horizontal Poisson noise line should be really close to the power at and above ~10 Hz.  2i.  For plotting purposes, we sometimes subtract the Poisson noise level from the power before plotting. Once we've done this and removed the noise, we can also plot the data in units of Power, instead of Power/Hz, by multiplying the power by the frequency. Recall that following the propagation of errors, you will need to multiply the error by the frequency as well, but not subtract the Poisson noise level there. Beautiful! This lets us see the components clearly above the noise and see their *relative* contributions to the power spectrum (and thus to the light curve).  Recap of what you learned in problem 2: You are now able to take a light curve, break it into appropriate segments using the given Good Time Intervals, compute the average power spectrum (without weird aliasing artefacts or a strong DC component), and plot it in such away that you can see the signals clearly.  Problem 3: It's pulsar timeWe are going to take these skills and now work on two different observations of the same source, the ultra-luminous X-ray pulsar Swift J0243.6+6124. The goal is for you to see how different harmonics in the pulse shape manifest in the power spectrum.  3a. Load the data and GTIUsing the files J0243-122_evt.fits and J0243-134_evt.fits, and the corresponding *x*_gti.fits.  3b. Apply a mask to remove energies below 0.5 keV and above 12 keV.  3c. Make the average power spectrum for each data file.Go through in the same way as 2d (using the `make_lc` function you already wrote). Re-bin and normalize (using the `rebin` function you already wrote). The spin period is 10 seconds, so I don't recommend using a segment length shorter than that (try 64 seconds). Since the period is quite long (for a pulsar), you can use a longer `dt`, like 1/8 seconds. Use the same segment length and dt for both data sets.<jupyter_code>seg_length = 64. # seconds
dt = 1./8.# seconds
n_bins = int(seg_length/dt) # Number of time bins in a segment of light curve<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08455385308697612, 'fraction_numerical': 0.01622352410995944, 'mean_word_length': 4.1102475532527345, 'pattern_counts': {'":': 0, '<': 81, '<?xml version=': 0, '>': 82, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '427257', 'license': 'MIT', 'n_tokens_mistral': 5627, 'n_tokens_neox': 5208, 'n_words': 2430, 'path': 'Session9/Day4/workbook_QPO.ipynb', 'repo_name': 'rmorgan10/LSSTC-DSFP-Sessions', 'size': '17761'}