1 files changed, 333 insertions, 0 deletions
diff --git a/gr-fft/lib/window.cc b/gr-fft/lib/window.cc
new file mode 100644
index 0000000000..5a28b3239d
--- /dev/null
+++ b/gr-fft/lib/window.cc
@@ -0,0 +1,333 @@
+/* -*- c++ -*- */
+/*
+ * Copyright 2002,2007,2008,2012,2013 Free Software Foundation, Inc.
+ *
+ * This file is part of GNU Radio
+ *
+ * GNU Radio is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3, or (at your option)
+ * any later version.
+ *
+ * GNU Radio is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with GNU Radio; see the file COPYING.  If not, write to
+ * the Free Software Foundation, Inc., 51 Franklin Street,
+ * Boston, MA 02110-1301, USA.
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <gnuradio/fft/window.h>
+#include <stdexcept>
+
+namespace gr {
+  namespace fft {
+
+#define IzeroEPSILON 1E-21               /* Max error acceptable in Izero */
+    
+    static double Izero(double x)
+    {
+      double sum, u, halfx, temp;
+      int n;
+      
+      sum = u = n = 1;
+      halfx = x/2.0;
+      do {
+        temp = halfx/(double)n;
+        n += 1;
+        temp *= temp;
+        u *= temp;
+        sum += u;
+      } while (u >= IzeroEPSILON*sum);
+      return(sum);
+    }
+
+    double midn(int ntaps)
+    {
+      return ntaps/2.0;
+    }
+
+    double midm1(int ntaps)
+    {
+      return (ntaps - 1.0)/2.0;
+    }
+
+    double midp1(int ntaps)
+    {
+      return (ntaps + 1.0)/2.0;
+    }
+
+    double freq(int ntaps)
+    {
+      return 2.0*M_PI/ntaps;
+    }
+
+    double rate(int ntaps)
+    {
+      return 1.0/(ntaps >> 1);
+    }
+
+    std::vector<float>
+    window::coswindow(int ntaps, float c0, float c1, float c2)
+    {
+      std::vector<float> taps(ntaps);
+      float M = static_cast<float>(ntaps - 1);
+
+      for(int n = 0; n < ntaps; n++)
+        taps[n] = c0 - c1*cos((2*M_PI*n)/M) + c2*cos((4*M_PI*n)/M);
+      return taps;
+    }
+
+    std::vector<float>
+    window::coswindow(int ntaps, float c0, float c1, float c2, float c3)
+    {
+      std::vector<float> taps(ntaps);
+      float M = static_cast<float>(ntaps - 1);
+
+      for(int n = 0; n < ntaps; n++)
+        taps[n] = c0 - c1*cos((2*M_PI*n)/M) + c2*cos((4*M_PI*n)/M) - c3*cos((6*M_PI*n)/M);
+      return taps;
+    }
+
+    std::vector<float>
+    window::coswindow(int ntaps, float c0, float c1, float c2, float c3, float c4)
+    {
+      std::vector<float> taps(ntaps);
+      float M = static_cast<float>(ntaps - 1);
+
+      for(int n = 0; n < ntaps; n++)
+        taps[n] = c0 - c1*cos((2*M_PI*n)/M) + c2*cos((4*M_PI*n)/M) - c3*cos((6*M_PI*n)/M) + c4*cos((8*M_PI*n)/M);
+      return taps;
+    }
+
+    std::vector<float>
+    window::rectangular(int ntaps)
+    {
+      std::vector<float> taps(ntaps);
+      for(int n = 0; n < ntaps; n++)
+        taps[n] = 1;
+      return taps;
+    }
+
+    std::vector<float>
+    window::hamming(int ntaps)
+    {
+      std::vector<float> taps(ntaps);
+      float M = static_cast<float>(ntaps - 1);
+
+      for(int n = 0; n < ntaps; n++)
+        taps[n] = 0.54 - 0.46 * cos((2 * M_PI * n) / M);
+      return taps;
+    }
+
+    std::vector<float>
+    window::hann(int ntaps)
+    {
+      std::vector<float> taps(ntaps);
+      float M = static_cast<float>(ntaps - 1);
+
+      for(int n = 0; n < ntaps; n++)
+        taps[n] = 0.5 - 0.5 * cos((2 * M_PI * n) / M);
+      return taps;
+    }
+
+    std::vector<float>
+    window::hanning(int ntaps)
+    {
+      return hann(ntaps);
+    }
+
+    std::vector<float>
+    window::blackman(int ntaps)
+    {
+      return coswindow(ntaps, 0.42, 0.5, 0.08);
+    }
+
+    std::vector<float>
+    window::blackman2(int ntaps)
+    {
+      return coswindow(ntaps, 0.34401, 0.49755, 0.15844);
+    }
+
+    std::vector<float>
+    window::blackman3(int ntaps)
+    {
+      return coswindow(ntaps, 0.21747, 0.45325, 0.28256, 0.04672);
+    }
+    
+    std::vector<float>
+    window::blackman4(int ntaps)
+    {
+      return coswindow(ntaps, 0.084037, 0.29145, 0.375696, 0.20762, 0.041194);
+    }
+
+    std::vector<float>
+    window::blackman_harris(int ntaps, int atten)
+    {
+      switch(atten) {
+      case(61): return coswindow(ntaps, 0.42323, 0.49755, 0.07922);
+      case(67): return coswindow(ntaps, 0.44959, 0.49364, 0.05677);
+      case(74): return coswindow(ntaps, 0.40271, 0.49703, 0.09392, 0.00183);
+      case(92): return coswindow(ntaps, 0.35875, 0.48829, 0.14128, 0.01168);
+      default:
+        throw std::out_of_range("window::blackman_harris: unknown attenuation value (must be 61, 67, 74, or 92)");
+      }
+    }
+
+    std::vector<float>
+    window::blackmanharris(int ntaps, int atten)
+    {
+      return blackman_harris(ntaps, atten);
+    }
+
+    std::vector<float>
+    window::nuttal(int ntaps)
+    {
+      return coswindow(ntaps, 0.3635819, 0.4891775, 0.1365995, 0.0106411);
+    }
+
+    std::vector<float>
+    window::blackman_nuttal(int ntaps)
+    {
+      return nuttal(ntaps);
+    }
+
+    std::vector<float>
+    window::nuttal_cfd(int ntaps)
+    {
+      return coswindow(ntaps, 0.355768, 0.487396, 0.144232, 0.012604);
+    }
+
+    std::vector<float>
+    window::flattop(int ntaps)
+    {
+      double scale = 4.63867;
+      return coswindow(ntaps, 1.0/scale, 1.93/scale, 1.29/scale, 0.388/scale, 0.028/scale);
+    }
+
+    std::vector<float>
+    window::kaiser(int ntaps, double beta)
+    {
+      if(beta < 0)
+        throw std::out_of_range("window::kaiser: beta must be >= 0");
+
+      std::vector<float> taps(ntaps);
+
+      double IBeta = 1.0/Izero(beta);
+      double inm1 = 1.0/((double)(ntaps-1));
+      double temp;
+
+      for(int i = 0; i <= ntaps; i++) {
+        temp = 2*i*inm1 - 1;
+        taps[i] = Izero(beta*sqrt(1.0-temp*temp)) * IBeta;
+      }
+      return taps;
+    }
+
+    std::vector<float>
+    window::bartlett(int ntaps)
+    {
+      std::vector<float> taps(ntaps);
+      float M = static_cast<float>(ntaps - 1);
+
+      for(int n = 0; n < ntaps/2; n++)
+        taps[n] = 2*n/M;
+      for(int n = ntaps/2; n < ntaps; n++)
+        taps[n] = 2 - 2*n/M;
+
+      return taps;
+    }
+
+    std::vector<float>
+    window::welch(int ntaps)
+    {
+      std::vector<float> taps(ntaps);
+      double m1 = midm1(ntaps);
+      double p1 = midp1(ntaps);
+      for(int i = 0; i < midn(ntaps)+1; i++) {
+        taps[i] = 1.0 - pow((i - m1) / p1, 2);
+        taps[ntaps-i-1] = taps[i];
+      }
+      return taps;
+    }
+
+    std::vector<float>
+    window::parzen(int ntaps)
+    {
+      std::vector<float> taps(ntaps);
+      double m1 = midm1(ntaps);
+      double m = midn(ntaps);
+      int i;
+      for(i = ntaps/4; i < 3*ntaps/4; i++) {
+        taps[i] = 1.0 - 6.0*pow((i-m1)/m, 2.0) * (1.0 - fabs(i-m1)/m);
+      }
+      for(; i < ntaps; i++) {
+        taps[i] = 2.0 * pow(1.0 - fabs(i-m1)/m, 3.0);
+        taps[ntaps-i-1] = taps[i];
+      }
+      return taps;
+    }
+
+    std::vector<float>
+    window::exponential(int ntaps, double d)
+    {
+      if(d < 0)
+        throw std::out_of_range("window::exponential: d must be >= 0");
+
+      double m1 = midm1(ntaps);
+      double p = 1.0 / (8.69*ntaps/(2.0*d));
+      std::vector<float> taps(ntaps);
+      for(int i = 0; i < midn(ntaps)+1; i++) {
+        taps[i] = exp(-fabs(i - m1)*p);
+        taps[ntaps-i-1] = taps[i];
+      }
+      return taps;
+    }
+
+    std::vector<float>
+    window::riemann(int ntaps)
+    {
+      double cx;
+      double sr1 = freq(ntaps);
+      double mid = midn(ntaps);
+      std::vector<float> taps(ntaps);
+      for(int i = 0; i < mid; i++) {
+        if(i == midn(ntaps)) {
+          taps[i] = 1.0;
+          taps[ntaps-i-1] = 1.0;
+        }
+        else {
+          cx = sr1*(i - mid);
+          taps[i] = sin(cx)/cx;
+          taps[ntaps-i-1] = sin(cx)/cx;
+        }
+      }
+      return taps;
+    }
+
+    std::vector<float>
+    window::build(win_type type, int ntaps, double beta)
+    {
+      switch (type) {
+      case WIN_RECTANGULAR: return rectangular(ntaps);
+      case WIN_HAMMING: return hamming(ntaps);
+      case WIN_HANN: return hann(ntaps);
+      case WIN_BLACKMAN: return blackman(ntaps);
+      case WIN_BLACKMAN_hARRIS: return blackman_harris(ntaps);
+      case WIN_KAISER: return kaiser(ntaps, beta);
+      case WIN_BARTLETT: return bartlett(ntaps);
+      case WIN_FLATTOP: return flattop(ntaps);
+      default:
+        throw std::out_of_range("window::build: type out of range");
+      }
+    }
+
+  } /* namespace fft */
+} /* namespace gr */