diff options
| author | Andy Walls <awalls@md.metrocast.net> | 2016-03-29 21:30:38 -0400 |
|---|---|---|
| committer | Andy Walls <awalls@md.metrocast.net> | 2016-03-29 21:30:38 -0400 |
| commit | 137ae3420b370de42cac2fcb7b1e212968f3a638 (patch) | |
| tree | 3988128f58bc6e32f8ae25ae39baaf6c0b7c4a6e /gr-analog/python/analog/fm_emph.py | |
| parent | 0a185537075ceb3a985f3a5cf198cb3b76c039a2 (diff) | |
gr-analog: Fix FM preemphasis filter and rework deemphasis filter
Add working filters designs for the FM preemphasis Tx filter.
Rework the FM deemphasis Rx filter as it was easier to rederive
the transfer function, than to determine if the one in use was correct.
Diffstat (limited to 'gr-analog/python/analog/fm_emph.py')
| -rw-r--r-- | gr-analog/python/analog/fm_emph.py | 267 |
1 files changed, 221 insertions, 46 deletions
diff --git a/gr-analog/python/analog/fm_emph.py b/gr-analog/python/analog/fm_emph.py index d2a38d4aff..e457f56dc0 100644 --- a/gr-analog/python/analog/fm_emph.py +++ b/gr-analog/python/analog/fm_emph.py @@ -21,17 +21,78 @@ from gnuradio import gr, filter import math +import cmath # -# 1 -# H(s) = ------- -# 1 + s +# An analog deemphasis filter: # -# tau is the RC time constant. -# critical frequency: w_p = 1/tau +# R +# o------/\/\/\/---+----o +# | +# = C +# | +# --- # -# We prewarp and use the bilinear z-transform to get our IIR coefficients. -# See "Digital Signal Processing: A Practical Approach" by Ifeachor and Jervis +# Has this transfer function: +# +# 1 1 +# ---- --- +# RC tau +# H(s) = ---------- = ---------- +# 1 1 +# s + ---- s + --- +# RC tau +# +# And has its -3 dB response, due to the pole, at +# +# |H(j w_c)|^2 = 1/2 => s = j w_c = j (1/(RC)) +# +# Historically, this corner frequency of analog audio deemphasis filters +# been specified by the RC time constant used, called tau. +# So w_c = 1/tau. +# +# FWIW, for standard tau values, some standard analog components would be: +# tau = 75 us = (50K)(1.5 nF) = (50 ohms)(1.5 uF) +# tau = 50 us = (50K)(1.0 nF) = (50 ohms)(1.0 uF) +# +# In specifying tau for this digital deemphasis filter, tau specifies +# the *digital* corner frequency, w_c, desired. +# +# The digital deemphasis filter design below, uses the +# "bilinear transformation" method of designing digital filters: +# +# 1. Convert digitial specifications into the analog domain, by prewarping +# digital frequency specifications into analog frequencies. +# +# w_a = (2/T)tan(wT/2) +# +# 2. Use an analog filter design technique to design the filter. +# +# 3. Use the bilinear transformation to convert the analog filter design to a +# digital filter design. +# +# H(z) = H(s)| +# s = (2/T)(1-z^-1)/(1+z^-1) +# +# +# w_ca 1 1 - (-1) z^-1 +# H(z) = ---- * ----------- * ----------------------- +# 2 fs -w_ca -w_ca +# 1 - ----- 1 + ----- +# 2 fs 2 fs +# 1 - ----------- z^-1 +# -w_ca +# 1 - ----- +# 2 fs +# +# We use this design technique, because it is an easy way to obtain a filter +# design with the -6 dB/octave roll-off required of the deemphasis filter. +# +# Jackson, Leland B., _Digital_Filters_and_Signal_Processing_Second_Edition_, +# Kluwer Academic Publishers, 1989, pp 201-212 +# +# Orfanidis, Sophocles J., _Introduction_to_Signal_Processing_, Prentice Hall, +# 1996, pp 573-583 # @@ -51,15 +112,24 @@ class fm_deemph(gr.hier_block2): gr.io_signature(1, 1, gr.sizeof_float), # Input signature gr.io_signature(1, 1, gr.sizeof_float)) # Output signature - w_p = 1/tau - w_pp = math.tan(w_p / (fs * 2)) # prewarped analog freq + # Digital corner frequency + w_c = 1.0 / tau + + # Prewarped analog corner frequency + w_ca = 2.0 * fs * math.tan(w_c / (2.0 * fs)) - a1 = (w_pp - 1)/(w_pp + 1) - b0 = w_pp/(1 + w_pp) - b1 = b0 + # Resulting digital pole, zero, and gain term from the bilinear + # transformation of H(s) = w_ca / (s + w_ca) to + # H(z) = b0 (1 - z1 z^-1)/(1 - p1 z^-1) + k = -w_ca / (2.0 * fs) + z1 = -1.0 + p1 = (1.0 + k) / (1.0 - k) + b0 = -k / (1.0 - k) - btaps = [b0, b1] - ataps = [1, a1] + btaps = [ b0 * 1.0, b0 * -z1 ] + ataps = [ 1.0, -p1 ] + + # Since H(s = 0) = 1.0, then H(z = 1) = 1.0 and has 0 dB gain at DC if 0: print "btaps =", btaps @@ -71,22 +141,13 @@ class fm_deemph(gr.hier_block2): self.connect(self, deemph, self) # -# 1 + s*t1 -# H(s) = ---------- -# 1 + s*t2 -# -# I think this is the right transfer function. -# +# An analog preemphasis filter, that flattens out again at the high end: # -# This fine ASCII rendition is based on Figure 5-15 -# in "Digital and Analog Communication Systems", Leon W. Couch II -# -# -# R1 +# C # +-----||------+ # | | # o------+ +-----+--------o -# | C1 | | +# | R1 | | # +----/\/\/\/--+ \ # / # \ R2 @@ -95,28 +156,94 @@ class fm_deemph(gr.hier_block2): # | # o--------------------------+--------o # -# f1 = 1/(2*pi*t1) = 1/(2*pi*R1*C) +# (This fine ASCII rendition is based on Figure 5-15 +# in "Digital and Analog Communication Systems", Leon W. Couch II) +# +# Has this transfer function: +# +# 1 +# s + --- +# R1C +# H(s) = ------------------ +# 1 R1 +# s + --- (1 + --) +# R1C R2 +# +# +# It has a corner due to the numerator, where the rise starts, at # -# 1 R1 + R2 -# f2 = ------- = ------------ -# 2*pi*t2 2*pi*R1*R2*C +# |Hn(j w_cl)|^2 = 2*|Hn(0)|^2 => s = j w_cl = j (1/(R1C)) # -# t1 is 75us in US, 50us in EUR -# f2 should be higher than our audio bandwidth. +# It has a corner due to the denominator, where it levels off again, at # +# |Hn(j w_ch)|^2 = 1/2*|Hd(0)|^2 => s = j w_ch = j (1/(R1C) * (1 + R1/R2)) # -# The Bode plot looks like this: +# Historically, the corner frequency of analog audio preemphasis filters +# been specified by the R1C time constant used, called tau. # +# So +# w_cl = 1/tau = 1/R1C; f_cl = 1/(2*pi*tau) = 1/(2*pi*R1*C) +# w_ch = 1/tau2 = (1+R1/R2)/R1C; f_ch = 1/(2*pi*tau2) = (1+R1/R2)/(2*pi*R1*C) # -# /---------------- -# / -# / <-- slope = 20dB/decade -# / -# -------------/ -# f1 f2 +# and note f_ch = f_cl * (1 + R1/R2). # -# We prewarp and use the bilinear z-transform to get our IIR coefficients. -# See "Digital Signal Processing: A Practical Approach" by Ifeachor and Jervis +# For broadcast FM audio, tau is 75us in the United States and 50us in Europe. +# f_ch should be higher than our digital audio bandwidth. +# +# The Bode plot looks like this: +# +# +# /---------------- +# / +# / <-- slope = 20dB/decade +# / +# -------------/ +# f_cl f_ch +# +# In specifying tau for this digital preemphasis filter, tau specifies +# the *digital* corner frequency, w_cl, desired. +# +# The digital preemphasis filter design below, uses the +# "bilinear transformation" method of designing digital filters: +# +# 1. Convert digitial specifications into the analog domain, by prewarping +# digital frequency specifications into analog frequencies. +# +# w_a = (2/T)tan(wT/2) +# +# 2. Use an analog filter design technique to design the filter. +# +# 3. Use the bilinear transformation to convert the analog filter design to a +# digital filter design. +# +# H(z) = H(s)| +# s = (2/T)(1-z^-1)/(1+z^-1) +# +# +# -w_cla +# 1 + ------ +# 2 fs +# 1 - ------------ z^-1 +# -w_cla -w_cla +# 1 - ------ 1 - ------ +# 2 fs 2 fs +# H(z) = ------------ * ----------------------- +# -w_cha -w_cha +# 1 - ------ 1 + ------ +# 2 fs 2 fs +# 1 - ------------ z^-1 +# -w_cha +# 1 - ------ +# 2 fs +# +# We use this design technique, because it is an easy way to obtain a filter +# design with the 6 dB/octave rise required of the premphasis filter. +# +# Jackson, Leland B., _Digital_Filters_and_Signal_Processing_Second_Edition_, +# Kluwer Academic Publishers, 1989, pp 201-212 +# +# Orfanidis, Sophocles J., _Introduction_to_Signal_Processing_, Prentice Hall, +# 1996, pp 573-583 # @@ -124,21 +251,69 @@ class fm_preemph(gr.hier_block2): """ FM Preemphasis IIR filter. """ - def __init__(self, fs, tau=75e-6): + def __init__(self, fs, tau=75e-6, fh=0.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) """ - gr.hier_block2.__init__(self, "fm_deemph", + 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 - # FIXME make this compute the right answer + if fh > 0.0 and fh < (fs / 2.0): + # Digital corner frequencies + w_cl = 1.0 / tau + w_ch = 2.0 * math.pi * fh + + # 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)) + + # 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) + + # 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))) + + btaps = [ g * b0 * 1.0, g * b0 * -z1 ] + ataps = [ 1.0, -p1 ] + + 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 = [1] - ataps = [1] + btaps = [ b0 * 1.0, b0 * -z1 ] + ataps = [ 1.0, -p1 ] if 0: print "btaps =", btaps |