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

   1 /* -*- c++ -*- */
   2 /*
   3  * Copyright 2002 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 2, 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 (double gain,
  55                      double sampling_freq,
  56                      double cutoff_freq,        // Hz center of transition band
  57                      double transition_width,   // Hz width of transition band
  58                      win_type window_type,
  59                      double beta)               // used only with Kaiser
  60 {
  61   sanity_check_1f (sampling_freq, cutoff_freq, transition_width);
  62
  63   int ntaps = compute_ntaps (sampling_freq, transition_width,
  64                              window_type, beta);
  65
  66   // construct the truncated ideal impulse response
  67   // [sin(x)/x for the low pass case]
  68
  69   vector<float> taps(ntaps);
  70   vector<float> w = window (window_type, ntaps, beta);
  71
  72   int M = (ntaps - 1) / 2;
  73   double fwT0 = 2 * M_PI * cutoff_freq / sampling_freq;
  74
  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   return taps;
  97 }
  98
  99 //
 100 //      === High Pass ===
 101 //
 102
 103 vector<float>
 104 gr_firdes::high_pass (double gain,
 105                       double sampling_freq,
 106                       double cutoff_freq,       // Hz center of transition band
 107                       double transition_width,  // Hz width of transition band
 108                       win_type window_type,
 109                       double beta)              // used only with Kaiser
 110 {
 111   sanity_check_1f (sampling_freq, cutoff_freq, transition_width);
 112
 113   int ntaps = compute_ntaps (sampling_freq, transition_width,
 114                              window_type, beta);
 115
 116   // construct the truncated ideal impulse response times the window function
 117
 118   vector<float> taps(ntaps);
 119   vector<float> w = window (window_type, ntaps, beta);
 120
 121   int M = (ntaps - 1) / 2;
 122   double fwT0 = 2 * M_PI * cutoff_freq / sampling_freq;
 123
 124   for (int n = -M; n <= M; n++){
 125     if (n == 0)
 126       taps[n + M] = (1 - (fwT0 / M_PI)) * w[n + M];
 127     else {
 128       // a little algebra gets this into the more familiar sin(x)/x form
 129       taps[n + M] =  -sin (n * fwT0) / (n * M_PI) * w[n + M];
 130     }
 131   }
 132
 133   // find the factor to normalize the gain, fmax.
 134   // For high-pass, gain @ fs/2 freq = 1.0
 135
 136   double fmax = taps[0 + M];
 137   for (int n = 1; n <= M; n++)
 138     fmax += 2 * taps[n + M] * cos (n * M_PI);
 139
 140   gain /= fmax; // normalize
 141
 142   for (int i = 0; i < ntaps; i++)
 143     taps[i] *= gain;
 144
 145   return taps;
 146 }
 147
 148 //
 149 //      === Band Pass ===
 150 //
 151
 152 vector<float>
 153 gr_firdes::band_pass (double gain,
 154                       double sampling_freq,
 155                       double low_cutoff_freq,   // Hz center of transition band
 156                       double high_cutoff_freq,  // Hz center of transition band
 157                       double transition_width,  // Hz width of transition band
 158                       win_type window_type,
 159                       double beta)              // used only with Kaiser
 160 {
 161   sanity_check_2f (sampling_freq,
 162                    low_cutoff_freq,
 163                    high_cutoff_freq, transition_width);
 164
 165   int ntaps = compute_ntaps (sampling_freq, transition_width,
 166                              window_type, beta);
 167
 168   // construct the truncated ideal impulse response times the window function
 169
 170   vector<float> taps(ntaps);
 171   vector<float> w = window (window_type, ntaps, beta);
 172
 173   int M = (ntaps - 1) / 2;
 174   double fwT0 = 2 * M_PI * low_cutoff_freq / sampling_freq;
 175   double fwT1 = 2 * M_PI * high_cutoff_freq / sampling_freq;
 176
 177   for (int n = -M; n <= M; n++){
 178     if (n == 0)
 179       taps[n + M] = (fwT1 - fwT0) / M_PI * w[n + M];
 180     else {
 181       taps[n + M] =  (sin (n * fwT1) - sin (n * fwT0)) / (n * M_PI) * w[n + M];
 182     }
 183   }
 184
 185   // find the factor to normalize the gain, fmax.
 186   // For band-pass, gain @ center freq = 1.0
 187
 188   double fmax = taps[0 + M];
 189   for (int n = 1; n <= M; n++)
 190     fmax += 2 * taps[n + M] * cos (n * (fwT0 + fwT1) * 0.5);
 191
 192   gain /= fmax; // normalize
 193
 194   for (int i = 0; i < ntaps; i++)
 195     taps[i] *= gain;
 196
 197   return taps;
 198 }
 199
 200 //
 201 //      === Complex Band Pass ===
 202 //
 203
 204 vector<gr_complex>
 205 gr_firdes::complex_band_pass (double gain,
 206                       double sampling_freq,
 207                       double low_cutoff_freq,   // Hz center of transition band
 208                       double high_cutoff_freq,  // Hz center of transition band
 209                       double transition_width,  // Hz width of transition band
 210                       win_type window_type,
 211                       double beta)              // used only with Kaiser
 212 {
 213   sanity_check_2f_c (sampling_freq,
 214                    low_cutoff_freq,
 215                    high_cutoff_freq, transition_width);
 216
 217   int ntaps = compute_ntaps (sampling_freq, transition_width,
 218                              window_type, beta);
 219
 220   // construct the truncated ideal impulse response times the window function
 221
 222   vector<gr_complex> taps(ntaps);
 223   vector<float> lptaps(ntaps);
 224   vector<float> w = window (window_type, ntaps, beta);
 225
 226   lptaps = low_pass(gain,sampling_freq,(high_cutoff_freq - low_cutoff_freq)/2,transition_width,window_type,beta);
 227
 228   gr_complex *optr = &taps[0];
 229   float *iptr = &lptaps[0];
 230   float freq = M_PI * (high_cutoff_freq + low_cutoff_freq)/sampling_freq;
 231   float phase=0;
 232   if (lptaps.size() & 01) {
 233       phase = - freq * ( lptaps.size() >> 1 );
 234   } else phase = - freq/2.0 * ((1 + 2*lptaps.size()) >> 1);
 235   for(unsigned int i=0;i<lptaps.size();i++) {
 236       *optr++ = gr_complex(*iptr * cos(phase),*iptr * sin(phase));
 237       iptr++, phase += freq;
 238   }
 239
 240   return taps;
 241 }
 242
 243
 244 //
 245 //      === Band Reject ===
 246 //
 247
 248 vector<float>
 249 gr_firdes::band_reject (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                         win_type window_type,
 255                         double beta)             // used only with Kaiser
 256 {
 257   sanity_check_2f (sampling_freq,
 258                    low_cutoff_freq,
 259                    high_cutoff_freq, transition_width);
 260
 261   int ntaps = compute_ntaps (sampling_freq, transition_width,
 262                              window_type, beta);
 263
 264   // construct the truncated ideal impulse response times the window function
 265
 266   vector<float> taps(ntaps);
 267   vector<float> w = window (window_type, ntaps, beta);
 268
 269   int M = (ntaps - 1) / 2;
 270   double fwT0 = 2 * M_PI * low_cutoff_freq / sampling_freq;
 271   double fwT1 = 2 * M_PI * high_cutoff_freq / sampling_freq;
 272
 273   for (int n = -M; n <= M; n++){
 274     if (n == 0)
 275       taps[n + M] = (1.0 + (fwT0 - fwT1)) / M_PI * w[n + M];
 276     else {
 277       taps[n + M] =  (sin (n * fwT0) - sin (n * fwT1)) / (n * M_PI) * w[n + M];
 278     }
 279   }
 280
 281   // find the factor to normalize the gain, fmax.
 282   // For band-reject, gain @ zero freq = 1.0
 283
 284   double fmax = taps[0 + M];
 285   for (int n = 1; n <= M; n++)
 286     fmax += 2 * taps[n + M];
 287
 288   gain /= fmax; // normalize
 289
 290   for (int i = 0; i < ntaps; i++)
 291     taps[i] *= gain;
 292
 293   return taps;
 294 }
 295
 296 //
 297 // Hilbert Transform
 298 //
 299
 300 vector<float>
 301 gr_firdes::hilbert (unsigned int ntaps,
 302                     win_type windowtype,
 303                     double beta)
 304 {
 305   if(!(ntaps & 1))
 306     throw std::out_of_range("Hilbert:  Must have odd number of taps");
 307
 308   vector<float> taps(ntaps);
 309   vector<float> w = window (windowtype, ntaps, beta);
 310   unsigned int h = (ntaps-1)/2;
 311   float gain=0;
 312   for (unsigned int i = 1; i <= h; i++)
 313     {
 314       if(i&1)
 315         {
 316           float x = 1/(float)i;
 317           taps[h+i] = x * w[h+i];
 318           taps[h-i] = -x * w[h-i];
 319           gain = taps[h+i] - gain;
 320         }
 321       else
 322         taps[h+i] = taps[h-i] = 0;
 323     }
 324   gain = 2 * fabs(gain);
 325   for ( unsigned int i = 0; i < ntaps; i++)
 326       taps[i] /= gain;
 327   return taps;
 328 }
 329
 330 //
 331 // Gaussian
 332 //
 333
 334 vector<float>
 335 gr_firdes::gaussian (double gain,
 336                      double spb,
 337                      double bt,
 338                      int ntaps)
 339 {
 340
 341   vector<float> taps(ntaps);
 342   double scale = 0;
 343   double dt = 1.0/spb;
 344   double s = 1.0/(sqrt(log(2)) / (2*M_PI*bt));
 345   double t0 = -0.5 * ntaps;
 346   double ts;
 347   for(int i=0;i<ntaps;i++)
 348     {
 349       t0++;
 350       ts = s*dt*t0;
 351       taps[i] = exp(-0.5*ts*ts);
 352       scale += taps[i];
 353     }
 354   for(int i=0;i<ntaps;i++)
 355     taps[i] = taps[i] / scale * gain;
 356   return taps;
 357 }
 358
 359
 360 //
 361 // Root Raised Cosine
 362 //
 363
 364 vector<float>
 365 gr_firdes::root_raised_cosine (double gain,
 366                                double sampling_freq,
 367                                double symbol_rate,
 368                                double alpha,
 369                                int ntaps)
 370 {
 371   ntaps |= 1;   // ensure that ntaps is odd
 372
 373   double spb = sampling_freq/symbol_rate; // samples per bit/symbol
 374   vector<float> taps(ntaps);
 375   double scale = 0;
 376   for(int i=0;i<ntaps;i++)
 377     {
 378       double x1,x2,x3,num,den;
 379       double xindx = i - ntaps/2;
 380       x1 = M_PI * xindx/spb;
 381       x2 = 4 * alpha * xindx / spb;
 382       x3 = x2*x2 - 1;
 383       if( fabs(x3) >= 0.000001 )  // Avoid Rounding errors...
 384         {
 385           if( i != ntaps/2 )
 386             num = cos((1+alpha)*x1) + sin((1-alpha)*x1)/(4*alpha*xindx/spb);
 387           else
 388             num = cos((1+alpha)*x1) + (1-alpha) * M_PI / (4*alpha);
 389           den = x3 * M_PI;
 390         }
 391       else
 392         {
 393           if(alpha==1)
 394             {
 395               taps[i] = -1;
 396               continue;
 397             }
 398           x3 = (1-alpha)*x1;
 399           x2 = (1+alpha)*x1;
 400           num = (sin(x2)*(1+alpha)*M_PI
 401                  - cos(x3)*((1-alpha)*M_PI*spb)/(4*alpha*xindx)
 402                  + sin(x3)*spb*spb/(4*alpha*xindx*xindx));
 403           den = -32 * M_PI * alpha * alpha * xindx/spb;
 404         }
 405       taps[i] = 4 * alpha * num / den;
 406       scale += taps[i];
 407     }
 408
 409   for(int i=0;i<ntaps;i++)
 410     taps[i] = taps[i] * gain / scale;
 411
 412   return taps;
 413 }
 414
 415 //
 416 //      === Utilities ===
 417 //
 418
 419 // delta_f / width_factor gives number of taps required.
 420 static const float width_factor[5] = {   // indexed by win_type
 421   3.3,                  // WIN_HAMMING
 422   3.1,                  // WIN_HANN
 423   5.5,                  // WIN_BLACKMAN
 424   2.0,                  // WIN_RECTANGULAR
 425   //5.0                   // WIN_KAISER  (guesstimate compromise)
 426   //2.0                   // WIN_KAISER  (guesstimate compromise)
 427   10.0                  // WIN_KAISER
 428 };
 429
 430 int
 431 gr_firdes::compute_ntaps (double sampling_freq,
 432                           double transition_width,
 433                           win_type window_type,
 434                           double beta)
 435 {
 436   // normalized transition width
 437   double delta_f = transition_width / sampling_freq;
 438
 439   // compute number of taps required for given transition width
 440   int ntaps = (int) (width_factor[window_type] / delta_f + 0.5);
 441   if ((ntaps & 1) == 0) // if even...
 442     ntaps++;            // ...make odd
 443
 444   return ntaps;
 445 }
 446
 447 double gr_firdes::bessi0(double x)
 448 {
 449         double ax,ans;
 450         double y;
 451
 452         ax=fabs(x);
 453         if (ax < 3.75)
 454         {
 455                 y=x/3.75;
 456                 y*=y;
 457                 ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
 458                         +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
 459         }
 460         else
 461         {
 462                 y=3.75/ax;
 463                 ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1
 464                         +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
 465                         +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
 466                         +y*0.392377e-2))))))));
 467         }
 468         return ans;
 469 }
 470 vector<float>
 471 gr_firdes::window (win_type type, int ntaps, double beta)
 472 {
 473   vector<float> taps(ntaps);
 474   int   M = ntaps - 1;          // filter order
 475
 476   switch (type){
 477   case WIN_RECTANGULAR:
 478     for (int n = 0; n < ntaps; n++)
 479       taps[n] = 1;
 480
 481   case WIN_HAMMING:
 482     for (int n = 0; n < ntaps; n++)
 483       taps[n] = 0.54 - 0.46 * cos ((2 * M_PI * n) / M);
 484     break;
 485
 486   case WIN_HANN:
 487     for (int n = 0; n < ntaps; n++)
 488       taps[n] = 0.5 - 0.5 * cos ((2 * M_PI * n) / M);
 489     break;
 490
 491   case WIN_BLACKMAN:
 492     for (int n = 0; n < ntaps; n++)
 493       taps[n] = 0.42 - 0.50 * cos ((2*M_PI * n) / (M-1)) - 0.08 * cos ((4*M_PI * n) / (M-1));
 494     break;
 495
 496 #if 0
 497   case WIN_KAISER:
 498     for (int n = 0; n < ntaps; n++)
 499         taps[n] = bessi0(beta*sqrt(1.0 - (4.0*n/(M*M))))/bessi0(beta);
 500     break;
 501 #else
 502
 503   case WIN_KAISER:
 504     {
 505       double IBeta = 1.0/Izero(beta);
 506       double inm1 = 1.0/((double)(ntaps));
 507       double temp;
 508       //fprintf(stderr, "IBeta = %g; inm1 = %g\n", IBeta, inm1);
 509
 510       for (int i=0; i<ntaps; i++) {
 511         temp = i * inm1;
 512         //fprintf(stderr, "temp = %g\n", temp);
 513         taps[i] = Izero(beta*sqrt(1.0-temp*temp)) * IBeta;
 514         //fprintf(stderr, "taps[%d] = %g\n", i, taps[i]);
 515       }
 516     }
 517     break;
 518
 519 #endif
 520   default:
 521     throw std::runtime_error ("not_implemented");
 522   }
 523
 524   return taps;
 525 }
 526
 527 void
 528 gr_firdes::sanity_check_1f (double sampling_freq,
 529                             double fa,                  // cutoff freq
 530                             double transition_width)
 531 {
 532   if (sampling_freq <= 0.0)
 533     throw std::out_of_range ("gr_firdes check failed: sampling_freq > 0");
 534
 535   if (fa <= 0.0 || fa > sampling_freq / 2)
 536     throw std::out_of_range ("gr_firdes check failed: 0 < fa <= sampling_freq / 2");
 537
 538   if (transition_width <= 0)
 539     throw std::out_of_range ("gr_dirdes check failed: transition_width > 0");
 540 }
 541
 542 void
 543 gr_firdes::sanity_check_2f (double sampling_freq,
 544                             double fa,                  // first cutoff freq
 545                             double fb,                  // second cutoff freq
 546                             double transition_width)
 547 {
 548   if (sampling_freq <= 0.0)
 549     throw std::out_of_range ("gr_firdes check failed: sampling_freq > 0");
 550
 551   if (fa <= 0.0 || fa > sampling_freq / 2)
 552     throw std::out_of_range ("gr_firdes check failed: 0 < fa <= sampling_freq / 2");
 553
 554   if (fb <= 0.0 || fb > sampling_freq / 2)
 555     throw std::out_of_range ("gr_firdes check failed: 0 < fb <= sampling_freq / 2");
 556
 557   if (fa > fb)
 558     throw std::out_of_range ("gr_firdes check failed: fa <= fb");
 559
 560   if (transition_width <= 0)
 561     throw std::out_of_range ("gr_firdes check failed: transition_width > 0");
 562 }
 563
 564 void
 565 gr_firdes::sanity_check_2f_c (double sampling_freq,
 566                             double fa,                  // first cutoff freq
 567                             double fb,                  // second cutoff freq
 568                             double transition_width)
 569 {
 570   if (sampling_freq <= 0.0)
 571     throw std::out_of_range ("gr_firdes check failed: sampling_freq > 0");
 572
 573   if (fa < -sampling_freq / 2 || fa > sampling_freq / 2)
 574     throw std::out_of_range ("gr_firdes check failed: 0 < fa <= sampling_freq / 2");
 575
 576   if (fb < -sampling_freq / 2  || fb > sampling_freq / 2)
 577     throw std::out_of_range ("gr_firdes check failed: 0 < fb <= sampling_freq / 2");
 578
 579   if (fa > fb)
 580     throw std::out_of_range ("gr_firdes check failed: fa <= fb");
 581
 582   if (transition_width <= 0)
 583     throw std::out_of_range ("gr_firdes check failed: transition_width > 0");
 584 }