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.37 2001/02/16 00:28:41 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* 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* 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* SignalFilter::s_aszDomainName[] = {
  86   {"frequency"},
  87   {"spatial"},
  88 };
  89
  90 const char* 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
 411   switch (filterID) {
 412   case FILTER_BANDLIMIT:
 413     if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
 414       q = 0.;
 415     else
 416       q = 1;
 417     break;
 418   case FILTER_ABS_BANDLIMIT:
 419     if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
 420       q = 0.;
 421     else
 422       q = au;
 423     break;
 424   case FILTER_TRIANGLE:
 425     if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
 426       q = 0;
 427     else
 428       q = 1 - au / bw;
 429     break;
 430   case FILTER_COSINE:
 431     if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
 432       q = 0;
 433     else
 434       q = cos(PI * u / bw);
 435     break;
 436   case FILTER_ABS_COSINE:
 437     if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
 438       q = 0;
 439     else
 440       q = au * cos(PI * u / bw);
 441     break;
 442   case FILTER_SINC:
 443     q = bw * sinc (PI * bw * u, 1.);
 444     break;
 445   case FILTER_ABS_SINC:
 446     q = au * bw * sinc (PI * bw * u, 1.);
 447     break;
 448   case FILTER_HANNING:
 449     param = 0.54;
 450     // follow through to G_HAMMING
 451   case FILTER_G_HAMMING:
 452     if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
 453       q = 0;
 454     else
 455       q = param + (1 - param) * cos (TWOPI * u / bw);
 456     break;
 457   case FILTER_ABS_HANNING:
 458     param = 0.54;
 459     // follow through to ABS_G_HAMMING
 460   case FILTER_ABS_G_HAMMING:
 461     if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
 462       q = 0;
 463     else
 464       q = au * (param + (1 - param) * cos(TWOPI * u / bw));
 465     break;
 466   default:
 467     q = 0;
 468     sys_error (ERR_WARNING, "Frequency response for filter %d not implemented [filter_frequency_response]", filterID);
 469     break;
 470   }
 471   return (q);
 472 }
 473
 474
 475
 476 /* NAME
 477  *   filter_spatial_response_analytic                   Calculate filter by analytic inverse fourier
 478  *                              transform of filters's frequency
 479  *                              response
 480  *
 481  * SYNOPSIS
 482  *   y = filter_spatial_response_analytic (filt_type, x, m_bw, param)
 483  *   double y                   Filter's response in spatial domain
 484  *   int filt_type              Type of filter (definitions in ct.h)
 485  *   double x                   Spatial position to evaluate filter
 486  *   double m_bw                        Bandwidth of window
 487  *   double param               General parameter for various filters
 488  */
 489
 490 double
 491 SignalFilter::spatialResponseAnalytic (double x) const
 492 {
 493   return spatialResponseAnalytic (m_idFilter, m_dBandwidth, x, m_dFilterParam);
 494 }
 495
 496 const bool
 497 SignalFilter::haveAnalyticSpatial (int filterID)
 498 {
 499   bool haveAnalytic = false;
 500
 501   switch (filterID) {
 502   case FILTER_BANDLIMIT:
 503   case FILTER_TRIANGLE:
 504   case FILTER_COSINE:
 505   case FILTER_G_HAMMING:
 506   case FILTER_HANNING:
 507   case FILTER_ABS_BANDLIMIT:
 508   case FILTER_ABS_COSINE:
 509   case FILTER_ABS_G_HAMMING:
 510   case FILTER_ABS_HANNING:
 511   case FILTER_SHEPP:
 512   case FILTER_SINC:
 513     haveAnalytic = true;
 514     break;
 515   default:
 516     break;
 517   }
 518
 519   return (haveAnalytic);
 520 }
 521
 522 double
 523 SignalFilter::spatialResponseAnalytic (int filterID, double bw, double x, double param)
 524 {
 525   double q, temp;
 526   double u = TWOPI * x;
 527   double w = bw / 2;
 528   double b = PI / bw;
 529   double b2 = TWOPI / bw;
 530
 531   switch (filterID) {
 532   case FILTER_BANDLIMIT:
 533     q = bw * sinc(u * w, 1.0);
 534     break;
 535   case FILTER_TRIANGLE:
 536     temp = sinc (u * w, 1.0);
 537     q = bw * temp * temp;
 538     break;
 539   case FILTER_COSINE:
 540     q = sinc(b-u,w) + sinc(b+u,w);
 541     break;
 542   case FILTER_HANNING:
 543     param = 0.54;
 544     // follow through to G_HAMMING
 545   case FILTER_G_HAMMING:
 546     q = 2 * param * sin(u*w)/u + (1-param) * (sinc(b2-u, w) + sinc(b2+u, w));
 547     break;
 548   case FILTER_ABS_BANDLIMIT:
 549     q = 2 * integral_abscos (u, w);
 550     break;
 551   case FILTER_ABS_COSINE:
 552     q = integral_abscos(b-u,w) + integral_abscos(b+u,w);
 553     break;
 554   case FILTER_ABS_HANNING:
 555     param = 0.54;
 556     // follow through to ABS_G_HAMMING
 557   case FILTER_ABS_G_HAMMING:
 558     q = 2 * param * integral_abscos(u,w) +
 559       (1-param)*(integral_abscos(u-b2,w)+integral_abscos(u+b2,w));
 560     break;
 561   case FILTER_SHEPP:
 562     if (fabs (u) < 1E-7)
 563       q = 4. / (PI * bw * bw);
 564     else
 565       q = fabs ((2 / bw) * sin (u * w)) * sinc (u * w, 1.) * sinc (u * w, 1.);
 566     break;
 567   case FILTER_SINC:
 568     if (fabs (x) < bw / 2)
 569       q = 1.;
 570     else
 571       q = 0.;
 572     break;
 573   case FILTER_ABS_SINC:
 574   default:
 575     sys_error (ERR_WARNING, "Analytic filter type %d not implemented [filter_spatial_response_analytic]", filterID);
 576     q = 0;
 577     break;
 578   }
 579
 580   return (q);
 581 }
 582
 583
 584
 585 // Functions that are inline in filter.h
 586
 587
 588 //  sinc                        Return sin(x)/x function
 589 //   v = sinc (x, mult)
 590 // Calculates sin(x * mult) / x;
 591
 592 //  integral_abscos     Returns integral of u*cos(u)
 593 //
 594 //   q = integral_abscos (u, w)
 595 //   double q                   Integral value
 596 //   double u                   Integration variable
 597 //   double w                   Upper integration boundary
 598 // Returns the value of integral of u*cos(u)*dV for V = 0 to w
 599
 600
 601