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