Improve design_fir_low_pass_filter_adjust(). Add filterbank_frequency_response().

applications/lofar2/model/rtdsp/firfilter.py

@@ -29,7 +29,7 @@
 
 import numpy as np
 from scipy import signal
-from .utilities import ceil_div, is_even, is_symmetrical
+from .utilities import c_rtol, c_atol, ceil_div, ceil_pow2, is_even, is_symmetrical
 from .fourier import fourier_interpolate, dtft, estimate_gain_at_frequency
 
 
@@ -215,6 +215,9 @@ def design_fir_low_pass_filter(method,
     if method == 'firls':
         NcoefsFirlsMax = 1024
         Q = ceil_div(Ncoefs, NcoefsFirlsMax)
+        # Band edges should be less than fNyquist = fs / 2
+        if fstop * Q >= fNyquist:
+            Q = int(fNyquist / fstop)
         N = Ncoefs // Q
         # If necessary -1 to make odd, because firls only supports odd number of FIR coefficients
         if is_even(N):
@@ -222,6 +225,9 @@ def design_fir_low_pass_filter(method,
     if method == 'remez':
         NcoefsRemezMax = 512
         Q = ceil_div(Ncoefs, NcoefsRemezMax)
+        # Band edges should be less than fNyquist = fs / 2
+        if fstop * Q >= fNyquist:
+            Q = int(fNyquist / fstop)
         N = Ncoefs // Q
 
     # Low pass transition band
@@ -276,6 +282,7 @@ def design_fir_low_pass_filter(method,
         print('. fpass               = %f' % fpass)
         print('. fstop               = %f' % fstop)
         print('. lpBW                = %f' % lpBW)
+        print('. transistionBW       = %f' % transistionBW)
         if fcutoff:
             print('. fcutoff             = %f' % fcutoff)
             print('. cutoffGain          = %f' % cutoffGain)
@@ -290,18 +297,25 @@ def design_fir_low_pass_filter(method,
 
 def design_fir_low_pass_filter_adjust(method,
                                       Ncoefs, fpass, fstop, fcutoff, cutoffGain=0.5, rippleWeights=[1, 1],
-                                      kaiserBeta=0, fs=1.0, resolution=1e-4, verbosity=1):
+                                      kaiserBeta=0, fs=1.0, resolution=1e-3, verbosity=1):
     """Derive FIR coefficients for low pass filter for exact fcutoff
 
     Uses design_fir_low_pass_filter() but achieves cutoffGain at fcutoff by adjusting fpass or fstop, and
     fpass < fcutoff < fstop.
+
+    Input:
+    . resolution: Resolution measure for frequency and for gain. It appears that at some fine resolution
+      even finer resolution does not yield more accurate cutoffGain. Possibly due to the method.
+    Return:
+    . h: LPF impulse response
     """
     # Determine whether to adjust fpass or fstop
     nofIterations = 1
-    nofIterationsMax = 1000
+    nofIterationsMax = 100
+    Ndtft = ceil_pow2(max(Ncoefs, 1 / resolution) * 10)
     h = design_fir_low_pass_filter(method, Ncoefs, fpass, fstop, rippleWeights=rippleWeights,
                                    kaiserBeta=kaiserBeta, fs=fs, verbosity=0)
-    _, fn, HF = dtft(h)
+    _, fn, HF = dtft(h, Ndtft)
     f = fn * fs
     fIndex, fValue, fGain = estimate_gain_at_frequency(f, HF, fcutoff)
     if fGain < cutoffGain:
@@ -312,9 +326,10 @@ def design_fir_low_pass_filter_adjust(method,
         adjustBW = fstop - fcutoff
 
     # Iterate to achieve cutoffGain at fcutoff
-    fRadix = 16
+    fRadix = 8
     fStep = adjustBW / fRadix
     gainResolution = cutoffGain * resolution
+    freqResolution = 0.1 * fs / Ndtft
     FP = fpass
     FS = fstop
     result = False
@@ -323,12 +338,15 @@ def design_fir_low_pass_filter_adjust(method,
         nofIterations += 1
         h = design_fir_low_pass_filter(method, Ncoefs, FP, FS, rippleWeights=rippleWeights,
                                        kaiserBeta=kaiserBeta, fs=fs, verbosity=0)
-        _, fn, HF = dtft(h)
-        fIndex, fValue, fGain = estimate_gain_at_frequency(fn, HF, fcutoff)
+        _, fn, HF = dtft(h, Ndtft)
+        fIndex, fValue, fGain = estimate_gain_at_frequency(f, HF, fcutoff)
         # print('DEBUG: FP, FS, fIndex, fValue, fGain = %s' % str((FP, FS, fIndex, fValue, fGain)))
         if np.abs(fGain - cutoffGain) < gainResolution:
             result = True
             break
+        if fStep < freqResolution:
+            result = True
+            break
         if nofIterations > nofIterationsMax:
             print('ERROR: Too many iterations.')
             break
@@ -337,44 +355,48 @@ def design_fir_low_pass_filter_adjust(method,
             if fGain < cutoffGain:
                 # Keep on increasing FP
                 FP += fStep
-                # print('DEBUG: Increase FP')
+                # print('DEBUG: Increase FP = %s' % str(FP))
                 # If necessary increase FS to keep fcutoff, this can occur when rippleWeights for stop band >
                 # rippleWeights for pass band
                 if FP >= fcutoff:
-                    # print('DEBUG: Increment FS')
                     FP = fpass
                     FS += fStep
+                    # print('DEBUG: Need to increment FS = %s' % str(FS))
             else:
                 # FP too high, one fStep back and increment FP with reduced fStep
                 FP -= fStep
                 fStep /= fRadix
                 FP += fStep
+                # print('DEBUG: FP too high, decrease fStep: FP, fStep = %s' % str((FP, fStep)))
         else:
             if fGain > cutoffGain:
                 # Keep on decreasing FS
                 FS -= fStep
-                # print('DEBUG: Decrease FS')
+                # print('DEBUG: Decrease FS = %s' % str(FS))
             else:
                 # FS too low, one fStep back and decrement FS with reduced fStep
                 FS += fStep
                 fStep /= fRadix
                 FS -= fStep
+                # print('DEBUG: FS too low, decrease fStep: FS, fStep = %s' % str((FS, fStep)))
 
     # Return result
     if result:
         if verbosity:
+            # Repeat design_fir_low_pass_filter() for logging
+            h = design_fir_low_pass_filter(method, Ncoefs, FP, FS, rippleWeights=rippleWeights,
+                                           kaiserBeta=kaiserBeta, fs=fs, verbosity=verbosity)
             print('> design_fir_low_pass_filter_adjust():')
             print('. nofIterations        = %d' % nofIterations)
-            print('. fpass, fcutoff, stop = %s' % str((fpass, fcutoff, fstop)))
-            print('. FP, fcutoff, FS      = %s' % str((FP, fcutoff, FS)))
-            print('. fcutoff / fpass      = %f' % (fcutoff / fpass))
-            print('. fstop / fcutoff      = %f' % (fstop / fcutoff))
             print('. FP / fpass           = %f' % (FP / fpass))
             print('. FS / fstop           = %f' % (FS / fstop))
-        return h
+            print('. fcutoff              = %f' % fcutoff)
+            print('. fGain                = %f' % fGain)
+            print('. fGain**2             = %f' % (fGain**2))
+            print('')
     else:
-        print('ERROR: Return None.')
-        return None
+        print('ERROR: Return iteration %d.' % nofIterations)
+    return h
 
 
 def prototype_fir_low_pass_filter(method='firls',
@@ -437,20 +459,40 @@ def prototype_fir_low_pass_filter(method='firls',
 # Filterbank
 ###############################################################################
 
-def filterbank_aggregate_transfer(HFprototype, Npoints):
+def filterbank_frequency_response(HFprototype, Npoints):
     """Aggregate frequency reponse of filterbank
 
     Sum of bandpass responses for all Npoints frequency bins in a filterbank.
 
-             N-1
-    HFbank = sum HF_k(w), with HF_k(w) = H0(w - 2pi/N k)
-             k=0
+                N-1
+    HFbank(w) = sum HF_k(w), with HF_k(w) = H0(w - 2pi/N k)
+                k=0
 
     and H0(w) is the prototype LPF [PROAKIS Eq. 10.9.2].
     """
     Lprototype = len(HFprototype)  # corresponds to 2 pi
-    Lbin = Lprototype / Npoints  # corresponds to 2 pi / N
+    Lbin = Lprototype // Npoints  # corresponds to 2 pi / N
     HFbank = np.zeros(Lprototype, dtype='cfloat')
     for bi in range(Npoints):
         HFbank += np.roll(HFprototype, bi * Lbin)
     return HFbank
+
+
+def filterbank_impulse_response(hprototype, Npoints):
+    """Impulse reponse of filterbank
+
+    Sum of bandpass impulse responses for all Npoints frequency bins in a
+    filterbank.
+
+               N-1
+    hbank[n] = sum hprototype[n] * exp(-j * 2 * np.pi / Npoints * k * n)
+               k=0
+    """
+    Lprototype = len(hprototype)
+    hbank = np.zeros(Lprototype, dtype='cfloat')
+    for k in range(Npoints):
+        hbank += hprototype * np.exp(-1.0j * 2 * np.pi / Npoints * k * np.arange(Lprototype))
+    # hbank.imag is zero
+    if not np.allclose(hbank.imag, np.zeros(Lprototype), rtol=c_rtol, atol=c_atol * Npoints):
+        print('ERROR: hbank.imag != 0 (max(abs) = %e)' % np.max(np.abs(hbank.imag)))
+    return hbank.real