diff --git a/applications/lofar2/model/rtdsp/firfilter.py b/applications/lofar2/model/rtdsp/firfilter.py index daa2d90926af16ce057eac64513a24a6189a7926..679b84149cafaf62adcf204562cb96cfdb8eb153 100644 --- a/applications/lofar2/model/rtdsp/firfilter.py +++ b/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,36 +353,28 @@ 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) - if 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)) - # 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 result + if result: + if 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 + else: + print('ERROR: Return None.') return None - return h def prototype_fir_low_pass_filter(method='firls',