<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Power spectrum example
This tutorial shows how to make and manipulate a power spectrum of two light curves using Stingray.
Step1:  1. Create a light curve
There are two ways to make Lightcurve objects. We'll show one way here. Check out "Lightcurve/Lightcurve\ tutorial.ipynb" for more examples.
Generate an array of relative timestamps that's 8 seconds long, with dt = 0.03125 s, and make two signals in units of counts. The signal is a sine wave with amplitude = 300 cts/s, frequency = 2 Hz, phase offset = 0 radians, and mean = 1000 cts/s. We then add Poisson noise to the light curve.
Step2:  Now let's turn noisy into a Lightcurve object.
Step3:  Here we plot it to see what it looks like.
Step4:  2. Pass the light curve to the Powerspectrum class to create a Powerspectrum object.
You can also specify the optional attribute norm if you wish to normalize power to squared fractional rms, Leahy, or squared absolute normalization. The default normalization is 'none'.
Step5:  Note that, in principle, the Powerspectrum object could have been initialized directly as
ps = Powerspectrum(lc, norm="leahy")
However, we recommend using this explicit syntax, for clarity. Equivalently, one can initialize a Powerspectrum object
Step6:  Properties
A Powerspectrum object has the following properties 
Step7:  We can plot the power as a function of Fourier frequency. Notice how there's a spike at our signal frequency of 2 Hz!
Step8:  You'll notice that the power spectrum is a bit noisy. This is because we're only using one segment of data. Let's try averaging together power spectra from multiple segments of data.
Averaged power spectrum example
You could use a long Lightcurve and have AveragedPowerspectrum chop it into specified segments, or give a list of Lightcurves where each segment of Lightcurve is the same length. We'll show the first way here.
1. Create a long light curve.
Generate an array of relative timestamps that's 1600 seconds long, and a signal in count units, with the same properties as the previous example. We then add Poisson noise and turn it into a Lightcurve object.
Step9:  2. Pass the light curve to the AveragedPowerspectrum class with a specified segment_size.
If the exposure (length) of the light curve cannot be divided by segment_size with a remainder of zero, the last incomplete segment is thrown out, to avoid signal artefacts. Here we're using 8 second segments.
Step10:  We can check how many segments were averaged together by printing the m attribute.
Step11:  AveragedPowerspectrum has the same properties as Powerspectrum, but with m $>$1.
Let's plot the averaged power spectrum!
Step12:  Now we'll show examples of all the things you can do with a Powerspectrum or AveragedPowerspectrum object using built-in stingray methods.
Normalizating the power spectrum
The three kinds of normalization are
Step13:  Re-binning a power spectrum in frequency
Typically, rebinning is done on an averaged, normalized power spectrum.
1. We can linearly re-bin a power spectrum
(although this is not done much in practice)
Step14:  2. And we can logarithmically/geometrically re-bin a power spectrum
In this re-binning, each bin size is 1+f times larger than the previous bin size, where f is user-specified and normally in the range 0.01-0.1. The default value is f=0.01.
Step15:  Like rebin, rebin_log returns a Powerspectrum or AveragedPowerspectrum object (depending on the input object)
Step16:  Power spectra of normal-distributed light curves
Starting with Stingray 0.3, we can also get Leahy-normalized power spectra of normally-distributed light curves.
Let us calculate such a light curve by subtracting the noise level and normalizing
<ASSISTANT_TASK:>
Python Code:
%load_ext autoreload
%autoreload 2
import numpy as np
from stingray import Lightcurve, Powerspectrum, AveragedPowerspectrum
import matplotlib.pyplot as plt
import matplotlib.font_manager as font_manager
%matplotlib inline
font_prop = font_manager.FontProperties(size=16)
Explanation: Power spectrum example
This tutorial shows how to make and manipulate a power spectrum of two light curves using Stingray.
End of explanation
dt = 0.03125  # seconds
exposure = 8.  # seconds
times = np.arange(0, exposure, dt)  # seconds
signal = 300 * np.sin(2.*np.pi*times/0.5) + 1000  # counts/s
noisy = np.random.poisson(signal*dt)  # counts
Explanation: 1. Create a light curve
There are two ways to make Lightcurve objects. We'll show one way here. Check out "Lightcurve/Lightcurve\ tutorial.ipynb" for more examples.
Generate an array of relative timestamps that's 8 seconds long, with dt = 0.03125 s, and make two signals in units of counts. The signal is a sine wave with amplitude = 300 cts/s, frequency = 2 Hz, phase offset = 0 radians, and mean = 1000 cts/s. We then add Poisson noise to the light curve.
End of explanation
lc = Lightcurve(times, noisy, dt=dt, skip_checks=True)
Explanation: Now let's turn noisy into a Lightcurve object.
End of explanation
fig, ax = plt.subplots(1,1,figsize=(10,6))
ax.plot(lc.time, lc.counts, lw=2, color='blue')
ax.set_xlabel("Time (s)", fontproperties=font_prop)
ax.set_ylabel("Counts (cts)", fontproperties=font_prop)
ax.tick_params(axis='x', labelsize=16)
ax.tick_params(axis='y', labelsize=16)
ax.tick_params(which='major', width=1.5, length=7)
ax.tick_params(which='minor', width=1.5, length=4)
plt.show()
Explanation: Here we plot it to see what it looks like.
End of explanation
ps = Powerspectrum.from_lightcurve(lc, norm="leahy")
print(ps)
Explanation: 2. Pass the light curve to the Powerspectrum class to create a Powerspectrum object.
You can also specify the optional attribute norm if you wish to normalize power to squared fractional rms, Leahy, or squared absolute normalization. The default normalization is 'none'.
End of explanation
print("\nSize of positive Fourier frequencies:", len(ps.freq))
print("Number of data points per segment:", ps.n)
Explanation: Note that, in principle, the Powerspectrum object could have been initialized directly as
ps = Powerspectrum(lc, norm="leahy")
However, we recommend using this explicit syntax, for clarity. Equivalently, one can initialize a Powerspectrum object:
from an EventList object as
bin_time = 0.1
   ps = Powerspectrum.from_events(events, dt=bin_time, norm="leahy")
   where the light curve, uniformly binned at 0.1 s, is created internally.
