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