From 515e40039e614ac86fb4c7be8dc4d067aa2bd9a3 Mon Sep 17 00:00:00 2001
From: Andy Walls <awalls@md.metrocast.net>
Date: Wed, 30 Mar 2016 09:15:45 -0400
Subject: gr-analog: Add safety and default for FM preemphasis filter
The FM preemphasis filter design now precludes the user from
inducing a pole on the unit circle at z = -1.0 and z = 1.0.
A pole at either of these locations makes the filter unstable and
useless: feeding back "+/-inf" into an IIR filter has no good
recovery.
Also provide a reasonable, maximally safe default of 0.925*fs/2.0
for the high frequency corner, fh. This keeps the slope of the
preemphasis filter looking reasonable sane in the whole band; at
least for tau=75e-6 and fs=48000.
---
gr-analog/python/analog/fm_emph.py | 79 +++++++++++++++-----------------------
1 file changed, 31 insertions(+), 48 deletions(-)
(limited to 'gr-analog/python/analog/fm_emph.py')
diff --git a/gr-analog/python/analog/fm_emph.py b/gr-analog/python/analog/fm_emph.py
index 7637743d91..bfa4742ace 100644
--- a/gr-analog/python/analog/fm_emph.py
+++ b/gr-analog/python/analog/fm_emph.py
@@ -251,69 +251,52 @@ class fm_preemph(gr.hier_block2):
"""
FM Preemphasis IIR filter.
"""
- def __init__(self, fs, tau=75e-6, fh=0.0):
+ def __init__(self, fs, tau=75e-6, fh=-1.0):
"""
Args:
fs: sampling frequency in Hz (float)
tau: Time constant in seconds (75us in US, 50us in EUR) (float)
- fh: High frequency at which to flatten out; 0.0 means none (float)
+ fh: High frequency at which to flatten out (< 0 means default of 0.925*fs/2.0) (float)
"""
gr.hier_block2.__init__(self, "fm_preemph",
gr.io_signature(1, 1, gr.sizeof_float), # Input signature
gr.io_signature(1, 1, gr.sizeof_float)) # Output signature
- if fh > 0.0 and fh < (fs / 2.0):
- # Digital corner frequencies
- w_cl = 1.0 / tau
- w_ch = 2.0 * math.pi * fh
+ # Set fh to something sensible, if needed.
+ # N.B. fh == fs/2.0 or fh == 0.0 results in a pole on the unit circle
+ # at z = -1.0 or z = 1.0 respectively. That makes the filter unstable
+ # and useless.
+ if fh <= 0.0 or fh >= fs/2.0:
+ fh = 0.925 * fs/2.0
- # Prewarped analog corner frequencies
- w_cla = 2.0 * fs * math.tan(w_cl / (2.0 * fs))
- w_cha = 2.0 * fs * math.tan(w_ch / (2.0 * fs))
+ # Digital corner frequencies
+ w_cl = 1.0 / tau
+ w_ch = 2.0 * math.pi * fh
- # Resulting digital pole, zero, and gain term from the bilinear
- # transformation of H(s) = (s + w_cla) / (s + w_cha) to
- # H(z) = b0 (1 - z1 z^-1)/(1 - p1 z^-1)
- kl = -w_cla / (2.0 * fs)
- kh = -w_cha / (2.0 * fs)
- z1 = (1.0 + kl) / (1.0 - kl)
- p1 = (1.0 + kh) / (1.0 - kh)
- b0 = (1.0 - kl) / (1.0 - kh)
+ # Prewarped analog corner frequencies
+ w_cla = 2.0 * fs * math.tan(w_cl / (2.0 * fs))
+ w_cha = 2.0 * fs * math.tan(w_ch / (2.0 * fs))
- # Since H(s = infinity) = 1.0, then H(z = -1) = 1.0 and
- # this filter has 0 dB gain at fs/2.0.
- # That isn't what users are going to expect, so adjust with a
- # gain, g, so that H(z = 1) = 1.0 for 0 dB gain at DC.
- w_0dB = 2.0 * math.pi * 0.0
- g = abs(1.0 - p1 * cmath.rect(1.0, -w_0dB)) \
- / (b0 * abs(1.0 - z1 * cmath.rect(1.0, -w_0dB)))
+ # Resulting digital pole, zero, and gain term from the bilinear
+ # transformation of H(s) = (s + w_cla) / (s + w_cha) to
+ # H(z) = b0 (1 - z1 z^-1)/(1 - p1 z^-1)
+ kl = -w_cla / (2.0 * fs)
+ kh = -w_cha / (2.0 * fs)
+ z1 = (1.0 + kl) / (1.0 - kl)
+ p1 = (1.0 + kh) / (1.0 - kh)
+ b0 = (1.0 - kl) / (1.0 - kh)
- btaps = [ g * b0 * 1.0, g * b0 * -z1 ]
- ataps = [ 1.0, -p1 ]
+ # Since H(s = infinity) = 1.0, then H(z = -1) = 1.0 and
+ # this filter has 0 dB gain at fs/2.0.
+ # That isn't what users are going to expect, so adjust with a
+ # gain, g, so that H(z = 1) = 1.0 for 0 dB gain at DC.
+ w_0dB = 2.0 * math.pi * 0.0
+ g = abs(1.0 - p1 * cmath.rect(1.0, -w_0dB)) \
+ / (b0 * abs(1.0 - z1 * cmath.rect(1.0, -w_0dB)))
- else:
- # Just use H(s) = (s + 1/RC)/(1/RC) as the transfer function
-
- # Digital corner frequencies
- w_cl = 1.0 / tau
-
- # Prewarped analog corner frequencies
- w_cla = 2.0 * fs * math.tan(w_cl / (2.0 * fs))
-
- # Resulting digital pole, zero, and gain term from the bilinear
- # transformation of H(s) = (s + w_cl)/w_cl to
- # H(z) = b0 (1 - z1 z^-1)/(1 - p1 z^-1)
- kl = -w_cla / (2.0 * fs)
- z1 = (1.0 + kl) / (1.0 - kl)
- p1 = -1.0
- b0 = (1.0 - kl) / -kl
-
- # Since H(s = 0) = 1.0, then H(z = 1) = 1.0 and
- # has 0 dB gain at DC
-
- btaps = [ b0 * 1.0, b0 * -z1 ]
- ataps = [ 1.0, -p1 ]
+ btaps = [ g * b0 * 1.0, g * b0 * -z1 ]
+ ataps = [ 1.0, -p1 ]
if 0:
print "btaps =", btaps
--
cgit v1.2.3