Commit a5864572 authored by Tammo Jan Dijkema's avatar Tammo Jan Dijkema
Browse files

Update with right figure

parent dad9d404
%% Cell type:code id: tags:
``` python
import numpy as np
```
%% Cell type:code id: tags:
``` python
import matplotlib.pyplot as plt
```
%% Cell type:code id: tags:
``` python
%matplotlib inline
```
%% Cell type:code id: tags:
``` python
origfile = "/Users/dijkema/Desktop/B0329+54.2016.11.18.1038.5min.dat"
```
%% Cell type:code id: tags:
``` python
origdata = np.fromfile(origfile, dtype='>u4')
```
%% Cell type:markdown id: tags:
The data consists of chunks of 256 integers:
* One 4-byte integer indicating a packet sequence number
* 255 4-byte integers with the data for 255 channels
The sequence numbers should be increasing (wrapping at the limit of what a 4-byte integer can hold). If a packet would be dropped, flagged data should be inserted. I have not implemented that yet.
%% Cell type:code id: tags:
``` python
if np.sum(origdata[256::256]-origdata[0:-256:256]-1)>0:
print("Missing packets, something will go wrong")
```
%% Cell type:code id: tags:
``` python
ntimes = len(origdata)//256 + 1
```
%% Cell type:code id: tags:
``` python
data = np.zeros((ntimes, 255), dtype='u4')
```
%% Cell type:code id: tags:
``` python
for t in range(ntimes-1):
data[t, :] = origdata[t*256+1:(t+1)*256]
```
%% Cell type:code id: tags:
``` python
data = np.transpose(data)
```
%% Cell type:markdown id: tags:
Flag data that has value larger than 1000 (this could be tweaked).
%% Cell type:code id: tags:
``` python
data_masked = np.ma.masked_greater(data, 1000)
```
%% Cell type:markdown id: tags:
Flag entire channels where there is too much bad data:
%% Cell type:code id: tags:
``` python
badchannels = (data.mean(axis=1)>300).reshape((255,1))
```
%% Cell type:code id: tags:
``` python
data_masked.mask = np.logical_or(data_masked.mask, badchannels)
```
%% Cell type:code id: tags:
``` python
# Make figure larger
# Make figure and labels larger
import matplotlib
matplotlib.rcParams['figure.figsize'] = [20,8]
matplotlib.rcParams['xtick.labelsize'] = 'x-large'
matplotlib.rcParams['ytick.labelsize'] = 'x-large'
matplotlib.rcParams['axes.labelsize'] = 'x-large'
```
%% Cell type:code id: tags:
``` python
# Duration of one sample
timeres = 0.00046811428
```
%% Cell type:code id: tags:
``` python
# Insert the period of the pulsar for some folding
period = 0.71448
offset = 1*period / timeres
```
%% Cell type:code id: tags:
``` python
# Plot exactly 2.5 seconds
maxsample = int(2.5/timeres)
```
%% Cell type:code id: tags:
``` python
foldmax = 200
period = 0.71451
offset = period / timeres
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
folded = data_masked[:, 0:maxsample].copy()
for j in range(1, foldmax+1):
folded += data_masked[:, int(j*offset):int(j*offset)+maxsample]
# Expand mask: if one of folded samples is masked, mask all
folded.mask = np.ma.mask_or(folded.mask, data_masked[:, int(j*offset):int(j*offset)+maxsample].mask)
im = ax.imshow(folded,
aspect='auto',
origin='lower',
vmin=0, vmax=np.percentile(folded[np.logical_not(folded.mask)],96)#,
# extent = [0*timeres, maxsample*timeres, 406.4, 441.12]
)
ax.set_title("Period: "+str(period))
ax.set_xlabel("Time (s)")
ax.set_ylabel("Channel (MHz)")
fig.colorbar(im)
fig.colorbar(im);
```
%% Output
<matplotlib.colorbar.Colorbar at 0x11abc75f8>
%% Cell type:markdown id: tags:
FFT one channel to determine the period of the pulsar (this does not work yet):
%% Cell type:code id: tags:
``` python
chan = 126
data_filled = np.ma.filled(data_masked, np.ma.mean(data_masked))
fft = np.abs(np.fft.fft(data_filled[chan,:]))
fft[0:5] = 0
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
im = plt.plot(fft[0:len(fft)//2])
ax.set_title("Channel: "+str(chan));
```
%% Output
%% Cell type:code id: tags:
``` python
```
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment