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

Add design_fir_low_pass_filter_adjust().

parent 97334f9f
No related branches found
No related tags found
1 merge request!411Resolve RTSD-271
...@@ -30,7 +30,7 @@ ...@@ -30,7 +30,7 @@
import numpy as np import numpy as np
from scipy import signal from scipy import signal
from .utilities import ceil_div, is_even, is_symmetrical from .utilities import ceil_div, is_even, is_symmetrical
from .fourier import fourier_interpolate from .fourier import fourier_interpolate, dtft, estimate_gain_at_frequency
############################################################################### ###############################################################################
...@@ -161,7 +161,7 @@ def ideal_low_pass_filter(Npoints, Npass, bandEdgeGain=1.0): ...@@ -161,7 +161,7 @@ def ideal_low_pass_filter(Npoints, Npass, bandEdgeGain=1.0):
def design_fir_low_pass_filter(method, def design_fir_low_pass_filter(method,
Ncoefs, fpass, fstop, fcutoff=0, cutoffGain=0.5, rippleWeights=[1, 1], Ncoefs, fpass, fstop, fcutoff=0, cutoffGain=0.5, rippleWeights=[1, 1],
kaiserBeta=0, fs=1.0, kaiserBeta=0, fs=1.0,
cutoffBW=0): cutoffBW=0, verbosity=True):
"""Derive FIR coefficients for low pass filter """Derive FIR coefficients for low pass filter
Use method 'firls' or 'remez', fs = 1.0 Use method 'firls' or 'remez', fs = 1.0
...@@ -196,6 +196,7 @@ def design_fir_low_pass_filter(method, ...@@ -196,6 +196,7 @@ def design_fir_low_pass_filter(method,
- kaiserBeta: When kaiserBeta > 0, then additionally apply a Kaiser window on FIR - kaiserBeta: When kaiserBeta > 0, then additionally apply a Kaiser window on FIR
coefficients coefficients
- cutoffBW: workaround bandwidth at fcutoff for firls, because firls does not support zero BW at fcutoff. - cutoffBW: workaround bandwidth at fcutoff for firls, because firls does not support zero BW at fcutoff.
- verbosity: when > 0 print() status, else no print()
Return: Return:
- h: FIR coefficients for requested Ncoefs - h: FIR coefficients for requested Ncoefs
""" """
...@@ -205,23 +206,22 @@ def design_fir_low_pass_filter(method, ...@@ -205,23 +206,22 @@ def design_fir_low_pass_filter(method,
# Initial FIR filter length N # Initial FIR filter length N
# . Use interpolation of factor Q shorter filter to ensure the FIR filter design converges # . Use interpolation of factor Q shorter filter to ensure the FIR filter design converges
# and to speed up calculation. N >> NFirlsMax = 1024 is not feasible. # and to speed up calculation. N >> NcoefsFirlsMax = 1024 is not feasible.
# . The passband ripple and stopband attenuation depend on the transition bandwidth # . The passband ripple and stopband attenuation depend on the transition bandwidth
# and the rippleWeights. Choose transition band width ~< fpass, to ensure the FIR filter # and the rippleWeights. Choose transition band width ~< fpass, to ensure the FIR filter
# design converges and improve the passband ripple and stopband attenuation. A too large # design converges and improve the passband ripple and stopband attenuation. A too large
# transition bandwidth also gives the design too much freedom and causes artifacts in # transition bandwidth also gives the design too much freedom and causes artifacts in
# the transition band. # the transition band.
# . scipy firls() defines fpass relative to fs, so can use fpass as cutoff frequency.
if method == 'firls': if method == 'firls':
NFirlsMax = 1024 NcoefsFirlsMax = 1024
Q = ceil_div(Ncoefs, NFirlsMax) Q = ceil_div(Ncoefs, NcoefsFirlsMax)
N = Ncoefs // Q N = Ncoefs // Q
# If necessary -1 to make odd, because firls only supports odd number of FIR coefficients # If necessary -1 to make odd, because firls only supports odd number of FIR coefficients
if is_even(N): if is_even(N):
N = N - 1 N = N - 1
if method == 'remez': if method == 'remez':
NRemezMax = 512 NcoefsRemezMax = 512
Q = ceil_div(Ncoefs, NRemezMax) Q = ceil_div(Ncoefs, NcoefsRemezMax)
N = Ncoefs // Q N = Ncoefs // Q
# Low pass transition band # Low pass transition band
...@@ -269,7 +269,8 @@ def design_fir_low_pass_filter(method, ...@@ -269,7 +269,8 @@ def design_fir_low_pass_filter(method,
HFfir = np.fft.fft(hFir) HFfir = np.fft.fft(hFir)
h = fourier_interpolate(HFfir, Ncoefs) h = fourier_interpolate(HFfir, Ncoefs)
print('> design_fir_low_pass_filter()') if verbosity:
print('> design_fir_low_pass_filter():')
print('. Method = %s' % method) print('. Method = %s' % method)
print('. Q = %f' % Q) print('. Q = %f' % Q)
print('. fpass = %f' % fpass) print('. fpass = %f' % fpass)
...@@ -287,6 +288,83 @@ def design_fir_low_pass_filter(method, ...@@ -287,6 +288,83 @@ def design_fir_low_pass_filter(method,
return h return h
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):
"""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.
"""
# Determine whether to adjust fpass or fstop
nofIterations = 1
nofIterationsMax = 100
h = design_fir_low_pass_filter(method, Ncoefs, fpass, fstop, rippleWeights=rippleWeights,
kaiserBeta=kaiserBeta, fs=fs, verbosity=0)
_, fn, HF = dtft(h)
f = fn * fs
fIndex, fValue, fGain = estimate_gain_at_frequency(f, HF, fcutoff)
if fGain < cutoffGain:
adjustBand = 0 # increase fpass
adjustBW = fcutoff - fpass
else:
adjustBand = 1 # decrease fstop
adjustBW = fstop - fcutoff
# Iterate to achieve cutoffGain at fcutoff
fRadix = 32
fStep = adjustBW / fRadix
gainResolution = cutoffGain * resolution
FP = fpass
FS = fstop
while True:
prevFP = FP
prevFS = FS
if adjustBand == 0:
if fGain < cutoffGain:
# Keep on increasing FP
FP += fStep
else:
# FP too high, one fStep back and increment FP with reduced fStep
FP -= fStep
fStep /= fRadix
FP += fStep
else:
if fGain > cutoffGain:
# Keep on decreasing FS
FS -= fStep
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 None
return h
def prototype_fir_low_pass_filter(method='firls', def prototype_fir_low_pass_filter(method='firls',
Npoints=1024, Ntaps=16, Ncoefs=1024*16, Npoints=1024, Ntaps=16, Ncoefs=1024*16,
hpFactor=0.9, transitionFactor=0.4, stopRippleFactor=1000000, kaiserBeta=0, fs=1.0): hpFactor=0.9, transitionFactor=0.4, stopRippleFactor=1000000, kaiserBeta=0, fs=1.0):
...@@ -344,12 +422,22 @@ def prototype_fir_low_pass_filter(method='firls', ...@@ -344,12 +422,22 @@ def prototype_fir_low_pass_filter(method='firls',
############################################################################### ###############################################################################
# Filterbank transfer function for prototype LPF # Filterbank
############################################################################### ###############################################################################
def derive_filterbank_transfer(HFprototype, Npoints): def filterbank_aggregate_transfer(HFprototype, Npoints):
Lprototype = len(HFprototype) """Aggregate frequency reponse of filterbank
Lbin = Lprototype / Npoints
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
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
HFbank = np.zeros(Lprototype, dtype='cfloat') HFbank = np.zeros(Lprototype, dtype='cfloat')
for bi in range(Npoints): for bi in range(Npoints):
HFbank += np.roll(HFprototype, bi * Lbin) HFbank += np.roll(HFprototype, bi * Lbin)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment