Merge pull request #56 from AUAS-Pulsar/dedisp_implementation

Dedisp implementation

Author: TVH7

dedisperse/__init__.py

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

dedisperse/dedisperse.py

@@ -1,19 +1,29 @@
 '''
-    Dedisperses data
+Dedisperses data
 '''
 # pylint: disable-msg=C0103
 import numpy as np
 
-def dedisperse(samples, dm):
+def dedisperse(samples, highest_x=None, max_delay=0, disperion_measure=None):
     '''
     This method performs dedispersion on the filterbank data
+    The maximum_delay specifies between the currently considered pulsar signal and the next pulsar
+    signal should be
+    The highest_x specifies the amount of intensities that are used for estimating the minimum
+    pulsar intensity
     '''
-    samples = np.asarray(samples)
-    # Distribute the DM over the amount of samples
-    delays_per_sample = np.round(np.linspace(dm, 0, samples.shape[1])).astype(int)
+
+    # Check if parameters contain a Dispersion Measure, if not, estimate one
+    if disperion_measure is None:
+        # Estimates the minimum for an intensity to be considered a pulsar
+        pulsar_intensity = find_estimation_intensity(samples, highest_x)
+        disperion_measure = find_dm(samples, pulsar_intensity, max_delay)
+
+    # Distribute the Dispersion Measure over the amount of samples
+    delays_per_sample = np.round(np.linspace(disperion_measure, 0, samples.shape[1])).astype(int)
 
     # Loop over the frequencies
-    for i, _ in enumerate(delays_per_sample):
+    for i in range(delays_per_sample.size):
 
         # Temporary array that is used to later delay the frequency
         temporary_samples = []
@@ -25,3 +35,80 @@ def dedisperse(samples, dm):
         samples[:, i] = np.roll(temporary_samples, delays_per_sample[i])
 
     return samples
+
+def find_dm(samples, pulsar_intensity, max_delay):
+    '''
+    This method attempts to find a dispersion measure
+    '''
+
+    # Loop through the samples to find a pulsar intensity to start calculating from
+    for s, sample in enumerate(samples[:, 0]):
+
+        # If the sample meets the minimum intensity, attempt to find a line continuing from
+        # this intensity
+        if sample > pulsar_intensity:
+            start_sample_index = s
+
+            # Attempt to find a line, line_coordinates contains first and last index of the pulsar
+            line_coordinates = find_line(samples, start_sample_index, max_delay, pulsar_intensity)
+
+            # If a line is found, calculate and return the dispersion measure
+            if line_coordinates is not None:
+                disperion_measure = line_coordinates[1] - line_coordinates[0]
+                return disperion_measure
+
+    return None
+
+
+def find_line(samples, start_sample_index, max_delay, pulsar_intensity):
+    '''
+    This method tries to find a line starting from the sample index given in the parameters
+    it stops if there is no intensity within the max_delay higher than the average_intensity
+    '''
+
+    previous_index = start_sample_index
+
+    failed_to_find_line = True
+
+    # Loop through the frequencies
+    for f in range(samples[1].size):
+
+        # Loop through previous intensity until the max delay is reached
+        for i, intensity in enumerate(samples[:, f][previous_index:previous_index + max_delay]):
+
+            # Skip the first frequency, since we start measuring from the first intensity
+            if f == 0:
+                failed_to_find_line = False
+                break
+
+            # If the intensity is higher than the pulsar_intensity, continue finding a signal
+            if intensity > pulsar_intensity:
+                previous_index = previous_index + i
+                failed_to_find_line = False
+                break
+
+        # If there is no line found, return None
+        if failed_to_find_line:
+            return None
+
+    # If all frequencies are looped, a continuous signal is found,
+    # so we return the first and last index of the line
+    return start_sample_index, previous_index
+
+def find_estimation_intensity(samples, highest_x):
+    '''
+    This method finds the average intensity for the highest x intensities
+    The average_intensity is considered a requirement for intensities to be considered a pulsar
+    '''
+
+    # Sum of all intensities
+    sum_intensities = 0
+
+    # Looks for the top x highest intensities in the samples and adds them up together
+    for sample in samples:
+        sum_intensities += np.sum(sorted(sample, reverse=True)[:highest_x])
+
+    # Calculates the average_intensity
+    average_intensity = (sum_intensities) / (samples.shape[0] * highest_x)
+
+    return average_intensity

examples/dedisperse_samples.py

@@ -16,43 +16,15 @@
 
 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")
 
 special_pspm.read_filterbank()
 
 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)
-
 # Dispersion Measure
-DM = 230
+DM = 235
 
 plt.subplot(2,1,1)
 data, extent = waterfall_plot(samples, frequencies)

examples/dedisperse_stream.py

@@ -0,0 +1,77 @@
+'''
+    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 clipping import clipping
+
+# Read filterbank data,
+
+# Standard file
+special_pspm = fb.Filterbank(filename = "./my_special_pspm.fil")
+highest_x=10
+max_delay=10
+
+special_pspm.read_filterbank()
+
+frequencies, samples = special_pspm.select_data()
+
+# Plot the original 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')
+
+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()
+
+# Dedisperse the samples
+samples = dedisperse.dedisperse(samples, highest_x, max_delay)
+
+# Plot the dedispersed 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()

pipeline/pipeline.py

@@ -112,12 +112,13 @@ def measure_methods(self, stopwatch, fil_data, freqs, DM, scale):
             Run and time all methods/modules
         """
         # clipping
+
         time_clipping = timer()
         _, _ = clipping.clipping(freqs, fil_data)
         stopwatch['time_clipping'] = timer() - time_clipping
         # dedisperse
         time_dedisp = timer()
-        fil_data = dedisperse.dedisperse(fil_data, DM)
+        fil_data = dedisperse.dedisperse(fil_data, disperion_measure=DM)
         stopwatch['time_dedisp'] = timer() - time_dedisp
         # timeseries
         time_t_series = timer()

tests/test_dedisperse.py

@@ -3,6 +3,7 @@
 """
 
 import unittest
+import numpy as np
 
 from .context import dedisperse # pylint: disable-msg=E0611
 
@@ -25,5 +26,5 @@ def test_dedisperse(self):
             expect moved frequencies per sample
         """
         disp_measure = 6
-        results = dedisperse.dedisperse(self.samples, disp_measure)
+        results = dedisperse.dedisperse(np.array(self.samples), disperion_measure=disp_measure)
         self.assertListEqual(list(results[len(self.samples)-1]), [10]*7)