[Commit-gnuradio] [gnuradio] 01/07: gr-analog: Fix FM preemphasis filter and rework deemphasis filter

[git]

This is an automated email from the git hooks/post-receive script.

jcorgan pushed a commit to branch master
in repository gnuradio.

commit 137ae3420b370de42cac2fcb7b1e212968f3a638
Author: Andy Walls <address@hidden>
Date:   Tue Mar 29 21:30:38 2016 -0400

    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.
---
 docs/exploring-gnuradio/fm_tx.grc   |   4 +
 gr-analog/examples/fmtest.py        |   3 +-
 gr-analog/grc/analog_fm_preemph.xml |   8 +-
 gr-analog/grc/analog_nbfm_tx.xml    |   8 ++
 gr-analog/grc/analog_wfm_tx.xml     |   7 +
 gr-analog/python/analog/fm_emph.py  | 267 +++++++++++++++++++++++++++++-------
 gr-analog/python/analog/nbfm_tx.py  |   5 +-
 gr-analog/python/analog/wfm_tx.py   |   5 +-
 gr-filter/examples/synth_to_chan.py |   2 +-
 gr-uhd/examples/python/fm_tx4.py    |   3 +-
 10 files changed, 258 insertions(+), 54 deletions(-)

diff --git a/docs/exploring-gnuradio/fm_tx.grc 
b/docs/exploring-gnuradio/fm_tx.grc
index 2f047bf..be8fe77 100644
--- a/docs/exploring-gnuradio/fm_tx.grc
+++ b/docs/exploring-gnuradio/fm_tx.grc
@@ -792,6 +792,10 @@
       <value>75e3</value>
     </param>
     <param>
+      <key>fh</key>
+      <value>0.0</value>
+    </param>
+    <param>
       <key>affinity</key>
       <value></value>
     </param>
diff --git a/gr-analog/examples/fmtest.py b/gr-analog/examples/fmtest.py
index 327da8e..22448e6 100755
--- a/gr-analog/examples/fmtest.py
+++ b/gr-analog/examples/fmtest.py
@@ -48,7 +48,8 @@ class fmtx(gr.hier_block2):
                                 gr.io_signature(1, 1, gr.sizeof_float),
                                 gr.io_signature(1, 1, gr.sizeof_gr_complex))
 
-        fmtx = analog.nbfm_tx(audio_rate, if_rate, max_dev=5e3, tau=75e-6)
+        fmtx = analog.nbfm_tx(audio_rate, if_rate, max_dev=5e3, tau=75e-6,
+                                                               fh=0.0)
 
         # Local oscillator
         lo = analog.sig_source_c(if_rate,            # sample rate
diff --git a/gr-analog/grc/analog_fm_preemph.xml 
b/gr-analog/grc/analog_fm_preemph.xml
index fb898b8..e5a8c35 100644
--- a/gr-analog/grc/analog_fm_preemph.xml
+++ b/gr-analog/grc/analog_fm_preemph.xml
@@ -8,7 +8,7 @@
        <name>FM Preemphasis</name>
        <key>analog_fm_preemph</key>
        <import>from gnuradio import analog</import>
-       <make>analog.fm_preemph(fs=$samp_rate, tau=$tau)</make>
+       <make>analog.fm_preemph(fs=$samp_rate, tau=$tau, fh=$fh)</make>
        <param>
                <name>Sample Rate</name>
                <key>samp_rate</key>
@@ -20,6 +20,12 @@
                <value>75e-6</value>
                <type>real</type>
        </param>
+       <param>
+               <name>High Corner Freq</name>
+               <key>fh</key>
+               <value>0.0</value>
+               <type>real</type>
+       </param>
        <sink>
                <name>in</name>
                <type>float</type>
diff --git a/gr-analog/grc/analog_nbfm_tx.xml b/gr-analog/grc/analog_nbfm_tx.xml
index ea13b8f..3e91e9c 100644
--- a/gr-analog/grc/analog_nbfm_tx.xml
+++ b/gr-analog/grc/analog_nbfm_tx.xml
@@ -13,6 +13,7 @@
        quad_rate=$quad_rate,
        tau=$tau,
        max_dev=$max_dev,
+       fh=$fh,
         )</make>
   <callback>set_max_deviation($max_dev)</callback>
 
@@ -42,6 +43,13 @@
     <type>real</type>
   </param>
 
+  <param>
+    <name>Preemphasis High Corner Freq</name>
+    <key>fh</key>
+    <value>0.0</value>
+    <type>real</type>
+  </param>
+
   <check>($quad_rate)%($audio_rate) == 0</check>
 
   <sink>
diff --git a/gr-analog/grc/analog_wfm_tx.xml b/gr-analog/grc/analog_wfm_tx.xml
index 89844f9..cd67d6a 100644
--- a/gr-analog/grc/analog_wfm_tx.xml
+++ b/gr-analog/grc/analog_wfm_tx.xml
@@ -13,6 +13,7 @@
        quad_rate=$quad_rate,
        tau=$tau,
        max_dev=$max_dev,
+       fh=$fh,
 )</make>
        <param>
                <name>Audio Rate</name>
@@ -36,6 +37,12 @@
                <value>75e3</value>
                <type>real</type>
        </param>
