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