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 x;
 202   int i;
 203   for (x = m_dFilterMin, i = 0; i < m_nFilterPoints; x += m_dFilterInc, i++)
 204     adFilter[i] = frequencyResponse (x);
 205 }
 206
 207
 208 void
 209 SignalFilter::createSpatialFilter (double* adFilter) const
 210 {
 211   if (m_idFilter == FILTER_SHEPP) {
 212     double a = 2 * m_dBandwidth;
 213     double c = - 4. / (a * a);
 214     int center = (m_nFilterPoints - 1) / 2;
 215     int sidelen = center;
 216     m_adFilter[center] = 4. / (a * a);
 217
 218     for (int i = 1; i <= sidelen; i++ )
 219       m_adFilter [center + i] = m_adFilter [center - i] = c / (4 * (i * i) - 1);
 220   } else {
 221     double x = m_dFilterMin;
 222     for (int i = 0; i < m_nFilterPoints; i++, x += m_dFilterInc) {
 223       if (haveAnalyticSpatial(m_idFilter))
 224         m_adFilter[i] = spatialResponseAnalytic (x);
 225       else
 226         m_adFilter[i] = spatialResponseCalc (x);
 227     }
 228   }
 229 }
 230
 231 int
 232 SignalFilter::convertFilterNameToID (const char *filterName)
 233 {
 234   int filterID = FILTER_INVALID;
 235
 236   for (int i = 0; i < s_iFilterCount; i++) {
 237     if (strcasecmp (filterName, s_aszFilterName[i]) == 0) {
 238       filterID = i;
 239       break;
 240     }
 241   }
 242
 243     return (filterID);
 244 }
 245
 246 const char *
 247 SignalFilter::convertFilterIDToName (const int filterID)
 248 {
 249   static const char *name = "";
 250
 251   if (filterID >= 0 && filterID < s_iFilterCount)
 252     return (s_aszFilterName [filterID]);
 253
 254   return (name);
 255 }
 256
 257 const char *
 258 SignalFilter::convertFilterIDToTitle (const int filterID)
 259 {
 260   static const char *title = "";
 261
 262   if (filterID >= 0 && filterID < s_iFilterCount) {
 263     return (s_aszFilterTitle [filterID]);
 264   }
 265   return (title);
 266 }
 267
 268 int
 269 SignalFilter::convertDomainNameToID (const char* const domainName)
 270 {
 271   int dID = DOMAIN_INVALID;
 272
 273   for (int i = 0; i < s_iDomainCount; i++) {
 274     if (strcasecmp (domainName, s_aszDomainName[i]) == 0) {
 275       dID = i;
 276       break;
 277     }
 278   }
 279   return (dID);
 280 }
 281
 282 const char *
 283 SignalFilter::convertDomainIDToName (const int domainID)
 284 {
 285   static const char *name = "";
 286
 287   if (domainID >= 0 && domainID < s_iDomainCount)
 288     return (s_aszDomainName [domainID]);
 289
 290   return (name);
 291 }
 292
 293 const char *
 294 SignalFilter::convertDomainIDToTitle (const int domainID)
 295 {
 296   static const char *title = "";
 297
 298   if (domainID >= 0 && domainID < s_iDomainCount)
 299     return (s_aszDomainTitle [domainID]);
 300
 301   return (title);
 302 }
 303
 304
 305 double
 306 SignalFilter::response (double x)
 307 {
 308   double response = 0;
 309
 310   if (m_idDomain == DOMAIN_SPATIAL)
 311     response = spatialResponse (m_idFilter, m_dBandwidth, x, m_dFilterParam);
 312   else if (m_idDomain == DOMAIN_FREQUENCY)
 313     response = frequencyResponse (m_idFilter, m_dBandwidth, x, m_dFilterParam);
 314
 315   return (response);
 316 }
 317
 318
 319 double
 320 SignalFilter::spatialResponse (int filterID, double bw, double x, double param)
 321 {
 322   if (haveAnalyticSpatial(filterID))
 323     return spatialResponseAnalytic (filterID, bw, x, param);
 324   else
 325     return spatialResponseCalc (filterID, bw, x, param, N_INTEGRAL);
 326 }
 327
 328 void
 329 SignalFilter::copyFilterData (double* pdFilter, const int iStart, const int nPoints) const
 330 {
 331   int iFirst = clamp (iStart, 0, m_nFilterPoints - 1);
 332   int iLast = clamp (iFirst + nPoints - 1, 0, m_nFilterPoints - 1);
 333
 334   for (int i = iFirst; i <= iLast; i++)
 335     pdFilter[i - iFirst] = m_adFilter[i];
 336 }
 337
 338 /* NAME
 339 *   filter_spatial_response_calc        Calculate filter by discrete inverse fourier
 340 *                                       transform of filters's frequency
 341 *                                       response
 342 *
 343 * SYNOPSIS
 344 *   y = filter_spatial_response_calc (filt_type, x, m_bw, param, n)
 345 *   double y                    Filter's response in spatial domain
 346 *   int filt_type               Type of filter (definitions in ct.h)
 347 *   double x                    Spatial position to evaluate filter
 348 *   double m_bw                 Bandwidth of window
 349 *   double param                General parameter for various filters
 350 *   int n                       Number of points to calculate integrations
 351 */
 352
 353 double
 354 SignalFilter::spatialResponseCalc (double x) const
 355 {
 356   return (spatialResponseCalc (m_idFilter, m_dBandwidth, x, m_dFilterParam, N_INTEGRAL));
 357 }
 358
 359 double
 360 SignalFilter::spatialResponseCalc (int filterID, double bw, double x, double param, int n)
 361 {
 362   double zmin, zmax;
 363
 364   if (filterID == FILTER_TRIANGLE) {
 365     zmin = 0;
 366     zmax = bw;
 367   } else {
 368     zmin = 0;
 369     zmax = bw / 2;
 370   }
 371   double zinc = (zmax - zmin) / (n - 1);
 372
 373   double z = zmin;
 374   double* q = new double [n];
 375   for (int i = 0; i < n; i++, z += zinc)
 376     q[i] = frequencyResponse (filterID, bw, z, param) * cos (TWOPI * z * x);
 377
 378   double y = 2 * integrateSimpson (zmin, zmax, q, n);
 379   delete q;
 380
 381   return (y);
 382 }
 383
 384
 385 /* NAME
 386 *    filter_frequency_response                  Return filter frequency response
 387 *
 388 * SYNOPSIS
 389 *    h = filter_frequency_response (filt_type, u, m_bw, param)
 390 *    double h                   Filters frequency response at u
 391 *    int filt_type              Type of filter
 392 *    double u                   Frequency to evaluate filter at
 393 *    double m_bw                        Bandwidth of filter
 394 *    double param               General input parameter for various filters
 395 */
 396
 397 double
 398 SignalFilter::frequencyResponse (double u) const
 399 {
 400   return frequencyResponse (m_idFilter, m_dBandwidth, u, m_dFilterParam);
 401 }
 402
 403
 404 double
 405 SignalFilter::frequencyResponse (int filterID, double bw, double u, double param)
 406 {
 407   double q;
 408   double au = fabs (u);
 409   double abw = fabs (bw);
 410
 411   switch (filterID) {
 412   case FILTER_BANDLIMIT:
 413     if (au >= (abw / 2) + F_EPSILON)
 414       q = 0.;
 415     else
 416       q = 1;
 417     break;
 418   case FILTER_ABS_BANDLIMIT:
 419     if (au >= (abw / 2) + F_EPSILON)
 420       q = 0.;
 421     else
 422       q = au;
 423     break;
 424   case FILTER_TRIANGLE:
 425     if (au >= (abw / 2) + F_EPSILON)
 426       q = 0;
 427     else
 428       q = 1 - au / abw;
 429     break;
 430   case FILTER_COSINE:
 431     if (au >= (abw / 2) + F_EPSILON)
 432       q = 0;
 433     else
 434       q = cos(PI * au / abw);
 435     break;
 436   case FILTER_ABS_COSINE:
 437     if (au >= (abw / 2) + F_EPSILON)
 438       q = 0;
 439     else
 440       q = au * cos(PI * au / abw);
 441     break;
 442   case FILTER_SINC:
 443     q = abw * sinc (PI * abw * au, 1.);
 444     break;
 445   case FILTER_ABS_SINC:
 446     if (au >= (abw / 2) + F_EPSILON)
 447       q = 0;
 448     else
 449       q = au * abw * sinc (PI * abw * au, 1.);
 450     break;
 451   case FILTER_HANNING:
 452     param = 0.5;
 453     // follow through to G_HAMMING
 454   case FILTER_G_HAMMING:
 455     if (au >= (abw / 2) + F_EPSILON)
 456       q = 0;
 457     else
 458       q = param + (1 - param) * cos (TWOPI * au / abw);
 459     break;
 460   case FILTER_ABS_HANNING:
 461     param = 0.5;
 462     // follow through to ABS_G_HAMMING
 463   case FILTER_ABS_G_HAMMING:
 464     if (au >= (abw / 2) + F_EPSILON)
 465       q = 0;
 466     else
 467       q = au * (param + (1 - param) * cos(TWOPI * au / abw));
 468     break;
 469   default:
 470     q = 0;
 471     sys_error (ERR_WARNING, "Frequency response for filter %d not implemented [filter_frequency_response]", filterID);
 472     break;
 473   }
 474
 475   return (q);
 476 }
 477
 478
 479
 480 /* NAME
 481 *   filter_spatial_response_analytic                    Calculate filter by analytic inverse fourier
 482 *                               transform of filters's frequency
 483 *                               response
 484 *
 485 * SYNOPSIS
 486 *   y = filter_spatial_response_analytic (filt_type, x, m_bw, param)
 487 *   double y                    Filter's response in spatial domain
 488 *   int filt_type               Type of filter (definitions in ct.h)
 489 *   double x                    Spatial position to evaluate filter
 490 *   double m_bw                 Bandwidth of window
 491 *   double param                General parameter for various filters
 492 */
 493
 494 double
 495 SignalFilter::spatialResponseAnalytic (double x) const
 496 {
 497   return spatialResponseAnalytic (m_idFilter, m_dBandwidth, x, m_dFilterParam);
 498 }
 499
 500 const bool
 501 SignalFilter::haveAnalyticSpatial (int filterID)
 502 {
 503   bool haveAnalytic = false;
 504
 505   switch (filterID) {
 506   case FILTER_BANDLIMIT:
 507   case FILTER_TRIANGLE:
 508   case FILTER_COSINE:
 509   case FILTER_G_HAMMING:
 510   case FILTER_HANNING:
 511   case FILTER_ABS_BANDLIMIT:
 512   case FILTER_ABS_COSINE:
 513   case FILTER_ABS_G_HAMMING:
 514   case FILTER_ABS_HANNING:
 515   case FILTER_SHEPP:
 516   case FILTER_SINC:
 517     haveAnalytic = true;
 518     break;
 519   default:
 520     break;
 521   }
 522
 523   return (haveAnalytic);
 524 }
 525
 526 double
 527 SignalFilter::spatialResponseAnalytic (int filterID, double bw, double x, double param)
 528 {
 529   double q, temp;
 530   double u = TWOPI * x;
 531   double w = bw / 2;
 532   double b = PI / bw;
 533   double b2 = TWOPI / bw;
 534
 535   switch (filterID) {
 536   case FILTER_BANDLIMIT:
 537     q = bw * sinc(u * w, 1.0);
 538     break;
 539   case FILTER_TRIANGLE:
 540     temp = sinc (u * w, 1.0);
 541     q = bw * temp * temp;
 542     break;
 543   case FILTER_COSINE:
 544     q = sinc(b-u,w) + sinc(b+u,w);
 545     break;
 546   case FILTER_HANNING:
 547     param = 0.5;
 548     // follow through to G_HAMMING
 549   case FILTER_G_HAMMING:
 550     q = 2 * param * sin(u*w)/u + (1-param) * (sinc(b2-u, w) + sinc(b2+u, w));
 551     break;
 552   case FILTER_ABS_BANDLIMIT:
 553     q = 2 * integral_abscos (u, w);
 554     break;
 555   case FILTER_ABS_COSINE:
 556     q = integral_abscos(b-u,w) + integral_abscos(b+u,w);
 557     break;
 558   case FILTER_ABS_HANNING:
 559     param = 0.5;
 560     // follow through to ABS_G_HAMMING
 561   case FILTER_ABS_G_HAMMING:
 562     q = 2 * param * integral_abscos(u,w) +
 563       (1-param)*(integral_abscos(u-b2,w)+integral_abscos(u+b2,w));
 564     break;
 565   case FILTER_SHEPP:
 566     if (fabs (u) < 1E-7)
 567       q = 4. / (PI * bw * bw);
 568     else
 569       q = fabs ((2 / bw) * sin (u * w)) * sinc (u * w, 1.) * sinc (u * w, 1.);
 570     break;
 571   case FILTER_SINC:
 572     if (fabs (x) < bw / 2)
 573       q = 1.;
 574     else
 575       q = 0.;
 576     break;
 577   case FILTER_ABS_SINC:
 578   default:
 579     sys_error (ERR_WARNING, "Analytic filter type %d not implemented [filter_spatial_response_analytic]", filterID);
 580     q = 0;
 581     break;
 582   }
 583
 584   return (q);
 585 }
 586
 587
 588
 589 // Functions that are inline in filter.h
 590
 591
 592 //  sinc                        Return sin(x)/x function
 593 //   v = sinc (x, mult)
 594 // Calculates sin(x * mult) / x;
 595
 596 //  integral_abscos     Returns integral of u*cos(u)
 597 //
 598 //   q = integral_abscos (u, w)
 599 //   double q                   Integral value
 600 //   double u                   Integration variable
 601 //   double w                   Upper integration boundary
 602 // Returns the value of integral of u*cos(u)*dV for V = 0 to w
 603
 604
 605