libctsim/procsignal.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: procsignal.cpp,v 1.2 2000/08/22 07:02:48 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 // FilterMethod ID/Names
  31 const int ProcessSignal::FILTER_METHOD_INVALID = -1;
  32 const int ProcessSignal::FILTER_METHOD_CONVOLUTION = 0;
  33 const int ProcessSignal::FILTER_METHOD_FOURIER = 1;
  34 const int ProcessSignal::FILTER_METHOD_FOURIER_TABLE = 2;
  35 const int ProcessSignal::FILTER_METHOD_FFT = 3;
  36 #if HAVE_FFTW
  37 const int ProcessSignal::FILTER_METHOD_FFTW = 4;
  38 const int ProcessSignal::FILTER_METHOD_RFFTW =5 ;
  39 #endif
  40 const char* ProcessSignal::s_aszFilterMethodName[] = {
  41   {"convolution"},
  42   {"fourier"},
  43   {"fouier_table"},
  44   {"fft"},
  45 #if HAVE_FFTW
  46   {"fftw"},
  47   {"rfftw"},
  48 #endif
  49 };
  50 const char* ProcessSignal::s_aszFilterMethodTitle[] = {
  51   {"Convolution"},
  52   {"Fourier"},
  53   {"Fouier Trigometric Table"},
  54   {"FFT"},
  55 #if HAVE_FFTW
  56   {"FFTW"},
  57   {"Real/Half-Complex FFTW"},
  58 #endif
  59 };
  60 const int ProcessSignal::s_iFilterMethodCount = sizeof(s_aszFilterMethodName) / sizeof(const char*);
  61
  62 // FilterGeneration ID/Names
  63 const int ProcessSignal::FILTER_GENERATION_INVALID = -1;
  64 const int ProcessSignal::FILTER_GENERATION_DIRECT = 0;
  65 const int ProcessSignal::FILTER_GENERATION_INVERSE_FOURIER = 1;
  66 const char* ProcessSignal::s_aszFilterGenerationName[] = {
  67   {"direct"},
  68   {"inverse_fourier"},
  69 };
  70 const char* ProcessSignal::s_aszFilterGenerationTitle[] = {
  71   {"Direct"},
  72   {"Inverse Fourier"},
  73 };
  74 const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGenerationName) / sizeof(const char*);
  75
  76
  77 // CLASS IDENTIFICATION
  78 //   ProcessSignal
  79 //
  80 ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMethodName, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const char* szDomainName, const char* szFilterGenerationName, int iZeropad = 0, int iPreinterpolationFactor = 1, int iTraceLevel = TRACE_NONE)
  81     : m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false)
  82 {
  83   m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName);
  84   if (m_idFilterMethod == FILTER_METHOD_INVALID) {
  85     m_fail = true;
  86     m_failMessage = "Invalid filter method name ";
  87     m_failMessage += szFilterMethodName;
  88     return;
  89   }
  90   m_idFilterGeneration = convertFilterGenerationNameToID (szFilterGenerationName);
  91   if (m_idFilterGeneration == FILTER_GENERATION_INVALID) {
  92     m_fail = true;
  93     m_failMessage = "Invalid frequency filter name ";
  94     m_failMessage += szFilterGenerationName;
  95     return;
  96   }
  97   m_idFilter = SignalFilter::convertFilterNameToID (szFilterName);
  98   if (m_idFilter == SignalFilter::FILTER_INVALID) {
  99     m_fail = true;
 100     m_failMessage = "Invalid Filter name ";
 101     m_failMessage += szFilterName;
 102     return;
 103   }
 104   m_idDomain = SignalFilter::convertDomainNameToID (szDomainName);
 105   if (m_idDomain == SignalFilter::DOMAIN_INVALID) {
 106     m_fail = true;
 107     m_failMessage = "Invalid domain name ";
 108     m_failMessage += szDomainName;
 109     return;
 110   }
 111
 112   init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel);
 113 }
 114
 115
 116 void
 117 ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const int idDomain, const int idFilterGeneration, const int iZeropad, const int iPreinterpolationFactor, int iTraceLevel)
 118 {
 119   m_idFilter = idFilter;
 120   m_idDomain = idDomain;
 121   m_idFilterMethod = idFilterMethod;
 122   m_idFilterGeneration = idFilterGeneration;
 123   if (m_idFilter == SignalFilter::FILTER_INVALID || m_idDomain == SignalFilter::DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID || m_idFilterGeneration == FILTER_GENERATION_INVALID) {
 124     m_fail = true;
 125     return;
 126   }
 127   m_traceLevel = iTraceLevel;
 128   m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod);
 129   m_nameFilterGeneration = convertFilterGenerationIDToName (m_idFilterGeneration);
 130   m_dBandwidth = dBandwidth;
 131   m_nSignalPoints = nSignalPoints;
 132   m_dSignalInc = dSignalIncrement;
 133   m_dFilterParam = dFilterParam;
 134   m_iZeropad = iZeropad;
 135   m_iPreinterpolationFactor = iPreinterpolationFactor;
 136
 137   if (m_idFilterMethod == FILTER_METHOD_FFT) {
 138 #if HAVE_FFTW
 139     m_idFilterMethod = FILTER_METHOD_RFFTW;
 140 #else
 141     m_fail = true;
 142     m_failMessage = "FFT not yet implemented";
 143     return;
 144 #endif
 145   }
 146
 147   bool m_bFrequencyFiltering = true;
 148   if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION)
 149     m_bFrequencyFiltering = false;
 150
 151   // Spatial-based filtering
 152   if (! m_bFrequencyFiltering) {
 153
 154     if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
 155         m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
 156         m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
 157         m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
 158         m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
 159         SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
 160         m_adFilter = new double[ m_nFilterPoints ];
 161         filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
 162     } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
 163         m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
 164         m_dFilterMin = -1. / (2 * m_dSignalInc);
 165         m_dFilterMax = 1. / (2 * m_dSignalInc);
 166         m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
 167         SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
 168         m_adFilter = new double[ m_nFilterPoints ];
 169         double adFrequencyFilter [m_nFilterPoints];
 170         filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints);
 171         if (m_traceLevel >= TRACE_PLOT) {
 172           SGPDriver sgpDriver ("Frequency Filter: Natural Order");
 173           SGP sgp (sgpDriver);
 174           EZPlot ezplot (sgp);
 175
 176           ezplot.ezset ("title Filter Response: Natural Order");
 177           ezplot.addCurve (adFrequencyFilter, m_nFilterPoints);
 178           ezplot.plot();
 179           cio_put_str ("Press any key to continue");
 180           cio_kb_getc ();
 181         }
 182
 183         shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
 184         if (m_traceLevel >= TRACE_PLOT) {
 185           SGPDriver sgpDriver ("Frequency Filter: Fourier Order");
 186           SGP sgp (sgpDriver);
 187           EZPlot ezplot (sgp);
 188
 189           ezplot.ezset ("title Filter Response: Fourier Order");
 190           ezplot.addCurve (adFrequencyFilter, m_nFilterPoints);
 191           ezplot.plot();
 192           cio_put_str ("Press any key to continue");
 193           cio_kb_getc ();
 194         }
 195         ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, -1);
 196         if (m_traceLevel >= TRACE_PLOT) {
 197           SGPDriver sgpDriver ("Inverse Fourier Frequency: Fourier Order");
 198           SGP sgp (sgpDriver);
 199           EZPlot ezplot (sgp);
 200
 201           ezplot.ezset ("title Inverse Fourier Frequency: Fourier Order");
 202           ezplot.addCurve (m_adFilter, m_nFilterPoints);
 203           ezplot.plot();
 204           cio_put_str ("Press any key to continue");
 205           cio_kb_getc ();
 206         }
 207         shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints);
 208         if (m_traceLevel >= TRACE_PLOT) {
 209           SGPDriver sgpDriver ("Inverse Fourier Frequency: Natural Order");
 210           SGP sgp (sgpDriver);
 211           EZPlot ezplot (sgp);
 212
 213           ezplot.ezset ("title Inverse Fourier Frequency: Natural Order");
 214           ezplot.addCurve (m_adFilter, m_nFilterPoints);
 215           ezplot.plot();
 216           cio_put_str ("Press any key to continue");
 217           cio_kb_getc ();
 218         }
 219         for (int i = 0; i < m_nFilterPoints; i++) {
 220             m_adFilter[i] /= m_dSignalInc;
 221         }
 222     }
 223   }
 224
 225   // Frequency-based filtering
 226   else if (m_bFrequencyFiltering) {
 227
 228     if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
 229       // calculate number of filter points with zeropadding
 230       m_nFilterPoints = m_nSignalPoints;
 231       if (m_iZeropad > 0) {
 232         double logBase2 = log(m_nFilterPoints) / log(2);
 233         int nextPowerOf2 = static_cast<int>(floor(logBase2));
 234         if (logBase2 != floor(logBase2))
 235           nextPowerOf2++;
 236         nextPowerOf2 += (m_iZeropad - 1);
 237         m_nFilterPoints = 1 << nextPowerOf2;
 238         if (m_traceLevel >= TRACE_TEXT)
 239         cout << "nFilterPoints = " << m_nFilterPoints << endl;
 240       }
 241       m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
 242
 243       if (m_nFilterPoints % 2) { // Odd
 244         m_dFilterMin = -1. / (2 * m_dSignalInc);
 245         m_dFilterMax = 1. / (2 * m_dSignalInc);
 246         m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
 247       } else { // Even
 248         m_dFilterMin = -1. / (2 * m_dSignalInc);
 249         m_dFilterMax = 1. / (2 * m_dSignalInc);
 250         m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints;
 251         m_dFilterMax -= m_dFilterInc;
 252       }
 253
 254       SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
 255       m_adFilter = new double [m_nFilterPoints];
 256       filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
 257       if (m_traceLevel >= TRACE_PLOT) {
 258         SGPDriver sgpDriver ("Frequency Filter: Natural Order");
 259         SGP sgp (sgpDriver);
 260         EZPlot ezplot (sgp);
 261
 262         ezplot.ezset ("title Filter Filter: Natural Order");
 263         ezplot.addCurve (m_adFilter, m_nFilterPoints);
 264           ezplot.plot();
 265           cio_put_str ("Press any key to continue");
 266           cio_kb_getc ();
 267       }
 268       shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
 269         if (m_traceLevel >= TRACE_PLOT) {
 270           SGPDriver sgpDriver ("Frequency Filter: Fourier Order");
 271           SGP sgp (sgpDriver);
 272           EZPlot ezplot (sgp);
 273
 274           ezplot.ezset ("title Filter Filter: Fourier Order");
 275           ezplot.addCurve (m_adFilter, m_nFilterPoints);
 276           ezplot.plot();
 277           cio_put_str ("Press any key to continue");
 278           cio_kb_getc ();
 279         }
 280     } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
 281       // calculate number of filter points with zeropadding
 282       int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
 283       m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
 284       m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
 285       m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1);
 286       m_nFilterPoints = nSpatialPoints;
 287       if (m_iZeropad > 0) {
 288         double logBase2 = log(nSpatialPoints) / log(2);
 289         int nextPowerOf2 = static_cast<int>(floor(logBase2));
 290         if (logBase2 != floor(logBase2))
 291           nextPowerOf2++;
 292         nextPowerOf2 += (m_iZeropad - 1);
 293         m_nFilterPoints = 1 << nextPowerOf2;
 294       }
 295       m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
 296       if (m_traceLevel >= TRACE_TEXT)
 297         cout << "nFilterPoints = " << m_nFilterPoints << endl;
 298       double adSpatialFilter [m_nFilterPoints];
 299       SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
 300       filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
 301       if (m_traceLevel >= TRACE_PLOT) {
 302         SGPDriver sgpDriver ("Spatial Filter: Natural Order");
 303         SGP sgp (sgpDriver);
 304         EZPlot ezplot (sgp);
 305
 306         ezplot.ezset ("title Spatial Filter: Natural Order");
 307         ezplot.addCurve (adSpatialFilter, nSpatialPoints);
 308         ezplot.plot();
 309         cio_put_str ("Press any key to continue");
 310         cio_kb_getc ();
 311       }
 312       for (int i = nSpatialPoints; i < m_nFilterPoints; i++)
 313         adSpatialFilter[i] = 0;
 314
 315       m_adFilter = new double [m_nFilterPoints];
 316       complex<double> acInverseFilter [m_nFilterPoints];
 317       finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, 1);
 318       for (int i = 0; i < m_nFilterPoints; i++)
 319         m_adFilter[i] = abs(acInverseFilter[i]) * m_dSignalInc;
 320       if (m_traceLevel >= TRACE_PLOT) {
 321         SGPDriver sgpDriver ("Spatial Filter: Inverse");
 322         SGP sgp (sgpDriver);
 323         EZPlot ezplot (sgp);
 324
 325         ezplot.ezset ("title Spatial Filter: Inverse");
 326         ezplot.addCurve (m_adFilter, m_nFilterPoints);
 327         ezplot.plot();
 328         cio_put_str ("Press any key to continue");
 329         cio_kb_getc ();
 330       }
 331     }
 332   }
 333
 334   // precalculate sin and cosine tables for fourier transform
 335   if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
 336     int nFourier = max(m_nFilterPoints,m_nOutputPoints) * max(m_nFilterPoints, m_nOutputPoints) + 1;
 337     double angleIncrement = (2. * PI) / m_nFilterPoints;
 338     m_adFourierCosTable = new double[ nFourier ];
 339     m_adFourierSinTable = new double[ nFourier ];
 340     double angle = 0;
 341     for (int i = 0; i < nFourier; i++) {
 342       m_adFourierCosTable[i] = cos (angle);
 343       m_adFourierSinTable[i] = sin (angle);
 344       angle += angleIncrement;
 345     }
 346   }
 347
 348 #if HAVE_FFTW
 349   if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
 350     for (int i = 0; i < m_nFilterPoints; i++)  //fftw uses unnormalized fft
 351       m_adFilter[i] /= m_nFilterPoints;
 352   }
 353
 354   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 355     m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
 356     m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
 357     m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
 358     m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
 359     for (int i = 0; i < m_nFilterPoints; i++)
 360       m_adRealFftInput[i] = 0;
 361   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 362     m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
 363     m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
 364     m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
 365     m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
 366     for (int i = 0; i < m_nFilterPoints; i++)
 367       m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
 368     for (int i = 0; i < m_nOutputPoints; i++)
 369       m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
 370   }
 371 #endif
 372
 373 }
 374
 375 ProcessSignal::~ProcessSignal (void)
 376 {
 377     delete [] m_adFourierSinTable;
 378     delete [] m_adFourierCosTable;
 379     delete [] m_adFilter;
 380
 381 #if HAVE_FFTW
 382     if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 383         fftw_destroy_plan(m_complexPlanForward);
 384         fftw_destroy_plan(m_complexPlanBackward);
 385         delete [] m_adComplexFftInput;
 386         delete [] m_adComplexFftSignal;
 387     }
 388     if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 389         rfftw_destroy_plan(m_realPlanForward);
 390         rfftw_destroy_plan(m_realPlanBackward);
 391         delete [] m_adRealFftInput;
 392         delete [] m_adRealFftSignal;
 393     }
 394 #endif
 395 }
 396
 397 int
 398 ProcessSignal::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 ProcessSignal::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 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
 424 {
 425   static const char *title = "";
 426
 427   if (fmID >= 0 && fmID < s_iFilterMethodCount)
 428       return (s_aszFilterMethodTitle [fmID]);
 429
 430   return (title);
 431 }
 432
 433
 434 int
 435 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
 436 {
 437   int fgID = FILTER_GENERATION_INVALID;
 438
 439   for (int i = 0; i < s_iFilterGenerationCount; i++)
 440    if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
 441       fgID = i;
 442       break;
 443     }
 444
 445   return (fgID);
 446 }
 447
 448 const char *
 449 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
 450 {
 451   static const char *name = "";
 452
 453   if (fgID >= 0 && fgID < s_iFilterGenerationCount)
 454       return (s_aszFilterGenerationName [fgID]);
 455
 456   return (name);
 457 }
 458
 459 const char *
 460 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
 461 {
 462   static const char *name = "";
 463
 464   if (fgID >= 0 && fgID < s_iFilterGenerationCount)
 465       return (s_aszFilterGenerationTitle [fgID]);
 466
 467   return (name);
 468 }
 469
 470 void
 471 ProcessSignal::filterSignal (const float input[], double output[]) const
 472 {
 473   if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
 474     for (int i = 0; i < m_nSignalPoints; i++)
 475       output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
 476   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
 477     double inputSignal[m_nFilterPoints];
 478     for (int i = 0; i < m_nSignalPoints; i++)
 479       inputSignal[i] = input[i];
 480     for (int i = m_nSignalPoints; i < m_nFilterPoints; i++)
 481       inputSignal[i] = 0;  // zeropad
 482     complex<double> fftSignal[m_nFilterPoints];
 483     finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);
 484     for (int i = 0; i < m_nFilterPoints; i++)
 485       fftSignal[i] *= m_adFilter[i];
 486     double inverseFourier[m_nFilterPoints];
 487     finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1);
 488     for (int i = 0; i < m_nSignalPoints; i++)
 489       output[i] = inverseFourier[i];
 490   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
 491     double inputSignal[m_nFilterPoints];
 492     for (int i = 0; i < m_nSignalPoints; i++)
 493       inputSignal[i] = input[i];
 494     for (int i = m_nSignalPoints; i < m_nFilterPoints; i++)
 495       inputSignal[i] = 0;  // zeropad
 496     complex<double> fftSignal[m_nFilterPoints];
 497     finiteFourierTransform (inputSignal, fftSignal, -1);
 498     for (int i = 0; i < m_nFilterPoints; i++)
 499       fftSignal[i] *= m_adFilter[i];
 500     double inverseFourier[m_nFilterPoints];
 501     finiteFourierTransform (fftSignal, inverseFourier, 1);
 502     for (int i = 0; i < m_nSignalPoints; i++)
 503       output[i] = inverseFourier[i];
 504   }
 505 #if HAVE_FFTW
 506   else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 507       for (int i = 0; i < m_nSignalPoints; i++)
 508           m_adRealFftInput[i] = input[i];
 509
 510       fftw_real fftOutput [ m_nFilterPoints ];
 511       rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
 512       for (int i = 0; i < m_nFilterPoints; i++)
 513           m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
 514       for (int i = m_nFilterPoints; i < m_nOutputPoints; i++)
 515         m_adRealFftSignal[i] = 0;
 516
 517       fftw_real ifftOutput [ m_nOutputPoints ];
 518       rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
 519       for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
 520           output[i] = ifftOutput[i];
 521   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 522       for (int i = 0; i < m_nSignalPoints; i++)
 523           m_adComplexFftInput[i].re = input[i];
 524
 525       fftw_complex fftOutput [ m_nFilterPoints ];
 526       fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
 527       for (int i = 0; i < m_nFilterPoints; i++) {
 528           m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
 529           m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
 530       }
 531       fftw_complex ifftOutput [ m_nOutputPoints ];
 532       fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
 533       for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
 534           output[i] = ifftOutput[i].re;
 535   }
 536 #endif
 537 }
 538
 539
 540 /* NAME
 541  *    convolve                  Discrete convolution of two functions
 542  *
 543  * SYNOPSIS
 544  *    r = convolve (f1, f2, dx, n, np, func_type)
 545  *    double r                  Convolved result
 546  *    double f1[], f2[]         Functions to be convolved
 547  *    double dx                 Difference between successive x values
 548  *    int n                     Array index to center convolution about
 549  *    int np                    Number of points in f1 array
 550  *    int func_type             EVEN or ODD or EVEN_AND_ODD function f2
 551  *
 552  * NOTES
 553  *    f1 is the projection data, its indices range from 0 to np - 1.
 554  *    The index for f2, the filter, ranges from -(np-1) to (np-1).
 555  *    There are 3 ways to handle the negative vertices of f2:
 556  *      1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
 557  *         All filters used in reconstruction are even.
 558  *      2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
 559  *      3. If f2 is both ODD AND EVEN, then we must store the value of f2
 560  *         for negative indices.  Since f2 must range from -(np-1) to (np-1),
 561  *         if we add (np - 1) to f2's array index, then f2's index will
 562  *         range from 0 to 2 * (np - 1), and the origin, x = 0, will be
 563  *         stored at f2[np-1].
 564  */
 565
 566 double
 567 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
 568 {
 569   double sum = 0.0;
 570
 571 #if UNOPTIMIZED_CONVOLUTION
 572   for (int i = 0; i < np; i++)
 573     sum += func[i] * m_adFilter[n - i + (np - 1)];
 574 #else
 575   double* f2 = m_adFilter + n + (np - 1);
 576   for (int i = 0; i < np; i++)
 577     sum += *func++ * *f2--;
 578 #endif
 579
 580   return (sum * dx);
 581 }
 582
 583
 584 double
 585 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
 586 {
 587   double sum = 0.0;
 588
 589 #if UNOPTIMIZED_CONVOLUTION
 590 for (int i = 0; i < np; i++)
 591   sum += func[i] * m_adFilter[n - i + (np - 1)];
 592 #else
 593 double* f2 = m_adFilter + n + (np - 1);
 594 for (int i = 0; i < np; i++)
 595   sum += *func++ * *f2--;
 596 #endif
 597
 598   return (sum * dx);
 599 }
 600
 601
 602 void
 603 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
 604 {
 605     complex<double> complexOutput[n];
 606
 607     finiteFourierTransform (input, complexOutput, n, direction);
 608     for (int i = 0; i < n; i++)
 609         output[i] = complexOutput[i].real();
 610 }
 611
 612 void
 613 ProcessSignal::finiteFourierTransform (const double input[], complex<double> output[], const int n, int direction)
 614 {
 615   if (direction < 0)
 616     direction = -1;
 617   else
 618     direction = 1;
 619
 620   double angleIncrement = direction * 2 * PI / n;
 621   for (int i = 0; i < n; i++) {
 622     double sumReal = 0;
 623     double sumImag = 0;
 624     for (int j = 0; j < n; j++) {
 625       double angle = i * j * angleIncrement;
 626       sumReal += input[j] * cos(angle);
 627       sumImag += input[j] * sin(angle);
 628     }
 629     if (direction < 0) {
 630       sumReal /= n;
 631       sumImag /= n;
 632     }
 633     output[i] = complex<double> (sumReal, sumImag);
 634   }
 635 }
 636
 637
 638 void
 639 ProcessSignal::finiteFourierTransform (const complex<double> input[], complex<double> output[], const int n, int direction)
 640 {
 641   if (direction < 0)
 642     direction = -1;
 643   else
 644     direction = 1;
 645
 646   double angleIncrement = direction * 2 * PI / n;
 647   for (int i = 0; i < n; i++) {
 648     complex<double> sum (0,0);
 649     for (int j = 0; j < n; j++) {
 650       double angle = i * j * angleIncrement;
 651       complex<double> exponentTerm (cos(angle), sin(angle));
 652       sum += input[j] * exponentTerm;
 653     }
 654     if (direction < 0) {
 655       sum /= n;
 656     }
 657     output[i] = sum;
 658   }
 659 }
 660
 661 void
 662 ProcessSignal::finiteFourierTransform (const complex<double> input[], double output[], const int n, int direction)
 663 {
 664   if (direction < 0)
 665     direction = -1;
 666   else
 667     direction = 1;
 668
 669   double angleIncrement = direction * 2 * PI / n;
 670   for (int i = 0; i < n; i++) {
 671       double sumReal = 0;
 672     for (int j = 0; j < n; j++) {
 673       double angle = i * j * angleIncrement;
 674       sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
 675     }
 676     if (direction < 0) {
 677       sumReal /= n;
 678     }
 679     output[i] = sumReal;
 680   }
 681 }
 682
 683 // Table-based routines
 684
 685 void
 686 ProcessSignal::finiteFourierTransform (const double input[], complex<double> output[], int direction) const
 687 {
 688   if (direction < 0)
 689     direction = -1;
 690   else
 691     direction = 1;
 692
 693   for (int i = 0; i < m_nFilterPoints; i++) {
 694     double sumReal = 0, sumImag = 0;
 695     for (int j = 0; j < m_nFilterPoints; j++) {
 696       int tableIndex = i * j;
 697       if (direction > 0) {
 698         sumReal += input[j] * m_adFourierCosTable[tableIndex];
 699         sumImag += input[j] * m_adFourierSinTable[tableIndex];
 700       } else {
 701         sumReal += input[j] * m_adFourierCosTable[tableIndex];
 702         sumImag -= input[j] * m_adFourierSinTable[tableIndex];
 703       }
 704     }
 705     if (direction < 0) {
 706       sumReal /= m_nFilterPoints;
 707       sumImag /= m_nFilterPoints;
 708     }
 709     output[i] = complex<double> (sumReal, sumImag);
 710   }
 711 }
 712
 713 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
 714 void
 715 ProcessSignal::finiteFourierTransform (const complex<double> input[], complex<double> output[], int direction) const
 716 {
 717   if (direction < 0)
 718     direction = -1;
 719   else
 720     direction = 1;
 721
 722   for (int i = 0; i < m_nFilterPoints; i++) {
 723     double sumReal = 0, sumImag = 0;
 724     for (int j = 0; j < m_nFilterPoints; j++) {
 725       int tableIndex = i * j;
 726       if (direction > 0) {
 727         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 728           - input[j].imag() * m_adFourierSinTable[tableIndex];
 729         sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
 730           + input[j].imag() * m_adFourierCosTable[tableIndex];
 731       } else {
 732         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 733           - input[j].imag() * -m_adFourierSinTable[tableIndex];
 734         sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
 735           + input[j].imag() * m_adFourierCosTable[tableIndex];
 736       }
 737     }
 738     if (direction < 0) {
 739       sumReal /= m_nFilterPoints;
 740       sumImag /= m_nFilterPoints;
 741     }
 742     output[i] = complex<double> (sumReal, sumImag);
 743   }
 744 }
 745
 746 void
 747 ProcessSignal::finiteFourierTransform (const complex<double> input[], double output[], int direction) const
 748 {
 749   if (direction < 0)
 750     direction = -1;
 751   else
 752     direction = 1;
 753
 754   for (int i = 0; i < m_nFilterPoints; i++) {
 755       double sumReal = 0;
 756     for (int j = 0; j < m_nFilterPoints; j++) {
 757       int tableIndex = i * j;
 758       if (direction > 0) {
 759         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 760           - input[j].imag() * m_adFourierSinTable[tableIndex];
 761       } else {
 762         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 763           - input[j].imag() * -m_adFourierSinTable[tableIndex];
 764       }
 765     }
 766     if (direction < 0) {
 767       sumReal /= m_nFilterPoints;
 768     }
 769     output[i] = sumReal;
 770   }
 771 }
 772
 773 // Odd Number of Points
 774 //   Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2
 775 //   Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1
 776 // Even Number of Points
 777 //   Natural Frequency Order: -n/2...-1,0,1...((n/2)-1)
 778 //   Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1
 779
 780 void
 781 ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n)
 782 {
 783   double* pdTemp = new double [n];
 784   if (n % 2) { // Odd
 785     int iHalfN = (n - 1) / 2;
 786
 787     pdTemp[0] = pdVector[iHalfN];
 788     for (int i = 0; i < iHalfN; i++)
 789       pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
 790     for (int i = 0; i < iHalfN; i++)
 791       pdTemp[i + iHalfN + 1] = pdVector[i];
 792   } else {     // Even
 793       int iHalfN = n / 2;
 794       pdTemp[0] = pdVector[iHalfN];
 795       for (int i = 0; i < iHalfN; i++)
 796         pdTemp[i + 1] = pdVector[i + iHalfN];
 797       for (int i = 0; i < iHalfN - 1; i++)
 798         pdTemp[i + iHalfN + 1] = pdVector[i];
 799   }
 800
 801   for (int i = 0; i < n; i++)
 802     pdVector[i] = pdTemp[i];
 803   delete pdTemp;
 804 }
 805
 806
 807 void
 808 ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n)
 809 {
 810   double* pdTemp = new double [n];
 811   if (n % 2) { // Odd
 812     int iHalfN = (n - 1) / 2;
 813
 814     pdTemp[iHalfN] = pdVector[0];
 815     for (int i = 0; i < iHalfN; i++)
 816       pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
 817     for (int i = 0; i < iHalfN; i++)
 818       pdTemp[i] = pdVector[i + iHalfN + 1];
 819   } else {     // Even
 820       int iHalfN = n / 2;
 821       pdTemp[iHalfN] = pdVector[0];
 822       for (int i = 0; i < iHalfN; i++)
 823         pdTemp[i] = pdVector[i + iHalfN];
 824       for (int i = 0; i < iHalfN - 1; i++)
 825         pdTemp[i + iHalfN + 1] = pdVector[i+1];
 826   }
 827
 828   for (int i = 0; i < n; i++)
 829     pdVector[i] = pdTemp[i];
 830   delete pdTemp;
 831 }
 832