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