Implemented dispersion measure guessing

Author: robertkleef

dedisperse/__init__.py

@@ -1,4 +1,4 @@
 """
     import file for dedisperse.py
 """
-from .dedisperse import dedisperse
+from .dedisperse import dedisperse, estimate_dm, find_initial_signal

dedisperse/dedisperse.py

@@ -4,11 +4,15 @@
 # pylint: disable-msg=C0103
 import numpy as np
 
-def dedisperse(samples, dm):
+def dedisperse(samples, dm=None):
     '''
     This method performs dedispersion on the filterbank data
     '''
-    samples = np.asarray(samples)
+
+    if dm is None:
+        print("Finding possible DM's")
+        dm = estimate_dm(samples)
+
     # Distribute the DM over the amount of samples
     delays_per_sample = np.round(np.linspace(dm, 0, samples.shape[1])).astype(int)
 
@@ -25,3 +29,65 @@ def dedisperse(samples, dm):
         samples[:, i] = np.roll(temporary_samples, delays_per_sample[i])
 
     return samples
+
+
+def estimate_dm(samples):
+    '''
+    This method attempts to find a dispersion measure
+    '''
+
+    # Tuple of frequency index and sample index
+    initial_signal_point = find_initial_signal(samples)
+    last_sample = initial_signal_point
+
+    for i, frequency in enumerate(samples[1]):
+        for j, data_point in enumerate(samples[:, i]):
+            #print(samples[:, i][previous_time_sample:].shape[0])
+            if(j > last_sample[1]):
+                if(data_point > 10):
+                    last_sample = i, j
+                    #print("At Frequency ", i, " found on Time sample ", j, " - ", data_point)
+                    break
+
+    highest_difference = 0
+    '''
+    for s, samples in enumerate(samples):
+        for i, data_point in enumerate(samples):
+            if(i > highest_difference):                
+                if(data_point > 10):
+                    print(s, i, " - ", data_point)                                        
+                    highest_difference = i
+                break
+    '''
+    estimated_dm = last_sample[1] - initial_signal_point[1]
+    print("Estimated DM is", estimated_dm)
+    return estimated_dm
+
+
+def find_initial_signal(samples):
+    '''
+    This method attempts to find a viable data point to start estimating a dispersion measure from
+    '''
+
+    # Tuple of frequency index and sample index
+    lowest_sample = 0, samples.shape[0]
+    
+    for i, frequency in enumerate(samples[1]):
+        for j, data_point in enumerate(samples[:, i]):
+            if(j < lowest_sample[1]):
+                if(data_point > 1):
+                    print("Initial signal found on freq, sample", i, j)
+                    return i, j
+                    '''
+                    print(lowest_sample, " ", j)
+                    print("At Frequency ", i, " found on Time sample ", j, " - ", data_point)
+                    lowest_sample = i, j
+                    break
+            if(lowest_sample[1] == 0):
+                print("", lowest_sample)
+                return i, j
+                    '''
+
+    print("NO INITIAL SIGNAL FOUND")
+    return None
+    

examples/dedisperse_samples.py

@@ -16,31 +16,6 @@
 
 from timeseries.timeseries import Timeseries
 
-def print_possible_dm(data):
-        '''
-        This method is for testing purposes and should probably be removed later on
-        '''
-        first_sample = 400
-
-        for i, data_point in enumerate(data[first_sample]):
-                if(data_point > 10):
-                        print("first sample ", i, " - ", data_point)
-                        break
-
-        highest_difference = 0
-
-        for s, samples in enumerate(data):
-                for i, data_point in enumerate(samples):
-                        if(i > highest_difference):
-                                
-                                if(data_point > 10):
-                                        print(s, i, " - ", data_point)                                        
-                                        highest_difference = i
-                                break
-
-        print(highest_difference)
-
-
 # Read filterbank data
 special_pspm = fb.Filterbank(filename = "../data/my_special_pspm.fil")
 
@@ -49,10 +24,10 @@ def print_possible_dm(data):
 frequencies, samples = special_pspm.select_data()
 
 # Use this if you have your own file with a clear pulsar signal, this method assumes all signals other than the pulsar are lower than 10
-#print_possible_dm(dispersed_samples)
+print_possible_dm(samples)
 
 # Dispersion Measure
-DM = 230
+DM = 240
 
 plt.subplot(2,1,1)
 data, extent = waterfall_plot(samples, frequencies)

examples/dedisperse_stream.py

@@ -0,0 +1,76 @@
+'''
+    Example that dedisperses filterbank data and plots it together with a timeseries plot
+'''
+# pylint: disable-all
+import os,sys,inspect
+currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
+parentdir = os.path.dirname(currentdir)
+sys.path.insert(0,parentdir)
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import filterbank.filterbank as fb
+import dedisperse as dedisperse
+from plot.static_waterfall import waterfall_plot
+
+from timeseries.timeseries import Timeseries
+
+from clipping import clipping
+
+# Read filterbank data
+special_pspm = fb.Filterbank(filename = "../data/my_uber_pspm.fil")
+
+special_pspm.read_filterbank()
+
+frequencies, samples = special_pspm.select_data()
+
+
+plt.subplot(2,1,1)
+data, extent = waterfall_plot(samples, frequencies)
+
+img = plt.imshow(data.T,
+                 aspect='auto',
+                 origin='lower',
+                 rasterized=True,
+                 interpolation='nearest',
+                 extent=extent,
+                 cmap='cubehelix')
+
+time_series = []
+
+for sample_set in samples:
+            summed_sample = np.sum(sample_set)
+            time_series.append(summed_sample)
+
+plt.subplot(2,1,2)
+plt.plot(time_series)
+plt.show()
+
+#clipped_samples = clipping(frequencies, samples)
+samples = dedisperse.dedisperse(samples)
+#samples = dedisperse.find_lowest_pulsar(samples)
+#samples = dedisperse.estimate_dm(samples)
+
+
+plt.subplot(2,1,1)
+data, extent = waterfall_plot(samples, frequencies)
+
+img = plt.imshow(data.T,
+                 aspect='auto',
+                 origin='lower',
+                 rasterized=True,
+                 interpolation='nearest',
+                 extent=extent,
+                 cmap='cubehelix')
+plt.colorbar()
+
+time_series = []
+
+for sample_set in samples:
+            summed_sample = np.sum(sample_set)
+            time_series.append(summed_sample)
+
+plt.subplot(2,1,2)
+plt.plot(time_series)
+plt.show()