+       <param>
+               <name>Preemphasis High Corner Freq</name>
+               <key>fh</key>
+               <value>0.0</value>
+               <type>real</type>
+       </param>
        <check>($quad_rate)%($audio_rate) == 0</check>
        <sink>
                <name>in</name>
diff --git a/gr-analog/python/analog/fm_emph.py 
b/gr-analog/python/analog/fm_emph.py
index d2a38d4..e457f56 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
diff --git a/gr-analog/python/analog/nbfm_tx.py 
b/gr-analog/python/analog/nbfm_tx.py
index ffd539e..cd11c2f 100644
--- a/gr-analog/python/analog/nbfm_tx.py
+++ b/gr-analog/python/analog/nbfm_tx.py
@@ -29,7 +29,7 @@ except ImportError:
     import analog_swig as analog
 
 class nbfm_tx(gr.hier_block2):
-    def __init__(self, audio_rate, quad_rate, tau=75e-6, max_dev=5e3):
+    def __init__(self, audio_rate, quad_rate, tau=75e-6, max_dev=5e3, fh=0.0):
         """
         Narrow Band FM Transmitter.
 
@@ -41,6 +41,7 @@ class nbfm_tx(gr.hier_block2):
             quad_rate: sample rate of output stream (integer)
             tau: preemphasis time constant (default 75e-6) (float)
             max_dev: maximum deviation in Hz (default 5e3) (float)
+            fh: high frequency at which to flatten preemphasis; 0.0 means none 
(float)
 
         quad_rate must be an integer multiple of audio_rate.
         """
@@ -71,7 +72,7 @@ class nbfm_tx(gr.hier_block2):
             #print "len(interp_taps) =", len(interp_taps)
             self.interpolator = filter.interp_fir_filter_fff (interp_factor, 
interp_taps)
 
-        self.preemph = fm_preemph(quad_rate, tau=tau)
+        self.preemph = fm_preemph(quad_rate, tau=tau, fh=fh)
 
         k = 2 * math.pi * max_dev / quad_rate
         self.modulator = analog.frequency_modulator_fc(k)
diff --git a/gr-analog/python/analog/wfm_tx.py 
b/gr-analog/python/analog/wfm_tx.py
index be66231..7363ccd 100644
--- a/gr-analog/python/analog/wfm_tx.py
+++ b/gr-analog/python/analog/wfm_tx.py
@@ -30,7 +30,7 @@ except ImportError:
     import analog_swig as analog
 
 class wfm_tx(gr.hier_block2):
-    def __init__(self, audio_rate, quad_rate, tau=75e-6, max_dev=75e3):
+    def __init__(self, audio_rate, quad_rate, tau=75e-6, max_dev=75e3, fh=0.0):
         """
         Wide Band FM Transmitter.
 
@@ -42,6 +42,7 @@ class wfm_tx(gr.hier_block2):
             quad_rate: sample rate of output stream (integer)
             tau: preemphasis time constant (default 75e-6) (float)
             max_dev: maximum deviation in Hz (default 75e3) (float)
+            fh: high frequency at which to flatten preemphasis; 0.0 means none 
(float)
 
         quad_rate must be an integer multiple of audio_rate.
         """
@@ -71,7 +72,7 @@ class wfm_tx(gr.hier_block2):
             print "len(interp_taps) =", len(interp_taps)
             self.interpolator = filter.interp_fir_filter_fff (interp_factor, 
interp_taps)
 
-        self.preemph = fm_preemph(quad_rate, tau=tau)
+        self.preemph = fm_preemph(quad_rate, tau=tau, fh=fh)
 
         k = 2 * math.pi * max_dev / quad_rate
         self.modulator = analog.frequency_modulator_fc (k)
diff --git a/gr-filter/examples/synth_to_chan.py 
b/gr-filter/examples/synth_to_chan.py
index 9e68202..7a295ea 100755
--- a/gr-filter/examples/synth_to_chan.py
+++ b/gr-filter/examples/synth_to_chan.py
@@ -54,7 +54,7 @@ def main():
     fmtx = list()
     for fi in freqs:
         s = analog.sig_source_f(fs, analog.GR_SIN_WAVE, fi, 1)
-        fm = analog.nbfm_tx(fs, 4*fs, max_dev=10000, tau=75e-6)
+        fm = analog.nbfm_tx(fs, 4*fs, max_dev=10000, tau=75e-6, fh=0.0)
         sigs.append(s)
         fmtx.append(fm)
 
diff --git a/gr-uhd/examples/python/fm_tx4.py b/gr-uhd/examples/python/fm_tx4.py
index fefa678..ffc74ab 100755
--- a/gr-uhd/examples/python/fm_tx4.py
+++ b/gr-uhd/examples/python/fm_tx4.py
@@ -63,7 +63,8 @@ class pipeline(gr.hier_block2):
             sys.exit(1)
 
         print audio_rate, if_rate
-        fmtx = analog.nbfm_tx(audio_rate, if_rate, max_dev=5e3, tau=75e-6)
+        fmtx = analog.nbfm_tx(audio_rate, if_rate, max_dev=5e3, tau=75e-6,
+                                                               fh=0.0)
 
         # Local oscillator
         lo = analog.sig_source_c(if_rate,            # sample rate