from a numpy array of times expressed in seconds, as
   bin_time = 0.1
   ps = Powerspectrum.from_events(times, dt=bin_time, gti=[[t0, t1], [t2, t3], ...], norm="leahy")
   where the light curve, uniformly binned at 0.1 s, is created internally, and the good time intervals (time interval where the instrument was collecting data nominally) are passed by hand.
from an iterable of light curves
   ps = Powerspectrum.from_lc_iter(lc_iterable, norm="leahy")
   where lc_iterable is any iterable of Lightcurve objects (list, tuple, generator, etc.)
Since the negative Fourier frequencies (and their associated powers) are discarded, the number of time bins per segment n is twice the length of freq and power.
End of explanation
print(ps.freq)
print(ps.power)
print(ps.df)
print(ps.m)
print(ps.n)
print(ps.nphots1)
Explanation: Properties
A Powerspectrum object has the following properties :
freq : Numpy array of mid-bin frequencies that the Fourier transform samples.
power : Numpy array of the power spectrum.
df : The frequency resolution.
m : The number of power spectra averaged together. For a Powerspectrum of a single segment, m=1.
n : The number of data points (time bins) in one segment of the light curve.
nphots1 : The total number of photons in the light curve.
norm : The normalization, one of leahy (Leahy et al. 1983), abs (absolute rms), frac (fractional rms), or none
End of explanation
fig, ax1 = plt.subplots(1,1,figsize=(9,6), sharex=True)
ax1.plot(ps.freq, ps.power, lw=2, color='blue')
ax1.set_ylabel("Frequency (Hz)", fontproperties=font_prop)
ax1.set_ylabel("Power (raw)", fontproperties=font_prop)
ax1.set_yscale('log')
ax1.tick_params(axis='x', labelsize=16)
ax1.tick_params(axis='y', labelsize=16)
ax1.tick_params(which='major', width=1.5, length=7)
ax1.tick_params(which='minor', width=1.5, length=4)
for axis in ['top', 'bottom', 'left', 'right']:
    ax1.spines[axis].set_linewidth(1.5)
