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.23 2000/07/31 14:48:35 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 #if LIMIT_BANDWIDTH_TRIAL
 326         if (i < m_nFilterPoints / 4 || i > (m_nFilterPoints * 3) / 4)
 327           m_vecFilter[i] = 0;
 328 #endif
 329       }
 330     } else {
 331       m_failMessage = "Illegal domain name ";
 332       m_failMessage += m_idDomain;
 333       m_fail = true;
 334     }
 335   }
 336 }
 337
 338 SignalFilter::~SignalFilter (void)
 339 {
 340     delete [] m_vecFilter;
 341     delete [] m_vecFourierSinTable;
 342     delete [] m_vecFourierCosTable;
 343
 344 #if HAVE_FFTW
 345     if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 346         fftw_destroy_plan(m_complexPlanForward);
 347         fftw_destroy_plan(m_complexPlanBackward);
 348         delete [] m_vecComplexFftInput;
 349         delete [] m_vecComplexFftSignal;
 350     }
 351     if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 352         rfftw_destroy_plan(m_realPlanForward);
 353         rfftw_destroy_plan(m_realPlanBackward);
 354         delete [] m_vecRealFftInput;
 355         delete [] m_vecRealFftSignal;
 356     }
 357 #endif
 358 }
 359
 360
 361 int
 362 SignalFilter::convertFilterNameToID (const char *filterName)
 363 {
 364   int filterID = FILTER_INVALID;
 365
 366   for (int i = 0; i < s_iFilterCount; i++)
 367     if (strcasecmp (filterName, s_aszFilterName[i]) == 0) {
 368       filterID = i;
 369       break;
 370     }
 371
 372   return (filterID);
 373 }
 374
 375 const char *
 376 SignalFilter::convertFilterIDToName (const int filterID)
 377 {
 378   static const char *name = "";
 379
 380   if (filterID >= 0 && filterID < s_iFilterCount)
 381       return (s_aszFilterName [filterID]);
 382
 383   return (name);
 384 }
 385
 386 const char *
 387 SignalFilter::convertFilterIDToTitle (const int filterID)
 388 {
 389   static const char *title = "";
 390
 391   if (filterID >= 0 && filterID < s_iFilterCount)
 392       return (s_aszFilterTitle [filterID]);
 393
 394   return (title);
 395 }
 396
 397 int
 398 SignalFilter::convertFilterMethodNameToID (const char* const filterMethodName)
 399 {
 400   int fmID = FILTER_METHOD_INVALID;
 401
 402   for (int i = 0; i < s_iFilterMethodCount; i++)
 403    if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
 404       fmID = i;
 405       break;
 406     }
 407
 408   return (fmID);
 409 }
 410
 411 const char *
 412 SignalFilter::convertFilterMethodIDToName (const int fmID)
 413 {
 414   static const char *name = "";
 415
 416   if (fmID >= 0 && fmID < s_iFilterMethodCount)
 417       return (s_aszFilterMethodName [fmID]);
 418
 419   return (name);
 420 }
 421
 422 const char *
 423 SignalFilter::convertFilterMethodIDToTitle (const int fmID)
 424 {
 425   static const char *title = "";
 426
 427   if (fmID >= 0 && fmID < s_iFilterMethodCount)
 428       return (s_aszFilterTitle [fmID]);
 429
 430   return (title);
 431 }
 432
 433 int
 434 SignalFilter::convertDomainNameToID (const char* const domainName)
 435 {
 436   int dID = DOMAIN_INVALID;
 437
 438   for (int i = 0; i < s_iDomainCount; i++)
 439    if (strcasecmp (domainName, s_aszDomainName[i]) == 0) {
 440       dID = i;
 441       break;
 442     }
 443
 444   return (dID);
 445 }
 446
 447 const char *
 448 SignalFilter::convertDomainIDToName (const int domainID)
 449 {
 450   static const char *name = "";
 451
 452   if (domainID >= 0 && domainID < s_iDomainCount)
 453       return (s_aszDomainName [domainID]);
 454
 455   return (name);
 456 }
 457
 458 const char *
 459 SignalFilter::convertDomainIDToTitle (const int domainID)
 460 {
 461   static const char *title = "";
 462
 463   if (domainID >= 0 && domainID < s_iDomainCount)
 464       return (s_aszDomainTitle [domainID]);
 465
 466   return (title);
 467 }
 468
 469 void
 470 SignalFilter::filterSignal (const float input[], double output[]) const
 471 {
 472   if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
 473     for (int i = 0; i < m_nSignalPoints; i++)
 474       output[i] = convolve (input, m_signalInc, i, m_nSignalPoints);
 475   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
 476     double inputSignal[m_nFilterPoints];
 477     for (int i = 0; i < m_nSignalPoints; i++)
 478       inputSignal[i] = input[i];
 479     for (int i = m_nSignalPoints; i < m_nFilterPoints; i++)
 480       inputSignal[i] = 0;  // zeropad
 481     complex<double> fftSignal[m_nFilterPoints];
 482     finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);
 483     for (int i = 0; i < m_nFilterPoints; i++)
 484       fftSignal[i] *= m_vecFilter[i];
 485     double inverseFourier[m_nFilterPoints];
 486     finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1);
 487     for (int i = 0; i < m_nSignalPoints; i++)
 488       output[i] = inverseFourier[i];
 489   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
 490     double inputSignal[m_nFilterPoints];
 491     for (int i = 0; i < m_nSignalPoints; i++)
 492       inputSignal[i] = input[i];
 493     for (int i = m_nSignalPoints; i < m_nFilterPoints; i++)
 494       inputSignal[i] = 0;  // zeropad
 495     complex<double> fftSignal[m_nFilterPoints];
 496     finiteFourierTransform (inputSignal, fftSignal, -1);
 497     for (int i = 0; i < m_nFilterPoints; i++)
 498       fftSignal[i] *= m_vecFilter[i];
 499     double inverseFourier[m_nFilterPoints];
 500     finiteFourierTransform (fftSignal, inverseFourier, 1);
 501     for (int i = 0; i < m_nSignalPoints; i++)
 502       output[i] = inverseFourier[i];
 503   }
 504 #if HAVE_FFTW
 505   else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 506       for (int i = 0; i < m_nSignalPoints; i++)
 507           m_vecRealFftInput[i] = input[i];
 508
 509       fftw_real fftOutput [ m_nFilterPoints ];
 510       rfftw_one (m_realPlanForward, m_vecRealFftInput, fftOutput);
 511       for (int i = 0; i < m_nFilterPoints; i++)
 512           m_vecRealFftSignal[i] = m_vecFilter[i] * fftOutput[i];
 513       for (int i = m_nFilterPoints; i < m_nOutputPoints; i++)
 514         m_vecRealFftSignal[i] = 0;
 515
 516       fftw_real ifftOutput [ m_nOutputPoints ];
 517       rfftw_one(m_realPlanBackward, m_vecRealFftSignal, ifftOutput);
 518       for (int i = 0; i < m_nSignalPoints * m_preinterpolationFactor; i++)
 519           output[i] = ifftOutput[i];
 520   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 521       for (int i = 0; i < m_nSignalPoints; i++)
 522           m_vecComplexFftInput[i].re = input[i];
 523
 524       fftw_complex fftOutput [ m_nFilterPoints ];
 525       fftw_one(m_complexPlanForward, m_vecComplexFftInput, fftOutput);
 526       for (int i = 0; i < m_nFilterPoints; i++) {
 527           m_vecComplexFftSignal[i].re = m_vecFilter[i] * fftOutput[i].re;
 528           m_vecComplexFftSignal[i].im = m_vecFilter[i] * fftOutput[i].im;
 529       }
 530       fftw_complex ifftOutput [ m_nOutputPoints ];
 531       fftw_one(m_complexPlanBackward, m_vecComplexFftSignal, ifftOutput);
 532       for (int i = 0; i < m_nSignalPoints * m_preinterpolationFactor; i++)
 533           output[i] = ifftOutput[i].re;
 534   }
 535 #endif
 536 }
 537
 538 double
 539 SignalFilter::response (double x)
 540 {
 541   double response = 0;
 542
 543   if (m_idDomain == DOMAIN_SPATIAL)
 544     response = spatialResponse (m_idFilter, m_bw, x, m_filterParam);
 545   else if (m_idDomain == DOMAIN_FREQUENCY)
 546     response = frequencyResponse (m_idFilter, m_bw, x, m_filterParam);
 547
 548   return (response);
 549 }
 550
 551
 552 double
 553 SignalFilter::spatialResponse (int filterID, double bw, double x, double param)
 554 {
 555   if (haveAnalyticSpatial(filterID))
 556     return spatialResponseAnalytic (filterID, bw, x, param);
 557   else
 558     return spatialResponseCalc (filterID, bw, x, param, N_INTEGRAL);
 559 }
 560
 561 /* NAME
 562  *   filter_spatial_response_calc       Calculate filter by discrete inverse fourier
 563  *                                      transform of filters's frequency
 564  *                                      response
 565  *
 566  * SYNOPSIS
 567  *   y = filter_spatial_response_calc (filt_type, x, m_bw, param, n)
 568  *   double y                   Filter's response in spatial domain
 569  *   int filt_type              Type of filter (definitions in ct.h)
 570  *   double x                   Spatial position to evaluate filter
 571  *   double m_bw                        Bandwidth of window
 572  *   double param               General parameter for various filters
 573  *   int n                      Number of points to calculate integrations
 574  */
 575
 576 double
 577 SignalFilter::spatialResponseCalc (double x, double param) const
 578 {
 579   return (spatialResponseCalc (m_idFilter, m_bw, x, param, N_INTEGRAL));
 580 }
 581
 582 double
 583 SignalFilter::spatialResponseCalc (int filterID, double bw, double x, double param, int n)
 584 {
 585   double zmin, zmax;
 586
 587   if (filterID == FILTER_TRIANGLE) {
 588     zmin = 0;
 589     zmax = bw;
 590   } else {
 591     zmin = 0;
 592     zmax = bw / 2;
 593   }
 594   double zinc = (zmax - zmin) / (n - 1);
 595
 596   double z = zmin;
 597   double q [n];
 598   for (int i = 0; i < n; i++, z += zinc)
 599     q[i] = frequencyResponse (filterID, bw, z, param) * cos (TWOPI * z * x);
 600
 601   double y = 2 * integrateSimpson (zmin, zmax, q, n);
 602
 603   return (y);
 604 }
 605
 606
 607 /* NAME
 608  *    filter_frequency_response                 Return filter frequency response
 609  *
 610  * SYNOPSIS
 611  *    h = filter_frequency_response (filt_type, u, m_bw, param)
 612  *    double h                  Filters frequency response at u
 613  *    int filt_type             Type of filter
 614  *    double u                  Frequency to evaluate filter at
 615  *    double m_bw                       Bandwidth of filter
 616  *    double param              General input parameter for various filters
 617  */
 618
 619 double
 620 SignalFilter::frequencyResponse (double u, double param) const
 621 {
 622   return frequencyResponse (m_idFilter, m_bw, u, param);
 623 }
 624
 625
 626 double
 627 SignalFilter::frequencyResponse (int filterID, double bw, double u, double param)
 628 {
 629   double q;
 630   double au = fabs (u);
 631
 632   switch (filterID) {
 633   case FILTER_BANDLIMIT:
 634     if (au >= bw / 2)
 635       q = 0.;
 636     else
 637       q = 1;
 638     break;
 639   case FILTER_ABS_BANDLIMIT:
 640     if (au >= bw / 2)
 641       q = 0.;
 642     else
 643       q = au;
 644     break;
 645   case FILTER_TRIANGLE:
 646     if (au >= bw)
 647       q = 0;
 648     else
 649       q = 1 - au / bw;
 650     break;
 651   case FILTER_COSINE:
 652     if (au >= bw / 2)
 653       q = 0;
 654     else
 655       q = cos(PI * u / bw);
 656     break;
 657   case FILTER_ABS_COSINE:
 658     if (au >= bw / 2)
 659       q = 0;
 660     else
 661       q = au * cos(PI * u / bw);
 662     break;
 663   case FILTER_SINC:
 664     q = bw * sinc (PI * bw * u, 1.);
 665     break;
 666   case FILTER_ABS_SINC:
 667     q = au * bw * sinc (PI * bw * u, 1.);
 668     break;
 669   case FILTER_G_HAMMING:
 670     if (au >= bw / 2)
 671       q = 0;
 672     else
 673       q = param + (1 - param) * cos (TWOPI * u / bw);
 674     break;
 675   case FILTER_ABS_G_HAMMING:
 676     if (au >= bw / 2)
 677       q = 0;
 678     else
 679       q = au * (param + (1 - param) * cos(TWOPI * u / bw));
 680     break;
 681   default:
 682     q = 0;
 683     sys_error (ERR_WARNING, "Frequency response for filter %d not implemented [filter_frequency_response]", filterID);
 684     break;
 685   }
 686   return (q);
 687 }
 688
 689
 690
 691 /* NAME
 692  *   filter_spatial_response_analytic                   Calculate filter by analytic inverse fourier
 693  *                              transform of filters's frequency
 694  *                              response
 695  *
 696  * SYNOPSIS
 697  *   y = filter_spatial_response_analytic (filt_type, x, m_bw, param)
 698  *   double y                   Filter's response in spatial domain
 699  *   int filt_type              Type of filter (definitions in ct.h)
 700  *   double x                   Spatial position to evaluate filter
 701  *   double m_bw                        Bandwidth of window
 702  *   double param               General parameter for various filters
 703  */
 704
 705 double
 706 SignalFilter::spatialResponseAnalytic (double x, double param) const
 707 {
 708   return spatialResponseAnalytic (m_idFilter, m_bw, x, param);
 709 }
 710
 711 const bool
 712 SignalFilter::haveAnalyticSpatial (int filterID)
 713 {
 714   bool haveAnalytic = false;
 715
 716   switch (filterID) {
 717   case FILTER_BANDLIMIT:
 718   case FILTER_TRIANGLE:
 719   case FILTER_COSINE:
 720   case FILTER_G_HAMMING:
 721   case FILTER_ABS_BANDLIMIT:
 722   case FILTER_ABS_COSINE:
 723   case FILTER_ABS_G_HAMMING:
 724   case FILTER_SHEPP:
 725   case FILTER_SINC:
 726     haveAnalytic = true;
 727     break;
 728   default:
 729     break;
 730   }
 731
 732   return (haveAnalytic);
 733 }
 734
 735 double
 736 SignalFilter::spatialResponseAnalytic (int filterID, double bw, double x, double param)
 737 {
 738   double q, temp;
 739   double u = TWOPI * x;
 740   double w = bw / 2;
 741   double b = PI / bw;
 742   double b2 = TWOPI / bw;
 743
 744   switch (filterID) {
 745   case FILTER_BANDLIMIT:
 746     q = bw * sinc(u * w, 1.0);
 747     break;
 748   case FILTER_TRIANGLE:
 749     temp = sinc (u * w, 1.0);
 750     q = bw * temp * temp;
 751     break;
 752   case FILTER_COSINE:
 753     q = sinc(b-u,w) + sinc(b+u,w);
 754     break;
 755   case FILTER_G_HAMMING:
 756     q = 2 * param * sin(u*w)/u + (1-param) * (sinc(b2-u, w) + sinc(b2+u, w));
 757     break;
 758   case FILTER_ABS_BANDLIMIT:
 759     q = 2 * integral_abscos (u, w);
 760     break;
 761   case FILTER_ABS_COSINE:
 762     q = integral_abscos(b-u,w) + integral_abscos(b+u,w);
 763     break;
 764   case FILTER_ABS_G_HAMMING:
 765     q = 2 * param * integral_abscos(u,w) +
 766       (1-param)*(integral_abscos(u-b2,w)+integral_abscos(u+b2,w));
 767     break;
 768   case FILTER_SHEPP:
 769     if (fabs (u) < 1E-7)
 770       q = 4. / (PI * bw * bw);
 771     else
 772       q = fabs ((2 / bw) * sin (u * w)) * sinc (u * w, 1.) * sinc (u * w, 1.);
 773     break;
 774   case FILTER_SINC:
 775     if (fabs (x) < bw / 2)
 776       q = 1.;
 777     else
 778       q = 0.;
 779     break;
 780   case FILTER_ABS_SINC:
 781   default:
 782     sys_error (ERR_WARNING, "Analytic filter type %d not implemented [filter_spatial_response_analytic]", filterID);
 783     q = 0;
 784     break;
 785   }
 786
 787   return (q);
 788 }
 789
 790
 791 /* NAME
 792  *   sinc                       Return sin(x)/x function
 793  *
 794  * SYNOPSIS
 795  *   v = sinc (x, mult)
 796  *   double v                   sinc value
 797  *   double x, mult
 798  *
 799  * DESCRIPTION
 800  *   v = sin(x * mult) / x;
 801  */
 802
 803
 804 /* NAME
 805  *   integral_abscos                    Returns integral of u*cos(u)
 806  *
 807  * SYNOPSIS
 808  *   q = integral_abscos (u, w)
 809  *   double q                   Integral value
 810  *   double u                   Integration variable
 811  *   double w                   Upper integration boundary
 812  *
 813  * DESCRIPTION
 814  *   Returns the value of integral of u*cos(u)*dV for V = 0 to w
 815  */
 816
 817 double
 818 SignalFilter::integral_abscos (double u, double w)
 819 {
 820   return (fabs (u) > F_EPSILON
 821      ? (cos(u * w) - 1) / (u * u) + w / u * sin (u * w)
 822      : (w * w / 2));
 823 }
 824
 825
 826 /* NAME
 827  *    convolve                  Discrete convolution of two functions
 828  *
 829  * SYNOPSIS
 830  *    r = convolve (f1, f2, dx, n, np, func_type)
 831  *    double r                  Convolved result
 832  *    double f1[], f2[]         Functions to be convolved
 833  *    double dx                 Difference between successive x values
 834  *    int n                     Array index to center convolution about
 835  *    int np                    Number of points in f1 array
 836  *    int func_type             EVEN or ODD or EVEN_AND_ODD function f2
 837  *
 838  * NOTES
 839  *    f1 is the projection data, its indices range from 0 to np - 1.
 840  *    The index for f2, the filter, ranges from -(np-1) to (np-1).
 841  *    There are 3 ways to handle the negative vertices of f2:
 842  *      1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
 843  *         All filters used in reconstruction are even.
 844  *      2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
 845  *      3. If f2 is both ODD AND EVEN, then we must store the value of f2
 846  *         for negative indices.  Since f2 must range from -(np-1) to (np-1),
 847  *         if we add (np - 1) to f2's array index, then f2's index will
 848  *         range from 0 to 2 * (np - 1), and the origin, x = 0, will be
 849  *         stored at f2[np-1].
 850  */
 851
 852 double
 853 SignalFilter::convolve (const double func[], const double dx, const int n, const int np) const
 854 {
 855   double sum = 0.0;
 856
 857 #if UNOPTIMIZED_CONVOLUTION
 858   for (int i = 0; i < np; i++)
 859     sum += func[i] * m_vecFilter[n - i + (np - 1)];
 860 #else
 861   double* f2 = m_vecFilter + n + (np - 1);
 862   for (int i = 0; i < np; i++)
 863     sum += *func++ * *f2--;
 864 #endif
 865
 866   return (sum * dx);
 867 }
 868
 869
 870 double
 871 SignalFilter::convolve (const float func[], const double dx, const int n, const int np) const
 872 {
 873   double sum = 0.0;
 874
 875 #if UNOPTIMIZED_CONVOLUTION
 876 for (int i = 0; i < np; i++)
 877   sum += func[i] * m_vecFilter[n - i + (np - 1)];
 878 #else
 879 double* f2 = m_vecFilter + n + (np - 1);
 880 for (int i = 0; i < np; i++)
 881   sum += *func++ * *f2--;
 882 #endif
 883
 884   return (sum * dx);
 885 }
 886
 887
 888 void
 889 SignalFilter::finiteFourierTransform (const double input[], complex<double> output[], const int n, int direction)
 890 {
 891   if (direction < 0)
 892     direction = -1;
 893   else
 894     direction = 1;
 895
 896   double angleIncrement = direction * 2 * PI / n;
 897   for (int i = 0; i < n; i++) {
 898     double sumReal = 0;
 899     double sumImag = 0;
 900     for (int j = 0; j < n; j++) {
 901       double angle = i * j * angleIncrement;
 902       sumReal += input[j] * cos(angle);
 903       sumImag += input[j] * sin(angle);
 904     }
 905     if (direction < 0) {
 906       sumReal /= n;
 907       sumImag /= n;
 908     }
 909     output[i] = complex<double> (sumReal, sumImag);
 910   }
 911 }
 912
 913
 914 void
 915 SignalFilter::finiteFourierTransform (const complex<double> input[], complex<double> output[], const int n, int direction)
 916 {
 917   if (direction < 0)
 918     direction = -1;
 919   else
 920     direction = 1;
 921
 922   double angleIncrement = direction * 2 * PI / n;
 923   for (int i = 0; i < n; i++) {
 924     complex<double> sum (0,0);
 925     for (int j = 0; j < n; j++) {
 926       double angle = i * j * angleIncrement;
 927       complex<double> exponentTerm (cos(angle), sin(angle));
 928       sum += input[j] * exponentTerm;
 929     }
 930     if (direction < 0) {
 931       sum /= n;
 932     }
 933     output[i] = sum;
 934   }
 935 }
 936
 937 void
 938 SignalFilter::finiteFourierTransform (const complex<double> input[], double output[], const int n, int direction)
 939 {
 940   if (direction < 0)
 941     direction = -1;
 942   else
 943     direction = 1;
 944
 945   double angleIncrement = direction * 2 * PI / n;
 946   for (int i = 0; i < n; i++) {
 947       double sumReal = 0;
 948     for (int j = 0; j < n; j++) {
 949       double angle = i * j * angleIncrement;
 950       sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
 951     }
 952     if (direction < 0) {
 953       sumReal /= n;
 954     }
 955     output[i] = sumReal;
 956   }
 957 }
 958
 959 void
 960 SignalFilter::finiteFourierTransform (const double input[], complex<double> output[], int direction) const
 961 {
 962   if (direction < 0)
 963     direction = -1;
 964   else
 965     direction = 1;
 966
 967   for (int i = 0; i < m_nFilterPoints; i++) {
 968     double sumReal = 0, sumImag = 0;
 969     for (int j = 0; j < m_nFilterPoints; j++) {
 970       int tableIndex = i * j;
 971       if (direction > 0) {
 972         sumReal += input[j] * m_vecFourierCosTable[tableIndex];
 973         sumImag += input[j] * m_vecFourierSinTable[tableIndex];
 974       } else {
 975         sumReal += input[j] * m_vecFourierCosTable[tableIndex];
 976         sumImag -= input[j] * m_vecFourierSinTable[tableIndex];
 977       }
 978     }
 979     if (direction < 0) {
 980       sumReal /= m_nFilterPoints;
 981       sumImag /= m_nFilterPoints;
 982     }
 983     output[i] = complex<double> (sumReal, sumImag);
 984   }
 985 }
 986
 987 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
 988 void
 989 SignalFilter::finiteFourierTransform (const complex<double> input[], complex<double> output[], int direction) const
 990 {
 991   if (direction < 0)
 992     direction = -1;
 993   else
 994     direction = 1;
 995
 996   for (int i = 0; i < m_nFilterPoints; i++) {
 997     double sumReal = 0, sumImag = 0;
 998     for (int j = 0; j < m_nFilterPoints; j++) {
 999       int tableIndex = i * j;
1000       if (direction > 0) {
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       } else {
1006         sumReal += input[j].real() * m_vecFourierCosTable[tableIndex]
1007           - input[j].imag() * -m_vecFourierSinTable[tableIndex];
1008         sumImag += input[j].real() * -m_vecFourierSinTable[tableIndex]
1009           + input[j].imag() * m_vecFourierCosTable[tableIndex];
1010       }
1011     }
1012     if (direction < 0) {
1013       sumReal /= m_nFilterPoints;
1014       sumImag /= m_nFilterPoints;
1015     }
1016     output[i] = complex<double> (sumReal, sumImag);
1017   }
1018 }
1019
1020 void
1021 SignalFilter::finiteFourierTransform (const complex<double> input[], double output[], int direction) const
1022 {
1023   if (direction < 0)
1024     direction = -1;
1025   else
1026     direction = 1;
1027
1028   for (int i = 0; i < m_nFilterPoints; i++) {
1029       double sumReal = 0;
1030     for (int j = 0; j < m_nFilterPoints; j++) {
1031       int tableIndex = i * j;
1032       if (direction > 0) {
1033         sumReal += input[j].real() * m_vecFourierCosTable[tableIndex]
1034           - input[j].imag() * m_vecFourierSinTable[tableIndex];
1035       } else {
1036         sumReal += input[j].real() * m_vecFourierCosTable[tableIndex]
1037           - input[j].imag() * -m_vecFourierSinTable[tableIndex];
1038       }
1039     }
1040     if (direction < 0) {
1041       sumReal /= m_nFilterPoints;
1042     }
1043     output[i] = sumReal;
1044   }
1045 }
1046
1047