From ccfdcc45736cfb103f2d09408f6c556aa67d2b04 Mon Sep 17 00:00:00 2001
From: Eric Kooistra <kooistra@astron.nl>
Date: Thu, 6 Jun 2024 08:28:28 +0200
Subject: [PATCH] Workaround for fcutoff with firls.
---
applications/lofar2/model/rtdsp/firfilter.py | 161 +++++++++++--------
1 file changed, 90 insertions(+), 71 deletions(-)
diff --git a/applications/lofar2/model/rtdsp/firfilter.py b/applications/lofar2/model/rtdsp/firfilter.py
index cbfb0a3915..beea36c2d2 100644
--- a/applications/lofar2/model/rtdsp/firfilter.py
+++ b/applications/lofar2/model/rtdsp/firfilter.py
@@ -156,65 +156,10 @@ def ideal_low_pass_filter(Npoints, Npass, bandEdgeGain=1.0):
return h, f, HF
-def prototype_fir_low_pass_filter(method='firls',
- Npoints=1024, Ntaps=16, Ncoefs=1024*16,
- hpFactor=0.9, transitionFactor=0.4, stopRippleFactor=1000000, kaiserBeta=1, fs=1.0):
- """Derive FIR coefficients for prototype low pass filter
-
- Use method 'firls' or 'remez'.
-
- Default use LPF specification for LOFAR subband filter. For subband filter
- BWbin = fs / Npoints and the number of subbands is Nsub = Npoints / 2.
-
- For multirate upsampling or or downsampling filter the Npoints is the
- multirate factor Nup in upsampling or Ndown in downsampling. The fpass =
- BWbin / 2 = fNyquist / Npoints, where fNyquist = fs / 2.
-
- Use hpFactor < 1.0 to increase attentuation at fpass, to reduce disturbance
- due to aliasing.
-
- The number of prototype FIR filter coefficients Ncoefs <= Npoints * Ntaps.
- For symmetrical-coefficient FIR filters the group delay is
- (Ncoefs - 1) / 2, so choose Ncoefs is odd to have integers number of
- samples output delay.
-
- Normalize DC gain to 1.0.
-
- Input:
- - method: FIR filter design method 'firls' or 'remez'
- - Npoints: Number of points of the DFT in the filterbank, or multirate
- factor Nup in upsampling or Ndown in downsampling
- - Ntaps: Number of taps per polyphase.
- - Ncoefs: Actual number of prototype FIR filter coefficients, any extra
- coefs are zero
- - hpFactor : Half power bandwidth of the filter relative to BWbin
- - transitionFactor: transition bandwidth factor relative to fpass
- - stopRippleFactor: stopband ripple factor relative to pass band ripple
- - kaiserBeta: When kaiserBeta > 0 then additionally apply a Kaiser window on FIR
- coefficients
- - fs: sample frequency, for logging
- Return:
- - h: FIR coefficients for requested Ncoefs
- """
- # LPF specification for subband filter
- BWbin = fs / Npoints # bandwidth of one bin (is subband frequency)
- # Use half power bandwidth factor to tune half power cutoff frequency of LPF, default 1.0
- BWpass = hpFactor * BWbin
- fpass = BWpass / 2 # pass band frequency of bin at DC: -fpass to +fpass
- transitionBW = transitionFactor * fpass # transition bandwidth
- fcutoff = 0 # not used, use hpFactor instead
- cutoffGain = 0.5
- fstop = fpass + transitionBW # stop band frequency
- rippleWeights = [1, stopRippleFactor]
-
- # Design subband filter
- h = design_fir_low_pass_filter(method, Ncoefs, fpass, fstop, fcutoff, cutoffGain, rippleWeights, kaiserBeta, fs)
- return h
-
-
def design_fir_low_pass_filter(method,
Ncoefs, fpass, fstop, fcutoff=0, cutoffGain=0.5, rippleWeights=[1, 1],
- kaiserBeta=0, fs=1.0):
+ kaiserBeta=0, fs=1.0,
+ cutoffBW=0):
"""Derive FIR coefficients for prototype low pass filter
Use method 'firls' or 'remez', fs = 1.0
@@ -225,6 +170,17 @@ def design_fir_low_pass_filter(method,
Normalize DC gain to 1.0.
+ Remark:
+ . The actual LPF gain at fcutoff will differ from specified cutoffGain by about 1%.
+
+ Workaround:
+ . Method 'firls' does not support zero BW at fcutoff. With single fcutoff point error occurs: 'bands must be
+ monotonically nondecreasing and have width > 0'. Work around: use a very narrow band cutoffBW at fcutoff.
+ The cutoffBW band has to be much smaller than the transition band transistionBW, to avoid a flattened gain
+ around fcutoff, but also not too small, to because then the cutoffGain at fcutoff is not achieved.
+ Default force cutoffBW = transistionBW / 10 for 'firls' when input cutoffBW=0. For 'remez' method the
+ cutoffBW=0 effectively, without work around, because 'remez' does support single fcutoff value specification.
+
Input:
- method: FIR filter design method 'firls' or 'remez'
- Ncoefs: Actual number of prototype FIR filter coefficients, any extra
@@ -237,11 +193,13 @@ def design_fir_low_pass_filter(method,
- rippleWeights: relative ripple factors for pass band, optional fcutoff, and stop band
- kaiserBeta: When kaiserBeta > 0, then additionally apply a Kaiser window on FIR
coefficients
+ - cutoffBW: workaround bandwidth at fcutoff for firls, because firls does not support zero BW at fcutoff.
Return:
- h: FIR coefficients for requested Ncoefs
"""
fNyquist = fs / 2
- fsub = fpass * 2 # double sided LPF band width, when used as prototype for a BPF
+ lpBW = fpass + fstop # low pass bandwidth is twice center of transition band
+ transistionBW = fstop - fpass
# Initial FIR filter length N
# . Use interpolation of factor Q shorter filter to ensure the FIR filter design converges
@@ -268,6 +226,7 @@ def design_fir_low_pass_filter(method,
f_pb = fpass * Q # pass band frequency
f_co = fcutoff * Q # optional cutoff frequency in transition band
f_sb = fstop * Q # stop band frequency
+ f_tb = transistionBW * Q
# Calculate initial FIR filter for length N
if method == 'firls':
@@ -275,7 +234,11 @@ def design_fir_low_pass_filter(method,
hFir = signal.firls(N, [0, f_pb, f_sb, fNyquist],
[1, 1, 0, 0], rippleWeights, fs=fs)
else:
- hFir = signal.firls(N, [0, f_pb, f_co, f_co, f_sb, fNyquist],
+ if cutoffBW == 0:
+ b_co = f_tb / 10 # default
+ else:
+ b_co = cutoffBW * Q # from input argument
+ hFir = signal.firls(N, [0, f_pb, f_co - b_co / 2, f_co + b_co / 2, f_sb, fNyquist],
[1, 1, cutoffGain, cutoffGain, 0, 0], rippleWeights, fs=fs)
if method == 'remez':
if not fcutoff:
@@ -305,18 +268,74 @@ def design_fir_low_pass_filter(method,
h = fourier_interpolate(HFfir, Ncoefs)
print('> design_fir_low_pass_filter()')
- print('. Method = %s' % method)
- print('. Q = %f' % Q)
- print('. fsub = fpass * 2 = %f' % fsub)
- print('. fpass = %f' % fpass)
+ print('. Method = %s' % method)
+ print('. Q = %f' % Q)
+ print('. fpass = %f' % fpass)
+ print('. fstop = %f' % fstop)
+ print('. lpBW = %f' % lpBW)
if fcutoff:
- print('. fcutoff = %f' % fcutoff)
- print('. cutoffGain = %f' % cutoffGain)
- print('. fstop = %f' % fstop)
- print('. fNyquist = %f' % fNyquist)
- print('. fs = %f' % fs)
- print('. Ncoefs = %d' % len(h))
- print('. DC sum = %f' % np.sum(h))
- print('. Symmetrical coefs = %s' % is_symmetrical(h))
+ print('. fcutoff = %f' % fcutoff)
+ print('. cutoffGain = %f' % cutoffGain)
+ print('. fNyquist = %f' % fNyquist)
+ print('. fs = %f' % fs)
+ print('. Ncoefs = %d' % len(h))
+ print('. DC sum = %f' % np.sum(h))
+ print('. Symmetrical coefs = %s' % is_symmetrical(h))
print('')
return h
+
+
+def prototype_fir_low_pass_filter(method='firls',
+ Npoints=1024, Ntaps=16, Ncoefs=1024*16,
+ hpFactor=0.9, transitionFactor=0.4, stopRippleFactor=1000000, kaiserBeta=0, fs=1.0):
+ """Derive FIR coefficients for prototype low pass filter
+
+ Use method 'firls' or 'remez'.
+
+ Default use LPF specification for LOFAR subband filter. For subband filter
+ BWbin = fs / Npoints and the number of subbands is Nsub = Npoints / 2.
+
+ For multirate upsampling or or downsampling filter the Npoints is the
+ multirate factor Nup in upsampling or Ndown in downsampling. The fpass =
+ BWbin / 2 = fNyquist / Npoints, where fNyquist = fs / 2.
+
+ Use hpFactor < 1.0 to increase attentuation at fpass, to reduce disturbance
+ due to aliasing.
+
+ The number of prototype FIR filter coefficients Ncoefs <= Npoints * Ntaps.
+ For symmetrical-coefficient FIR filters the group delay is
+ (Ncoefs - 1) / 2, so choose Ncoefs is odd to have integers number of
+ samples output delay.
+
+ Normalize DC gain to 1.0.
+
+ Input:
+ - method: FIR filter design method 'firls' or 'remez'
+ - Npoints: Number of points of the DFT in the filterbank, or multirate
+ factor Nup in upsampling or Ndown in downsampling
+ - Ntaps: Number of taps per polyphase.
+ - Ncoefs: Actual number of prototype FIR filter coefficients, any extra
+ coefs are zero
+ - hpFactor : Half power bandwidth of the filter relative to BWbin
+ - transitionFactor: transition bandwidth factor relative to fpass
+ - stopRippleFactor: stopband ripple factor relative to pass band ripple
+ - kaiserBeta: When kaiserBeta > 0 then additionally apply a Kaiser window on FIR
+ coefficients
+ - fs: sample frequency, for logging
+ Return:
+ - h: FIR coefficients for requested Ncoefs
+ """
+ # LPF specification for subband filter
+ BWbin = fs / Npoints # bandwidth of one bin (is subband frequency)
+ # Use half power bandwidth factor to tune half power cutoff frequency of LPF, default 1.0
+ BWpass = hpFactor * BWbin
+ fpass = BWpass / 2 # pass band frequency of bin at DC: -fpass to +fpass
+ transitionBW = transitionFactor * fpass # transition bandwidth
+ fcutoff = 0 # not used, use hpFactor instead
+ cutoffGain = 0.5
+ fstop = fpass + transitionBW # stop band frequency
+ rippleWeights = [1, stopRippleFactor]
+
+ # Design subband filter
+ h = design_fir_low_pass_filter(method, Ncoefs, fpass, fstop, fcutoff, cutoffGain, rippleWeights, kaiserBeta, fs)
+ return h
--
GitLab