In [4]:
import numpy as np
from stingray import Lightcurve
from stingray.crosscorrelation import CrossCorrelation
import matplotlib.pyplot as plt
import matplotlib.font_manager as font_manager
%matplotlib inline
font_prop = font_manager.FontProperties(size=16)
There are two ways to create a Lightcurve.
1) Using an array of time stamps and an array of counts.
2) From the Photon Arrival times.
In this example, Lightcurve is created using arrays of time stamps and counts.
Generate an array of relative timestamps that's 10 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 of pi/2 radians, and mean = 1000 cts/s. We then add Poisson noise to the light curve.
In [5]:
dt = 0.03125 # seconds
exposure = 10. # seconds
freq = 1 # Hz
times = np.arange(0, exposure, dt) # seconds
signal_1 = 300 * np.sin(2.*np.pi*freq*times) + 1000 # counts/s
signal_2 = 300 * np.sin(2.*np.pi*freq*times + np.pi/2) + 1000 # counts/s
noisy_1 = np.random.poisson(signal_1*dt) # counts
noisy_2 = np.random.poisson(signal_2*dt) # counts
Now let's turn noisy_1 and noisy_2 into Lightcurve objects. This way we have two Lightcurves to calculate CrossCorrelation.
In [6]:
lc1 = Lightcurve(times, noisy_1)
lc2 = Lightcurve(times, noisy_2)
len(lc1)
Out[6]:
In [7]:
fig, ax = plt.subplots(1,1,figsize=(10,6))
ax.plot(lc1.time, lc1.counts, lw=2, color='blue')
ax.plot(lc1.time, lc2.counts, lw=2, color='red')
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()
In [8]:
cr = CrossCorrelation(lc1, lc2)
Now, Cross Correlation values are stored in attribute corr, which is called below.
In [9]:
cr.corr[:10]
Out[9]:
In [10]:
# Time Resolution for Cross Correlation is same as that of each of the Lightcurves
cr.dt
Out[10]:
In [11]:
cr.plot(labels = ['Time Lags (seconds)','Correlation'])
Out[11]:
Given the Phase offset of pi/2 between two lightcurves created above, and freq=1 Hz, time_shift should be close to 0.25 sec. Small error is due to time resolution.
In [12]:
cr.time_shift #seconds
Out[12]:
You can also specify an optional argument on modes of cross-correlation.
There are three modes : 1) same 2) valid 3) full
Visit following ink on more details on mode : https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.signal.correlate.html
Default mode is 'same' and it gives output equal to the size of larger lightcurve and is most common in astronomy. You can see mode of your CrossCorrelation by calling mode attribute on the object.
In [13]:
cr.mode
Out[13]:
The number of data points in corr and largest lightcurve are same in this mode.
In [14]:
cr.n
Out[14]:
Creating CrossCorrelation with full mode now using same data as above.
In [15]:
cr1 = CrossCorrelation(lc1, lc2, mode = 'full')
In [16]:
cr1.plot()
Out[16]:
In [17]:
cr1.mode
Out[17]:
Full mode does a full cross-correlation.
In [18]:
cr1.n
Out[18]:
You can also create CrossCorrelation Object by using Cross Correlation data. This can be useful in some cases when you have correlation data and want to calculate time shift for max. correlation. You need to specify time resolution for correlation(default value of 1.0 seconds is used otherwise).
In [19]:
cs = CrossCorrelation()
cs.corr = np.array([ 660, 1790, 3026, 4019, 5164, 6647, 8105, 7023, 6012, 5162])
time_shift, time_lags, n = cs.cal_timeshift(dt=0.5)
In [20]:
time_shift
Out[20]:
In [21]:
cs.plot( ['Time Lags (seconds)','Correlation'])
Out[21]:
I will be using same lightcurves as in the example above but with much longer duration and shorter lags.
Both Lightcurves are chosen to be more or less same with a certain phase shift to demonstrate Correlation in a better way.
Again Generating two signals this time without poission noise so that time lag can be demonstrated. For noisy lightcurves, accurate calculation requires interpolation.
In [22]:
dt = 0.0001 # seconds
exposure = 50. # seconds
freq = 1 # Hz
times = np.arange(0, exposure, dt) # seconds
signal_1 = 300 * np.sin(2.*np.pi*freq*times) + 1000 * dt # counts/s
signal_2 = 200 * np.sin(2.*np.pi*freq*times + np.pi/2) + 900 * dt # counts/s
Converting noisy signals into Lightcurves.
In [23]:
lc1 = Lightcurve(times, signal_1)
lc2 = Lightcurve(times, signal_2)
len(lc1)
Out[23]:
In [24]:
fig, ax = plt.subplots(1,1,figsize=(10,6))
ax.plot(lc1.time, lc1.counts, lw=2, color='blue')
ax.plot(lc1.time, lc2.counts, lw=2, color='red')
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()
Now, creating CrossCorrelation Object by passing lc1 and lc2 into the constructor.
In [25]:
cs = CrossCorrelation(lc1, lc2)
print('Done')
In [26]:
cs.corr[:50]
Out[26]:
In [27]:
# Time Resolution for Cross Correlation is same as that of each of the Lightcurves
cs.dt
Out[27]:
In [28]:
cs.plot( ['Time Lags (seconds)','Correlation'])
Out[28]:
In [29]:
cs.time_shift #seconds
Out[29]:
time_shift is very close to 0.25 sec, in this case.
Stingray has also separate class for AutoCorrelation. AutoCorrealtion is similar to crosscorrelation but involves only One Lightcurve.i.e. Correlation of Lightcurve with itself.
AutoCorrelation is part of stingray.crosscorrelation module. Following line imports AutoCorrelation.
In [30]:
from stingray.crosscorrelation import AutoCorrelation
To create AutoCorrelation object, simply pass lightcurve into AutoCorrelation Constructor.
Using same Lighrcurve created above to demonstrate AutoCorrelation.
In [31]:
lc = lc1
In [32]:
ac = AutoCorrelation(lc)
ac.n
Out[32]:
In [33]:
ac.corr[:10]
Out[33]:
In [34]:
ac.time_lags
Out[34]:
time_Shift for AutoCorrelation is always zero. Since signals are maximally correlated at zero lag.
In [35]:
ac.time_shift
Out[35]:
In [36]:
ac.plot()
Out[36]:
Another example is demonstrated using a Lightcurve with Poisson Noise.
In [37]:
dt = 0.001 # seconds
exposure = 20. # seconds
freq = 1 # Hz
times = np.arange(0, exposure, dt) # seconds
signal_1 = 300 * np.sin(2.*np.pi*freq*times) + 1000 # counts/s
noisy_1 = np.random.poisson(signal_1*dt) # counts
lc = Lightcurve(times, noisy_1)
AutoCorrelation also supports {full,same,valid} modes similar to CrossCorrelation
In [38]:
ac = AutoCorrelation(lc, mode = 'full')
In [39]:
ac.corr
Out[39]:
In [40]:
ac.time_lags
Out[40]:
In [41]:
ac.time_shift
Out[41]:
In [42]:
ac.plot()
Out[42]: