libctsim/filter.cpp

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