Add unit_circle_loops_phasor_arr(), time_shift_phasor_arr().

applications/lofar2/model/rtdsp/multirate.py

@@ -612,13 +612,35 @@ def resample(x, Nup, Ndown, coefs, verify=False, verbosity=1):  # interpolate an
 # Single bandpass channel up and down sampling and up and down conversion
 ###############################################################################
 
-def phasor_arr(k, Ndft, sign):
-    """Return array of phasors: exp(+-j 2pi k / Ndft) for k in 0 : Ndft - 1
+def unit_circle_loops_phasor_arr(k, N, sign):
+    """Return array of N phasors on k loops along the unit circle.
+
+    Polyphase dependent phase offsets for bin k, when N = Ndown = Ndft [HARRIS
+    Eq 6.8]. For k = 1 this yields the roots of unity.
+
+    Return:
+    . pArr: exp(+-j 2pi k / N) for k in 0 : N - 1
     """
     if sign == 'positive':
-        return np.array([np.exp(2j * np.pi * p * k / Ndft) for p in range(Ndft)])
+        pArr = np.array([np.exp(2j * np.pi * p * k / N) for p in range(N)])
     else:  # 'negative'
-        return np.array([np.exp(-2j * np.pi * p * k / Ndft) for p in range(Ndft)])
+        pArr = np.array([np.exp(-2j * np.pi * p * k / N) for p in range(N)])
+    return pArr
+
+
+def time_shift_phasor_arr(k, Ndown, Ndft, Msamples):
+    """Return array of Msamples phasors in time to compensate for oversampling
+    time shift.
+
+    The time shift due to downsampling causes a frequency component of
+    k * Ndown / Ndft. With oversampling Ndown < Ndft, and then after
+    downsampling there remains a frequency offset [HARRIS Eq 9.3].
+
+    Return:
+    . mArr = exp(2j pi k * Ndown / Ndft * m),  for m in 0 : Msamples - 1
+    """
+    mArr = np.exp(-2j * np.pi * k * Ndown / Ndft * np.arange(Msamples))
+    return mArr
 
 
 def maximal_downsample_bpf(x, Ndown, k, coefs, verbosity=1):
@@ -649,9 +671,9 @@ def maximal_downsample_bpf(x, Ndown, k, coefs, verbosity=1):
     # Phase rotate per polyphase for bin k, due to delay line at branch inputs
     # [HARRIS Eq 6.8]
     polyYC = np.zeros((Ndown, Nxp), dtype='cfloat')
-    phasors = phasor_arr(k, Ndown, 'positive')
+    kPhasors = unit_circle_loops_phasor_arr(k, Ndown, 'positive')
     for p in range(Ndown):
-        polyYC[p] = polyY[p] * phasors[p]  # row = row * scalar
+        polyYC[p] = polyY[p] * kPhasors[p]  # row = row * scalar
 
     # Sum the branch outputs to get single downsampled and downconverted output
     # complex baseband value yc.
@@ -707,15 +729,19 @@ def non_maximal_downsample_bpf(x, Ndown, k, Ndft, coefs, verbosity=1):
 
     # PFS with Ndft polyphases
     pfs = PolyPhaseFirFilterStructure(Ndft, coefs)
-    phasors = phasor_arr(k, Ndft, 'positive')
+    kPhasors = unit_circle_loops_phasor_arr(k, Ndft, 'positive')  # [HARRIS Eq 6.8]
+
+    # Oversampling time shift compensation via frequency dependent phase shift
+    tPhasors = time_shift_phasor_arr(k, Ndown, Ndft, Nblocks)  # [HARRIS Eq 9.3]
 
     for b in range(Nblocks):
         # Filter block
         inData = xBlocks[:, b]
         pfsData = pfs.filter_block(inData, flipped=True)
-        # Phase rotate polyphases for bin k [HARRIS Eq 6.8]
-        pfsBinData = pfsData * phasors
-        # Sum the polyphases to get single downsampled and downconverted output value
+        # Phase rotate polyphases for bin k
+        pfsBinData = pfsData * kPhasors * tPhasors[b]  # [HARRIS Eq 6.8, 9.3]
+        # Sum the polyphases to get single downsampled and downconverted output
+        # value
         yc[b] = np.sum(pfsBinData)
 
     if verbosity: