irfft: Restore (true) asymmetric apodization, move ramp into apodize()

Author: stuart-cls

orangecontrib/spectroscopy/irfft.py

@@ -56,15 +56,14 @@ def find_zpd(ifg, peak_search):
         raise NotImplementedError
 
 def wing_size(ifg, zpd):
-    """Calculate negative and positive wing size (including zpd point)"""
+    """Calculate negative and positive wing size (including zpd point in both)"""
     wing_n = zpd + 1
     wing_p = ifg.shape[-1] - zpd
     return wing_n, wing_p
 
-def apodize(ifg, zpd, apod_func):
+def apodize(ifg, zpd, apod_func, apod_asym):
     """
-    Perform apodization of asymmetric interferogram using selected apodization
-    function
+    Perform apodization of interferogram using selected apodization function
 
     Args:
         ifg (np.array): interferogram array (1D or 2D row-wise)
@@ -74,6 +73,8 @@ def apodize(ifg, zpd, apod_func):
                 <ApodFunc.BLACKMAN_HARRIS_3: 1> : Blackman-Harris (3-term)
                 <ApodFunc.BLACKMAN_HARRIS_4: 2> : Blackman-Harris (4-term)
                 <ApodFunc.BLACKMAN_NUTTALL: 3>  : Blackman-Nuttall (Eric Peach implementation)
+        apod_asym (bool): Apply asymmetric apodization for asymmetric interferometers (s-SNOM, DFTS)
+                          Otherwise, truncate symmetric apodization function as required and apply ramp
 
     Returns:
         ifg_apod (np.array): apodized interferogram(s)
@@ -87,8 +88,6 @@ def apodize(ifg, zpd, apod_func):
     # correcting zpd from 0-based index
     ifg_N = ifg.shape[-1]
     wing_n, wing_p = wing_size(ifg, zpd)
-    # Use larger wing to define apodization window width
-    M = max(wing_n, wing_p)
 
     # Generate apodization function
     if apod_func == ApodFunc.BLACKMAN_HARRIS_3:
@@ -118,17 +117,45 @@ def apodize(ifg, zpd, apod_func):
     else:
         raise ValueError(f"Invalid apodization function {apod_func}")
 
-    n = np.arange(1-M, M, 1)
-    Bs = A0 \
-         + A1 * np.cos(np.pi*n/(M-1)) \
-         + A2 * np.cos(np.pi*2*n/(M-1)) \
-         + A3 * np.cos(np.pi*3*n/(M-1))
-
-    # Select appropriate subsection of apodization function
-    if wing_n > wing_p:
-        Bs = Bs[:ifg_N]
+    if apod_asym:
+        # Scale window to each wing length (asymmetric interferometer)
+        n_n = np.arange(wing_n)
+        n_p = np.arange(wing_p)
+        Bs_n = A0 \
+               + A1 * np.cos(np.pi*n_n/wing_n) \
+               + A2 * np.cos(np.pi*2*n_n/wing_n) \
+               + A3 * np.cos(np.pi*3*n_n/wing_n)
+        Bs_p = A0 \
+               + A1 * np.cos(np.pi*n_p/wing_p) \
+               + A2 * np.cos(np.pi*2*n_p/wing_p) \
+               + A3 * np.cos(np.pi*3*n_p/wing_p)
+        # Flip Bs_n and overlap at zpd (Bs_n[0] == Bs_p[0])
+        Bs = np.hstack((Bs_n[::-1], Bs_p[1:]))
     else:
-        Bs = Bs[2 * M - 1 - ifg_N:]
+        # Use larger wing to define apodization window width (symmetric interferometer)
+        M = max(wing_n, wing_p)
+        n = np.arange(1-M, M, 1)
+        Bs = A0 \
+             + A1 * np.cos(np.pi*n/(M-1)) \
+             + A2 * np.cos(np.pi*2*n/(M-1)) \
+             + A3 * np.cos(np.pi*3*n/(M-1))
+
+        # Select appropriate subsection of apodization function
+        if wing_n > wing_p:
+            Bs = Bs[:ifg_N]
+        else:
+            Bs = Bs[2 * M - 1 - ifg_N:]
+
+        # Add ramp to avoid uneven counting of truncated section
+        # Use smaller wing to define symmetric subset size
+        # Skip if fully symmetric (phase data)
+        wing_min = min(wing_n, wing_p)
+        ramp = np.linspace(0, 1, 2 * wing_min - 1)
+        if wing_n < wing_p:
+            Bs[..., :ramp.shape[0]] *= ramp
+        elif wing_n > wing_p:
+            Bs[..., -ramp.shape[0]:] *= ramp[::-1]
+
     # Apodize the sampled Interferogram
     try:
         ifg_apod = ifg * Bs
@@ -231,14 +258,13 @@ def __call__(self, ifg, zpd=None, phase=None):
         else:
             # Select phase interferogram as copy
             # Note that L is now the zpd index
-            Ixs = ifg[self.zpd - L : self.zpd + L].copy()
-            ifg = apodize(ifg, self.zpd, self.apod_func)
-            ifg = ramp(ifg, self.zpd)
+            Ixs = ifg[self.zpd - L : self.zpd + 1 + L].copy()
+            ifg = apodize(ifg, self.zpd, self.apod_func, self.apod_asym)
             ifg = zero_fill(ifg, self.zff)
             ifg = np.hstack((ifg[self.zpd:], ifg[0:self.zpd]))
             Nzff = ifg.shape[0]
 
-            Ixs = apodize(Ixs, L, self.apod_func)
+            Ixs = apodize(Ixs, L, self.apod_func, self.apod_asym)
             # Zero-fill Ixs to same size as ifg (instead of interpolating later)
             Ixs = _zero_fill_pad(Ixs, Nzff - Ixs.shape[0])
             Ixs = np.hstack((Ixs[L:], Ixs[0:L]))
@@ -306,7 +332,7 @@ def __call__(self, ifg, zpd=None, phase=None):
         # Calculate phase on interferogram of specified size 2*L
         L = self.phase_ifg_size(ifg.shape[1])
         if L == 0: # Use full ifg for phase #TODO multi is this code tested
-            ifg = apodize(ifg, self.zpd, self.apod_func)
+            ifg = apodize(ifg, self.zpd, self.apod_func, self.apod_asym)
             ifg = zero_fill(ifg, self.zff)
             # Rotate the Complete IFG so that the centerburst is at edges.
             ifg = np.hstack((ifg[self.zpd:], ifg[0:self.zpd]))
@@ -318,13 +344,12 @@ def __call__(self, ifg, zpd=None, phase=None):
             # Select phase interferogram as copy
             # Note that L is now the zpd index
             Ixs = ifg[:, self.zpd - L : self.zpd + L].copy()
-            ifg = apodize(ifg, self.zpd, self.apod_func)
-            ifg = ramp(ifg, self.zpd)
+            ifg = apodize(ifg, self.zpd, self.apod_func, self.apod_asym)
             ifg = zero_fill(ifg, self.zff)
             ifg = np.hstack((ifg[:, self.zpd:], ifg[:, 0:self.zpd]))
             Nzff = ifg.shape[1]
 
-            Ixs = apodize(Ixs, L, self.apod_func)
+            Ixs = apodize(Ixs, L, self.apod_func, self.apod_asym)
             # Zero-fill Ixs to same size as ifg (instead of interpolating later)
             Ixs = _zero_fill_pad(Ixs, Nzff - Ixs.shape[1])
             Ixs = np.hstack((Ixs[:, L:], Ixs[:, 0:L]))
@@ -357,7 +382,7 @@ def __call__(self, ifg, zpd=None, phase=None):
         else:
             self.zpd = find_zpd(ifg, self.peak_search)
 
-        ifg = apodize(ifg, self.zpd, self.apod_func)
+        ifg = apodize(ifg, self.zpd, self.apod_func, self.apod_asym)
         ifg = zero_fill(ifg, self.zff)
         # Rotate the Complete IFG so that the centerburst is at edges.
         ifg = np.hstack((ifg[self.zpd:], ifg[0:self.zpd]))

orangecontrib/spectroscopy/tests/test_irfft.py

@@ -7,7 +7,7 @@
 
 from orangecontrib.spectroscopy.irfft import (IRFFT, zero_fill, PhaseCorrection,
                                               find_zpd, PeakSearch, ApodFunc,
-                                              MultiIRFFT, apodize, ramp,
+                                              MultiIRFFT, apodize,
                                               )
 
 dx = 1.0 / 15797.337544 / 2.0
@@ -26,6 +26,7 @@ def setUp(self):
             'odd_trunc': self.ifg_seq_ref.X[0],
             'even_trunc_reverse': self.ifg_seq_ref.X[0][1:][::-1],
             'odd_trunc_reverse': self.ifg_seq_ref.X[0][::-1],
+            'fake_asym': self.ifg_seq_ref.X[0],
         }
 
     def test_zero_fill(self):
@@ -97,6 +98,7 @@ def test_agilent_fft_ab(self):
         dx_ag = (1 / 1.57980039e+04 / 2) * 4
         fft = IRFFT(dx=dx_ag,
                     apod_func=ApodFunc.BLACKMAN_HARRIS_4,
+                    apod_asym=False,
                     zff=1,
                     phase_res=None,
                     phase_corr=PhaseCorrection.MERTZ,
@@ -133,6 +135,7 @@ def test_multi_ab(self):
         dx_ag = (1 / 1.57980039e+04 / 2) * 4
         fft = MultiIRFFT(dx=dx_ag,
                          apod_func=ApodFunc.BLACKMAN_HARRIS_4,
+                         apod_asym=False,
                          zff=1,
                          phase_res=None,
                          phase_corr=PhaseCorrection.MERTZ,
@@ -153,18 +156,24 @@ def test_multi_ab(self):
 
     def test_apodization(self):
         for apod_func in ApodFunc:
+            if apod_func == ApodFunc.BOXCAR:
+                ifg = self.ifgs['even_sym']
+                zpd = find_zpd(ifg, PeakSearch.ABSOLUTE)
+                out = apodize(ifg, zpd, apod_func, False)
+                np.testing.assert_array_equal(ifg, out)
+                continue
             for k, ifg in self.ifgs.items():
-                with self.subTest(apod_func=apod_func.name, ifg=k):
+                if 'asym' in k:
+                    # Test asymmetric ifgs
+                    apod_asym = True
+                    ramp_value = 1
+                else:
+                    # Test symmetric and truncated ifgs
+                    apod_asym = False
+                    ramp_value = 0.5
+                with self.subTest(apod_func=apod_func.name, ifg=k, asym=apod_asym):
                     zpd = find_zpd(ifg, PeakSearch.ABSOLUTE)
-                    out = apodize(ifg, zpd, apod_func)
+                    out = apodize(ifg, zpd, apod_func, apod_asym)
                     # Apodization should not change value at zpd
-                    self.assertAlmostEqual(ifg[zpd], out[zpd])
-
-    def test_ramp(self):
-        ifg = self.ifgs['odd_trunc']
-        zpd = find_zpd(ifg, PeakSearch.ABSOLUTE)
-        ramp_fwd = ramp(ifg, zpd)
-        ifg_rev = self.ifgs['odd_trunc_reverse']
-        zpd_rev = find_zpd(ifg_rev, PeakSearch.ABSOLUTE)
-        ramp_rev = ramp(ifg_rev, zpd_rev)
-        np.testing.assert_array_equal(ramp_fwd, ramp_rev[::-1])
+                    # Ramp should reduce zpd value by half
+                    self.assertAlmostEqual(ifg[zpd] * ramp_value, out[zpd])