Merge pull request #49 from hfmark/multitap

update mtspec -> multitaper

.github/workflows/tests.yml

@@ -46,11 +46,11 @@ jobs:
         with:
           python-version: ${{ matrix.python }}
           environment-file: .github/test_conda_env.yml
-      - name: install mtspec
+      - name: install multitaper
         if: ${{ matrix.add_package == 'features' }}
         continue-on-error: ${{ matrix.continue-on-error }}
         run: |
-          conda install mtspec
+          conda install multitaper
       - name: install other optional dependencies
         if: ${{ matrix.add_package == 'features' }}
         continue-on-error: ${{ matrix.continue-on-error }}

rf/__init__.py

@@ -55,7 +55,7 @@
 
     * ObsPy_ and some of its dependencies,
     * cartopy, geographiclib, shapely,
-    * mtspec_ for multitaper deconvolution,
+    * multitaper_ for multitaper deconvolution,
     * toeplitz_ for faster time domain deconvolution (optional),
     * obspyh5_ for hdf5 file support (optional),
     * tqdm for progress bar in batch processing (optional).
@@ -75,7 +75,7 @@
 Here are some instructions to install rf into a fresh conda environment::
 
     conda config --add channels conda-forge
-    conda create -n rfenv obspy cartopy geographiclib shapely h5py mtspec tqdm fortran-compiler
+    conda create -n rfenv obspy cartopy geographiclib shapely h5py multitaper tqdm fortran-compiler
     conda activate rfenv
     pip install obspyh5 toeplitz rf
     rf-runtests
@@ -259,7 +259,7 @@
 .. _pip: http://www.pip-installer.org/
 .. _obspyh5: https://github.com/trichter/obspyh5/
 .. _toeplitz: https://github.com/trichter/toeplitz/
-.. _mtspec: https://github.com/krischer/mtspec
+.. _multitaper: https://github.com/gaprieto/multitaper
 
 .. _notebook1: http://nbviewer.jupyter.org/github/trichter/notebooks/blob/master/receiver_function_minimal_example.ipynb
 .. _notebook2: http://nbviewer.jupyter.org/github/trichter/notebooks/blob/master/receiver_function_profile_chile.ipynb

rf/deconvolve.py

@@ -558,9 +558,9 @@ def deconv_multitaper(rsp, src, nse, sampling_rate, tshift, gauss=0.5,
     :return: (list of) array(s) with deconvolution(s)
     """
     try:
-        import mtspec as mt
+        import multitaper as mt
     except ImportError as ex:
-        msg = 'mtspec package is needed for multitaper deconvolution'
+        msg = 'multitaper package is needed for multitaper deconvolution'
         raise ImportError(msg) from ex
 
     # check src trace length < rsp trace length:
@@ -570,10 +570,10 @@ def deconv_multitaper(rsp, src, nse, sampling_rate, tshift, gauss=0.5,
     nft = len(rsp[0])  # length of final arrays
     dt = 1./sampling_rate  # sample spacing
 
-    # calculate multitapers with mtspec
+    # calculate multitapers with German's multitaper package
     ntap = int(round(T/dt))  # #points in each taper
-
-    tap, el, _ = mt.multitaper.dpss(ntap, tband, K); tap = tap.T
+    slepians = mt.MTSpec(np.arange(ntap),nw=tband,kspec=K)
+    tap = slepians.vn.T; el = slepians.lamb
     if K >= 2:
         tap[1] = -tap[1]  # adjust to match sign convention
     if K >= 3:
@@ -583,7 +583,7 @@ def deconv_multitaper(rsp, src, nse, sampling_rate, tshift, gauss=0.5,
 
     ofac = 1 / (1 - olap)  # calculate overlap factor
 
-    sfft = np.zeros((K,nft), dtype=np.complex)  # array for holding freq-domain source estimate
+    sfft = np.zeros((K,nft), dtype=complex)  # array for holding freq-domain source estimate
     src = detrend(src, type='constant')  # demean and detrend source
     src = detrend(src, type='linear')
 
@@ -612,8 +612,8 @@ def deconv_multitaper(rsp, src, nse, sampling_rate, tshift, gauss=0.5,
     RF_out = np.zeros((ncomp, nft))
     for c in range(ncomp):
         dat = rsp[c]; nos = nse[c]  # pick out component
-        dfft = np.zeros((K,nft), dtype=np.complex)  # arrays for storing freq-domain estimates
-        nfft = np.zeros((K,nft), dtype=np.complex)
+        dfft = np.zeros((K,nft), dtype=complex)  # arrays for storing freq-domain estimates
+        nfft = np.zeros((K,nft), dtype=complex)
 
         # window, taper, and transform the noise
         nos = detrend(nos, type='constant')

setup.py

@@ -30,7 +30,7 @@ def find_version(*paths):
 
 EXTRAS_REQUIRE = {
     'doc': ['sphinx', 'alabaster'],  # and decorator, obspy
-    'deconv_multitaper': ['mtspec'],
+    'deconv_multitaper': ['multitaper'],
     'h5': ['obspyh5>=0.3'],
     'toeplitz': ['toeplitz'],
     'batch': ['tqdm']