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