CrossCorrelation¶
This Tutorial is intended to give a demostration of How to make a CrossCorrelation Object in Stingray Library.
[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)
CrossCorrelation Example¶
1. Create two light curves¶
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.
[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.
[6]:
lc1 = Lightcurve(times, noisy_1)
lc2 = Lightcurve(times, noisy_2)
len(lc1)
[6]:
320
[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()
2. Create a CrossCorrelation Object from two Light curves created above¶
To create a CrossCorrelation Object from LightCurves, simply pass both Lightvurves created above into the CrossCorrelation.
[8]:
cr = CrossCorrelation(lc1, lc2)
Now, Cross Correlation values are stored in attribute corr, which is called below.
[9]:
cr.corr[:10]
[9]:
array([ 201.553125 , 1412.10121094, 2828.54304688, 3948.95050781,
5370.02359375, 5750.04355469, 6222.50101563, 6664.92722656,
5969.0503125 , 6770.80464844])
[10]:
# Time Resolution for Cross Correlation is same as that of each of the Lightcurves
cr.dt
[10]:
0.03125
3. Plot Cross Correlation for Different lags¶
To visulaize correlation for different values of time lags, simply call plot function on cs.
[11]:
cr.plot(labels = ['Time Lags (seconds)','Correlation'])
[11]:
<module 'matplotlib.pyplot' from 'C:\\Users\\Haroon Rashid\\Anaconda3\\lib\\site-packages\\matplotlib\\pyplot.py'>
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.
[12]:
cr.time_shift #seconds
[12]:
0.26645768025078276
Modes of Correlation¶
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.
[13]:
cr.mode
[13]:
'same'
The number of data points in corr and largest lightcurve are same in this mode.
[14]:
cr.n
[14]:
320
Creating CrossCorrelation with full mode now using same data as above.
[15]:
cr1 = CrossCorrelation(lc1, lc2, mode = 'full')
[16]:
cr1.plot()
[16]:
<module 'matplotlib.pyplot' from 'C:\\Users\\Haroon Rashid\\Anaconda3\\lib\\site-packages\\matplotlib\\pyplot.py'>
[17]:
cr1.mode
[17]:
'full'
Full mode does a full cross-correlation.
[18]:
cr1.n
[18]:
639
Another Example¶
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).
[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)
[20]:
time_shift
[20]:
0.83333333333333348
[21]:
cs.plot( ['Time Lags (seconds)','Correlation'])
[21]:
<module 'matplotlib.pyplot' from 'C:\\Users\\Haroon Rashid\\Anaconda3\\lib\\site-packages\\matplotlib\\pyplot.py'>
Yet another Example with longer Lightcurve¶
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.
[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.
[23]:
lc1 = Lightcurve(times, signal_1)
lc2 = Lightcurve(times, signal_2)
len(lc1)
[23]:
500000
[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.
[25]:
cs = CrossCorrelation(lc1, lc2)
print('Done')
Done
[26]:
cs.corr[:50]
[26]:
array([ 2.86241768e-05, 4.71238867e+06, 9.42481318e+06,
1.41372717e+07, 1.88497623e+07, 2.35622831e+07,
2.82748324e+07, 3.29874082e+07, 3.77000087e+07,
4.24126319e+07, 4.71252762e+07, 5.18379395e+07,
5.65506201e+07, 6.12633160e+07, 6.59760255e+07,
7.06887466e+07, 7.54014775e+07, 8.01142163e+07,
8.48269612e+07, 8.95397103e+07, 9.42524618e+07,
9.89652137e+07, 1.03677964e+08, 1.08390712e+08,
1.13103454e+08, 1.17816189e+08, 1.22528916e+08,
1.27241631e+08, 1.31954335e+08, 1.36667023e+08,
1.41379696e+08, 1.46092350e+08, 1.50804985e+08,
1.55517598e+08, 1.60230186e+08, 1.64942750e+08,
1.69655286e+08, 1.74367792e+08, 1.79080268e+08,
1.83792710e+08, 1.88505118e+08, 1.93217489e+08,
1.97929821e+08, 2.02642113e+08, 2.07354363e+08,
2.12066568e+08, 2.16778727e+08, 2.21490839e+08,
2.26202900e+08, 2.30914910e+08])
[27]:
# Time Resolution for Cross Correlation is same as that of each of the Lightcurves
cs.dt
[27]:
9.9999999999766942e-05
[28]:
cs.plot( ['Time Lags (seconds)','Correlation'])
[28]:
<module 'matplotlib.pyplot' from 'C:\\Users\\Haroon Rashid\\Anaconda3\\lib\\site-packages\\matplotlib\\pyplot.py'>
[29]:
cs.time_shift #seconds
[29]:
0.2495504991004161
time_shift is very close to 0.25 sec, in this case.
AutoCorrelation¶
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.
[30]:
from stingray.crosscorrelation import AutoCorrelation
To create AutoCorrelation object, simply pass lightcurve into AutoCorrelation Constructor. Using same Lighrcurve created above to demonstrate AutoCorrelation.
[31]:
lc = lc1
[32]:
ac = AutoCorrelation(lc)
ac.n
[32]:
500000
[33]:
ac.corr[:10]
[33]:
array([ 1.12500000e+10, 1.12499978e+10, 1.12499911e+10,
1.12499800e+10, 1.12499645e+10, 1.12499445e+10,
1.12499201e+10, 1.12498912e+10, 1.12498579e+10,
1.12498201e+10])
[34]:
ac.time_lags
[34]:
array([-25. , -24.9999, -24.9998, ..., 24.9998, 24.9999, 25. ])
time_Shift for AutoCorrelation is always zero. Since signals are maximally correlated at zero lag.
[35]:
ac.time_shift
[35]:
5.0000099997734535e-05
[36]:
ac.plot()
[36]:
<module 'matplotlib.pyplot' from 'C:\\Users\\Haroon Rashid\\Anaconda3\\lib\\site-packages\\matplotlib\\pyplot.py'>
Another Example¶
Another example is demonstrated using a Lightcurve with Poisson Noise.
[37]:
from stingray.crosscorrelation import AutoCorrelation
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
[38]:
ac = AutoCorrelation(lc, mode = 'full')
[39]:
ac.corr
[39]:
array([-0.00487599, -0.00485198, -0.99992797, ..., -0.99992797,
-0.00485198, -0.00487599])
[40]:
ac.time_lags
[40]:
array([-19.999, -19.998, -19.997, ..., 19.997, 19.998, 19.999])
[41]:
ac.time_shift
[41]:
0.0
[42]:
ac.plot()
[42]:
<module 'matplotlib.pyplot' from 'C:\\Users\\Haroon Rashid\\Anaconda3\\lib\\site-packages\\matplotlib\\pyplot.py'>