plt.show()
Explanation: We can plot the power as a function of Fourier frequency. Notice how there's a spike at our signal frequency of 2 Hz!
End of explanation
long_dt = 0.03125  # seconds
long_exposure = 1600.  # seconds
long_times = np.arange(0, long_exposure, long_dt)  # seconds
# In count rate units here
long_signal = 300 * np.sin(2.*np.pi*long_times/0.5) + 1000
# Multiply by dt to get count units, then add Poisson noise
long_noisy = np.random.poisson(long_signal*dt)
long_lc = Lightcurve(long_times, long_noisy, dt=long_dt, skip_checks=True)
fig, ax = plt.subplots(1,1,figsize=(10,6))
ax.plot(long_lc.time, long_lc.counts, lw=2, color='blue')
ax.set_xlim(0,20)
ax.set_xlabel("Time (s)", fontproperties=font_prop)
ax.set_ylabel("Counts (cts)", fontproperties=font_prop)
ax.tick_params(axis='x', labelsize=16)
ax.tick_params(axis='y', labelsize=16)
ax.tick_params(which='major', width=1.5, length=7)
ax.tick_params(which='minor', width=1.5, length=4)
plt.show()
Explanation: You'll notice that the power spectrum is a bit noisy. This is because we're only using one segment of data. Let's try averaging together power spectra from multiple segments of data.
Averaged power spectrum example
You could use a long Lightcurve and have AveragedPowerspectrum chop it into specified segments, or give a list of Lightcurves where each segment of Lightcurve is the same length. We'll show the first way here.
1. Create a long light curve.
Generate an array of relative timestamps that's 1600 seconds long, and a signal in count units, with the same properties as the previous example. We then add Poisson noise and turn it into a Lightcurve object.
End of explanation
avg_ps = AveragedPowerspectrum.from_lightcurve(long_lc, 8., norm="leahy")
Explanation: 2. Pass the light curve to the AveragedPowerspectrum class with a specified segment_size.
If the exposure (length) of the light curve cannot be divided by segment_size with a remainder of zero, the last incomplete segment is thrown out, to avoid signal artefacts. Here we're using 8 second segments.
End of explanation
print("Number of segments: %d" % avg_ps.m)
Explanation: We can check how many segments were averaged together by printing the m attribute.
End of explanation
fig, ax1 = plt.subplots(1,1,figsize=(9,6))
ax1.plot(avg_ps.freq, avg_ps.power, lw=2, color='blue')
ax1.set_xlabel("Frequency (Hz)", fontproperties=font_prop)
ax1.set_ylabel("Power (raw)", fontproperties=font_prop)
ax1.set_yscale('log')
ax1.tick_params(axis='x', labelsize=16)
ax1.tick_params(axis='y', labelsize=16)
ax1.tick_params(which='major', width=1.5, length=7)
ax1.tick_params(which='minor', width=1.5, length=4)
for axis in ['top', 'bottom', 'left', 'right']:
    ax1.spines[axis].set_linewidth(1.5)
plt.show()
Explanation: AveragedPowerspectrum has the same properties as Powerspectrum, but with m $>$1.
Let's plot the averaged power spectrum!
End of explanation
avg_ps_leahy = AveragedPowerspectrum.from_lightcurve(long_lc, 8, norm='leahy')
avg_ps_frac = AveragedPowerspectrum.from_lightcurve(long_lc, 8., norm='frac')
avg_ps_abs = AveragedPowerspectrum.from_lightcurve(long_lc, 8., norm='abs')
fig, [ax1, ax2, ax3] = plt.subplots(3,1,figsize=(6,12))
ax1.plot(avg_ps_leahy.freq, avg_ps_leahy.power, lw=2, color='black')
ax1.set_xlabel("Frequency (Hz)", fontproperties=font_prop)
ax1.set_ylabel("Power (Leahy)", fontproperties=font_prop)
ax1.set_yscale('log')
ax1.tick_params(axis='x', labelsize=14)
ax1.tick_params(axis='y', labelsize=14)
ax1.tick_params(which='major', width=1.5, length=7)
ax1.tick_params(which='minor', width=1.5, length=4)
ax1.set_title("Leahy norm.", fontproperties=font_prop)
    
ax2.plot(avg_ps_frac.freq, avg_ps_frac.power, lw=2, color='black')
ax2.set_xlabel("Frequency (Hz)", fontproperties=font_prop)
ax2.set_ylabel("Power (rms)", fontproperties=font_prop)
ax2.tick_params(axis='x', labelsize=14)
ax2.tick_params(axis='y', labelsize=14)
ax2.set_yscale('log')
ax2.tick_params(which='major', width=1.5, length=7)
ax2.tick_params(which='minor', width=1.5, length=4)
ax2.set_title("Fractional rms-squared norm.", fontproperties=font_prop)
ax3.plot(avg_ps_abs.freq, avg_ps_abs.power, lw=2, color='black')
ax3.set_xlabel("Frequency (Hz)", fontproperties=font_prop)
ax3.set_ylabel("Power (abs)", fontproperties=font_prop)
ax3.tick_params(axis='x', labelsize=14)
ax3.tick_params(axis='y', labelsize=14)
ax3.set_yscale('log')
ax3.tick_params(which='major', width=1.5, length=7)
ax3.tick_params(which='minor', width=1.5, length=4)
ax3.set_title("Absolute rms-squared norm.", fontproperties=font_prop)
for axis in ['top', 'bottom', 'left', 'right']:
    ax1.spines[axis].set_linewidth(1.5)
    ax2.spines[axis].set_linewidth(1.5)
    ax3.spines[axis].set_linewidth(1.5)
plt.tight_layout()
plt.show()
Explanation: Now we'll show examples of all the things you can do with a Powerspectrum or AveragedPowerspectrum object using built-in stingray methods.
Normalizating the power spectrum
The three kinds of normalization are:
* leahy: Leahy normalization. Makes the Poisson noise level $= 2$. See Leahy et al. 1983, ApJ, 266, 160L. 
* frac: Fractional rms-squared normalization, also known as rms normalization. Makes the Poisson noise level $= 2 / meanrate$. See Belloni & Hasinger 1990, A&A, 227, L33, and Miyamoto et al. 1992, ApJ, 391, L21.
* abs: Absolute rms-squared normalization, also known as absolute normalization. Makes the Poisson noise level $= 2 \times meanrate$. See insert citation.
* none: No normalization applied. This is the default.
End of explanation
print("DF before:", avg_ps.df)
# Both of the following ways are allowed syntax:
# lin_rb_ps = Powerspectrum.rebin(avg_ps, 0.25, method='mean')
lin_rb_ps = avg_ps.rebin(0.25, method='mean')
print("DF after:", lin_rb_ps.df)
Explanation: Re-binning a power spectrum in frequency
Typically, rebinning is done on an averaged, normalized power spectrum.
1. We can linearly re-bin a power spectrum
(although this is not done much in practice)
End of explanation
# Both of the following ways are allowed syntax:
# log_rb_ps, log_rb_freq, binning = Powerspectrum.rebin_log(avg_ps, f=0.02)
log_rb_ps = ps.rebin_log(f=0.02)
Explanation: 2. And we can logarithmically/geometrically re-bin a power spectrum
In this re-binning, each bin size is 1+f times larger than the previous bin size, where f is user-specified and normally in the range 0.01-0.1. The default value is f=0.01.
End of explanation
print(type(lin_rb_ps))
Explanation: Like rebin, rebin_log returns a Powerspectrum or AveragedPowerspectrum object (depending on the input object):
End of explanation
long_norm = (long_noisy - long_noisy.mean()) / long_noisy.max()
err = np.sqrt(long_noisy.mean()) / long_noisy.max()
long_lc_gauss = Lightcurve(long_times, long_norm, err=np.zeros_like(long_norm) + err, dt=long_dt, skip_checks=True, err_dist='gauss')
fig, ax = plt.subplots(1,1,figsize=(10, 6))
ax.plot(long_lc.time, long_lc.counts, lw=2, color='blue', label='Original light curve')
ax.plot(long_lc_gauss.time, long_lc_gauss.counts, lw=2, color='red', label='Normalized light curve')
ax.set_xlim(0,20)
ax.set_xlabel("Time (s)", fontproperties=font_prop)
ax.set_ylabel("Counts (cts)", fontproperties=font_prop)
ax.tick_params(axis='x', labelsize=16)
ax.tick_params(axis='y', labelsize=16)
ax.tick_params(which='major', width=1.5, length=7)
ax.tick_params(which='minor', width=1.5, length=4)
plt.legend()
plt.show()
avg_ps_gauss_leahy = AveragedPowerspectrum.from_lightcurve(long_lc_gauss, 8, norm='leahy')
avg_ps_gauss_frac = AveragedPowerspectrum.from_lightcurve(long_lc_gauss, 8., norm='frac')
avg_ps_gauss_abs = AveragedPowerspectrum.from_lightcurve(long_lc_gauss, 8., norm='abs')
fig, [ax1, ax2, ax3] = plt.subplots(3,1,figsize=(6,12))
ax1.plot(avg_ps_leahy.freq, avg_ps_leahy.power, lw=2, color='black')
ax1.plot(avg_ps_gauss_leahy.freq, avg_ps_gauss_leahy.power, lw=2, color='red', zorder=10)
ax1.set_xlabel("Frequency (Hz)", fontproperties=font_prop)
ax1.set_ylabel("Power (Leahy)", fontproperties=font_prop)
ax1.set_yscale('log')
ax1.tick_params(axis='x', labelsize=14)
ax1.tick_params(axis='y', labelsize=14)
ax1.tick_params(which='major', width=1.5, length=7)
ax1.tick_params(which='minor', width=1.5, length=4)
ax1.set_title("Leahy norm.", fontproperties=font_prop)
  
