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