Merge branch 'L2SDP-260' into 'master'

Resolve L2SDP-260

Closes L2SDP-260

See merge request !405

applications/lofar2/model/pfb_os/README.txt

@@ -13,7 +13,7 @@ Author: Eric Kooistra, nov 2023
   can not exactly reproduce the actual LOFAR1 coefficients, therefore these
   are loaded from a file Coeffs16384Kaiser-quant.dat
 
-* Try low pass filter design methods using windowed sync, firls, remez [3]
+* Try low pass filter design methods using windowed sync, firls, remez [4]
   The windowed sync method, firls leased squares method and remez method all
   yield comparable results, but firls and remez perform slightly better near
   the transition band. The firls and remez functions from scipy.signal use
@@ -27,8 +27,20 @@ Author: Eric Kooistra, nov 2023
 
 [1] dsp_study_erko.txt, summary of DSP books
 [2] pyfda, dsp, at https://github.com/chipmuenk
-[3] Try FIR filter design methods
-    * dsp.py import for Python jupyter notebooks
-    * filter_design_firls.ipynb
-    * filter_design_remez.ipynb
-    * filter_design_windowed_sync.ipynb
+[3] dsp.py import for Python jupyter notebooks
+[4] python jupyter notebooks
+    * Try FIR filter design methods
+      - rectangular_window_and_ideal_lpf.ipynb
+      - filter_design_windowed_sync.ipynb
+      - filter_design_firls.ipynb
+      - filter_design_remez.ipynb
+    * Try Hilbert transform
+      - hilbert_transform_design.ipynb
+      - hilbert_transform_application.ipynb
+    * Try IIR filter design methods
+      - iir_filter.ipynb
+    * Try multirate processing
+      - up_down_sampling.ipynb
+      - narrowband_noise_generator.ipynb
+    * Try polyphase filterbanks
+      - one_pfb.ipynb

applications/lofar2/model/pfb_os/dsp.py

