git.gag.com Git - debian/gnuradio/blob - gnuradio-core/src/lib/general/gr_firdes.cc

   1 /* -*- c++ -*- */
   2 /*
   3  * Copyright 2002,2007,2008 Free Software Foundation, Inc.
   4  *
   5  * This file is part of GNU Radio
   6  *
   7  * GNU Radio is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU General Public License as published by
   9  * the Free Software Foundation; either version 3, or (at your option)
  10  * any later version.
  11  *
  12  * GNU Radio is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with GNU Radio; see the file COPYING.  If not, write to
  19  * the Free Software Foundation, Inc., 51 Franklin Street,
  20  * Boston, MA 02110-1301, USA.
  21  */
  22
  23 #include <gr_firdes.h>
  24 #include <stdexcept>
  25
  26
  27 using std::vector;
  28
  29 #define IzeroEPSILON 1E-21               /* Max error acceptable in Izero */
  30
  31 static double Izero(double x)
  32 {
  33    double sum, u, halfx, temp;
  34    int n;
  35
  36    sum = u = n = 1;
  37    halfx = x/2.0;
  38    do {
  39       temp = halfx/(double)n;
  40       n += 1;
  41       temp *= temp;
  42       u *= temp;
  43       sum += u;
  44       } while (u >= IzeroEPSILON*sum);
  45    return(sum);
  46 }
  47
  48
  49 //
  50 //      === Low Pass ===
  51 //
  52
  53 vector<float>
  54 gr_firdes::low_pass_2(double gain,
  55                       double sampling_freq,    // Hz
  56                       double cutoff_freq,      // Hz BEGINNING of transition band
  57                       double transition_width, // Hz width of transition band
  58                       double attenuation_dB,   // attenuation dB
  59                       win_type window_type,
  60                       double beta)             // used only with Kaiser
  61 {
  62   sanity_check_1f (sampling_freq, cutoff_freq, transition_width);
  63
  64   int ntaps = compute_ntaps_windes (sampling_freq, transition_width,
  65                                     attenuation_dB);
  66
  67   // construct the truncated ideal impulse response
  68   // [sin(x)/x for the low pass case]
  69
  70   vector<float> taps(ntaps);
  71   vector<float> w = window (window_type, ntaps, beta);
  72
  73   int M = (ntaps - 1) / 2;
  74   double fwT0 = 2 * M_PI * cutoff_freq / sampling_freq;
  75   for (int n = -M; n <= M; n++){
  76     if (n == 0)
  77       taps[n + M] = fwT0 / M_PI * w[n + M];
  78     else {
  79       // a little algebra gets this into the more familiar sin(x)/x form
  80       taps[n + M] =  sin (n * fwT0) / (n * M_PI) * w[n + M];
  81     }
  82   }
  83
  84   // find the factor to normalize the gain, fmax.
  85   // For low-pass, gain @ zero freq = 1.0
  86
  87   double fmax = taps[0 + M];
  88   for (int n = 1; n <= M; n++)
  89     fmax += 2 * taps[n + M];
  90
  91   gain /= fmax; // normalize
  92
  93   for (int i = 0; i < ntaps; i++)
  94     taps[i] *= gain;
  95
  96
  97   return taps;
  98 }
  99
 100 vector<float>
 101 gr_firdes::low_pass (double gain,
 102                      double sampling_freq,
 103                      double cutoff_freq,        // Hz center of transition band
 104                      double transition_width,   // Hz width of transition band
 105                      win_type window_type,
 106                      double beta)               // used only with Kaiser
 107 {
 108   sanity_check_1f (sampling_freq, cutoff_freq, transition_width);
 109
 110   int ntaps = compute_ntaps (sampling_freq, transition_width,
 111                              window_type, beta);
 112
 113   // construct the truncated ideal impulse response
 114   // [sin(x)/x for the low pass case]
 115
 116   vector<float> taps(ntaps);
 117   vector<float> w = window (window_type, ntaps, beta);
 118
 119   int M = (ntaps - 1) / 2;
 120   double fwT0 = 2 * M_PI * cutoff_freq / sampling_freq;
 121
 122   for (int n = -M; n <= M; n++){
 123     if (n == 0)
 124       taps[n + M] = fwT0 / M_PI * w[n + M];
 125     else {
 126       // a little algebra gets this into the more familiar sin(x)/x form
 127       taps[n + M] =  sin (n * fwT0) / (n * M_PI) * w[n + M];
 128     }
 129   }
 130
 131   // find the factor to normalize the gain, fmax.
 132   // For low-pass, gain @ zero freq = 1.0
 133
 134   double fmax = taps[0 + M];
 135   for (int n = 1; n <= M; n++)
 136     fmax += 2 * taps[n + M];
 137
 138   gain /= fmax; // normalize
 139
 140   for (int i = 0; i < ntaps; i++)
 141     taps[i] *= gain;
 142
 143   return taps;
 144 }
 145
 146
 147 //
 148 //      === High Pass ===
 149 //
 150
 151 vector<float>
 152 gr_firdes::high_pass_2 (double gain,
 153                       double sampling_freq,
 154                       double cutoff_freq,       // Hz center of transition band
 155                       double transition_width,  // Hz width of transition band
 156                       double attenuation_dB,   // attenuation dB
 157                       win_type window_type,
 158                       double beta)              // used only with Kaiser
 159 {
 160   sanity_check_1f (sampling_freq, cutoff_freq, transition_width);
 161
 162   int ntaps = compute_ntaps_windes (sampling_freq, transition_width,
 163                                     attenuation_dB);
 164
 165   // construct the truncated ideal impulse response times the window function
 166
 167   vector<float> taps(ntaps);
 168   vector<float> w = window (window_type, ntaps, beta);
 169
 170   int M = (ntaps - 1) / 2;
 171   double fwT0 = 2 * M_PI * cutoff_freq / sampling_freq;
 172
 173   for (int n = -M; n <= M; n++){
 174     if (n == 0)
 175       taps[n + M] = (1 - (fwT0 / M_PI)) * w[n + M];
 176     else {
 177       // a little algebra gets this into the more familiar sin(x)/x form
 178       taps[n + M] =  -sin (n * fwT0) / (n * M_PI) * w[n + M];
 179     }
 180   }
 181
 182   // find the factor to normalize the gain, fmax.
 183   // For high-pass, gain @ fs/2 freq = 1.0
 184
 185   double fmax = taps[0 + M];
 186   for (int n = 1; n <= M; n++)
 187     fmax += 2 * taps[n + M] * cos (n * M_PI);
 188
 189   gain /= fmax; // normalize
 190
 191   for (int i = 0; i < ntaps; i++)
 192     taps[i] *= gain;
 193
 194
 195   return taps;
 196 }
 197
 198
 199 vector<float>
 200 gr_firdes::high_pass (double gain,
 201                       double sampling_freq,
 202                       double cutoff_freq,       // Hz center of transition band
 203                       double transition_width,  // Hz width of transition band
 204                       win_type window_type,
 205                       double beta)              // used only with Kaiser
 206 {
 207   sanity_check_1f (sampling_freq, cutoff_freq, transition_width);
 208
 209   int ntaps = compute_ntaps (sampling_freq, transition_width,
 210                              window_type, beta);
 211
 212   // construct the truncated ideal impulse response times the window function
 213
 214   vector<float> taps(ntaps);
 215   vector<float> w = window (window_type, ntaps, beta);
 216
 217   int M = (ntaps - 1) / 2;
 218   double fwT0 = 2 * M_PI * cutoff_freq / sampling_freq;
 219
 220   for (int n = -M; n <= M; n++){
 221     if (n == 0)
 222       taps[n + M] = (1 - (fwT0 / M_PI)) * w[n + M];
 223     else {
 224       // a little algebra gets this into the more familiar sin(x)/x form
 225       taps[n + M] =  -sin (n * fwT0) / (n * M_PI) * w[n + M];
 226     }
 227   }
 228
 229   // find the factor to normalize the gain, fmax.
 230   // For high-pass, gain @ fs/2 freq = 1.0
 231
 232   double fmax = taps[0 + M];
 233   for (int n = 1; n <= M; n++)
 234     fmax += 2 * taps[n + M] * cos (n * M_PI);
 235
 236   gain /= fmax; // normalize
 237
 238   for (int i = 0; i < ntaps; i++)
 239     taps[i] *= gain;
 240
 241   return taps;
 242 }
 243
 244 //
 245 //      === Band Pass ===
 246 //
 247
 248 vector<float>
 249 gr_firdes::band_pass_2 (double gain,
 250                       double sampling_freq,
 251                       double low_cutoff_freq,   // Hz center of transition band
 252                       double high_cutoff_freq,  // Hz center of transition band
 253                       double transition_width,  // Hz width of transition band
 254                       double attenuation_dB,   // attenuation dB
 255                       win_type window_type,
 256                       double beta)              // used only with Kaiser
 257 {
 258   sanity_check_2f (sampling_freq,
 259                    low_cutoff_freq,
 260                    high_cutoff_freq, transition_width);
 261
 262   int ntaps = compute_ntaps_windes (sampling_freq, transition_width,
 263                                     attenuation_dB);
 264
 265   vector<float> taps(ntaps);
 266   vector<float> w = window (window_type, ntaps, beta);
 267
 268   int M = (ntaps - 1) / 2;
 269   double fwT0 = 2 * M_PI * low_cutoff_freq / sampling_freq;
 270   double fwT1 = 2 * M_PI * high_cutoff_freq / sampling_freq;
 271
 272   for (int n = -M; n <= M; n++){
 273     if (n == 0)
 274       taps[n + M] = (fwT1 - fwT0) / M_PI * w[n + M];
 275     else {
 276       taps[n + M] =  (sin (n * fwT1) - sin (n * fwT0)) / (n * M_PI) * w[n + M];
 277     }
 278   }
 279
 280   // find the factor to normalize the gain, fmax.
 281   // For band-pass, gain @ center freq = 1.0
 282
 283   double fmax = taps[0 + M];
 284   for (int n = 1; n <= M; n++)
 285     fmax += 2 * taps[n + M] * cos (n * (fwT0 + fwT1) * 0.5);
 286
 287   gain /= fmax; // normalize
 288
 289   for (int i = 0; i < ntaps; i++)
 290     taps[i] *= gain;
 291
 292   return taps;
 293 }
 294
 295
 296 vector<float>
 297 gr_firdes::band_pass (double gain,
 298                       double sampling_freq,
 299                       double low_cutoff_freq,   // Hz center of transition band
 300                       double high_cutoff_freq,  // Hz center of transition band
 301                       double transition_width,  // Hz width of transition band
 302                       win_type window_type,
 303                       double beta)              // used only with Kaiser
 304 {
 305   sanity_check_2f (sampling_freq,
 306                    low_cutoff_freq,
 307                    high_cutoff_freq, transition_width);
 308
 309   int ntaps = compute_ntaps (sampling_freq, transition_width,
 310                              window_type, beta);
 311
 312   // construct the truncated ideal impulse response times the window function
 313
 314   vector<float> taps(ntaps);
 315   vector<float> w = window (window_type, ntaps, beta);
 316
 317   int M = (ntaps - 1) / 2;
 318   double fwT0 = 2 * M_PI * low_cutoff_freq / sampling_freq;
 319   double fwT1 = 2 * M_PI * high_cutoff_freq / sampling_freq;
 320
 321   for (int n = -M; n <= M; n++){
 322     if (n == 0)
 323       taps[n + M] = (fwT1 - fwT0) / M_PI * w[n + M];
 324     else {
 325       taps[n + M] =  (sin (n * fwT1) - sin (n * fwT0)) / (n * M_PI) * w[n + M];
 326     }
 327   }
 328
 329   // find the factor to normalize the gain, fmax.
 330   // For band-pass, gain @ center freq = 1.0
 331
 332   double fmax = taps[0 + M];
 333   for (int n = 1; n <= M; n++)
 334     fmax += 2 * taps[n + M] * cos (n * (fwT0 + fwT1) * 0.5);
 335
 336   gain /= fmax; // normalize
 337
 338   for (int i = 0; i < ntaps; i++)
 339     taps[i] *= gain;
 340
 341   return taps;
 342 }
 343
 344 //
 345 //      === Complex Band Pass ===
 346 //
 347
 348 vector<gr_complex>
 349 gr_firdes::complex_band_pass_2 (double gain,
 350                       double sampling_freq,
 351                       double low_cutoff_freq,   // Hz center of transition band
 352                       double high_cutoff_freq,  // Hz center of transition band
 353                       double transition_width,  // Hz width of transition band
 354                       double attenuation_dB,      // attenuation dB
 355                       win_type window_type,
 356                       double beta)              // used only with Kaiser
 357 {
 358   sanity_check_2f_c (sampling_freq,
 359                    low_cutoff_freq,
 360                    high_cutoff_freq, transition_width);
 361
 362   int ntaps = compute_ntaps_windes (sampling_freq, transition_width,
 363                                     attenuation_dB);
 364
 365
 366
 367   vector<gr_complex> taps(ntaps);
 368   vector<float> lptaps(ntaps);
 369   vector<float> w = window (window_type, ntaps, beta);
 370
 371   lptaps = low_pass_2(gain,sampling_freq,(high_cutoff_freq - low_cutoff_freq)/2,transition_width,attenuation_dB,window_type,beta);
 372
 373   gr_complex *optr = &taps[0];
 374   float *iptr = &lptaps[0];
 375   float freq = M_PI * (high_cutoff_freq + low_cutoff_freq)/sampling_freq;
 376   float phase=0;
 377   if (lptaps.size() & 01) {
 378       phase = - freq * ( lptaps.size() >> 1 );
 379   } else phase = - freq/2.0 * ((1 + 2*lptaps.size()) >> 1);
 380   for(unsigned int i=0;i<lptaps.size();i++) {
 381       *optr++ = gr_complex(*iptr * cos(phase),*iptr * sin(phase));
 382       iptr++, phase += freq;
 383   }
 384
 385   return taps;
 386 }
 387
 388
 389 vector<gr_complex>
 390 gr_firdes::complex_band_pass (double gain,
 391                       double sampling_freq,
 392                       double low_cutoff_freq,   // Hz center of transition band
 393                       double high_cutoff_freq,  // Hz center of transition band
 394                       double transition_width,  // Hz width of transition band
 395                       win_type window_type,
 396                       double beta)              // used only with Kaiser
 397 {
 398   sanity_check_2f_c (sampling_freq,
 399                    low_cutoff_freq,
 400                    high_cutoff_freq, transition_width);
 401
 402   int ntaps = compute_ntaps (sampling_freq, transition_width,
 403                              window_type, beta);
 404
 405   // construct the truncated ideal impulse response times the window function
 406
 407   vector<gr_complex> taps(ntaps);
 408   vector<float> lptaps(ntaps);
 409   vector<float> w = window (window_type, ntaps, beta);
 410
 411   lptaps = low_pass(gain,sampling_freq,(high_cutoff_freq - low_cutoff_freq)/2,transition_width,window_type,beta);
 412
 413   gr_complex *optr = &taps[0];
 414   float *iptr = &lptaps[0];
 415   float freq = M_PI * (high_cutoff_freq + low_cutoff_freq)/sampling_freq;
 416   float phase=0;
 417   if (lptaps.size() & 01) {
 418       phase = - freq * ( lptaps.size() >> 1 );
 419   } else phase = - freq/2.0 * ((1 + 2*lptaps.size()) >> 1);
 420   for(unsigned int i=0;i<lptaps.size();i++) {
 421       *optr++ = gr_complex(*iptr * cos(phase),*iptr * sin(phase));
 422       iptr++, phase += freq;
 423   }
 424
 425   return taps;
 426 }
 427
 428 //
 429 //      === Band Reject ===
 430 //
 431
 432 vector<float>
 433 gr_firdes::band_reject_2 (double gain,
 434                         double sampling_freq,
 435                         double low_cutoff_freq,  // Hz center of transition band
 436                         double high_cutoff_freq, // Hz center of transition band
 437                         double transition_width, // Hz width of transition band
 438                         double attenuation_dB,   // attenuation dB
 439                         win_type window_type,
 440                         double beta)             // used only with Kaiser
 441 {
 442   sanity_check_2f (sampling_freq,
 443                    low_cutoff_freq,
 444                    high_cutoff_freq, transition_width);
 445
 446   int ntaps = compute_ntaps_windes (sampling_freq, transition_width,
 447                                     attenuation_dB);
 448
 449   // construct the truncated ideal impulse response times the window function
 450
 451   vector<float> taps(ntaps);
 452   vector<float> w = window (window_type, ntaps, beta);
 453
 454   int M = (ntaps - 1) / 2;
 455   double fwT0 = 2 * M_PI * low_cutoff_freq / sampling_freq;
 456   double fwT1 = 2 * M_PI * high_cutoff_freq / sampling_freq;
 457
 458   for (int n = -M; n <= M; n++){
 459     if (n == 0)
 460       taps[n + M] = 1.0 + ((fwT0 - fwT1) / M_PI * w[n + M]);
 461     else {
 462       taps[n + M] =  (sin (n * fwT0) - sin (n * fwT1)) / (n * M_PI) * w[n + M];
 463     }
 464   }
 465
 466   // find the factor to normalize the gain, fmax.
 467   // For band-reject, gain @ zero freq = 1.0
 468
 469   double fmax = taps[0 + M];
 470   for (int n = 1; n <= M; n++)
 471     fmax += 2 * taps[n + M];
 472
 473   gain /= fmax; // normalize
 474
 475   for (int i = 0; i < ntaps; i++)
 476     taps[i] *= gain;
 477
 478   return taps;
 479 }
 480
 481 vector<float>
 482 gr_firdes::band_reject (double gain,
 483                         double sampling_freq,
 484                         double low_cutoff_freq,  // Hz center of transition band
 485                         double high_cutoff_freq, // Hz center of transition band
 486                         double transition_width, // Hz width of transition band
 487                         win_type window_type,
 488                         double beta)             // used only with Kaiser
 489 {
 490   sanity_check_2f (sampling_freq,
 491                    low_cutoff_freq,
 492                    high_cutoff_freq, transition_width);
 493
 494   int ntaps = compute_ntaps (sampling_freq, transition_width,
 495                              window_type, beta);
 496
 497   // construct the truncated ideal impulse response times the window function
 498
 499   vector<float> taps(ntaps);
 500   vector<float> w = window (window_type, ntaps, beta);
 501
 502   int M = (ntaps - 1) / 2;
 503   double fwT0 = 2 * M_PI * low_cutoff_freq / sampling_freq;
 504   double fwT1 = 2 * M_PI * high_cutoff_freq / sampling_freq;
 505
 506   for (int n = -M; n <= M; n++){
 507     if (n == 0)
 508       taps[n + M] = 1.0 + ((fwT0 - fwT1) / M_PI * w[n + M]);
 509     else {
 510       taps[n + M] =  (sin (n * fwT0) - sin (n * fwT1)) / (n * M_PI) * w[n + M];
 511     }
 512   }
 513
 514   // find the factor to normalize the gain, fmax.
 515   // For band-reject, gain @ zero freq = 1.0
 516
 517   double fmax = taps[0 + M];
 518   for (int n = 1; n <= M; n++)
 519     fmax += 2 * taps[n + M];
 520
 521   gain /= fmax; // normalize
 522
 523   for (int i = 0; i < ntaps; i++)
 524     taps[i] *= gain;
 525
 526   return taps;
 527 }
 528
 529 //
 530 // Hilbert Transform
 531 //
 532
 533 vector<float>
 534 gr_firdes::hilbert (unsigned int ntaps,
 535                     win_type windowtype,
 536                     double beta)
 537 {
 538   if(!(ntaps & 1))
 539     throw std::out_of_range("Hilbert:  Must have odd number of taps");
 540
 541   vector<float> taps(ntaps);
 542   vector<float> w = window (windowtype, ntaps, beta);
 543   unsigned int h = (ntaps-1)/2;
 544   float gain=0;
 545   for (unsigned int i = 1; i <= h; i++)
 546     {
 547       if(i&1)
 548         {
 549           float x = 1/(float)i;
 550           taps[h+i] = x * w[h+i];
 551           taps[h-i] = -x * w[h-i];
 552           gain = taps[h+i] - gain;
 553         }
 554       else
 555         taps[h+i] = taps[h-i] = 0;
 556     }
 557   gain = 2 * fabs(gain);
 558   for ( unsigned int i = 0; i < ntaps; i++)
 559       taps[i] /= gain;
 560   return taps;
 561 }
 562
 563 //
 564 // Gaussian
 565 //
 566
 567 vector<float>
 568 gr_firdes::gaussian (double gain,
 569                      double spb,
 570                      double bt,
 571                      int ntaps)
 572 {
 573
 574   vector<float> taps(ntaps);
 575   double scale = 0;
 576   double dt = 1.0/spb;
 577   double s = 1.0/(sqrt(log(2)) / (2*M_PI*bt));
 578   double t0 = -0.5 * ntaps;
 579   double ts;
 580   for(int i=0;i<ntaps;i++)
 581     {
 582       t0++;
 583       ts = s*dt*t0;
 584       taps[i] = exp(-0.5*ts*ts);
 585       scale += taps[i];
 586     }
 587   for(int i=0;i<ntaps;i++)
 588     taps[i] = taps[i] / scale * gain;
 589   return taps;
 590 }
 591
 592
 593 //
 594 // Root Raised Cosine
 595 //
 596
 597 vector<float>
 598 gr_firdes::root_raised_cosine (double gain,
 599                                double sampling_freq,
 600                                double symbol_rate,
 601                                double alpha,
 602                                int ntaps)
 603 {
 604   ntaps |= 1;   // ensure that ntaps is odd
 605
 606   double spb = sampling_freq/symbol_rate; // samples per bit/symbol
 607   vector<float> taps(ntaps);
 608   double scale = 0;
 609   for(int i=0;i<ntaps;i++)
 610     {
 611       double x1,x2,x3,num,den;
 612       double xindx = i - ntaps/2;
 613       x1 = M_PI * xindx/spb;
 614       x2 = 4 * alpha * xindx / spb;
 615       x3 = x2*x2 - 1;
 616       if( fabs(x3) >= 0.000001 )  // Avoid Rounding errors...
 617         {
 618           if( i != ntaps/2 )
 619             num = cos((1+alpha)*x1) + sin((1-alpha)*x1)/(4*alpha*xindx/spb);
 620           else
 621             num = cos((1+alpha)*x1) + (1-alpha) * M_PI / (4*alpha);
 622           den = x3 * M_PI;
 623         }
 624       else
 625         {
 626           if(alpha==1)
 627             {
 628               taps[i] = -1;
 629               continue;
 630             }
 631           x3 = (1-alpha)*x1;
 632           x2 = (1+alpha)*x1;
 633           num = (sin(x2)*(1+alpha)*M_PI
 634                  - cos(x3)*((1-alpha)*M_PI*spb)/(4*alpha*xindx)
 635                  + sin(x3)*spb*spb/(4*alpha*xindx*xindx));
 636           den = -32 * M_PI * alpha * alpha * xindx/spb;
 637         }
 638       taps[i] = 4 * alpha * num / den;
 639       scale += taps[i];
 640     }
 641
 642   for(int i=0;i<ntaps;i++)
 643     taps[i] = taps[i] * gain / scale;
 644
 645   return taps;
 646 }
 647
 648 //
 649 //      === Utilities ===
 650 //
 651
 652 // delta_f / width_factor gives number of taps required.
 653 static const float width_factor[5] = {   // indexed by win_type
 654   3.3,                  // WIN_HAMMING
 655   3.1,                  // WIN_HANN
 656   5.5,                  // WIN_BLACKMAN
 657   2.0,                  // WIN_RECTANGULAR
 658   //5.0                   // WIN_KAISER  (guesstimate compromise)
 659   //2.0                   // WIN_KAISER  (guesstimate compromise)
 660   10.0                  // WIN_KAISER
 661 };
 662
 663 int
 664 gr_firdes::compute_ntaps_windes(double sampling_freq,
 665                                  double transition_width, // this is frequency, not relative frequency
 666                                  double attenuation_dB)
 667 {
 668   // Based on formula from Multirate Signal Processing for
 669   // Communications Systems, fredric j harris
 670   int ntaps = (int)(attenuation_dB*sampling_freq/(22.0*transition_width));
 671   if ((ntaps & 1) == 0) // if even...
 672     ntaps++;            // ...make odd
 673   return ntaps;
 674 }
 675
 676 int
 677 gr_firdes::compute_ntaps (double sampling_freq,
 678                           double transition_width,
 679                           win_type window_type,
 680                           double beta)
 681 {
 682   // normalized transition width
 683   double delta_f = transition_width / sampling_freq;
 684
 685   // compute number of taps required for given transition width
 686   int ntaps = (int) (width_factor[window_type] / delta_f + 0.5);
 687   if ((ntaps & 1) == 0) // if even...
 688     ntaps++;            // ...make odd
 689
 690   return ntaps;
 691 }
 692
 693 double gr_firdes::bessi0(double x)
 694 {
 695         double ax,ans;
 696         double y;
 697
 698         ax=fabs(x);
 699         if (ax < 3.75)
 700         {
 701                 y=x/3.75;
 702                 y*=y;
 703                 ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
 704                         +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
 705         }
 706         else
 707         {
 708                 y=3.75/ax;
 709                 ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1
 710                         +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
 711                         +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
 712                         +y*0.392377e-2))))))));
 713         }
 714         return ans;
 715 }
 716 vector<float>
 717 gr_firdes::window (win_type type, int ntaps, double beta)
 718 {
 719   vector<float> taps(ntaps);
 720   int   M = ntaps - 1;          // filter order
 721
 722   switch (type){
 723   case WIN_RECTANGULAR:
 724     for (int n = 0; n < ntaps; n++)
 725       taps[n] = 1;
 726
 727   case WIN_HAMMING:
 728     for (int n = 0; n < ntaps; n++)
 729       taps[n] = 0.54 - 0.46 * cos ((2 * M_PI * n) / M);
 730     break;
 731
 732   case WIN_HANN:
 733     for (int n = 0; n < ntaps; n++)
 734       taps[n] = 0.5 - 0.5 * cos ((2 * M_PI * n) / M);
 735     break;
 736
 737   case WIN_BLACKMAN:
 738     for (int n = 0; n < ntaps; n++)
 739       taps[n] = 0.42 - 0.50 * cos ((2*M_PI * n) / (M-1)) - 0.08 * cos ((4*M_PI * n) / (M-1));
 740     break;
 741
 742   case WIN_BLACKMAN_hARRIS:
 743     for (int n = -ntaps/2; n < ntaps/2; n++)
 744       taps[n+ntaps/2] = 0.35875 + 0.48829*cos((2*M_PI * n) / (float)M) +
 745         0.14128*cos((4*M_PI * n) / (float)M) + 0.01168*cos((6*M_PI * n) / (float)M);
 746     break;
 747
 748 #if 0
 749   case WIN_KAISER:
 750     for (int n = 0; n < ntaps; n++)
 751         taps[n] = bessi0(beta*sqrt(1.0 - (4.0*n/(M*M))))/bessi0(beta);
 752     break;
 753 #else
 754
 755   case WIN_KAISER:
 756     {
 757       double IBeta = 1.0/Izero(beta);
 758       double inm1 = 1.0/((double)(ntaps));
 759       double temp;
 760       //fprintf(stderr, "IBeta = %g; inm1 = %g\n", IBeta, inm1);
 761
 762       for (int i=0; i<ntaps; i++) {
 763         temp = i * inm1;
 764         //fprintf(stderr, "temp = %g\n", temp);
 765         taps[i] = Izero(beta*sqrt(1.0-temp*temp)) * IBeta;
 766         //fprintf(stderr, "taps[%d] = %g\n", i, taps[i]);
 767       }
 768     }
 769     break;
 770
 771 #endif
 772   default:
 773     throw std::out_of_range ("gr_firdes:window: type out of range");
 774   }
 775
 776   return taps;
 777 }
 778
 779 void
 780 gr_firdes::sanity_check_1f (double sampling_freq,
 781                             double fa,                  // cutoff freq
 782                             double transition_width)
 783 {
 784   if (sampling_freq <= 0.0)
 785     throw std::out_of_range ("gr_firdes check failed: sampling_freq > 0");
 786
 787   if (fa <= 0.0 || fa > sampling_freq / 2)
 788     throw std::out_of_range ("gr_firdes check failed: 0 < fa <= sampling_freq / 2");
 789
 790   if (transition_width <= 0)
 791     throw std::out_of_range ("gr_dirdes check failed: transition_width > 0");
 792 }
 793
 794 void
 795 gr_firdes::sanity_check_2f (double sampling_freq,
 796                             double fa,                  // first cutoff freq
 797                             double fb,                  // second cutoff freq
 798                             double transition_width)
 799 {
 800   if (sampling_freq <= 0.0)
 801     throw std::out_of_range ("gr_firdes check failed: sampling_freq > 0");
 802
 803   if (fa <= 0.0 || fa > sampling_freq / 2)
 804     throw std::out_of_range ("gr_firdes check failed: 0 < fa <= sampling_freq / 2");
 805
 806   if (fb <= 0.0 || fb > sampling_freq / 2)
 807     throw std::out_of_range ("gr_firdes check failed: 0 < fb <= sampling_freq / 2");
 808
 809   if (fa > fb)
 810     throw std::out_of_range ("gr_firdes check failed: fa <= fb");
 811
 812   if (transition_width <= 0)
 813     throw std::out_of_range ("gr_firdes check failed: transition_width > 0");
 814 }
 815
 816 void
 817 gr_firdes::sanity_check_2f_c (double sampling_freq,
 818                             double fa,                  // first cutoff freq
 819                             double fb,                  // second cutoff freq
 820                             double transition_width)
 821 {
 822   if (sampling_freq <= 0.0)
 823     throw std::out_of_range ("gr_firdes check failed: sampling_freq > 0");
 824
 825   if (fa < -sampling_freq / 2 || fa > sampling_freq / 2)
 826     throw std::out_of_range ("gr_firdes check failed: 0 < fa <= sampling_freq / 2");
 827
 828   if (fb < -sampling_freq / 2  || fb > sampling_freq / 2)
 829     throw std::out_of_range ("gr_firdes check failed: 0 < fb <= sampling_freq / 2");
 830
 831   if (fa > fb)
 832     throw std::out_of_range ("gr_firdes check failed: fa <= fb");
 833
 834   if (transition_width <= 0)
 835     throw std::out_of_range ("gr_firdes check failed: transition_width > 0");
 836 }