Skip to content
Snippets Groups Projects
Commit ccfdcc45 authored by Eric Kooistra's avatar Eric Kooistra
Browse files

Workaround for fcutoff with firls.

parent f413e509
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
applications/lofar2/model/rtdsp/firfilter.py
+ 90
− 71
View file @ ccfdcc45
......@@ -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:
......@@ -307,12 +270,12 @@ def design_fir_low_pass_filter(method,
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('. 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))
......@@ -320,3 +283,59 @@ def design_fir_low_pass_filter(method,
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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment