Repository: UCBerkeleySETI/blimpy (language: python)

Issue:
Add functionality to perform dedoppler
This required subclassing numpy array, and adding a dedoppler function


Discussion / hints:
@telegraphic What do you think about this solution to dedoppler (https://github.com/UCBerkeleySETI/blimpy/issues/92)
# [Codecov](https://codecov.io/gh/UCBerkeleySETI/blimpy/pull/220?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=UCBerkeleySETI) Report
> Merging [#220](https://codecov.io/gh/UCBerkeleySETI/blimpy/pull/220?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=UCBerkeleySETI) (d808cf4) into [master](https://codecov.io/gh/UCBerkeleySETI/blimpy/commit/314d3619e2d01332d58e1f1a06e999bea00352e8?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=UCBerkeleySETI) (314d361) will **decrease** coverage by `0.03%`.
> The diff coverage is `73.68%`.

[![Impacted file tree graph](https://codecov.io/gh/UCBerkeleySETI/blimpy/pull/220/graphs/tree.svg?width=650&height=150&src=pr&token=qw1I8qrUk2&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=UCBerkeleySETI)](https://codecov.io/gh/UCBerkeleySETI/blimpy/pull/220?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=UCBerkeleySETI)

```diff
@@            Coverage Diff             @@
##           master     #220      +/-   ##
==========================================
- Coverage   78.65%   78.62%   -0.04%     
==========================================
  Files          60       61       +1     
  Lines        2999     3018      +19     
==========================================
+ Hits         2359     2373      +14     
- Misses        640      645       +5     
```


| [Impacted Files](https://codecov.io/gh/UCBerkeleySETI/blimpy/pull/220?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=UCBerkeleySETI) | Coverage Δ | |
|---|---|---|
| [blimpy/io/waterfall\_blob.py](https://codecov.io/gh/UCBerkeleySETI/blimpy/pull/220/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=UCBerkeleySETI#diff-YmxpbXB5L2lvL3dhdGVyZmFsbF9ibG9iLnB5) | `68.75% <68.75%> (ø)` | |
| [blimpy/io/fil\_reader.py](https://codecov.io/gh/UCBerkeleySETI/blimpy/pull/220/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=UCBerkeleySETI#diff-YmxpbXB5L2lvL2ZpbF9yZWFkZXIucHk=) | `74.81% <100.00%> (+0.57%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/UCBerkeleySETI/blimpy/pull/220?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=UCBerkeleySETI).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=UCBerkeleySETI)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/UCBerkeleySETI/blimpy/pull/220?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=UCBerkeleySETI). Last update [314d361...d808cf4](https://codecov.io/gh/UCBerkeleySETI/blimpy/pull/220?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=UCBerkeleySETI). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=UCBerkeleySETI).

To test:
```
from astropy import units as u
import setigen as stg
import matplotlib.pyplot as plt

frame = stg.Frame(fchans=1024*u.pixel,
                  tchans=32*u.pixel,
                  df=2.7939677238464355*u.Hz,
                  dt=18.253611008*u.s,
                  fch1=6095.214842353016*u.MHz)
noise = frame.add_noise(x_mean=10, noise_type='chi2')
signal = frame.add_signal(stg.constant_path(f_start=frame.get_frequency(index=200),
                                            drift_rate=2*u.Hz/u.s),
                          stg.constant_t_profile(level=frame.get_intensity(snr=30)),
                          stg.gaussian_f_profile(width=40*u.Hz),
                          stg.constant_bp_profile(level=1))

frame.save_fil("test_fil.fil")
```

The above will create a `test_fil.fil` file, and then:

```
from blimpy import Waterfall
import matplotlib.pyplot as plt

wf = Waterfall("test_fil.fil")
f,d = wf.grab_data()
plt.figure(1)
plt.imshow(d, interpolation='nearest', aspect='auto')

plt.figure(2)
d.dedoppler(2) #2 Hz/s
plt.imshow(d, interpolation='nearest', aspect='auto')
plt.show()
```
Hey @wfarah, this is a very clever method! However, I am going to suggest that we follow a similar approach to the plotting, which uses `six.create_bound_method(plot_spectrum, self)`, e.g. https://github.com/UCBerkeleySETI/blimpy/blob/314d3619e2d01332d58e1f1a06e999bea00352e8/blimpy/waterfall.py#L130

So we would:
* create a new directory `signal_processing` (or similar), with `__init__.py` in it.
* create a `dedoppler.py` file in `signal_processing` directory
* add the code with `wf` as the first argument (essentially wf becomes `self`)

I favor this approach as we can then add more methods in `signal_processing`, e.g. dedispersion, without making the Waterfall codebase massive. @texadactyl and @bbrzycki may want to add their 2c here too! 

Code is then:

```python
   def dedoppler(wf, drift_rate):
        """ 
        Simple de-doppler code for a filterbank 
        Parameters:
        ----------
        drift_rate : float
            drift_rate [Hz/s]
        """
        tsamp   = wf.header['tsamp']
        chan_bw = wf.header['foff'] * 1e6 #foff is in MHz
        n_roll = (drift_rate * tsamp) / chan_bw
        for ii in range(wf.shape[0]):
            wf.data[ii] = np.roll(wf.data[ii], -int(n_roll*ii))
```


This code is of course the simplest possible, we could attempt to zero-pad or trim the array and update `fch1` and `nchans` accordingly. 
I see 2 issues:

1. Only Filterbank files are addressed but not HDF5 files.  The suggestion by @telegraphic solves that by waiting until the data matrix is formed in the Waterfall object in addition to avoid further complications in the existing code base.
2. One problem remains with both approaches:  What to do about files with a data matrix that will not fit in RAM?

Number 2 might be approached the same way that the new `stix` utility did for Waterfall plotting: divide and conquer.  But, how/where to store the updated data matrix?  original file?  new file?
I used @wfarah 's Filterbank generation code with ascending=True and @telegraphic 's dedoppler code and got a few plots.  No changes to the rest of blimpy.  Looks like a light sabre to me!  (:
![test_fil](https://user-images.githubusercontent.com/11318756/128576264-f6b4200a-20ec-4a1a-89b3-5116aff80aa5.png)



Source code:
[signal_processing.zip](https://github.com/UCBerkeleySETI/blimpy/files/6948150/signal_processing.zip)


Fix patch:
diff --git a/VERSION-HISTORY.md b/VERSION-HISTORY.md
index 7765d76..4d6fe8f 100644
--- a/VERSION-HISTORY.md
+++ b/VERSION-HISTORY.md
@@ -3,6 +3,7 @@ This file is a version history of turbo_seti amendments, beginning with version
 <br>
 |    Date    | Version | Contents |
 | :--: | :--: | :-- |
+| 2021-08-07 | 2.0.21 | New signal_processing source file, "dedoppler.py" (discussion in PR #220).  |
 | 2021-08-06 | 2.0.20 | New utility, "stix" (Issue #221).  |
 | 2021-07-28 | 2.0.19 | Update fil2h5 to handle YUGE data matrixes.  |
 | 2021-07-17 | 2.0.18 | Get rid of numpy "RuntimeWarning: Mean of empty slice" messages (Issue #212).  |
diff --git a/blimpy/signal_processing/__init__.py b/blimpy/signal_processing/__init__.py
new file mode 100644
index 0000000..0f449e1
--- /dev/null
+++ b/blimpy/signal_processing/__init__.py
@@ -0,0 +1,2 @@
+from .dedoppler import dedoppler_1
+
diff --git a/blimpy/signal_processing/dedoppler.py b/blimpy/signal_processing/dedoppler.py
new file mode 100644
index 0000000..d2e07a6
--- /dev/null
+++ b/blimpy/signal_processing/dedoppler.py
@@ -0,0 +1,31 @@
+r""" De-doppler signal processing functions """
+
+
+import numpy as np
+
+
+def dedoppler_1(wf, drift_rate):
+    """
+    Simple de-doppler code for a Filterbank or HDF5 file.
+    Parameters:
+    ----------
+    wf : object
+        Blimpy Waterfall object, previously instantiated with a loaded data matrix.
+    drift_rate : float
+        Signal drift rate over time [Hz/s]
+    """
+
+    # Get the time sampling interval in seconda.
+    tsamp = wf.header['tsamp']
+
+    # Get the fine channel bandwidth in Hz.
+    chan_bw = wf.header['foff'] * 1e6
+
+    # Compute the number of numpy rolls to perform.
+    n_roll = (drift_rate * tsamp) / chan_bw
+
+    # For each time-row,
+    #      roll all of the data power values in each fine channel frequency column
+    #      given by -(n_roll * row number).
+    for ii in range(wf.data.shape[0]):
+        wf.data[ii][0][:] = np.roll(wf.data[ii][0][:], -int(n_roll * ii))
diff --git a/requirements.txt b/requirements.txt
index 111c62b..661583f 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -10,3 +10,4 @@ hdf5plugin
 pandas
 Pillow
 psutil
+setigen
diff --git a/setup.py b/setup.py
index 5757d0d..86302fa 100644
--- a/setup.py
+++ b/setup.py
@@ -3,7 +3,7 @@
 """
 from setuptools import setup, find_packages
 
-__version__ = '2.0.20'
+__version__ = '2.0.21'
 
 with open("README.md", "r") as fh:
     long_description = fh.read()


Test patch:
diff --git a/tests/test_dedoppler.py b/tests/test_dedoppler.py
new file mode 100644
index 0000000..aed81ed
--- /dev/null
+++ b/tests/test_dedoppler.py
@@ -0,0 +1,100 @@
+from os.path import dirname
+import numpy as np
+from astropy import units as u
+import setigen as stg
+import matplotlib.pyplot as plt
+from blimpy.signal_processing.dedoppler import dedoppler_1
+from blimpy import Waterfall
+from blimpy.stax import plot_waterfall, sort2
+from tests.data import voyager_fil
+
+
+PLOT_DIR = dirname(voyager_fil)
+FIL_FILE = PLOT_DIR + "/test_dedoppler.fil"
+PNG_FILE = PLOT_DIR + "/test_dedoppler.png"
+
+
+# Plotting constants
+fontsize = 16
+font_dict = {"family" : "DejaVu Sans", "size" : fontsize}
+N_PLOTS = 6
+
+
+def plotter(counter, drift_rate):
+    wf = Waterfall(FIL_FILE)
+    dedoppler_1(wf, drift_rate)
+    wf.header["source_name"] = "Dedoppler at D.R. {} Hz".format(drift_rate)
+    plt.subplot(N_PLOTS, 1, counter)
+    return plot_waterfall(wf)
+
+
+def test_dedoppler_1():
+    
+    # Generate the Filterbank file.
+    print("test_dedoppler_1: Creating Filterbank file {}".format(FIL_FILE))
+    frame = stg.Frame(fchans=1024*u.pixel,
+                      tchans=32*u.pixel,
+                      df=2.7939677238464355*u.Hz,
+                      dt=18.253611008*u.s,
+                      fch1=6095.214842353016*u.MHz,
+                      ascending=True)
+    frame.add_noise(x_mean=10, noise_type='chi2')
+    frame.add_signal(stg.constant_path(f_start=frame.get_frequency(index=200),
+                                       drift_rate=2*u.Hz/u.s),
+                                       stg.constant_t_profile(level=frame.get_intensity(snr=30)),
+                                       stg.gaussian_f_profile(width=40*u.Hz),
+                                       stg.constant_bp_profile(level=1))    
+    frame.save_fil(FIL_FILE)
+    
+    # Load Filterban file.
+    print("test_dedoppler_1: Loading Filterbank file {}".format(FIL_FILE))
+    wf = Waterfall(FIL_FILE)
+    freqs = wf.get_freqs()
+    the_lowest, the_highest = sort2(freqs[0], freqs[-1])
+    the_midpoint = np.abs(the_lowest + the_highest) / 2.
+
+    # Initialise plotting.
+    print("test_dedoppler_1: Plotting to file {}".format(PNG_FILE))
+    fig = plt.subplots(N_PLOTS, sharex=True, sharey=True, figsize=(10, 2 * N_PLOTS))
+    wf.header["source_name"] = "Initial Data"
+
+    # Plot 1.
+    plt.subplot(N_PLOTS, 1, 1)
+    this_plot = plot_waterfall(wf)
+    
+    # Plot #2.
+    this_plot = plotter(2, 1.0)
+    
+    # Plot #3.
+    this_plot = plotter(3, 1.5)
+    
+    # Plot #4.
+    this_plot = plotter(4, 2.2)
+    
+    # Plot #5.
+    this_plot = plotter(5, 2.7)
+    
+    # Plot #6.
+    this_plot = plotter(6, 3.6)
+    
+    # Finish up plots.
+    plt.xticks(np.linspace(the_lowest, the_highest, num=4), ["","","",""])
+    factor = 1e6
+    units = "Hz"
+    xloc = np.linspace(the_lowest, the_highest, 5)
+    xticks = [round(loc_freq) for loc_freq in (xloc - the_midpoint) * factor]
+    if np.max(xticks) > 1000:
+        xticks = [xt / 1000 for xt in xticks]
+        units = "kHz"
+    plt.xticks(xloc, xticks)
+    plt.xlabel("Relative Frequency [%s] from %f MHz" % (units, the_midpoint), fontdict=font_dict)
+    cax = fig[0].add_axes([0.94, 0.11, 0.03, 0.77])
+    fig[0].colorbar(this_plot, cax=cax, label="Normalized Power (Arbitrary Units)")
+    plt.subplots_adjust(hspace=0, wspace=0)
+
+    plt.savefig(PNG_FILE, dpi=200, bbox_inches="tight")
+    print("test_dedoppler_1: End")
+
+
+if __name__ == "__main__":
+    test_dedoppler_1()