gr-fft: adding additional window options to win_type enum

Adding Nuttall, Welch, Parzen, Exponential, Riemann, Gaussian, and Tukey window options to the win_type enum for use with the various filter design functions

Signed-off-by: Jacob Gilbert <jacob.gilbert@protonmail.com>

3 files changed, 112 insertions, 33 deletions

diff --git a/gr-fft/include/gnuradio/fft/window.h b/gr-fft/include/gnuradio/fft/window.h
index 6739f9a0b..8ef85d758 100644
--- a/gr-fft/include/gnuradio/fft/window.h
+++ b/gr-fft/include/gnuradio/fft/window.h
@@ -22,10 +22,14 @@ namespace fft {
 class FFT_API window
 {
 public:
+    // illegal value for any window that requires a parameter
+    static constexpr double INVALID_WIN_PARAM = -1;
+
     enum win_type {
         WIN_NONE = -1,       //!< don't use a window
         WIN_HAMMING = 0,     //!< Hamming window; max attenuation 53 dB
         WIN_HANN = 1,        //!< Hann window; max attenuation 44 dB
+        WIN_HANNING = 1,     //!< alias to WIN_HANN
         WIN_BLACKMAN = 2,    //!< Blackman window; max attenuation 74 dB
         WIN_RECTANGULAR = 3, //!< Basic rectangular window; max attenuation 21 dB
         WIN_KAISER = 4, //!< Kaiser window; max attenuation see window::max_attenuation
@@ -34,30 +38,52 @@ public:
             5,            //!< alias to WIN_BLACKMAN_hARRIS for capitalization consistency
         WIN_BARTLETT = 6, //!< Barlett (triangular) window; max attenuation 26 dB
         WIN_FLATTOP = 7,  //!< flat top window; useful in FFTs; max attenuation 93 dB
+        WIN_NUTTALL = 8,  //!< Nuttall window; max attenuation 114 dB
+        WIN_BLACKMAN_NUTTALL = 8, //!< Nuttall window; max attenuation 114 dB
+        WIN_NUTTALL_CFD =
+            9, //!< Nuttall continuous-first-derivative window; max attenuation 112 dB
+        WIN_WELCH = 10,  //!< Welch window; max attenuation 31 dB
+        WIN_PARZEN = 11, //!< Parzen window; max attenuation 56 dB
+        WIN_EXPONENTIAL =
+            12, //!< Exponential window; max attenuation see window::max_attenuation
+        WIN_RIEMANN = 13, //!< Riemann window; max attenuation 39 dB
+        WIN_GAUSSIAN =
+            14,         //!< Gaussian window; max attenuation see window::max_attenuation
+        WIN_TUKEY = 15, //!< Tukey window; max attenuation see window::max_attenuation
     };
 
     /*!
      * \brief Given a window::win_type, this tells you the maximum
-     * attenuation you can expect.
+     * attenuation (really the maximum approximation error) you can expect.
      *
      * \details
-     * For most windows, this is a set value. For the Kaiser window,
-     * the attenuation is based on the value of beta. The actual
-     * relationship is a piece-wise exponential relationship to
-     * calculate beta from the desired attenuation and can be found
-     * on page 542 of Oppenheim and Schafer (Discrete-Time Signal
-     * Processing, 3rd edition). To simplify this function to solve
-     * for A given beta, we use a linear form that is exact for
-     * attenuation >= 50 dB.
+     * For most windows, this is a set value. For the Kaiser, Exponential, Gaussian, and
+     * Tukey windows, the attenuation is based on the value of a provided parameter.
+     *
+     * For the Kaiser window the actual relationship is a piece-wise exponential
+     * relationship to calculate beta from the desired attenuation and can be found on
+     * page 542 of Oppenheim and Schafer (Discrete-Time Signal Processing, 3rd edition).
+     * To simplify this function to solve for A given beta, we use a linear form that is
+     * exact for attenuation >= 50 dB. For an attenuation of 50 dB, beta = 4.55; for an
+     * attenuation of 70 dB, beta = 6.76.
+     *
+     * Exponential attenuation is complicated to measure due to the irregular error ripple
+     * structure, but it ranges between 23 and 26 dB depending on the decay factor; 26 dB
+     * is a good bound.
      *
-     * For an attenuation of 50 dB, beta = 4.55.
+     * The Gaussian window should not be used for window based filter construction;
+     * instead there is a dedicated gaussian filter construction fuction. There is no
+     * meaningful way to measure approximation error 'delta' as shown in Fig 7.23 of
+     * Oppenheim and Schafer (Discrete-Time Signal Processing, 3rd edition).
      *
-     * For an attenuation of 70 dB, beta = 6.76.
+     * Tukey windows provide attenuation that varies non-linearily between Rectangular (21
+     * dB) and Hann (44 dB) windows.
      *
      * \param type The window::win_type enumeration of the window type.
-     * \param beta Beta value only used for the Kaiser window.
+     * \param param Parameter value used for Kaiser (beta), Exponential (d), Gaussian
+     * (sigma) and Tukey (alpha) window creation.
      */
-    static double max_attenuation(win_type type, double beta = 6.76);
+    static double max_attenuation(win_type type, double param = INVALID_WIN_PARAM);
 
     /*!
      * \brief Helper function to build cosine-based windows. 3-coefficient version.
@@ -330,11 +356,11 @@ public:
      *
      * \param type a gr::fft::win_type index for the type of window.
      * \param ntaps Number of coefficients in the window.
-     * \param beta Used only for building Kaiser windows.
+     * \param param Parameter value used for Kaiser (beta), Exponential (d), Gaussian (sigma) and Tukey (alpha) window creation.
      * \param normalize If true, return a window with unit power
      */
     static std::vector<float>
-    build(win_type type, int ntaps, double beta = 6.76, const bool normalize = false);
+    build(win_type type, int ntaps, double param = INVALID_WIN_PARAM, const bool normalize = false);
 };
 
 } /* namespace fft */
