From 925839ffdfabecc5a38f7e0ede0a41877f7d5032 Mon Sep 17 00:00:00 2001
From: Eric Kooistra <kooistra@astron.nl>
Date: Wed, 17 Jan 2024 11:10:27 +0100
Subject: [PATCH] Add iir_biquad_butter_coefficients(), for biquads used in
TBB.
---
applications/lofar2/model/pfb_os/dsp.py | 41 ++++++++++++++++++++++---
1 file changed, 37 insertions(+), 4 deletions(-)
diff --git a/applications/lofar2/model/pfb_os/dsp.py b/applications/lofar2/model/pfb_os/dsp.py
index e3a8a068fb..e9e99dac0b 100644
--- a/applications/lofar2/model/pfb_os/dsp.py
+++ b/applications/lofar2/model/pfb_os/dsp.py
@@ -117,7 +117,7 @@ def impulse_at_zero_crossing(x):
###############################################################################
-# Filter design
+# FIR Filter design
###############################################################################
def nof_taps_kaiser_window(fs, fpass, fstop, atten_db):
@@ -157,7 +157,7 @@ def ideal_low_pass_filter(Npoints, Npass, bandEdgeGain=1.0):
# Magnitude frequency reponse
HF = np.zeros([Npoints])
HF[0] = bandEdgeGain
- HF[1 : Npass - 1] = 1
+ HF[1: Npass - 1] = 1
HF[Npass - 1] = bandEdgeGain
# Zero center HF to make it even
HF = np.roll(HF, -(Npass // 2))
@@ -222,6 +222,39 @@ def fourier_interpolate(HFfilter, Ncoefs):
return hInterpolated.real
+###############################################################################
+# IIR Filter design
+###############################################################################
+
+def iir_biquad_butter_coefficients(f0, BW, band='bandstop', fs=1):
+ """Calculate b0, b1, b2, a0, a1, a2 for biquad IIR filter section.
+
+ Bilinear transform of Butterworth filters. From
+ https://webaudio.github.io/Audio-EQ-Cookbook/Audio-EQ-Cookbook.txt
+
+ For single biquad bandpass alternatively use signal.iirpeak()
+ For single biquad bandstop alternatively use signal.iirnotch()
+ """
+ w0 = 2 * np.pi * f0 / fs
+ Q = f0 / BW
+ si = np.sin(w0)
+ cs = np.cos(w0)
+ alpha = np.sin(w0) / (2 * Q)
+ a1 = -2 * cs
+ if band == 'lowpass':
+ b0 = 1 - cs
+ b = [b0 / 2, b0, b0 / 2]
+ elif band == 'highpass':
+ b0 = 1 + cs
+ b = [b0 / 2, -b0, b0 / 2]
+ elif band == 'bandstop':
+ b = [1, a1, 1]
+ elif band == 'bandpass':
+ b = [si / 2, 0, -si / 2]
+ a = [1 + alpha, a1, 1 - alpha]
+ return b, a
+
+
###############################################################################
# Hilbert transform filter
###############################################################################
@@ -475,7 +508,7 @@ def plot_spectra(f, HF, fs=1.0, fLim=None, dbLim=None):
plt.grid(True)
-def plot_power_spectrum(f, HF, fmt='r', fs=1.0, fLim=None, dbLim=None, showRoll=False):
+def plot_power_spectrum(f, HF, fmt='r', fs=1.0, fLim=None, dbLim=None):
"""Plot power spectrum
Input:
@@ -485,7 +518,7 @@ def plot_power_spectrum(f, HF, fmt='r', fs=1.0, fLim=None, dbLim=None, showRoll=
. fs: sample frequency in Hz, scale f by fs, fs >= 1
"""
if fs > 1:
- flabel = 'Frequency [fs / %d]' % fs
+ flabel = 'Frequency'
else:
flabel = 'Frequency [fs]'
--
GitLab