Correct design_fir_low_pass_filter_adjust() for large stop band weights.

applications/lofar2/model/rtdsp/firfilter.py

@@ -161,7 +161,7 @@ def ideal_low_pass_filter(Npoints, Npass, bandEdgeGain=1.0):
 def design_fir_low_pass_filter(method,
                                Ncoefs, fpass, fstop, fcutoff=0, cutoffGain=0.5, rippleWeights=[1, 1],
                                kaiserBeta=0, fs=1.0,
-                               cutoffBW=0, verbosity=True):
+                               cutoffBW=0, verbosity=1):
     """Derive FIR coefficients for low pass filter
 
     Use method 'firls' or 'remez', fs = 1.0
@@ -267,7 +267,7 @@ def design_fir_low_pass_filter(method,
         # Use Fourier interpolation to create final FIR filter coefs when
         # Q > 1 or when N != Ncoefs
         HFfir = np.fft.fft(hFir)
-        h = fourier_interpolate(HFfir, Ncoefs)
+        h = fourier_interpolate(HFfir, Ncoefs, verbosity=0)
 
     if verbosity:
         print('> design_fir_low_pass_filter():')
@@ -298,7 +298,7 @@ def design_fir_low_pass_filter_adjust(method,
     """
     # Determine whether to adjust fpass or fstop
     nofIterations = 1
-    nofIterationsMax = 100
+    nofIterationsMax = 1000
     h = design_fir_low_pass_filter(method, Ncoefs, fpass, fstop, rippleWeights=rippleWeights,
                                    kaiserBeta=kaiserBeta, fs=fs, verbosity=0)
     _, fn, HF = dtft(h)
@@ -312,18 +312,38 @@ def design_fir_low_pass_filter_adjust(method,
         adjustBW = fstop - fcutoff
 
     # Iterate to achieve cutoffGain at fcutoff
-    fRadix = 32
+    fRadix = 16
     fStep = adjustBW / fRadix
     gainResolution = cutoffGain * resolution
     FP = fpass
     FS = fstop
+    result = False
     while True:
-        prevFP = FP
-        prevFS = FS
+        # Calculate iteration
+        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)
+        # print('DEBUG: FP, FS, fIndex, fValue, fGain = %s' % str((FP, FS, fIndex, fValue, fGain)))
+        if np.abs(fGain - cutoffGain) < gainResolution:
+            result = True
+            break
+        if nofIterations > nofIterationsMax:
+            print('ERROR: Too many iterations.')
+            break
+        # Prepare next iteration
         if adjustBand == 0:
             if fGain < cutoffGain:
                 # Keep on increasing FP
                 FP += fStep
+                # print('DEBUG: Increase 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
             else:
                 # FP too high, one fStep back and increment FP with reduced fStep
                 FP -= fStep
@@ -333,21 +353,15 @@ def design_fir_low_pass_filter_adjust(method,
             if fGain > cutoffGain:
                 # Keep on decreasing FS
                 FS -= fStep
+                # print('DEBUG: Decrease FS')
             else:
                 # FS too low, one fStep back and decrement FS with reduced fStep
                 FS += fStep
                 fStep /= fRadix
                 FS -= fStep
-        if np.abs(fGain - cutoffGain) < gainResolution:
-            break
-        if nofIterations > nofIterationsMax:
-            break
-        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(f, HF, fcutoff)
 
+    # Return result
+    if result:
         if verbosity:
             print('> design_fir_low_pass_filter_adjust():')
             print('. nofIterations        = %d' % nofIterations)
@@ -357,12 +371,10 @@ def design_fir_low_pass_filter_adjust(method,
             print('. fstop / fcutoff      = %f' % (fstop / fcutoff))
             print('. FP / fpass           = %f' % (FP / fpass))
             print('. FS / fstop           = %f' % (FS / fstop))
-        # Rerun with same parameters to get logging
-        h = design_fir_low_pass_filter(method, Ncoefs, prevFP, prevFS, rippleWeights=rippleWeights,
-                                       kaiserBeta=kaiserBeta, fs=fs)
-    if nofIterations > nofIterationsMax:
-        return None
         return h
+    else:
+        print('ERROR: Return None.')
+        return None
 
 
 def prototype_fir_low_pass_filter(method='firls',