libctsim/filter.cpp

   1 /*****************************************************************************
   2 ** File IDENTIFICATION
   3 **
   4 **     Name:                   filter.cpp
   5 **     Purpose:                Routines for signal-procesing filters
   6 **     Progammer:              Kevin Rosenberg
   7 **     Date Started:           Aug 1984
   8 **
   9 **  This is part of the CTSim program
  10 **  Copyright (c) 1983-2000 Kevin Rosenberg
  11 **
  12 **  $Id: filter.cpp,v 1.36 2001/02/11 04:56:37 kevin Exp $
  13 **
  14 **  This program is free software; you can redistribute it and/or modify
  15 **  it under the terms of the GNU General Public License (version 2) as
  16 **  published by the Free Software Foundation.
  17 **
  18 **  This program is distributed in the hope that it will be useful,
  19 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  20 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  21 **  GNU General Public License for more details.
  22 **
  23 **  You should have received a copy of the GNU General Public License
  24 **  along with this program; if not, write to the Free Software
  25 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  26 ******************************************************************************/
  27
  28 #include "ct.h"
  29
  30 int SignalFilter::N_INTEGRAL=500;  //static member
  31
  32 const int SignalFilter::FILTER_INVALID = -1 ;
  33 const int SignalFilter::FILTER_ABS_BANDLIMIT = 0;       // filter times |x|
  34 const int SignalFilter::FILTER_ABS_G_HAMMING = 1;
  35 const int SignalFilter::FILTER_ABS_COSINE = 2;
  36 const int SignalFilter::FILTER_ABS_SINC = 3;
  37 const int SignalFilter::FILTER_SHEPP = 4;
  38 const int SignalFilter::FILTER_BANDLIMIT = 5;
  39 const int SignalFilter::FILTER_SINC = 6;
  40 const int SignalFilter::FILTER_G_HAMMING = 7;
  41 const int SignalFilter::FILTER_COSINE = 8;
  42 const int SignalFilter::FILTER_TRIANGLE = 9;
  43
  44 const int SignalFilter::s_iReconstructFilterCount = 5;
  45
  46 const char* SignalFilter::s_aszFilterName[] = {
  47   {"abs_bandlimit"},
  48   {"abs_hamming"},
  49   {"abs_cosine"},
  50   {"abs_sinc"},
  51   {"shepp"},
  52   {"bandlimit"},
  53   {"sinc"},
  54   {"hamming"},
  55   {"cosine"},
  56   {"triangle"},
  57 };
  58
  59 const char* SignalFilter::s_aszFilterTitle[] = {
  60   {"Abs(w) * Bandlimit"},
  61   {"Abs(w) * Hamming"},
  62   {"Abs(w) * Cosine"},
  63   {"Abs(w) * Sinc"},
  64   {"Shepp"},
  65   {"Bandlimit"},
  66   {"Sinc"},
  67   {"Hamming"},
  68   {"Cosine"},
  69   {"Triangle"},
  70 };
  71
  72 const int SignalFilter::s_iFilterCount = sizeof(s_aszFilterName) / sizeof(const char*);
  73
  74
  75 const int SignalFilter::DOMAIN_INVALID = -1;
  76 const int SignalFilter::DOMAIN_FREQUENCY = 0;
  77 const int SignalFilter::DOMAIN_SPATIAL = 1;
  78
  79 const char* SignalFilter::s_aszDomainName[] = {
  80   {"frequency"},
  81   {"spatial"},
  82 };
  83
  84 const char* SignalFilter::s_aszDomainTitle[] = {
  85   {"Frequency"},
  86   {"Spatial"},
  87 };
  88
  89 const int SignalFilter::s_iDomainCount = sizeof(s_aszDomainName) / sizeof(const char*);
  90
  91
  92 /* NAME
  93  *   SignalFilter::SignalFilter     Construct a signal
  94  *
  95  * SYNOPSIS
  96  *   f = SignalFilter (filt_type, bw, filterMin, filterMax, n, param, domain, analytic)
  97  *   double f           Generated filter vector
  98  *   int filt_type      Type of filter wanted
  99  *   double bw          Bandwidth of filter
 100  *   double filterMin, filterMax        Filter limits
 101  *   int nFilterPoints  Number of points in signal
 102  *   double param       General input parameter to filters
 103  *   int domain         FREQUENCY or SPATIAL domain wanted
 104  */
 105
 106 SignalFilter::SignalFilter (const char* szFilterName, double dFilterMinimum, double dFilterMaximum, int nFilterPoints, double dBandwidth, double dFilterParam, const char* szDomainName)
 107   : m_adFilter(NULL), m_fail(false)
 108 {
 109   m_idFilter = convertFilterNameToID (szFilterName);
 110   if (m_idFilter == FILTER_INVALID) {
 111     m_fail = true;
 112     m_failMessage = "Invalid Filter name ";
 113     m_failMessage += szFilterName;
 114     return;
 115   }
 116   m_idDomain = convertDomainNameToID (szDomainName);
 117   if (m_idDomain == DOMAIN_INVALID) {
 118     m_fail = true;
 119     m_failMessage = "Invalid domain name ";
 120     m_failMessage += szDomainName;
 121     return;
 122   }
 123   init (m_idFilter, dFilterMinimum, dFilterMaximum, nFilterPoints, dBandwidth, dFilterParam, m_idDomain);
 124 }
 125
 126 SignalFilter::SignalFilter (const int idFilter, double dFilterMinimum, double dFilterMaximum, int nFilterPoints, double dBandwidth, double dFilterParam, const int idDomain)
 127   : m_adFilter(NULL), m_fail(false)
 128 {
 129   init (idFilter, dFilterMinimum, dFilterMaximum, nFilterPoints, dBandwidth, dFilterParam, idDomain);
 130 }
 131
 132 SignalFilter::SignalFilter (const char* szFilterName, const char* szDomainName, double dBandwidth, double dFilterParam)
 133   : m_adFilter(NULL), m_fail(false)
 134 {
 135   m_nFilterPoints = 0;
 136   m_dBandwidth = dBandwidth;
 137   m_dFilterParam = dFilterParam;
 138   m_idFilter = convertFilterNameToID (szFilterName);
 139   if (m_idFilter == FILTER_INVALID) {
 140     m_fail = true;
 141     m_failMessage = "Invalid Filter name ";
 142     m_failMessage += szFilterName;
 143     return;
 144   }
 145   m_idDomain = convertDomainNameToID (szDomainName);
 146   if (m_idDomain == DOMAIN_INVALID) {
 147     m_fail = true;
 148     m_failMessage = "Invalid domain name ";
 149     m_failMessage += szDomainName;
 150     return;
 151   }
 152 }
 153
 154 void
 155 SignalFilter::init (const int idFilter, double dFilterMinimum, double dFilterMaximum, int nFilterPoints, double dBandwidth, double dFilterParam, const int idDomain)
 156 {
 157   m_idFilter = idFilter;
 158   m_idDomain = idDomain;
 159   if (m_idFilter == FILTER_INVALID || m_idDomain == DOMAIN_INVALID) {
 160     m_fail = true;
 161     return;
 162   }
 163   if (nFilterPoints < 2) {
 164     m_fail = true;
 165     m_failMessage = "Number of filter points ";
 166     m_failMessage += nFilterPoints;
 167     m_failMessage = " less than 2";
 168     return;
 169   }
 170
 171   m_nameFilter = convertFilterIDToName (m_idFilter);
 172   m_nameDomain = convertDomainIDToName (m_idDomain);
 173   m_nFilterPoints = nFilterPoints;
 174   m_dFilterParam = dFilterParam;
 175   m_dBandwidth = dBandwidth;
 176   m_dFilterMin = dFilterMinimum;
 177   m_dFilterMax = dFilterMaximum;
 178
 179   m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
 180   m_adFilter = new double [m_nFilterPoints];
 181
 182   if (m_idDomain == DOMAIN_FREQUENCY)
 183       createFrequencyFilter (m_adFilter);
 184   else if (m_idDomain == DOMAIN_SPATIAL)
 185       createSpatialFilter (m_adFilter);
 186 }
 187
 188
 189 SignalFilter::~SignalFilter (void)
 190 {
 191     delete [] m_adFilter;
 192 }
 193
 194 void
 195 SignalFilter::createFrequencyFilter (double* adFilter) const
 196 {
 197   double x;
 198   int i;
 199   for (x = m_dFilterMin, i = 0; i < m_nFilterPoints; x += m_dFilterInc, i++)
 200     adFilter[i] = frequencyResponse (x);
 201 }
 202
 203
 204 void
 205 SignalFilter::createSpatialFilter (double* adFilter) const
 206 {
 207   if (m_idFilter == FILTER_SHEPP) {
 208     double a = 2 * m_dBandwidth;
 209     double c = - 4. / (a * a);
 210     int center = (m_nFilterPoints - 1) / 2;
 211     int sidelen = center;
 212     m_adFilter[center] = 4. / (a * a);
 213
 214     for (int i = 1; i <= sidelen; i++ )
 215       m_adFilter [center + i] = m_adFilter [center - i] = c / (4 * (i * i) - 1);
 216   } else {
 217     double x = m_dFilterMin;
 218     for (int i = 0; i < m_nFilterPoints; i++, x += m_dFilterInc) {
 219       if (haveAnalyticSpatial(m_idFilter))
 220         m_adFilter[i] = spatialResponseAnalytic (x);
 221       else
 222         m_adFilter[i] = spatialResponseCalc (x);
 223     }
 224   }
 225 }
 226
 227 int
 228 SignalFilter::convertFilterNameToID (const char *filterName)
 229 {
 230   int filterID = FILTER_INVALID;
 231
 232   for (int i = 0; i < s_iFilterCount; i++)
 233     if (strcasecmp (filterName, s_aszFilterName[i]) == 0) {
 234       filterID = i;
 235       break;
 236     }
 237
 238   return (filterID);
 239 }
 240
 241 const char *
 242 SignalFilter::convertFilterIDToName (const int filterID)
 243 {
 244   static const char *name = "";
 245
 246   if (filterID >= 0 && filterID < s_iFilterCount)
 247       return (s_aszFilterName [filterID]);
 248
 249   return (name);
 250 }
 251
 252 const char *
 253 SignalFilter::convertFilterIDToTitle (const int filterID)
 254 {
 255   static const char *title = "";
 256
 257   if (filterID >= 0 && filterID < s_iFilterCount)
 258       return (s_aszFilterTitle [filterID]);
 259
 260   return (title);
 261 }
 262
 263 int
 264 SignalFilter::convertDomainNameToID (const char* const domainName)
 265 {
 266   int dID = DOMAIN_INVALID;
 267
 268   for (int i = 0; i < s_iDomainCount; i++)
 269    if (strcasecmp (domainName, s_aszDomainName[i]) == 0) {
 270       dID = i;
 271       break;
 272     }
 273
 274   return (dID);
 275 }
 276
 277 const char *
 278 SignalFilter::convertDomainIDToName (const int domainID)
 279 {
 280   static const char *name = "";
 281
 282   if (domainID >= 0 && domainID < s_iDomainCount)
 283       return (s_aszDomainName [domainID]);
 284
 285   return (name);
 286 }
 287
 288 const char *
 289 SignalFilter::convertDomainIDToTitle (const int domainID)
 290 {
 291   static const char *title = "";
 292
 293   if (domainID >= 0 && domainID < s_iDomainCount)
 294       return (s_aszDomainTitle [domainID]);
 295
 296   return (title);
 297 }
 298
 299
 300 double
 301 SignalFilter::response (double x)
 302 {
 303   double response = 0;
 304
 305   if (m_idDomain == DOMAIN_SPATIAL)
 306     response = spatialResponse (m_idFilter, m_dBandwidth, x, m_dFilterParam);
 307   else if (m_idDomain == DOMAIN_FREQUENCY)
 308     response = frequencyResponse (m_idFilter, m_dBandwidth, x, m_dFilterParam);
 309
 310   return (response);
 311 }
 312
 313
 314 double
 315 SignalFilter::spatialResponse (int filterID, double bw, double x, double param)
 316 {
 317   if (haveAnalyticSpatial(filterID))
 318     return spatialResponseAnalytic (filterID, bw, x, param);
 319   else
 320     return spatialResponseCalc (filterID, bw, x, param, N_INTEGRAL);
 321 }
 322
 323 void
 324 SignalFilter::copyFilterData (double* pdFilter, const int iStart, const int nPoints) const
 325 {
 326     int iFirst = clamp (iStart, 0, m_nFilterPoints - 1);
 327     int iLast = clamp (iFirst + nPoints - 1, 0, m_nFilterPoints - 1);
 328
 329     for (int i = iFirst; i <= iLast; i++)
 330         pdFilter[i - iFirst] = m_adFilter[i];
 331 }
 332
 333 /* NAME
 334  *   filter_spatial_response_calc       Calculate filter by discrete inverse fourier
 335  *                                      transform of filters's frequency
 336  *                                      response
 337  *
 338  * SYNOPSIS
 339  *   y = filter_spatial_response_calc (filt_type, x, m_bw, param, n)
 340  *   double y                   Filter's response in spatial domain
 341  *   int filt_type              Type of filter (definitions in ct.h)
 342  *   double x                   Spatial position to evaluate filter
 343  *   double m_bw                        Bandwidth of window
 344  *   double param               General parameter for various filters
 345  *   int n                      Number of points to calculate integrations
 346  */
 347
 348 double
 349 SignalFilter::spatialResponseCalc (double x) const
 350 {
 351   return (spatialResponseCalc (m_idFilter, m_dBandwidth, x, m_dFilterParam, N_INTEGRAL));
 352 }
 353
 354 double
 355 SignalFilter::spatialResponseCalc (int filterID, double bw, double x, double param, int n)
 356 {
 357   double zmin, zmax;
 358
 359   if (filterID == FILTER_TRIANGLE) {
 360     zmin = 0;
 361     zmax = bw;
 362   } else {
 363     zmin = 0;
 364     zmax = bw / 2;
 365   }
 366   double zinc = (zmax - zmin) / (n - 1);
 367
 368   double z = zmin;
 369   double* q = new double [n];
 370   for (int i = 0; i < n; i++, z += zinc)
 371     q[i] = frequencyResponse (filterID, bw, z, param) * cos (TWOPI * z * x);
 372
 373   double y = 2 * integrateSimpson (zmin, zmax, q, n);
 374   delete q;
 375
 376   return (y);
 377 }
 378
 379
 380 /* NAME
 381  *    filter_frequency_response                 Return filter frequency response
 382  *
 383  * SYNOPSIS
 384  *    h = filter_frequency_response (filt_type, u, m_bw, param)
 385  *    double h                  Filters frequency response at u
 386  *    int filt_type             Type of filter
 387  *    double u                  Frequency to evaluate filter at
 388  *    double m_bw                       Bandwidth of filter
 389  *    double param              General input parameter for various filters
 390  */
 391
 392 double
 393 SignalFilter::frequencyResponse (double u) const
 394 {
 395   return frequencyResponse (m_idFilter, m_dBandwidth, u, m_dFilterParam);
 396 }
 397
 398
 399 double
 400 SignalFilter::frequencyResponse (int filterID, double bw, double u, double param)
 401 {
 402   double q;
 403   double au = fabs (u);
 404
 405   switch (filterID) {
 406   case FILTER_BANDLIMIT:
 407     if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
 408       q = 0.;
 409     else
 410       q = 1;
 411     break;
 412   case FILTER_ABS_BANDLIMIT:
 413     if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
 414       q = 0.;
 415     else
 416       q = au;
 417     break;
 418   case FILTER_TRIANGLE:
 419     if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
 420       q = 0;
 421     else
 422       q = 1 - au / bw;
 423     break;
 424   case FILTER_COSINE:
 425     if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
 426       q = 0;
 427     else
 428       q = cos(PI * u / bw);
 429     break;
 430   case FILTER_ABS_COSINE:
 431     if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
 432       q = 0;
 433     else
 434       q = au * cos(PI * u / bw);
 435     break;
 436   case FILTER_SINC:
 437     q = bw * sinc (PI * bw * u, 1.);
 438     break;
 439   case FILTER_ABS_SINC:
 440     q = au * bw * sinc (PI * bw * u, 1.);
 441     break;
 442   case FILTER_G_HAMMING:
 443     if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
 444       q = 0;
 445     else
 446       q = param + (1 - param) * cos (TWOPI * u / bw);
 447     break;
 448   case FILTER_ABS_G_HAMMING:
 449     if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
 450       q = 0;
 451     else
 452       q = au * (param + (1 - param) * cos(TWOPI * u / bw));
 453     break;
 454   default:
 455     q = 0;
 456     sys_error (ERR_WARNING, "Frequency response for filter %d not implemented [filter_frequency_response]", filterID);
 457     break;
 458   }
 459   return (q);
 460 }
 461
 462
 463
 464 /* NAME
 465  *   filter_spatial_response_analytic                   Calculate filter by analytic inverse fourier
 466  *                              transform of filters's frequency
 467  *                              response
 468  *
 469  * SYNOPSIS
 470  *   y = filter_spatial_response_analytic (filt_type, x, m_bw, param)
 471  *   double y                   Filter's response in spatial domain
 472  *   int filt_type              Type of filter (definitions in ct.h)
 473  *   double x                   Spatial position to evaluate filter
 474  *   double m_bw                        Bandwidth of window
 475  *   double param               General parameter for various filters
 476  */
 477
 478 double
 479 SignalFilter::spatialResponseAnalytic (double x) const
 480 {
 481   return spatialResponseAnalytic (m_idFilter, m_dBandwidth, x, m_dFilterParam);
 482 }
 483
 484 const bool
 485 SignalFilter::haveAnalyticSpatial (int filterID)
 486 {
 487   bool haveAnalytic = false;
 488
 489   switch (filterID) {
 490   case FILTER_BANDLIMIT:
 491   case FILTER_TRIANGLE:
 492   case FILTER_COSINE:
 493   case FILTER_G_HAMMING:
 494   case FILTER_ABS_BANDLIMIT:
 495   case FILTER_ABS_COSINE:
 496   case FILTER_ABS_G_HAMMING:
 497   case FILTER_SHEPP:
 498   case FILTER_SINC:
 499     haveAnalytic = true;
 500     break;
 501   default:
 502     break;
 503   }
 504
 505   return (haveAnalytic);
 506 }
 507
 508 double
 509 SignalFilter::spatialResponseAnalytic (int filterID, double bw, double x, double param)
 510 {
 511   double q, temp;
 512   double u = TWOPI * x;
 513   double w = bw / 2;
 514   double b = PI / bw;
 515   double b2 = TWOPI / bw;
 516
 517   switch (filterID) {
 518   case FILTER_BANDLIMIT:
 519     q = bw * sinc(u * w, 1.0);
 520     break;
 521   case FILTER_TRIANGLE:
 522     temp = sinc (u * w, 1.0);
 523     q = bw * temp * temp;
 524     break;
 525   case FILTER_COSINE:
 526     q = sinc(b-u,w) + sinc(b+u,w);
 527     break;
 528   case FILTER_G_HAMMING:
 529     q = 2 * param * sin(u*w)/u + (1-param) * (sinc(b2-u, w) + sinc(b2+u, w));
 530     break;
 531   case FILTER_ABS_BANDLIMIT:
 532     q = 2 * integral_abscos (u, w);
 533     break;
 534   case FILTER_ABS_COSINE:
 535     q = integral_abscos(b-u,w) + integral_abscos(b+u,w);
 536     break;
 537   case FILTER_ABS_G_HAMMING:
 538     q = 2 * param * integral_abscos(u,w) +
 539       (1-param)*(integral_abscos(u-b2,w)+integral_abscos(u+b2,w));
 540     break;
 541   case FILTER_SHEPP:
 542     if (fabs (u) < 1E-7)
 543       q = 4. / (PI * bw * bw);
 544     else
 545       q = fabs ((2 / bw) * sin (u * w)) * sinc (u * w, 1.) * sinc (u * w, 1.);
 546     break;
 547   case FILTER_SINC:
 548     if (fabs (x) < bw / 2)
 549       q = 1.;
 550     else
 551       q = 0.;
 552     break;
 553   case FILTER_ABS_SINC:
 554   default:
 555     sys_error (ERR_WARNING, "Analytic filter type %d not implemented [filter_spatial_response_analytic]", filterID);
 556     q = 0;
 557     break;
 558   }
 559
 560   return (q);
 561 }
 562
 563
 564
 565 // Functions that are inline in filter.h
 566
 567
 568 //  sinc                        Return sin(x)/x function
 569 //   v = sinc (x, mult)
 570 // Calculates sin(x * mult) / x;
 571
 572 //  integral_abscos     Returns integral of u*cos(u)
 573 //
 574 //   q = integral_abscos (u, w)
 575 //   double q                   Integral value
 576 //   double u                   Integration variable
 577 //   double w                   Upper integration boundary
 578 // Returns the value of integral of u*cos(u)*dV for V = 0 to w
 579
 580
 581