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.1 2000/08/19 23:00:05 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   {"Direct Fourier"},
  53   {"Fouier Trigometric Table Lookout"},
  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)
  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);
 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)
 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 = TRACE_NONE;
 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;
 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 = m_nSignalPoints;
 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         double adInverseFilter [m_nFilterPoints];
 171         filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints);
 172         shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
 173         ProcessSignal::finiteFourierTransform (adFrequencyFilter, adInverseFilter, m_nFilterPoints, 1);
 174         for (int i = 0; i < m_nFilterPoints; i++)
 175             m_adFilter [i] = adInverseFilter[i];
 176     }
 177   }
 178
 179   // Frequency-based filtering
 180   else if (m_bFrequencyFiltering) {
 181
 182     // calculate number of filter points with zeropadding
 183     m_nFilterPoints = m_nSignalPoints;
 184     if (m_iZeropad > 0) {
 185       double logBase2 = log(m_nSignalPoints) / log(2);
 186       int nextPowerOf2 = static_cast<int>(floor(logBase2));
 187       if (logBase2 != floor(logBase2))
 188         nextPowerOf2++;
 189       nextPowerOf2 += (m_iZeropad - 1);
 190       m_nFilterPoints = 1 << nextPowerOf2;
 191       if (m_traceLevel >= TRACE_TEXT)
 192         cout << "nFilterPoints = " << m_nFilterPoints << endl;
 193     }
 194     m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
 195
 196     if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
 197         m_dFilterMin = -1. / (2 * m_dSignalInc);
 198         m_dFilterMax = 1. / (2 * m_dSignalInc);
 199         m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
 200         SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
 201         m_adFilter = new double [m_nFilterPoints];
 202         filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
 203         shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
 204     } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
 205         m_nFilterPoints = 2 * m_nSignalPoints - 1;
 206         m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
 207         m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
 208         m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
 209         double adSpatialFilter [m_nFilterPoints];
 210         double adInverseFilter [m_nFilterPoints];
 211         SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
 212         filter.copyFilterData (adSpatialFilter, 0, m_nFilterPoints);
 213         m_adFilter = new double [m_nFilterPoints];
 214         finiteFourierTransform (adSpatialFilter, adInverseFilter, m_nFilterPoints, -1);
 215         for (int i = 0; i < m_nFilterPoints; i++)
 216             m_adFilter [i] = adInverseFilter[i];
 217       }
 218     }
 219
 220   // precalculate sin and cosine tables for fourier transform
 221   if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
 222     int nFourier = max(m_nFilterPoints,m_nOutputPoints) * max(m_nFilterPoints, m_nOutputPoints) + 1;
 223     double angleIncrement = (2. * PI) / m_nFilterPoints;
 224     m_adFourierCosTable = new double[ nFourier ];
 225     m_adFourierSinTable = new double[ nFourier ];
 226     double angle = 0;
 227     for (int i = 0; i < nFourier; i++) {
 228       m_adFourierCosTable[i] = cos (angle);
 229       m_adFourierSinTable[i] = sin (angle);
 230       angle += angleIncrement;
 231     }
 232   }
 233
 234 #if HAVE_FFTW
 235   if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
 236     for (int i = 0; i < m_nFilterPoints; i++)  //fftw uses unnormalized fft
 237       m_adFilter[i] /= m_nFilterPoints;
 238   }
 239
 240   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 241     m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
 242     m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
 243     m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
 244     m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
 245     for (int i = 0; i < m_nFilterPoints; i++)
 246       m_adRealFftInput[i] = 0;
 247   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 248     m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
 249     m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
 250     m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
 251     m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
 252     for (int i = 0; i < m_nFilterPoints; i++)
 253       m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
 254     for (int i = 0; i < m_nOutputPoints; i++)
 255       m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
 256   }
 257 #endif
 258
 259 }
 260
 261 ProcessSignal::~ProcessSignal (void)
 262 {
 263     delete [] m_adFourierSinTable;
 264     delete [] m_adFourierCosTable;
 265
 266 #if HAVE_FFTW
 267     if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 268         fftw_destroy_plan(m_complexPlanForward);
 269         fftw_destroy_plan(m_complexPlanBackward);
 270         delete [] m_adComplexFftInput;
 271         delete [] m_adComplexFftSignal;
 272     }
 273     if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 274         rfftw_destroy_plan(m_realPlanForward);
 275         rfftw_destroy_plan(m_realPlanBackward);
 276         delete [] m_adRealFftInput;
 277         delete [] m_adRealFftSignal;
 278     }
 279 #endif
 280 }
 281
 282 int
 283 ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
 284 {
 285   int fmID = FILTER_METHOD_INVALID;
 286
 287   for (int i = 0; i < s_iFilterMethodCount; i++)
 288    if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
 289       fmID = i;
 290       break;
 291     }
 292
 293   return (fmID);
 294 }
 295
 296 const char *
 297 ProcessSignal::convertFilterMethodIDToName (const int fmID)
 298 {
 299   static const char *name = "";
 300
 301   if (fmID >= 0 && fmID < s_iFilterMethodCount)
 302       return (s_aszFilterMethodName [fmID]);
 303
 304   return (name);
 305 }
 306
 307 const char *
 308 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
 309 {
 310   static const char *title = "";
 311
 312   if (fmID >= 0 && fmID < s_iFilterMethodCount)
 313       return (s_aszFilterMethodTitle [fmID]);
 314
 315   return (title);
 316 }
 317
 318
 319 int
 320 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
 321 {
 322   int fgID = FILTER_GENERATION_INVALID;
 323
 324   for (int i = 0; i < s_iFilterGenerationCount; i++)
 325    if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
 326       fgID = i;
 327       break;
 328     }
 329
 330   return (fgID);
 331 }
 332
 333 const char *
 334 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
 335 {
 336   static const char *name = "";
 337
 338   if (fgID >= 0 && fgID < s_iFilterGenerationCount)
 339       return (s_aszFilterGenerationName [fgID]);
 340
 341   return (name);
 342 }
 343
 344 const char *
 345 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
 346 {
 347   static const char *name = "";
 348
 349   if (fgID >= 0 && fgID < s_iFilterGenerationCount)
 350       return (s_aszFilterGenerationTitle [fgID]);
 351
 352   return (name);
 353 }
 354
 355 void
 356 ProcessSignal::filterSignal (const float input[], double output[]) const
 357 {
 358   if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
 359     for (int i = 0; i < m_nSignalPoints; i++)
 360       output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
 361   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
 362     double inputSignal[m_nFilterPoints];
 363     for (int i = 0; i < m_nSignalPoints; i++)
 364       inputSignal[i] = input[i];
 365     for (int i = m_nSignalPoints; i < m_nFilterPoints; i++)
 366       inputSignal[i] = 0;  // zeropad
 367     complex<double> fftSignal[m_nFilterPoints];
 368     finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);
 369     for (int i = 0; i < m_nFilterPoints; i++)
 370       fftSignal[i] *= m_adFilter[i];
 371     double inverseFourier[m_nFilterPoints];
 372     finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1);
 373     for (int i = 0; i < m_nSignalPoints; i++)
 374       output[i] = inverseFourier[i];
 375   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
 376     double inputSignal[m_nFilterPoints];
 377     for (int i = 0; i < m_nSignalPoints; i++)
 378       inputSignal[i] = input[i];
 379     for (int i = m_nSignalPoints; i < m_nFilterPoints; i++)
 380       inputSignal[i] = 0;  // zeropad
 381     complex<double> fftSignal[m_nFilterPoints];
 382     finiteFourierTransform (inputSignal, fftSignal, -1);
 383     for (int i = 0; i < m_nFilterPoints; i++)
 384       fftSignal[i] *= m_adFilter[i];
 385     double inverseFourier[m_nFilterPoints];
 386     finiteFourierTransform (fftSignal, inverseFourier, 1);
 387     for (int i = 0; i < m_nSignalPoints; i++)
 388       output[i] = inverseFourier[i];
 389   }
 390 #if HAVE_FFTW
 391   else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 392       for (int i = 0; i < m_nSignalPoints; i++)
 393           m_adRealFftInput[i] = input[i];
 394
 395       fftw_real fftOutput [ m_nFilterPoints ];
 396       rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
 397       for (int i = 0; i < m_nFilterPoints; i++)
 398           m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
 399       for (int i = m_nFilterPoints; i < m_nOutputPoints; i++)
 400         m_adRealFftSignal[i] = 0;
 401
 402       fftw_real ifftOutput [ m_nOutputPoints ];
 403       rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
 404       for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
 405           output[i] = ifftOutput[i];
 406   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 407       for (int i = 0; i < m_nSignalPoints; i++)
 408           m_adComplexFftInput[i].re = input[i];
 409
 410       fftw_complex fftOutput [ m_nFilterPoints ];
 411       fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
 412       for (int i = 0; i < m_nFilterPoints; i++) {
 413           m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
 414           m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
 415       }
 416       fftw_complex ifftOutput [ m_nOutputPoints ];
 417       fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
 418       for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
 419           output[i] = ifftOutput[i].re;
 420   }
 421 #endif
 422 }
 423
 424
 425 /* NAME
 426  *    convolve                  Discrete convolution of two functions
 427  *
 428  * SYNOPSIS
 429  *    r = convolve (f1, f2, dx, n, np, func_type)
 430  *    double r                  Convolved result
 431  *    double f1[], f2[]         Functions to be convolved
 432  *    double dx                 Difference between successive x values
 433  *    int n                     Array index to center convolution about
 434  *    int np                    Number of points in f1 array
 435  *    int func_type             EVEN or ODD or EVEN_AND_ODD function f2
 436  *
 437  * NOTES
 438  *    f1 is the projection data, its indices range from 0 to np - 1.
 439  *    The index for f2, the filter, ranges from -(np-1) to (np-1).
 440  *    There are 3 ways to handle the negative vertices of f2:
 441  *      1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
 442  *         All filters used in reconstruction are even.
 443  *      2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
 444  *      3. If f2 is both ODD AND EVEN, then we must store the value of f2
 445  *         for negative indices.  Since f2 must range from -(np-1) to (np-1),
 446  *         if we add (np - 1) to f2's array index, then f2's index will
 447  *         range from 0 to 2 * (np - 1), and the origin, x = 0, will be
 448  *         stored at f2[np-1].
 449  */
 450
 451 double
 452 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
 453 {
 454   double sum = 0.0;
 455
 456 #if UNOPTIMIZED_CONVOLUTION
 457   for (int i = 0; i < np; i++)
 458     sum += func[i] * m_adFilter[n - i + (np - 1)];
 459 #else
 460   double* f2 = m_adFilter + n + (np - 1);
 461   for (int i = 0; i < np; i++)
 462     sum += *func++ * *f2--;
 463 #endif
 464
 465   return (sum * dx);
 466 }
 467
 468
 469 double
 470 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
 471 {
 472   double sum = 0.0;
 473
 474 #if UNOPTIMIZED_CONVOLUTION
 475 for (int i = 0; i < np; i++)
 476   sum += func[i] * m_adFilter[n - i + (np - 1)];
 477 #else
 478 double* f2 = m_adFilter + n + (np - 1);
 479 for (int i = 0; i < np; i++)
 480   sum += *func++ * *f2--;
 481 #endif
 482
 483   return (sum * dx);
 484 }
 485
 486
 487 void
 488 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
 489 {
 490     complex<double> complexOutput[n];
 491
 492     finiteFourierTransform (input, complexOutput, n, direction);
 493     for (int i = 0; i < n; i++)
 494         output[i] = abs(complexOutput[n]);
 495 }
 496
 497 void
 498 ProcessSignal::finiteFourierTransform (const double input[], complex<double> output[], const int n, int direction)
 499 {
 500   if (direction < 0)
 501     direction = -1;
 502   else
 503     direction = 1;
 504
 505   double angleIncrement = direction * 2 * PI / n;
 506   for (int i = 0; i < n; i++) {
 507     double sumReal = 0;
 508     double sumImag = 0;
 509     for (int j = 0; j < n; j++) {
 510       double angle = i * j * angleIncrement;
 511       sumReal += input[j] * cos(angle);
 512       sumImag += input[j] * sin(angle);
 513     }
 514     if (direction < 0) {
 515       sumReal /= n;
 516       sumImag /= n;
 517     }
 518     output[i] = complex<double> (sumReal, sumImag);
 519   }
 520 }
 521
 522
 523 void
 524 ProcessSignal::finiteFourierTransform (const complex<double> input[], complex<double> output[], const int n, int direction)
 525 {
 526   if (direction < 0)
 527     direction = -1;
 528   else
 529     direction = 1;
 530
 531   double angleIncrement = direction * 2 * PI / n;
 532   for (int i = 0; i < n; i++) {
 533     complex<double> sum (0,0);
 534     for (int j = 0; j < n; j++) {
 535       double angle = i * j * angleIncrement;
 536       complex<double> exponentTerm (cos(angle), sin(angle));
 537       sum += input[j] * exponentTerm;
 538     }
 539     if (direction < 0) {
 540       sum /= n;
 541     }
 542     output[i] = sum;
 543   }
 544 }
 545
 546 void
 547 ProcessSignal::finiteFourierTransform (const complex<double> input[], double output[], const int n, int direction)
 548 {
 549   if (direction < 0)
 550     direction = -1;
 551   else
 552     direction = 1;
 553
 554   double angleIncrement = direction * 2 * PI / n;
 555   for (int i = 0; i < n; i++) {
 556       double sumReal = 0;
 557     for (int j = 0; j < n; j++) {
 558       double angle = i * j * angleIncrement;
 559       sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
 560     }
 561     if (direction < 0) {
 562       sumReal /= n;
 563     }
 564     output[i] = sumReal;
 565   }
 566 }
 567
 568 // Table-based routines
 569
 570 void
 571 ProcessSignal::finiteFourierTransform (const double input[], complex<double> output[], int direction) const
 572 {
 573   if (direction < 0)
 574     direction = -1;
 575   else
 576     direction = 1;
 577
 578   for (int i = 0; i < m_nFilterPoints; i++) {
 579     double sumReal = 0, sumImag = 0;
 580     for (int j = 0; j < m_nFilterPoints; j++) {
 581       int tableIndex = i * j;
 582       if (direction > 0) {
 583         sumReal += input[j] * m_adFourierCosTable[tableIndex];
 584         sumImag += input[j] * m_adFourierSinTable[tableIndex];
 585       } else {
 586         sumReal += input[j] * m_adFourierCosTable[tableIndex];
 587         sumImag -= input[j] * m_adFourierSinTable[tableIndex];
 588       }
 589     }
 590     if (direction < 0) {
 591       sumReal /= m_nFilterPoints;
 592       sumImag /= m_nFilterPoints;
 593     }
 594     output[i] = complex<double> (sumReal, sumImag);
 595   }
 596 }
 597
 598 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
 599 void
 600 ProcessSignal::finiteFourierTransform (const complex<double> input[], complex<double> output[], int direction) const
 601 {
 602   if (direction < 0)
 603     direction = -1;
 604   else
 605     direction = 1;
 606
 607   for (int i = 0; i < m_nFilterPoints; i++) {
 608     double sumReal = 0, sumImag = 0;
 609     for (int j = 0; j < m_nFilterPoints; j++) {
 610       int tableIndex = i * j;
 611       if (direction > 0) {
 612         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 613           - input[j].imag() * m_adFourierSinTable[tableIndex];
 614         sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
 615           + input[j].imag() * m_adFourierCosTable[tableIndex];
 616       } else {
 617         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 618           - input[j].imag() * -m_adFourierSinTable[tableIndex];
 619         sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
 620           + input[j].imag() * m_adFourierCosTable[tableIndex];
 621       }
 622     }
 623     if (direction < 0) {
 624       sumReal /= m_nFilterPoints;
 625       sumImag /= m_nFilterPoints;
 626     }
 627     output[i] = complex<double> (sumReal, sumImag);
 628   }
 629 }
 630
 631 void
 632 ProcessSignal::finiteFourierTransform (const complex<double> input[], double output[], int direction) const
 633 {
 634   if (direction < 0)
 635     direction = -1;
 636   else
 637     direction = 1;
 638
 639   for (int i = 0; i < m_nFilterPoints; i++) {
 640       double sumReal = 0;
 641     for (int j = 0; j < m_nFilterPoints; j++) {
 642       int tableIndex = i * j;
 643       if (direction > 0) {
 644         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 645           - input[j].imag() * m_adFourierSinTable[tableIndex];
 646       } else {
 647         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 648           - input[j].imag() * -m_adFourierSinTable[tableIndex];
 649       }
 650     }
 651     if (direction < 0) {
 652       sumReal /= m_nFilterPoints;
 653     }
 654     output[i] = sumReal;
 655   }
 656 }
 657
 658 // Odd Number of Points
 659 //   Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2
 660 //   Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1
 661 // Even Number of Points
 662 //   Natural Frequency Order: -n/2...-1,0,1...((n/2)-1)
 663 //   Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1
 664
 665 void
 666 ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n)
 667 {
 668   double* pdTemp = new double [n];
 669   if (n % 2) { // Odd
 670     int iHalfN = (n - 1) / 2;
 671
 672     pdTemp[0] = pdVector[iHalfN];
 673     for (int i = 1; i <= iHalfN; i++)
 674       pdTemp[i] = pdVector[i+iHalfN];
 675     for (int i = iHalfN+1; i < n; i++)
 676       pdTemp[i] = pdVector[i-iHalfN];
 677   } else {     // Even
 678       int iHalfN = n / 2;
 679       pdTemp[0] = pdVector[iHalfN];
 680   }
 681
 682   for (int i = 0; i < n; i++)
 683     pdVector[i] = pdTemp[i];
 684   delete pdTemp;
 685 }
 686
 687
 688 void
 689 ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n)
 690 {
 691   double* pdTemp = new double [n];
 692   if (n % 2) { // Odd
 693     int iHalfN = (n - 1) / 2;
 694
 695     pdTemp[iHalfN] = pdVector[0];
 696     for (int i = 1; i <= iHalfN; i++)
 697       pdTemp[i] = pdVector[i+iHalfN];
 698     for (int i = iHalfN+1; i < n; i++)
 699       pdTemp[i] = pdVector[i-iHalfN];
 700   } else {     // Even
 701       int iHalfN = n / 2;
 702       pdTemp[iHalfN] = pdVector[0];
 703   }
 704
 705   for (int i = 0; i < n; i++)
 706     pdVector[i] = pdTemp[i];
 707   delete pdTemp;
 708 }
 709