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