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.22 2000/07/23 01:49:03 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_SINC = 1;
  35 const int SignalFilter::FILTER_ABS_G_HAMMING = 2;
  36 const int SignalFilter::FILTER_ABS_COSINE = 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_sinc"},
  47   {"abs_hamming"},
  48   {"abs_cosine"},
  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) * Sinc"},
  60   {"Abs(w) * Hamming"},
  61   {"Abs(w) * Cosine"},
  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 const int SignalFilter::FILTER_METHOD_INVALID = -1;
  73 const int SignalFilter::FILTER_METHOD_CONVOLUTION = 0;
  74 const int SignalFilter::FILTER_METHOD_FOURIER = 1;
  75 const int SignalFilter::FILTER_METHOD_FOURIER_TABLE = 2;
  76 const int SignalFilter::FILTER_METHOD_FFT = 3;
  77 #if HAVE_FFTW
  78 const int SignalFilter::FILTER_METHOD_FFTW = 4;
  79 const int SignalFilter::FILTER_METHOD_RFFTW =5 ;
  80 #endif
  81
  82 const char* SignalFilter::s_aszFilterMethodName[] = {
  83   {"convolution"},
  84   {"fourier"},
  85   {"fouier_table"},
  86   {"fft"},
  87 #if HAVE_FFTW
  88   {"fftw"},
  89   {"rfftw"},
  90 #endif
  91 };
  92
  93 const char* SignalFilter::s_aszFilterMethodTitle[] = {
  94   {"Convolution"},
  95   {"Direct Fourier"},
  96   {"Fouier Trigometric Table Lookout"},
  97   {"FFT"},
  98 #if HAVE_FFTW
  99   {"FFTW"},
 100   {"Real/Half-Complex FFTW"},
 101 #endif
 102 };
 103
 104 const int SignalFilter::s_iFilterMethodCount = sizeof(s_aszFilterMethodName) / sizeof(const char*);
 105
 106
 107 const int SignalFilter::DOMAIN_INVALID = -1;
 108 const int SignalFilter::DOMAIN_FREQUENCY = 0;
 109 const int SignalFilter::DOMAIN_SPATIAL = 1;
 110
 111 const char* SignalFilter::s_aszDomainName[] = {
 112   {"frequency"},
 113   {"spatial"},
 114 };
 115
 116 const char* SignalFilter::s_aszDomainTitle[] = {
 117   {"Frequency"},
 118   {"Spatial"},
 119 };
 120
 121 const int SignalFilter::s_iDomainCount = sizeof(s_aszDomainName) / sizeof(const char*);
 122
 123
 124 /* NAME
 125  *   SignalFilter::SignalFilter     Construct a signal
 126  *
 127  * SYNOPSIS
 128  *   f = SignalFilter (filt_type, bw, filterMin, filterMax, n, param, domain, analytic)
 129  *   double f           Generated filter vector
 130  *   int filt_type      Type of filter wanted
 131  *   double bw          Bandwidth of filter
 132  *   double filterMin, filterMax        Filter limits
 133  *   int nSignalPoints  Number of points in signal
 134  *   double param       General input parameter to filters
 135  *   int domain         FREQUENCY or SPATIAL domain wanted
 136  */
 137
 138 SignalFilter::SignalFilter (const char* filterName, const char* filterMethodName, double bw, double signalIncrement, int nSignalPoints, double param, const char* domainName, int zeropad = 0, int preinterpolationFactor = 1)
 139   : m_vecFilter(NULL), m_vecFourierCosTable(NULL), m_vecFourierSinTable(NULL), m_fail(false)
 140 {
 141   m_idFilter = convertFilterNameToID (filterName);
 142   if (m_idFilter == FILTER_INVALID) {
 143     m_fail = true;
 144     m_failMessage = "Invalid Filter name ";
 145     m_failMessage += filterName;
 146     return;
 147   }
 148   m_idFilterMethod = convertFilterMethodNameToID (filterMethodName);
 149   if (m_idFilterMethod == FILTER_METHOD_INVALID) {
 150     m_fail = true;
 151     m_failMessage = "Invalid filter method name ";
 152     m_failMessage += filterMethodName;
 153     return;
 154   }
 155   m_idDomain = convertDomainNameToID (domainName);
 156   if (m_idDomain == DOMAIN_INVALID) {
 157     m_fail = true;
 158     m_failMessage = "Invalid domain name ";
 159     m_failMessage += domainName;
 160     return;
 161   }
 162   init (m_idFilter, m_idFilterMethod, bw, signalIncrement, nSignalPoints, param, m_idDomain, zeropad, preinterpolationFactor);
 163 }
 164
 165 SignalFilter::SignalFilter (const int filterID, const int filterMethodID, double bw, double signalIncrement, int nSignalPoints, double param, const int domainID, int zeropad = 0, int preinterpolationFactor = 1)
 166   : m_vecFilter(NULL), m_vecFourierCosTable(NULL), m_vecFourierSinTable(NULL), m_fail(false)
 167 {
 168   init (filterID, filterMethodID, bw, signalIncrement, nSignalPoints, param, domainID, zeropad, preinterpolationFactor);
 169 }
 170
 171 SignalFilter::SignalFilter (const char* filterName, const char* domainName, double bw, double param)
 172   : m_vecFilter(NULL), m_vecFourierCosTable(NULL), m_vecFourierSinTable(NULL), m_fail(false)
 173 {
 174   m_bw = bw;
 175   m_nSignalPoints = 0;
 176   m_nFilterPoints = 0;
 177   m_filterParam = param;
 178   m_idFilter = convertFilterNameToID (filterName);
 179   if (m_idFilter == FILTER_INVALID) {
 180     m_fail = true;
 181     m_failMessage = "Invalid Filter name ";
 182     m_failMessage += filterName;
 183     return;
 184   }
 185   m_idDomain = convertDomainNameToID (domainName);
 186   if (m_idDomain == DOMAIN_INVALID) {
 187     m_fail = true;
 188     m_failMessage = "Invalid domain name ";
 189     m_failMessage += domainName;
 190     return;
 191   }
 192 }
 193
 194 void
 195 SignalFilter::init (const int filterID, const int filterMethodID, double bw, double signalIncrement, int nSignalPoints, double filterParam, const int domainID, int zeropad, int preinterpolationFactor)
 196 {
 197   m_bw = bw;
 198   m_idFilter = filterID;
 199   m_idDomain = domainID;
 200   m_idFilterMethod = filterMethodID;
 201   if (m_idFilter == FILTER_INVALID || m_idDomain == DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID) {
 202     m_fail = true;
 203     return;
 204   }
 205   m_traceLevel = TRACE_NONE;
 206   m_nameFilter = convertFilterIDToName (m_idFilter);
 207   m_nameDomain = convertDomainIDToName (m_idDomain);
 208   m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod);
 209   m_nSignalPoints = nSignalPoints;
 210   m_signalInc = signalIncrement;
 211   m_filterParam = filterParam;
 212   m_zeropad = zeropad;
 213   m_preinterpolationFactor = preinterpolationFactor;
 214
 215   m_vecFourierCosTable = NULL;
 216   m_vecFourierSinTable = NULL;
 217   m_vecFilter = NULL;
 218
 219   if (m_idFilterMethod == FILTER_METHOD_FFT) {
 220 #if HAVE_FFTW
 221     m_idFilterMethod = FILTER_METHOD_RFFTW;
 222 #else
 223     m_fail = true;
 224     m_failMessage = "FFT not yet implemented";
 225     return;
 226 #endif
 227   }
 228
 229   if (m_idFilterMethod == FILTER_METHOD_FOURIER || m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE || m_idFilterMethod == FILTER_METHOD_FFT
 230 #if HAVE_FFTW
 231       || m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW
 232 #endif
 233       ) {
 234     m_nFilterPoints = m_nSignalPoints;
 235     if (m_zeropad > 0) {
 236       double logBase2 = log(m_nSignalPoints) / log(2);
 237       int nextPowerOf2 = static_cast<int>(floor(logBase2));
 238       if (logBase2 != floor(logBase2))
 239         nextPowerOf2++;
 240       nextPowerOf2 += (m_zeropad - 1);
 241       m_nFilterPoints = 1 << nextPowerOf2;
 242       if (m_traceLevel >= TRACE_TEXT)
 243         cout << "nFilterPoints = " << m_nFilterPoints << endl;
 244     }
 245     m_nOutputPoints = m_nFilterPoints * m_preinterpolationFactor;
 246     m_filterMin = -1. / (2 * m_signalInc);
 247     m_filterMax = 1. / (2 * m_signalInc);
 248     m_filterInc = (m_filterMax - m_filterMin) / m_nFilterPoints;
 249     m_vecFilter = new double [m_nFilterPoints];
 250     int halfFilter = m_nFilterPoints / 2;
 251     for (int i = 0; i <= halfFilter; i++)
 252         m_vecFilter[i] = static_cast<double>(i) / halfFilter/ (2. * m_signalInc);
 253     for (int i = 1; i <= halfFilter; i++)
 254         m_vecFilter[m_nFilterPoints - i] = static_cast<double>(i) / halfFilter / (2. * m_signalInc);
 255   }
 256
 257   // precalculate sin and cosine tables for fourier transform
 258   if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
 259       int nFourier = max(m_nFilterPoints,m_nOutputPoints) * max(m_nFilterPoints, m_nOutputPoints) + 1;
 260       double angleIncrement = (2. * PI) / m_nFilterPoints;
 261       m_vecFourierCosTable = new double[ nFourier ];
 262       m_vecFourierSinTable = new double[ nFourier ];
 263       double angle = 0;
 264       for (int i = 0; i < nFourier; i++) {
 265           m_vecFourierCosTable[i] = cos (angle);
 266           m_vecFourierSinTable[i] = sin (angle);
 267           angle += angleIncrement;
 268       }
 269   }
 270
 271 #if HAVE_FFTW
 272   if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
 273       for (int i = 0; i < m_nFilterPoints; i++)  //fftw uses unnormalized fft
 274         m_vecFilter[i] /= m_nFilterPoints;
 275   }
 276
 277   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 278       m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
 279       m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
 280       m_vecRealFftInput = new fftw_real [ m_nFilterPoints ];
 281       m_vecRealFftSignal = new fftw_real [ m_nOutputPoints ];
 282       for (int i = 0; i < m_nFilterPoints; i++)
 283           m_vecRealFftInput[i] = 0;
 284   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 285        m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
 286       m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
 287       m_vecComplexFftInput = new fftw_complex [ m_nFilterPoints ];
 288       m_vecComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
 289       for (int i = 0; i < m_nFilterPoints; i++)
 290           m_vecComplexFftInput[i].re = m_vecComplexFftInput[i].im = 0;
 291       for (int i = 0; i < m_nOutputPoints; i++)
 292           m_vecComplexFftSignal[i].re = m_vecComplexFftSignal[i].im = 0;
 293   }
 294 #endif
 295
 296  if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
 297     m_nFilterPoints = 2 * m_nSignalPoints - 1;
 298     m_filterMin = -m_signalInc * (m_nSignalPoints - 1);
 299     m_filterMax = m_signalInc * (m_nSignalPoints - 1);
 300     m_filterInc = (m_filterMax - m_filterMin) / (m_nFilterPoints - 1);
 301     m_vecFilter = new double[ m_nFilterPoints ];
 302
 303     if (m_idFilter == FILTER_SHEPP) {
 304       double a = 2 * m_bw;
 305       double c = - 4. / (a * a);
 306       int center = (m_nFilterPoints - 1) / 2;
 307       int sidelen = center;
 308       m_vecFilter[center] = 4. / (a * a);
 309
 310       for (int i = 1; i <= sidelen; i++ )
 311         m_vecFilter [center + i] = m_vecFilter [center - i] = c / (4 * (i * i) - 1);
 312     } else if (m_idDomain == DOMAIN_FREQUENCY) {
 313       double x;
 314       int i;
 315       for (x = m_filterMin, i = 0; i < m_nFilterPoints; x += m_filterInc, i++)
 316         m_vecFilter[i] = frequencyResponse (x, m_filterParam);
 317     } else if (m_idDomain == DOMAIN_SPATIAL) {
 318       double x;
 319       int i;
 320       for (x = m_filterMin, i = 0; i < m_nFilterPoints; x += m_filterInc, i++)
 321         if (haveAnalyticSpatial(m_idFilter))
 322           m_vecFilter[i] = spatialResponseAnalytic (x, m_filterParam);
 323         else
 324           m_vecFilter[i] = spatialResponseCalc (x, m_filterParam);
 325     } else {
 326       m_failMessage = "Illegal domain name ";
 327       m_failMessage += m_idDomain;
 328       m_fail = true;
 329     }
 330   }
 331 }
 332
 333 SignalFilter::~SignalFilter (void)
 334 {
 335     delete [] m_vecFilter;
 336     delete [] m_vecFourierSinTable;
 337     delete [] m_vecFourierCosTable;
 338
 339 #if HAVE_FFTW
 340     if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 341         fftw_destroy_plan(m_complexPlanForward);
 342         fftw_destroy_plan(m_complexPlanBackward);
 343         delete [] m_vecComplexFftInput;
 344         delete [] m_vecComplexFftSignal;
 345     }
 346     if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 347         rfftw_destroy_plan(m_realPlanForward);
 348         rfftw_destroy_plan(m_realPlanBackward);
 349         delete [] m_vecRealFftInput;
 350         delete [] m_vecRealFftSignal;
 351     }
 352 #endif
 353 }
 354
 355
 356 int
 357 SignalFilter::convertFilterNameToID (const char *filterName)
 358 {
 359   int filterID = FILTER_INVALID;
 360
 361   for (int i = 0; i < s_iFilterCount; i++)
 362     if (strcasecmp (filterName, s_aszFilterName[i]) == 0) {
 363       filterID = i;
 364       break;
 365     }
 366
 367   return (filterID);
 368 }
 369
 370 const char *
 371 SignalFilter::convertFilterIDToName (const int filterID)
 372 {
 373   static const char *name = "";
 374
 375   if (filterID >= 0 && filterID < s_iFilterCount)
 376       return (s_aszFilterName [filterID]);
 377
 378   return (name);
 379 }
 380
 381 const char *
 382 SignalFilter::convertFilterIDToTitle (const int filterID)
 383 {
 384   static const char *title = "";
 385
 386   if (filterID >= 0 && filterID < s_iFilterCount)
 387       return (s_aszFilterTitle [filterID]);
 388
 389   return (title);
 390 }
 391
 392 int
 393 SignalFilter::convertFilterMethodNameToID (const char* const filterMethodName)
 394 {
 395   int fmID = FILTER_METHOD_INVALID;
 396
 397   for (int i = 0; i < s_iFilterMethodCount; i++)
 398    if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
 399       fmID = i;
 400       break;
 401     }
 402
 403   return (fmID);
 404 }
 405
 406 const char *
 407 SignalFilter::convertFilterMethodIDToName (const int fmID)
 408 {
 409   static const char *name = "";
 410
 411   if (fmID >= 0 && fmID < s_iFilterMethodCount)
 412       return (s_aszFilterMethodName [fmID]);
 413
 414   return (name);
 415 }
 416
 417 const char *
 418 SignalFilter::convertFilterMethodIDToTitle (const int fmID)
 419 {
 420   static const char *title = "";
 421
 422   if (fmID >= 0 && fmID < s_iFilterMethodCount)
 423       return (s_aszFilterTitle [fmID]);
 424
 425   return (title);
 426 }
 427
 428 int
 429 SignalFilter::convertDomainNameToID (const char* const domainName)
 430 {
 431   int dID = DOMAIN_INVALID;
 432
 433   for (int i = 0; i < s_iDomainCount; i++)
 434    if (strcasecmp (domainName, s_aszDomainName[i]) == 0) {
 435       dID = i;
 436       break;
 437     }
 438
 439   return (dID);
 440 }
 441
 442 const char *
 443 SignalFilter::convertDomainIDToName (const int domainID)
 444 {
 445   static const char *name = "";
 446
 447   if (domainID >= 0 && domainID < s_iDomainCount)
 448       return (s_aszDomainName [domainID]);
 449
 450   return (name);
 451 }
 452
 453 const char *
 454 SignalFilter::convertDomainIDToTitle (const int domainID)
 455 {
 456   static const char *title = "";
 457
 458   if (domainID >= 0 && domainID < s_iDomainCount)
 459       return (s_aszDomainTitle [domainID]);
 460
 461   return (title);
 462 }
 463
 464 void
 465 SignalFilter::filterSignal (const float input[], double output[]) const
 466 {
 467   if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
 468     for (int i = 0; i < m_nSignalPoints; i++)
 469       output[i] = convolve (input, m_signalInc, i, m_nSignalPoints);
 470   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
 471     double inputSignal[m_nFilterPoints];
 472     for (int i = 0; i < m_nSignalPoints; i++)
 473       inputSignal[i] = input[i];
 474     for (int i = m_nSignalPoints; i < m_nFilterPoints; i++)
 475       inputSignal[i] = 0;  // zeropad
 476     complex<double> fftSignal[m_nFilterPoints];
 477     finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);
 478     for (int i = 0; i < m_nFilterPoints; i++)
 479       fftSignal[i] *= m_vecFilter[i];
 480     double inverseFourier[m_nFilterPoints];
 481     finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1);
 482     for (int i = 0; i < m_nSignalPoints; i++)
 483       output[i] = inverseFourier[i];
 484   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
 485     double inputSignal[m_nFilterPoints];
 486     for (int i = 0; i < m_nSignalPoints; i++)
 487       inputSignal[i] = input[i];
 488     for (int i = m_nSignalPoints; i < m_nFilterPoints; i++)
 489       inputSignal[i] = 0;  // zeropad
 490     complex<double> fftSignal[m_nFilterPoints];
 491     finiteFourierTransform (inputSignal, fftSignal, -1);
 492     for (int i = 0; i < m_nFilterPoints; i++)
 493       fftSignal[i] *= m_vecFilter[i];
 494     double inverseFourier[m_nFilterPoints];
 495     finiteFourierTransform (fftSignal, inverseFourier, 1);
 496     for (int i = 0; i < m_nSignalPoints; i++)
 497       output[i] = inverseFourier[i];
 498   }
 499 #if HAVE_FFTW
 500   else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 501       for (int i = 0; i < m_nSignalPoints; i++)
 502           m_vecRealFftInput[i] = input[i];
 503
 504       fftw_real fftOutput [ m_nFilterPoints ];
 505       rfftw_one (m_realPlanForward, m_vecRealFftInput, fftOutput);
 506       for (int i = 0; i < m_nFilterPoints; i++)
 507           m_vecRealFftSignal[i] = m_vecFilter[i] * fftOutput[i];
 508       for (int i = m_nFilterPoints; i < m_nOutputPoints; i++)
 509         m_vecRealFftSignal[i] = 0;
 510
 511       fftw_real ifftOutput [ m_nOutputPoints ];
 512       rfftw_one(m_realPlanBackward, m_vecRealFftSignal, ifftOutput);
 513       for (int i = 0; i < m_nSignalPoints * m_preinterpolationFactor; i++)
 514           output[i] = ifftOutput[i];
 515   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 516       for (int i = 0; i < m_nSignalPoints; i++)
 517           m_vecComplexFftInput[i].re = input[i];
 518
 519       fftw_complex fftOutput [ m_nFilterPoints ];
 520       fftw_one(m_complexPlanForward, m_vecComplexFftInput, fftOutput);
 521       for (int i = 0; i < m_nFilterPoints; i++) {
 522           m_vecComplexFftSignal[i].re = m_vecFilter[i] * fftOutput[i].re;
 523           m_vecComplexFftSignal[i].im = m_vecFilter[i] * fftOutput[i].im;
 524       }
 525       fftw_complex ifftOutput [ m_nOutputPoints ];
 526       fftw_one(m_complexPlanBackward, m_vecComplexFftSignal, ifftOutput);
 527       for (int i = 0; i < m_nSignalPoints * m_preinterpolationFactor; i++)
 528           output[i] = ifftOutput[i].re;
 529   }
 530 #endif
 531 }
 532
 533 double
 534 SignalFilter::response (double x)
 535 {
 536   double response = 0;
 537
 538   if (m_idDomain == DOMAIN_SPATIAL)
 539     response = spatialResponse (m_idFilter, m_bw, x, m_filterParam);
 540   else if (m_idDomain == DOMAIN_FREQUENCY)
 541     response = frequencyResponse (m_idFilter, m_bw, x, m_filterParam);
 542
 543   return (response);
 544 }
 545
 546
 547 double
 548 SignalFilter::spatialResponse (int filterID, double bw, double x, double param)
 549 {
 550   if (haveAnalyticSpatial(filterID))
 551     return spatialResponseAnalytic (filterID, bw, x, param);
 552   else
 553     return spatialResponseCalc (filterID, bw, x, param, N_INTEGRAL);
 554 }
 555
 556 /* NAME
 557  *   filter_spatial_response_calc       Calculate filter by discrete inverse fourier
 558  *                                      transform of filters's frequency
 559  *                                      response
 560  *
 561  * SYNOPSIS
 562  *   y = filter_spatial_response_calc (filt_type, x, m_bw, param, n)
 563  *   double y                   Filter's response in spatial domain
 564  *   int filt_type              Type of filter (definitions in ct.h)
 565  *   double x                   Spatial position to evaluate filter
 566  *   double m_bw                        Bandwidth of window
 567  *   double param               General parameter for various filters
 568  *   int n                      Number of points to calculate integrations
 569  */
 570
 571 double
 572 SignalFilter::spatialResponseCalc (double x, double param) const
 573 {
 574   return (spatialResponseCalc (m_idFilter, m_bw, x, param, N_INTEGRAL));
 575 }
 576
 577 double
 578 SignalFilter::spatialResponseCalc (int filterID, double bw, double x, double param, int n)
 579 {
 580   double zmin, zmax;
 581
 582   if (filterID == FILTER_TRIANGLE) {
 583     zmin = 0;
 584     zmax = bw;
 585   } else {
 586     zmin = 0;
 587     zmax = bw / 2;
 588   }
 589   double zinc = (zmax - zmin) / (n - 1);
 590
 591   double z = zmin;
 592   double q [n];
 593   for (int i = 0; i < n; i++, z += zinc)
 594     q[i] = frequencyResponse (filterID, bw, z, param) * cos (TWOPI * z * x);
 595
 596   double y = 2 * integrateSimpson (zmin, zmax, q, n);
 597
 598   return (y);
 599 }
 600
 601
 602 /* NAME
 603  *    filter_frequency_response                 Return filter frequency response
 604  *
 605  * SYNOPSIS
 606  *    h = filter_frequency_response (filt_type, u, m_bw, param)
 607  *    double h                  Filters frequency response at u
 608  *    int filt_type             Type of filter
 609  *    double u                  Frequency to evaluate filter at
 610  *    double m_bw                       Bandwidth of filter
 611  *    double param              General input parameter for various filters
 612  */
 613
 614 double
 615 SignalFilter::frequencyResponse (double u, double param) const
 616 {
 617   return frequencyResponse (m_idFilter, m_bw, u, param);
 618 }
 619
 620
 621 double
 622 SignalFilter::frequencyResponse (int filterID, double bw, double u, double param)
 623 {
 624   double q;
 625   double au = fabs (u);
 626
 627   switch (filterID) {
 628   case FILTER_BANDLIMIT:
 629     if (au >= bw / 2)
 630       q = 0.;
 631     else
 632       q = 1;
 633     break;
 634   case FILTER_ABS_BANDLIMIT:
 635     if (au >= bw / 2)
 636       q = 0.;
 637     else
 638       q = au;
 639     break;
 640   case FILTER_TRIANGLE:
 641     if (au >= bw)
 642       q = 0;
 643     else
 644       q = 1 - au / bw;
 645     break;
 646   case FILTER_COSINE:
 647     if (au >= bw / 2)
 648       q = 0;
 649     else
 650       q = cos(PI * u / bw);
 651     break;
 652   case FILTER_ABS_COSINE:
 653     if (au >= bw / 2)
 654       q = 0;
 655     else
 656       q = au * cos(PI * u / bw);
 657     break;
 658   case FILTER_SINC:
 659     q = bw * sinc (PI * bw * u, 1.);
 660     break;
 661   case FILTER_ABS_SINC:
 662     q = au * bw * sinc (PI * bw * u, 1.);
 663     break;
 664   case FILTER_G_HAMMING:
 665     if (au >= bw / 2)
 666       q = 0;
 667     else
 668       q = param + (1 - param) * cos (TWOPI * u / bw);
 669     break;
 670   case FILTER_ABS_G_HAMMING:
 671     if (au >= bw / 2)
 672       q = 0;
 673     else
 674       q = au * (param + (1 - param) * cos(TWOPI * u / bw));
 675     break;
 676   default:
 677     q = 0;
 678     sys_error (ERR_WARNING, "Frequency response for filter %d not implemented [filter_frequency_response]", filterID);
 679     break;
 680   }
 681   return (q);
 682 }
 683
 684
 685
 686 /* NAME
 687  *   filter_spatial_response_analytic                   Calculate filter by analytic inverse fourier
 688  *                              transform of filters's frequency
 689  *                              response
 690  *
 691  * SYNOPSIS
 692  *   y = filter_spatial_response_analytic (filt_type, x, m_bw, param)
 693  *   double y                   Filter's response in spatial domain
 694  *   int filt_type              Type of filter (definitions in ct.h)
 695  *   double x                   Spatial position to evaluate filter
 696  *   double m_bw                        Bandwidth of window
 697  *   double param               General parameter for various filters
 698  */
 699
 700 double
 701 SignalFilter::spatialResponseAnalytic (double x, double param) const
 702 {
 703   return spatialResponseAnalytic (m_idFilter, m_bw, x, param);
 704 }
 705
 706 const bool
 707 SignalFilter::haveAnalyticSpatial (int filterID)
 708 {
 709   bool haveAnalytic = false;
 710
 711   switch (filterID) {
 712   case FILTER_BANDLIMIT:
 713   case FILTER_TRIANGLE:
 714   case FILTER_COSINE:
 715   case FILTER_G_HAMMING:
 716   case FILTER_ABS_BANDLIMIT:
 717   case FILTER_ABS_COSINE:
 718   case FILTER_ABS_G_HAMMING:
 719   case FILTER_SHEPP:
 720   case FILTER_SINC:
 721     haveAnalytic = true;
 722     break;
 723   default:
 724     break;
 725   }
 726
 727   return (haveAnalytic);
 728 }
 729
 730 double
 731 SignalFilter::spatialResponseAnalytic (int filterID, double bw, double x, double param)
 732 {
 733   double q, temp;
 734   double u = TWOPI * x;
 735   double w = bw / 2;
 736   double b = PI / bw;
 737   double b2 = TWOPI / bw;
 738
 739   switch (filterID) {
 740   case FILTER_BANDLIMIT:
 741     q = bw * sinc(u * w, 1.0);
 742     break;
 743   case FILTER_TRIANGLE:
 744     temp = sinc (u * w, 1.0);
 745     q = bw * temp * temp;
 746     break;
 747   case FILTER_COSINE:
 748     q = sinc(b-u,w) + sinc(b+u,w);
 749     break;
 750   case FILTER_G_HAMMING:
 751     q = 2 * param * sin(u*w)/u + (1-param) * (sinc(b2-u, w) + sinc(b2+u, w));
 752     break;
 753   case FILTER_ABS_BANDLIMIT:
 754     q = 2 * integral_abscos (u, w);
 755     break;
 756   case FILTER_ABS_COSINE:
 757     q = integral_abscos(b-u,w) + integral_abscos(b+u,w);
 758     break;
 759   case FILTER_ABS_G_HAMMING:
 760     q = 2 * param * integral_abscos(u,w) +
 761       (1-param)*(integral_abscos(u-b2,w)+integral_abscos(u+b2,w));
 762     break;
 763   case FILTER_SHEPP:
 764     if (fabs (u) < 1E-7)
 765       q = 4. / (PI * bw * bw);
 766     else
 767       q = fabs ((2 / bw) * sin (u * w)) * sinc (u * w, 1.) * sinc (u * w, 1.);
 768     break;
 769   case FILTER_SINC:
 770     if (fabs (x) < bw / 2)
 771       q = 1.;
 772     else
 773       q = 0.;
 774     break;
 775   case FILTER_ABS_SINC:
 776   default:
 777     sys_error (ERR_WARNING, "Analytic filter type %d not implemented [filter_spatial_response_analytic]", filterID);
 778     q = 0;
 779     break;
 780   }
 781
 782   return (q);
 783 }
 784
 785
 786 /* NAME
 787  *   sinc                       Return sin(x)/x function
 788  *
 789  * SYNOPSIS
 790  *   v = sinc (x, mult)
 791  *   double v                   sinc value
 792  *   double x, mult
 793  *
 794  * DESCRIPTION
 795  *   v = sin(x * mult) / x;
 796  */
 797
 798
 799 /* NAME
 800  *   integral_abscos                    Returns integral of u*cos(u)
 801  *
 802  * SYNOPSIS
 803  *   q = integral_abscos (u, w)
 804  *   double q                   Integral value
 805  *   double u                   Integration variable
 806  *   double w                   Upper integration boundary
 807  *
 808  * DESCRIPTION
 809  *   Returns the value of integral of u*cos(u)*dV for V = 0 to w
 810  */
 811
 812 double
 813 SignalFilter::integral_abscos (double u, double w)
 814 {
 815   return (fabs (u) > F_EPSILON
 816      ? (cos(u * w) - 1) / (u * u) + w / u * sin (u * w)
 817      : (w * w / 2));
 818 }
 819
 820
 821 /* NAME
 822  *    convolve                  Discrete convolution of two functions
 823  *
 824  * SYNOPSIS
 825  *    r = convolve (f1, f2, dx, n, np, func_type)
 826  *    double r                  Convolved result
 827  *    double f1[], f2[]         Functions to be convolved
 828  *    double dx                 Difference between successive x values
 829  *    int n                     Array index to center convolution about
 830  *    int np                    Number of points in f1 array
 831  *    int func_type             EVEN or ODD or EVEN_AND_ODD function f2
 832  *
 833  * NOTES
 834  *    f1 is the projection data, its indices range from 0 to np - 1.
 835  *    The index for f2, the filter, ranges from -(np-1) to (np-1).
 836  *    There are 3 ways to handle the negative vertices of f2:
 837  *      1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
 838  *         All filters used in reconstruction are even.
 839  *      2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
 840  *      3. If f2 is both ODD AND EVEN, then we must store the value of f2
 841  *         for negative indices.  Since f2 must range from -(np-1) to (np-1),
 842  *         if we add (np - 1) to f2's array index, then f2's index will
 843  *         range from 0 to 2 * (np - 1), and the origin, x = 0, will be
 844  *         stored at f2[np-1].
 845  */
 846
 847 double
 848 SignalFilter::convolve (const double func[], const double dx, const int n, const int np) const
 849 {
 850   double sum = 0.0;
 851
 852 #if UNOPTIMIZED_CONVOLUTION
 853   for (int i = 0; i < np; i++)
 854     sum += func[i] * m_vecFilter[n - i + (np - 1)];
 855 #else
 856   double* f2 = m_vecFilter + n + (np - 1);
 857   for (int i = 0; i < np; i++)
 858     sum += *func++ * *f2--;
 859 #endif
 860
 861   return (sum * dx);
 862 }
 863
 864
 865 double
 866 SignalFilter::convolve (const float func[], const double dx, const int n, const int np) const
 867 {
 868   double sum = 0.0;
 869
 870 #if UNOPTIMIZED_CONVOLUTION
 871 for (int i = 0; i < np; i++)
 872   sum += func[i] * m_vecFilter[n - i + (np - 1)];
 873 #else
 874 double* f2 = m_vecFilter + n + (np - 1);
 875 for (int i = 0; i < np; i++)
 876   sum += *func++ * *f2--;
 877 #endif
 878
 879   return (sum * dx);
 880 }
 881
 882
 883 void
 884 SignalFilter::finiteFourierTransform (const double input[], complex<double> output[], const int n, int direction)
 885 {
 886   if (direction < 0)
 887     direction = -1;
 888   else
 889     direction = 1;
 890
 891   double angleIncrement = direction * 2 * PI / n;
 892   for (int i = 0; i < n; i++) {
 893     double sumReal = 0;
 894     double sumImag = 0;
 895     for (int j = 0; j < n; j++) {
 896       double angle = i * j * angleIncrement;
 897       sumReal += input[j] * cos(angle);
 898       sumImag += input[j] * sin(angle);
 899     }
 900     if (direction < 0) {
 901       sumReal /= n;
 902       sumImag /= n;
 903     }
 904     output[i] = complex<double> (sumReal, sumImag);
 905   }
 906 }
 907
 908
 909 void
 910 SignalFilter::finiteFourierTransform (const complex<double> input[], complex<double> output[], const int n, int direction)
 911 {
 912   if (direction < 0)
 913     direction = -1;
 914   else
 915     direction = 1;
 916
 917   double angleIncrement = direction * 2 * PI / n;
 918   for (int i = 0; i < n; i++) {
 919     complex<double> sum (0,0);
 920     for (int j = 0; j < n; j++) {
 921       double angle = i * j * angleIncrement;
 922       complex<double> exponentTerm (cos(angle), sin(angle));
 923       sum += input[j] * exponentTerm;
 924     }
 925     if (direction < 0) {
 926       sum /= n;
 927     }
 928     output[i] = sum;
 929   }
 930 }
 931
 932 void
 933 SignalFilter::finiteFourierTransform (const complex<double> input[], double output[], const int n, int direction)
 934 {
 935   if (direction < 0)
 936     direction = -1;
 937   else
 938     direction = 1;
 939
 940   double angleIncrement = direction * 2 * PI / n;
 941   for (int i = 0; i < n; i++) {
 942       double sumReal = 0;
 943     for (int j = 0; j < n; j++) {
 944       double angle = i * j * angleIncrement;
 945       sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
 946     }
 947     if (direction < 0) {
 948       sumReal /= n;
 949     }
 950     output[i] = sumReal;
 951   }
 952 }
 953
 954 void
 955 SignalFilter::finiteFourierTransform (const double input[], complex<double> output[], int direction) const
 956 {
 957   if (direction < 0)
 958     direction = -1;
 959   else
 960     direction = 1;
 961
 962   for (int i = 0; i < m_nFilterPoints; i++) {
 963     double sumReal = 0, sumImag = 0;
 964     for (int j = 0; j < m_nFilterPoints; j++) {
 965       int tableIndex = i * j;
 966       if (direction > 0) {
 967         sumReal += input[j] * m_vecFourierCosTable[tableIndex];
 968         sumImag += input[j] * m_vecFourierSinTable[tableIndex];
 969       } else {
 970         sumReal += input[j] * m_vecFourierCosTable[tableIndex];
 971         sumImag -= input[j] * m_vecFourierSinTable[tableIndex];
 972       }
 973     }
 974     if (direction < 0) {
 975       sumReal /= m_nFilterPoints;
 976       sumImag /= m_nFilterPoints;
 977     }
 978     output[i] = complex<double> (sumReal, sumImag);
 979   }
 980 }
 981
 982 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
 983 void
 984 SignalFilter::finiteFourierTransform (const complex<double> input[], complex<double> output[], int direction) const
 985 {
 986   if (direction < 0)
 987     direction = -1;
 988   else
 989     direction = 1;
 990
 991   for (int i = 0; i < m_nFilterPoints; i++) {
 992     double sumReal = 0, sumImag = 0;
 993     for (int j = 0; j < m_nFilterPoints; j++) {
 994       int tableIndex = i * j;
 995       if (direction > 0) {
 996         sumReal += input[j].real() * m_vecFourierCosTable[tableIndex]
 997           - input[j].imag() * m_vecFourierSinTable[tableIndex];
 998         sumImag += input[j].real() * m_vecFourierSinTable[tableIndex]
 999           + input[j].imag() * m_vecFourierCosTable[tableIndex];
1000       } else {
1001         sumReal += input[j].real() * m_vecFourierCosTable[tableIndex]
1002           - input[j].imag() * -m_vecFourierSinTable[tableIndex];
1003         sumImag += input[j].real() * -m_vecFourierSinTable[tableIndex]
1004           + input[j].imag() * m_vecFourierCosTable[tableIndex];
1005       }
1006     }
1007     if (direction < 0) {
1008       sumReal /= m_nFilterPoints;
1009       sumImag /= m_nFilterPoints;
1010     }
1011     output[i] = complex<double> (sumReal, sumImag);
1012   }
1013 }
1014
1015 void
1016 SignalFilter::finiteFourierTransform (const complex<double> input[], double output[], int direction) const
1017 {
1018   if (direction < 0)
1019     direction = -1;
1020   else
1021     direction = 1;
1022
1023   for (int i = 0; i < m_nFilterPoints; i++) {
1024       double sumReal = 0;
1025     for (int j = 0; j < m_nFilterPoints; j++) {
1026       int tableIndex = i * j;
1027       if (direction > 0) {
1028         sumReal += input[j].real() * m_vecFourierCosTable[tableIndex]
1029           - input[j].imag() * m_vecFourierSinTable[tableIndex];
1030       } else {
1031         sumReal += input[j].real() * m_vecFourierCosTable[tableIndex]
1032           - input[j].imag() * -m_vecFourierSinTable[tableIndex];
1033       }
1034     }
1035     if (direction < 0) {
1036       sumReal /= m_nFilterPoints;
1037     }
1038     output[i] = sumReal;
1039   }
1040 }
1041
1042