Improve description of upsample(). Add verify argument.

applications/lofar2/model/pfb_os/dsp.py

@@ -169,7 +169,8 @@ def ideal_low_pass_filter(Npoints, Npass, bandEdgeGain=1.0):
     return h, f, HF
 
 
-def prototype_fir_low_pass_filter(Npoints=1024, Ntaps=16, Ncoefs=1024*16, hpFactor=0.9, transitionFactor=0.4, stopRippleFactor=1000000, beta=1):
+def prototype_fir_low_pass_filter(
+        Npoints=1024, Ntaps=16, Ncoefs=1024*16, hpFactor=0.9, transitionFactor=0.4, stopRippleFactor=1000000, beta=1):
     """Derive FIR coefficients for prototype low pass filter using firls
 
     Default use LPF specification for LOFAR subband filter. For subband filter
@@ -536,28 +537,45 @@ class PolyPhaseFirFilterStructure:
         return outData
 
 
-def upsample(x, Nup, coefs):  # interpolate
+def upsample(x, Nup, coefs, verify=False):  # interpolate
     """Upsample x by factor I = Nup
 
+    Input:
+    . x: Input signal x[n]
+    . Nup: upsample (interpolation) factor
+    . coefs: FIR filter coefficients for antialiasing LPF
+    . verify: when True then verify that output y is the same when calculated directly or when calculatged using the
+              polyphase implementation.
+    Return:
+    . y: Upsampled output signal y[m]
+
     Assumptions:
     . x[n] = 0 for n < 0
     . no y[m] for m < 0
-    . insert I - 1 zeros after each x, so len(y) = I * len(x)
+    . insert I - 1 zeros after each x[n] to get v[m], so len(v) = len(y) = I * len(x)
     . use coefs as anti aliasing filter, must be LPF with BW < fNyquist / I
-    . len(coefs) should be multiple of Nup
+    . len(coefs) typically is multiple of Nup. If shorter, then the coefs are extended with zeros.
+
+    Remarks:
+    . The input sample period is ts and the output sample period of the upsampled (= interpolated signal) is tsUp =
+      ts / I
+    . The group delay is (Ncoefs - 1) / 2 * tsUp. With odd Ncoefs and symmetrical coefs to have linear phase, the
+      group delay is an integer number of tsUp periods.
+
+    Procedure:
 
-    x[n]:   0           1           2           3  ...
+      x[n]:   0           1           2           3 --> time n with unit Ts
 
       v[m] = x[m / I], for m = 0, +-I, +-2I, ...
            = 0, else
-    v[m]:   0  1  2  3  4  5  6  7  8  9 10 11 12  ...
+      v[m]:   0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 --> time m with unit Ts / I
+                                                  |  |  |  |
+      h[k]:     11 10  9  8  7  6  5  4  3  2  1  0  |  |  | --> coefs for m = 12
+                   11 10  9  8  7  6  5  4  3  2  1  0  |  | --> coefs for m = 13
+                      11 10  9  8  7  6  5  4  3  2  1  0  | --> coefs for m = 14
+                         11 10  9  8  7  6  5  4  3  2  1  0 --> coefs for m = 15
 
-    h[k]:     11 10  9  8  7  6  5  4  3  2  1  0  coefs
-                 11 10  9  8  7  6  5  4  3  2  1  0
-                    11 10  9  8  7  6  5  4  3  2  1  0
-                          11 10  9  8  7  6  5  4  3  2  1  0
-
-    y[m]:   0  1  2  3  4  5  6  7  8  9 10 11 12  ...
+      y[m]:   0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 --> time m with unit Ts / I
 
             Calucate y[0, 1, 2, 3] at x[0]
               y[0] = h[0] x[0]
@@ -577,7 +595,7 @@ def upsample(x, Nup, coefs):  # interpolate
               y[10] = h[2] x[2] + h[6] x[1] + h[10] x[0]
               y[11] = h[3] x[2] + h[7] x[1] + h[11] x[0]
 
-          Calucate y[12, 13, 14, 15] at x[3]
+            Calucate y[12, 13, 14, 15] at x[3], see '|' markers between v[m] (is zero padded x[n]) and h[k] above
               y[12] = h[0] x[3] + h[4] x[2] + h[ 8] x[1]
               y[13] = h[1] x[3] + h[5] x[2] + h[ 9] x[1]
               y[14] = h[2] x[3] + h[6] x[2] + h[10] x[1]
@@ -591,22 +609,28 @@ def upsample(x, Nup, coefs):  # interpolate
               y[n * I + 3] = lfilter(h[3, 7,11], [1], x)
     """
     Nx = len(x)
+    Ncoefs = len(coefs)
     a = [1.0]
-    if False:
-        # Inefficient implementation with multiply by zero values
-        xZeros = np.zeros((Nup, Nx))
-        xZeros[0] = x
-        xZeros = xZeros.T.reshape(1, Nup * Nx)[0]
-        y = signal.lfilter(coefs, a, xZeros)
-    else:
+
     # Polyphase implementation to avoid multiply by zero values
-        Ntaps = len(coefs) // Nup
+    Ntaps = Ncoefs // Nup
     polyCoefs = coefs.reshape(Ntaps, Nup).T
     polyY = np.zeros((Nup, Nx))
     # Filter x per polyphase
     for p in range(Nup):
         polyY[p] = signal.lfilter(polyCoefs[p], a, x)
     y = polyY.T.reshape(1, Nup * Nx)[0]
+
+    if verify:
+        # Inefficient implementation with multiply by zero values
+        xZeros = np.zeros((Nup, Nx))
+        xZeros[0] = x
+        xZeros = xZeros.T.reshape(1, Nup * Nx)[0]
+        yVerify = signal.lfilter(coefs, a, xZeros)
+        if np.all(y == yVerify):
+            print('ERROR: wrong upsample result')
+        else:
+            print('PASSED: correct upsample result')
     return y