ax2.plot(avg_ps_frac.freq, avg_ps_frac.power, lw=2, color='black')
ax2.plot(avg_ps_gauss_frac.freq, avg_ps_gauss_frac.power, lw=2, color='red')
ax2.set_xlabel("Frequency (Hz)", fontproperties=font_prop)
ax2.set_ylabel("Power (rms)", fontproperties=font_prop)
ax2.tick_params(axis='x', labelsize=14)
ax2.tick_params(axis='y', labelsize=14)
ax2.set_yscale('log')
ax2.tick_params(which='major', width=1.5, length=7)
ax2.tick_params(which='minor', width=1.5, length=4)
ax2.set_title("Fractional rms-squared norm.", fontproperties=font_prop)
ax3.plot(avg_ps_abs.freq, avg_ps_abs.power, lw=2, color='black')
ax3.plot(avg_ps_gauss_abs.freq, avg_ps_gauss_abs.power, lw=2, color='red')
ax3.set_xlabel("Frequency (Hz)", fontproperties=font_prop)
ax3.set_ylabel("Power (abs)", fontproperties=font_prop)
ax3.tick_params(axis='x', labelsize=14)
ax3.tick_params(axis='y', labelsize=14)
ax3.set_yscale('log')
ax3.tick_params(which='major', width=1.5, length=7)
ax3.tick_params(which='minor', width=1.5, length=4)
ax3.set_title("Absolute rms-squared norm.", fontproperties=font_prop)
for axis in ['top', 'bottom', 'left', 'right']:
    ax1.spines[axis].set_linewidth(1.5)
    ax2.spines[axis].set_linewidth(1.5)
    ax3.spines[axis].set_linewidth(1.5)
plt.tight_layout()
plt.show()
Explanation: Power spectra of normal-distributed light curves
Starting with Stingray 0.3, we can also get Leahy-normalized power spectra of normally-distributed light curves.
Let us calculate such a light curve by subtracting the noise level and normalizing
End of explanation
<END_TASK>