Remove structure bunton, instead integrate it as variant that is always verified.

applications/lofar2/model/rtdsp/dftfilterbank.py

@@ -38,7 +38,6 @@
 
 import numpy as np
 from sys import exit
-from .utilities import ceil_div
 from .multirate import fold, \
                        PolyPhaseFirFilterStructure, \
                        polyphase_data_for_downsampling_whole_x
@@ -70,6 +69,7 @@ def analysis_dft_filterbank(x, Ndown, Ndft, coefs, structure, commutator=None, v
       theorem [CROCHIERE Fig 7.21]:
         x[n] <= DFT => X(k), then
         x[(n - l) % K] <= DFT => X(k) exp(-j 2pi k l / K)
+    - Note: fold = roll -1 then flip = flip then roll +1
 
     Input:
     . x: Input signal x[n]
@@ -85,72 +85,37 @@ def analysis_dft_filterbank(x, Ndown, Ndft, coefs, structure, commutator=None, v
     . Yc: Downsampled and downconverted output signal Yc[k, m], m = n // M for
           all Ndft bins k. Complex baseband signals.
     """
-    Ros = Ndft / Ndown
-
-    # Prepend x with Ndown - 1 zeros, and represent x in Nblocks of Ndown
-    # samples
-    Nzeros = Ndown - 1
-    xBlocks, xData, Nx, Nblocks = polyphase_data_for_downsampling_whole_x(x, Ndown, Nzeros)
-
-    # Prepare output
-    Yc = np.zeros((Ndft, Nblocks), dtype='cfloat')
 
     # Oversampling analysis DFT filterbank
-    # . note: fold = roll -1 then flip = flip then roll +1
     if structure == 'wola':
-        # >>> Use WOLA structure
-        # Use PFS with Ndft polyphases and shift in xData in overlapping
-        # blocks of Ncoefs samples, to apply coefs as window. The coefs in
-        # the pfs are already in flipped order conform [CROCHIERE Eq 7.70].
+        # >>> Use WOLA structure [CROCHIERE Fig 7.19]
+        # Apply the coefs as a window conform [CROCHIERE Eq 7.70].
+        # For the pfs the coefs are already in the correct order
         pfs = PolyPhaseFirFilterStructure(Ndft, coefs)
-        # print(pfs.coefs)
-        # print(pfs.polyCoefs)
+        # For the delay line the coefs still need to be flipped. This flip is
+        # omitted in reconstruct.m
+        wCoefs = np.flip(coefs)
 
-        # Prepend xData with pfs.Ncoefs - Ndown zeros to align with 'pfs'
-        # implementation
-        xData = np.concatenate((np.zeros(pfs.Ncoefs - Ndown), xData))
+        # Get parameters
+        Ncoefs = pfs.Ncoefs
+        Ntaps = pfs.Ntaps
+        Ndelays = pfs.Ndelays
 
-        # PFB loop per Ndown input xData samples
-        for b in range(Nblocks):
-            # Window
-            tRange = np.arange(pfs.Ncoefs) + b * Ndown
-            tData = xData[tRange]
-            pfs.shift_in_data(np.flip(tData))
-            if b < 3:
-                # print(b)
-                # print(tData)
-                # print(pfs.polyDelays)
-                pass
-            zPoly = pfs.polyDelays * pfs.polyCoefs
-            # Sum pfs.Ntaps number of blocks
-            pfsData = np.sum(zPoly, axis=1)
-            # Time shift
-            tShift = b * Ndown  # roll is modulo Ndft
-            # DFT
-            pfsShifted = np.roll(pfsData, -tShift)
-            # fold = roll -1 then flip
-            pfsShifted = fold(pfsShifted)
-            Yc[:, b] = np.fft.fft(pfsShifted)
-            if b < 3:
-                # print(b, pfs.map_to_delay_line(pfs.polyDelays))
-                # print(b, np.flip(pfs.map_to_delay_line(zPoly)))
-                # print(b, np.flip(pfsData))
-                # print(b, Yc[:, b])
-                pass
-
-    elif structure == 'bunton':
-        # >>> Use WOLA structure as in reconstruct.m by John Bunton
-        Ncoefs = len(coefs)
-        Ntaps = ceil_div(Ncoefs, Ndft)  # taps, columns
-        Ndelays = Ntaps * Ndft  # zero padded coefs
-        # Prepend xData with pfs.Ncoefs - Ndown zeros to align with 'pfs'
-        # implementation
-        xData = np.concatenate((np.zeros(Ncoefs - Ndown), xData))
-        # Use coefs as window
-        wCoefs = coefs
-        wCoefs = np.flip(coefs)
-        # Time shift offset, due to roll -1 in fold() ?, to align with LO in
-        # single channel reference.
+        # Determine Nblocks of Ndown samples in x
+        _, _, _, Nblocks = polyphase_data_for_downsampling_whole_x(x, Ndown, Nzeros=0)
+
+        # Prepend x with Ncoefs - 1 zeros to have first down sampled sample at
+        # m = 0 start at n = 0, and to fit Nblocks for WOLA
+        Nzeros = pfs.Ncoefs - 1
+        xData = np.concatenate((np.zeros(Nzeros), x))
+
+        # Prepare output
+        Yc = np.zeros((Ndft, Nblocks), dtype='cfloat')
+
+        # Need time shift offset to align with LO in single channel reference.
+        # This tOffset = 1 together with the flip of pfsData makes a fold(),
+        # so that the use of DFT in WOLA analysis is equivalent to using
+        # IDFT as in LO downconverter and PFS analysis.
         tOffset = 1
         # PFB loop per Ndown input xData samples
         for b in range(Nblocks):
@@ -158,22 +123,56 @@ def analysis_dft_filterbank(x, Ndown, Ndft, coefs, structure, commutator=None, v
             tShift = b * Ndown
             # Apply window
             tRange = np.arange(Ncoefs) + tShift
-            tData = xData[tRange]
-            zLine = np.zeros(Ndelays, dtype='float')
-            zLine[0:Ncoefs] = tData * wCoefs
+            tData = xData[tRange]  # data for whole window
+
+            # 1) Variant using PFS
+            # Shifting in the (old and new) data for the whole window is
+            # equivalent to shifting in only the new data
+            if 0:
+                # load whole window data
+                pfs.shift_in_data(np.flip(tData))
+            else:
+                # shift in new data block
+                tBlock = tData[Ncoefs - Ndown:]
+                pfs.shift_in_data(np.flip(tBlock))
+            zPoly = pfs.polyDelays * pfs.polyCoefs
             # Sum Ntaps
-            sumLine = np.zeros(Ndft, dtype='cfloat')
+            pfsData = np.sum(zPoly, axis=1)
+            # Flip polyphases 0:Ndft to match time order for DFT input
+            pfsData = np.flip(pfsData)
+
+            # 2) Variant using a delay line like in reconstruct.m
+            dLine = np.zeros(Ndelays, dtype='float')
+            # Need to load data at end of delay line to sum the taps correctly
+            # in case Ncoefs < Ndelays
+            dLine[Ndelays - Ncoefs : Ndelays] = tData * wCoefs
+            # Sum Ntaps
+            sumLine = np.zeros(Ndft, dtype='float')
             for t in range(Ntaps):
-                sumLine += zLine[np.arange(Ndft) + t * Ndft]
+                sumLine += dLine[np.arange(Ndft) + t * Ndft]
+
+            # Verify that 1) PFS and 2) delay line yield same result
+            # . The two variants differ in structure, but yield same result.
+            if not np.allclose(sumLine, pfsData):
+                exit('ERROR: wrong analysis WOLA')
+
             # Fixed time reference, roll is modulo Ndft
-            sumShifted = np.roll(sumLine, tOffset + tShift)
+            pfsShifted = np.roll(pfsData, tOffset + tShift)
             # DFT
-            Yc[:, b] = np.fft.fft(sumShifted)
+            Yc[:, b] = np.fft.fft(pfsShifted)
 
     elif structure == 'pfs':
-        # >>> Use PFS
         # PFS with Ndft polyphases and shift in xBlocks of Ndown samples
         pfs = PolyPhaseFirFilterStructure(Ndft, coefs)
+
+        # Prepend x with Ndown - 1 zeros, and represent x in Nblocks of Ndown
+        # samples
+        Nzeros = Ndown - 1
+        xBlocks, xData, Nx, Nblocks = polyphase_data_for_downsampling_whole_x(x, Ndown, Nzeros)
+
+        # Prepare output
+        Yc = np.zeros((Ndft, Nblocks), dtype='cfloat')
+
         # PFB loop per Ndown input xData samples
         for b in range(Nblocks):
             # Filter block, the inData blocks are already flipped in
@@ -183,8 +182,8 @@ def analysis_dft_filterbank(x, Ndown, Ndft, coefs, structure, commutator=None, v
 
             # Time (m) and bin (k) dependend phase shift W_K^(kmM) by circular
             # time shift of DFT input
-            tShift = b * Ndown  # roll is modulo Ndft
-            pfsShifted = np.roll(pfsData, -tShift)
+            tShift = b * Ndown
+            pfsShifted = np.roll(pfsData, -tShift)  # roll is modulo Ndft
 
             # For 'ccw' apply IDFT, for 'cw' apply DFT(fold()), due to fold()
             # both yield identical result
@@ -202,9 +201,9 @@ def analysis_dft_filterbank(x, Ndown, Ndft, coefs, structure, commutator=None, v
         exit('ERROR: invalid structure ' + str(structure))
 
     if verbosity:
+        Ros = Ndft / Ndown
         print('> analysis_dft_filterbank():')
         print('  . len(x)     =', str(len(x)))
-        print('  . Nx         =', str(Nx))
         print('  . Nblocks    =', str(Nblocks))
         print('  . Ros        =', str(Ros))
         print('  . Ndown      =', str(Ndown))
@@ -230,8 +229,8 @@ def synthesis_dft_filterbank(Xbase, Nup, Ndft, coefs, structure, commutator=None
     . Ndft: DFT size, number of polyphases in PFS FIR filter
     . coefs: prototype LPF FIR filter coefficients for anti aliasing and
       interpolating BPF
-    . structure: 'wola' or 'pfs'
-    . commutator: 'ccw' to use DFT or 'cw' to use IDFT
+    . structure: 'wola'
+    . commutator: 'cw' to use DFT or 'ccw' to use IDFT
     - verbosity: when > 0 print() status, else no print()
     Return:
     . y: Upsampled and upconverted output signal y[n].
@@ -248,14 +247,24 @@ def synthesis_dft_filterbank(Xbase, Nup, Ndft, coefs, structure, commutator=None
     # Adjust LO phase for group delay of analysis LPF and synthesis LPF in
     # series. This is not needed in practise, but is needed to be able to
     # exactly compare reconstructed signal with original input time series
-    # signal. Assume coefs has group delay (len(coefs) - 1) / 2.
+    # signal. Assume analysis coefs and synthesis coefs have linear phase and
+    # same length (but not necessarily same coefs), then the total group
+    # delay is (len(coefs) - 1) / 2 + (len(coefs) - 1) / 2.
     hPairGroupDelay = len(coefs) - 1
 
-    # Oversampling analysis DFT filterbank
+    # Oversampling synthesis DFT filterbank
     if structure == 'wola':
-        # >>> Use WOLA structure
-        # Use PFS with Ndft polyphases for delay line and window
+        # >>> Use WOLA structure [CROCHIERE Fig 7.20]
+        # Apply the coefs as a window conform [CROCHIERE Eq 7.76].
+        # Use PFS with Ndft polyphases for delay line and coefs window
         pfs = PolyPhaseFirFilterStructure(Ndft, coefs)
+        # Use coefs as window
+        wCoefs = coefs
+
+        # Get parameters
+        Ncoefs = pfs.Ncoefs
+        Ntaps = pfs.Ntaps
+        Ndelays = pfs.Ndelays  # zero padded coefs
 
         # PFB loop per Nup output y samples
         for b in range(Nblocks):
@@ -270,51 +279,37 @@ def synthesis_dft_filterbank(Xbase, Nup, Ndft, coefs, structure, commutator=None
 
             # Apply filter group delay (hPairGroupDelay), and time (m) and bin
             # (k) dependend phase shift W_K^(kmM) by circular time shift of
-            # DFT output
-            tShift = b * Nup - hPairGroupDelay  # roll is modulo Ndft
-            xShifted = np.roll(xTime, -tShift)
+            # DFT output. To get from fixed time reference to sliding time
+            # reference.
+            tShift = b * Nup - hPairGroupDelay
+            xShifted = np.roll(xTime, -tShift)  # roll is modulo Ndft
 
-            # Load xTime at each tap in polyDelays
+            # 1) Variant using PFS.
+            # Load xTime at Ntaps in polyDelays
             pfs.parallel_load(xShifted)
-
             # Apply FIR coefs per delay element
             zPoly = Nup * pfs.polyDelays * pfs.polyCoefs
             zLine = pfs.map_to_delay_line(zPoly)
+            zLine = zLine[0:Ncoefs]
 
-            # Overlap add weigthed input to the output
-            tRange = np.arange(pfs.Ndelays) + b * Nup
-            y[tRange] += zLine
-            if b < 3:
-                # print(b, Xbase[:, b])
-                # print(b, xTime / Ndft)
-                # print(b, zLine / Ndft  / Nup)
-                # print(b, y[tRange] / Ndft  / Nup)
-                pass
-
-    elif structure == 'bunton':
-        # >>> Use WOLA structure as in reconstruct.m by John Bunton
-        Ncoefs = len(coefs)
-        Ntaps = ceil_div(Ncoefs, Ndft)  # taps, columns
-        Ndelays = Ntaps * Ndft  # zero padded coefs
-        # Use coefs as window
-        wCoefs = coefs
-        # wCoefs = np.flip(coefs)
-        for b in range(Nblocks):
-            # IDFT
-            xTime = Ndft * np.real(np.fft.ifft(Xbase[:, b]))
-            # Sliding time reference
-            tShift = b * Nup - hPairGroupDelay  # roll is modulo Ndft
-            xShifted = np.roll(xTime, -tShift)
+            # 2) Variant using a delay line like in reconstruct.m
             # Copy data at Ntaps
-            zLine = np.zeros(Ndelays, dtype='float')
+            dLine = np.zeros(Ndelays, dtype='float')
             for t in range(Ntaps):
-                zLine[np.arange(Ndft) + t * Ndft] = xShifted
+                dLine[np.arange(Ndft) + t * Ndft] = xShifted
             # Apply window
             yLine = np.zeros(Ncoefs, dtype='float')
-            yLine = Nup * zLine[0:Ncoefs] * wCoefs
+            yLine = Nup * dLine[0:Ncoefs] * wCoefs
+
+            # Verify that 1) PFS and 2) delay line yield same result
+            # . Both structures are in fact identical, because both use
+            #   parallel load of the xShifted data.
+            if not np.allclose(zLine, yLine):
+                exit('ERROR: wrong sythesis WOLA')
+
             # Overlap add weigthed input to the output
-            tRange = np.arange(Ncoefs) + b * Nup
-            y[tRange] += yLine
+            tRange = np.arange(pfs.Ncoefs) + b * Nup
+            y[tRange] += zLine
 
     else:
         # No 'pfs' for fractional oversampled synthesis, because the