diff --git a/gr-fft/lib/window.cc b/gr-fft/lib/window.cc
index faebdd4f5..98f0b897e 100644
--- a/gr-fft/lib/window.cc
+++ b/gr-fft/lib/window.cc
@@ -50,7 +50,7 @@ double freq(int ntaps) { return 2.0 * GR_M_PI / ntaps; }
 
 double rate(int ntaps) { return 1.0 / (ntaps >> 1); }
 
-double window::max_attenuation(win_type type, double beta)
+double window::max_attenuation(win_type type, double param)
 {
     switch (type) {
     case (WIN_HAMMING):
@@ -66,7 +66,7 @@ double window::max_attenuation(win_type type, double beta)
         return 21;
         break;
     case (WIN_KAISER):
-        return (beta / 0.1102 + 8.7);
+        return (param / 0.1102 + 8.7);
         break;
     case (WIN_BLACKMAN_hARRIS):
         return 92;
@@ -77,6 +77,33 @@ double window::max_attenuation(win_type type, double beta)
     case (WIN_FLATTOP):
         return 93;
         break;
+    case WIN_NUTTALL:
+        return 114;
+    case WIN_NUTTALL_CFD:
+        return 112;
+    case WIN_WELCH:
+        return 31;
+    case WIN_PARZEN:
+        return 56;
+    case WIN_EXPONENTIAL:
+        // varies slightly depending on the decay factor, but this is a safe return value
+        return 26;
+    case WIN_RIEMANN:
+        return 39;
+    case WIN_GAUSSIAN:
+        // not meaningful for gaussian windows, but return something reasonable
+        return 100;
+    case WIN_TUKEY:
+        // low end is a rectangular window, attenuation exponentialy approaches Hann
+        // piecewise linear estimate, determined empirically via curve fitting, median
+        // error is less than 0.5dB and maximum error is 2.5dB; the returned value will
+        // never be less than expected attenuation to ensure that window designed filters
+        // are never below expected quality.
+        if (param > 0.9)
+            return ((param - 0.9) * 135 + 30.5);
+        else if (param > 0.7)
+            return ((param - 0.6) * 20 + 24);
+        return (param * 5 + 21);
     default:
         throw std::out_of_range("window::max_attenuation: unknown window type provided.");
     }
@@ -358,12 +385,12 @@ std::vector<float> window::gaussian(int ntaps, float sigma)
 }
 
 std::vector<float>