@@ -230,9 +230,11 @@ def prototype_fir_low_pass_filter(
     w_tb = transitionFactor * fpass * Q  # transition bandwidth
     f_sb = f_pb + w_tb  # stop band frequency
     weight = [1, stopRippleFactor]  # weight pass band ripple versus stop band ripple
-    print('f_pb = ', str(f_pb))
-    print('w_tb = ', str(w_tb))
-    print('f_sb = ', str(f_sb))
+    print('> prototype_fir_low_pass_filter()')
+    print('Pass band specification')
+    print('. f_pb = ', str(f_pb))
+    print('. w_tb = ', str(w_tb))
+    print('. f_sb = ', str(f_sb))
     hFirls = signal.firls(N, [0, f_pb, f_sb, 0.5], [1, 1, 0, 0], weight, fs=1.0)
 
     # Apply Kaiser window with beta = 1 like in pfs_coeff_final.m, this improves the
@@ -258,7 +260,7 @@ def prototype_fir_low_pass_filter(
     print('. Symmetrical coefs = %s' % is_symmetrical(hFirls))
 
     if len(hFirls) == Ncoefs:
-        return hFirls
+        h = hFirls
     else:
         # Use Fourier interpolation to create final FIR filter coefs when
         # Q > 1 or when Ncoefs is even
@@ -269,7 +271,9 @@ def prototype_fir_low_pass_filter(
         print('. Ncoefs            = %d' % len(hInterpolated))
         print('. DC sum            = %f' % np.sum(hInterpolated))
         print('. Symmetrical coefs = %s' % is_symmetrical(hInterpolated))
-        return hInterpolated
+        h = hInterpolated
+    print('')
+    return h
 
 
 def fourier_interpolate(HFfilter, Ncoefs):
@@ -498,9 +502,7 @@ class PolyPhaseFirFilterStructure:
         #                [ 1,  5,  9],
         #                [ 0,  4,  8]]), Nphases = 4 rows (axis=0)
         #                                Ntaps = 3 columns (axis=1)
-        self.polyCoefs = self.poly_coeffs()
-        if commutator == 'down':
-            self.polyCoefs = np.flip(self.polyCoefs, axis=0)
+        self.polyCoefs = self.poly_coeffs(commutator)
 
         #             oldest sample[0]
         #             newest sample [15]
@@ -514,25 +516,28 @@ class PolyPhaseFirFilterStructure:
         #                 [11,  7,  3]])
         self.polyDelays = np.zeros((Nphases, self.Ntaps))
 
-    def poly_coeffs(self):
+    def poly_coeffs(self, commutator):
         """Polyphase structure for FIR filter coefficients coefs in Nphases
 
         Input:
         . coefs = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
+        . commutator: 'up' or 'down'
         Return:
-        . polyCoefs = [[ 0,  4,  8],   # p = 0
-                       [ 1,  5,  9],   # p = 1
-                       [ 2,  6, 10],   # p = 2
-                       [ 3,  7, 11]])  # p = 3, Nphases = 4 rows (axis=0)
+        . commutator: 'up', see title Figure in [HARRIS 6], Figure 10.22 and
+          10.23 in [LYONS] seem to have commutator arrows in wrong direction.
+          polyCoefs = [[ 0,  4,  8],   # p = 0, H0(z)
+                       [ 1,  5,  9],   # p = 1, H1(z)
+                       [ 2,  6, 10],   # p = 2, H2(z)
+                       [ 3,  7, 11]])  # p = 3, H3(z),
+                                         Nphases = 4 rows (axis=0)
                                          Ntaps = 3 columns (axis=1)
         If Ncoefs < Ntaps * Nphases, then pad polyCoefs with zeros.
         """
-        Ncoefs = self.Ncoefs
-        Nphases = self.Nphases
-        Ntaps = self.Ntaps
-        polyCoefs = np.zeros(Nphases * Ntaps)
-        polyCoefs[0:Ncoefs] = self.coefs
-        polyCoefs = polyCoefs.reshape(Ntaps, Nphases).T
+        polyCoefs = np.zeros(self.Nphases * self.Ntaps)
+        polyCoefs[0 : self.Ncoefs] = self.coefs
+        polyCoefs = polyCoefs.reshape(self.Ntaps, self.Nphases).T
+        if commutator == 'down':
+            polyCoefs = np.flip(polyCoefs, axis=0)
         return polyCoefs
 
     def filter_block(self, inData):
@@ -546,7 +551,6 @@ class PolyPhaseFirFilterStructure:
                    oldest data and outData[Nphases-1] is newest data, so time n
                    increments, like with inData
         """
-        Ntaps = self.Ntaps
         # Shift delay memory by one block
         self.polyDelays = np.roll(self.polyDelays, 1, axis=1)
         # Shift in inData block at column 0
@@ -633,7 +637,6 @@ def upsample(x, Nup, coefs, verify=False):  # interpolate
               y[n * I + 3] = lfilter(h[3, 7,11], [1], x)   # p = 3
     """
     Nx = len(x)
-    Ncoefs = len(coefs)
     a = [1.0]
 
     # Polyphase implementation to avoid multiply by zero values
@@ -667,10 +670,11 @@ def upsample(x, Nup, coefs, verify=False):  # interpolate
         else:
             print('ERROR: wrong upsample result')
 
-    print('Upsampling:')
+    print('> upsample():')
     print('. Nup = ' + str(Nup))
     print('. Nx = ' + str(Nx))
     print('. len(y) = ' + str(len(y)))
+    print('')
     return y
 
 
@@ -755,7 +759,6 @@ def downsample(x, Ndown, coefs, verify=False):  # decimate
     Nxp = (Ndown - 1 + len(x)) // Ndown
     # Used number of samples from x
     Nx = 1 + Ndown * (Nxp - 1)
-    Ncoefs = len(coefs)
     a = [1.0]
 
     # Polyphase implementation to avoid calculating values that are removed
@@ -793,12 +796,13 @@ def downsample(x, Ndown, coefs, verify=False):  # decimate
         else:
             print('ERROR: wrong downsample result')
 
-    print('Downsampling:')
+    print('> downsample():')
     print('. len(x) = ' + str(len(x)))
     print('. Ndown = ' + str(Ndown))
     print('. Nx = ' + str(Nx))
     print('. Nxp = ' + str(Nxp))
     print('. len(y) = ' + str(len(y)))  # = Nxp
+    print('')
     return y
 
 
@@ -849,7 +853,6 @@ def resample(x, Nup, Ndown, coefs, verify=False):  # interpolate and decimate by
     """
     Nx = len(x)
     Nv = Nup * Nx
-    Ncoefs = len(coefs)
     a = [1.0]
     Nyp = (Nv + Ndown - 1) // Ndown   # Number of samples from v, per poly phase in polyY
     Ny = 1 + Ndown * (Nyp - 1)  # Used number of samples from v
@@ -883,13 +886,14 @@ def resample(x, Nup, Ndown, coefs, verify=False):  # interpolate and decimate by
         else:
             print('ERROR: wrong resample result')
 
-    print('Resampling:')
+    print('> resample():')
     print('. len(x) = ' + str(len(x)))
     print('. Nx = ' + str(Nx))
     print('. len(v) = ' + str(len(v)))
     print('. Ny = ' + str(Ny))
     print('. Nyp = ' + str(Nyp))
     print('. len(y) = ' + str(len(y)))
+    print('')
     return y
 
 
@@ -952,6 +956,7 @@ def plot_iir_filter_analysis(b, a, fs=1, whole=False, Ntime=100, step=False, log
         z, p, k = dsp_fpga_lib.zplane(b, a)  # no plot, only calculate z, p, k
 
     if log:
+        print('plot_iir_filter_analysis()')
         print('Zeros, poles and gain from b, a coefficients:')
         if len(z) > 0:
             print('. zeros:')
@@ -967,6 +972,7 @@ def plot_iir_filter_analysis(b, a, fs=1, whole=False, Ntime=100, step=False, log
         print('Coefficients back from z, p, k:')
         print('  b = %s' % str(np.poly(z)))
         print('  a = %s' % str(np.poly(p) / k))
+        print('')
 
     # Plot transfer function H(f), is H(z) for z = exp(j w), so along the unit circle
     # . 0 Hz at 1 + 0j, fs / 4 at 0 + 1j, fNyquist = fs / 2 at -1 + 0j
@@ -1072,10 +1078,7 @@ def plot_power_spectrum(f, HF, fmt='r', fs=1.0, fLim=None, dbLim=None):
     . fmt: curve format string
     . fs: sample frequency in Hz, scale f by fs, fs >= 1
     """
-    if fs > 1:
-        flabel = 'Frequency'
-    else:
-        flabel = 'Frequency [fs]'
+    flabel = 'Frequency [fs = %f]' % fs
 
     plt.plot(f * fs, pow_db(HF), fmt)
     plt.title('Power spectrum')

applications/lofar2/model/pfb_os/dsp_study_erko.txt

@@ -597,15 +597,17 @@
   then W_N can be used directly in equation and matrix:
 
   W_N = exp(-j 2pi/N) is primitive Nth root of unity
-  W_N^k = exp(-j 2pi k / N)
-  W_N^kn = exp(-j 2pi k / N * n) = exp(-j w_k * t_n)
+  W_N^k = exp(-j 2pi/N k)
+  W_N^(kn) = exp(-j 2pi/N k n) = exp(-j w_k * t_n)
   . N is number of time samples and number of frequency point
-  . w_k = k 2pi fs / N
+  . w_k = 2pi/N k fs
   . t_n = n Ts
   . fs / N is the frequency sampling interval in Hz
   . Ts = 1 / fs is time sampling interval in seconds
   . ws = 2pi fs is the frequency sampling interval in rad/s
 
+  W_N^(nk) = W_N^(nk % N), because W_N^N = exp(-j2pi) = 1
+
 - Normalized frequency axis [LYONS 3.5, 3.13.4]:
   . fs = 1, ws = 2pi fs = 2pi
   .  -fs/2, 0, fs/2 [Hz]
@@ -647,10 +649,10 @@
                     N-1
     X(w_k) = X(k) = sum x(n) W_N^kn
                     n=0      exp(-j w_k t_n)
-                             exp(-j 2pi k n / N)
+                             exp(-j 2pi/N k n)
       with:
       . N is number of time samples and number of frequency point
-      . w_k = k 2pi fs / N
+      . w_k = 2pi/N k fs
       . t_n = n Ts
 
   Inverse DFT:
@@ -676,13 +678,13 @@
     |...   |   |...                     ...               | |...   |
     |X(N-1)|   |1  W_N^(N-1) W_N^2(N-1) ... W_N^(N-1)(N-1)| |x(N-1)|
 
-  DFT matrix is Square matrix, because N is number of time samples in x(n) and
-  number of frequency points in X(k).
+  The DFT matrix WN is a square N x N matrix, because N is number of time
+  samples in x(n) and number of frequency points in X(k).
 
   IDFT:
-  xN = WN^-1 XN = 1/N conj(WN) XN, so
+    xN = WN^-1 XN = 1/N conj(WN) XN,
 
-  WN conj(WN) = N IN, where IN i N x N identity matrix
+       so: WN conj(WN) = N IN, where IN is the N x N identity matrix.
 
 - Real input: X(k) = conj(X(N - k))
 
@@ -694,7 +696,7 @@
   . Zero padding N to Ndtft increases the interpolation density, but does not
     increase the frequency resolution in the ability to resolve, to distinguish
     between closely space features in the spectrum. The frequency resolution
-    can only be increased by increasing N
+    can only be increased by increasing N.
 - N point DFT of rectangular function yields aliased sinc (= Dirchlet kernel):
 
   x(n) = 1 for K samples