-window::build(win_type type, int ntaps, double beta, const bool normalize)
+window::build(win_type type, int ntaps, double param, const bool normalize)
 {
     // If we want a normalized window, we get a non-normalized one first, then
     // normalize it here:
     if (normalize) {
-        auto win = build(type, ntaps, beta, false);
+        auto win = build(type, ntaps, param, false);
         const double pwr_acc = // sum(win**2) / len(win)
             std::accumulate(win.cbegin(),
                             win.cend(),
@@ -390,11 +417,27 @@ window::build(win_type type, int ntaps, double beta, const bool normalize)
     case WIN_BLACKMAN_hARRIS:
         return blackman_harris(ntaps);
     case WIN_KAISER:
-        return kaiser(ntaps, beta);
+        return kaiser(ntaps, param);
     case WIN_BARTLETT:
         return bartlett(ntaps);
     case WIN_FLATTOP:
         return flattop(ntaps);
+    case WIN_NUTTALL:
+        return nuttall(ntaps);
+    case WIN_NUTTALL_CFD:
+        return nuttall_cfd(ntaps);
+    case WIN_WELCH:
+        return welch(ntaps);
+    case WIN_PARZEN:
+        return parzen(ntaps);
+    case WIN_EXPONENTIAL:
+        return exponential(ntaps, param);
+    case WIN_RIEMANN:
+        return riemann(ntaps);
+    case WIN_GAUSSIAN:
+        return gaussian(ntaps, param);
+    case WIN_TUKEY:
+        return tukey(ntaps, param);
     default:
         throw std::out_of_range("window::build: type out of range");
     }
diff --git a/gr-fft/python/fft/bindings/window_python.cc b/gr-fft/python/fft/bindings/window_python.cc
index 3e9cdeec5..cc0fb6d58 100644
--- a/gr-fft/python/fft/bindings/window_python.cc
+++ b/gr-fft/python/fft/bindings/window_python.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright 2020 Free Software Foundation, Inc.
+ * Copyright 2020,2021 Free Software Foundation, Inc.
  *
  * This file is part of GNU Radio
  *
@@ -14,7 +14,7 @@
 /* BINDTOOL_GEN_AUTOMATIC(0)                                                       */
 /* BINDTOOL_USE_PYGCCXML(0)                                                        */
 /* BINDTOOL_HEADER_FILE(window.h)                                                  */
-/* BINDTOOL_HEADER_FILE_HASH(22de6d8875628eec777952b4902a09e9)                     */
+/* BINDTOOL_HEADER_FILE_HASH(de72e082a5bc1eeed7c4b3221025eb02)                     */
 /***********************************************************************************/
 
 #include <pybind11/complex.h>
@@ -34,15 +34,25 @@ void bind_window(py::module& m)
     py::class_<window, std::shared_ptr<window>> window_class(m, "window", D(window));
 
     py::enum_<gr::fft::window::win_type>(window_class, "win_type")
-        .value("WIN_HAMMING", gr::fft::window::WIN_HAMMING)                 // 0
-        .value("WIN_HANN", gr::fft::window::WIN_HANN)                       // 1
-        .value("WIN_BLACKMAN", gr::fft::window::WIN_BLACKMAN)               // 2
-        .value("WIN_RECTANGULAR", gr::fft::window::WIN_RECTANGULAR)         // 3
-        .value("WIN_KAISER", gr::fft::window::WIN_KAISER)                   // 4
-        .value("WIN_BLACKMAN_hARRIS", gr::fft::window::WIN_BLACKMAN_hARRIS) // 5
-        .value("WIN_BLACKMAN_HARRIS", gr::fft::window::WIN_BLACKMAN_HARRIS) // 5
-        .value("WIN_BARTLETT", gr::fft::window::WIN_BARTLETT)               // 6
-        .value("WIN_FLATTOP", gr::fft::window::WIN_FLATTOP)                 // 7
+        .value("WIN_HAMMING", gr::fft::window::WIN_HAMMING)                     // 0
+        .value("WIN_HANN", gr::fft::window::WIN_HANN)                           // 1
+        .value("WIN_HANNING", gr::fft::window::WIN_HANNING)                     // 1
+        .value("WIN_BLACKMAN", gr::fft::window::WIN_BLACKMAN)                   // 2
+        .value("WIN_RECTANGULAR", gr::fft::window::WIN_RECTANGULAR)             // 3
+        .value("WIN_KAISER", gr::fft::window::WIN_KAISER)                       // 4
+        .value("WIN_BLACKMAN_hARRIS", gr::fft::window::WIN_BLACKMAN_hARRIS)     // 5
+        .value("WIN_BLACKMAN_HARRIS", gr::fft::window::WIN_BLACKMAN_HARRIS)     // 5
+        .value("WIN_BARTLETT", gr::fft::window::WIN_BARTLETT)                   // 6
+        .value("WIN_FLATTOP", gr::fft::window::WIN_FLATTOP)                     // 7
+        .value("WIN_NUTTALL", gr::fft::window::WIN_NUTTALL)                     // 8
+        .value("WIN_BLACKMAN_NUTTALL", gr::fft::window::WIN_BLACKMAN_NUTTALL)   // 8
+        .value("WIN_NUTTALL_CFD", gr::fft::window::WIN_NUTTALL_CFD)             // 9
+        .value("WIN_WELCH", gr::fft::window::WIN_WELCH)                         // 10
+        .value("WIN_PARZEN", gr::fft::window::WIN_PARZEN)                       // 11
+        .value("WIN_EXPONENTIAL", gr::fft::window::WIN_EXPONENTIAL)             // 12
+        .value("WIN_RIEMANN", gr::fft::window::WIN_RIEMANN)                     // 13
+        .value("WIN_GAUSSIAN", gr::fft::window::WIN_GAUSSIAN)                   // 14
+        .value("WIN_TUKEY", gr::fft::window::WIN_TUKEY)                         // 15
         .export_values();
 
     py::implicitly_convertible<int, gr::fft::window::win_type>();
@@ -51,7 +61,7 @@ void bind_window(py::module& m)
         .def_static("max_attenuation",
                     &window::max_attenuation,
                     py::arg("type"),
-                    py::arg("beta") = 6.7599999999999998,
+                    py::arg("param") = 6.7599999999999998,
                     D(window, max_attenuation))
 
 
@@ -188,7 +198,7 @@ void bind_window(py::module& m)
                     &window::build,
                     py::arg("type"),
                     py::arg("ntaps"),
-                    py::arg("beta") = 6.76,
+                    py::arg("param") = 6.76,
                     py::arg("normalize") = false,
                     D(window, build))