r138: *** empty log message ***

[ctsim.git] / libctsim / filter.cpp

diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp
index 6a344d42a6273957a8d847c264bcdebd438c0357..65163f3185b0be2de8beb7342ad643f0701ff875 100644 (file)
--- a/libctsim/filter.cpp
+++ b/libctsim/filter.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: filter.cpp,v 1.14 2000/07/06 18:37:24 kevin Exp $
+**  $Id: filter.cpp,v 1.15 2000/07/07 15:30:59 kevin Exp $
 **
 **  This program is free software; you can redistribute it and/or modify
 **  it under the terms of the GNU General Public License (version 2) as
@@ -129,10 +129,21 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
   m_vecFourierSinTable = NULL;
   m_vecFilter = NULL;
 
-  if (m_idFilterMethod == FILTER_METHOD_FFT)
-    m_idFilterMethod = FILTER_METHOD_FFTW;
+  if (m_idFilterMethod == FILTER_METHOD_FFT) {
+#if HAVE_FFTW
+    m_idFilterMethod = FILTER_METHOD_RFFTW;
+#else
+    m_fail = true;
+    m_failMessage = "FFT not yet implemented";
+    return;
+#endif
+  }
 
-  if (m_idFilterMethod == FILTER_METHOD_FOURIER || FILTER_METHOD_FOURIER_TABLE || m_idFilterMethod == FILTER_METHOD_FFT || m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
+  if (m_idFilterMethod == FILTER_METHOD_FOURIER || FILTER_METHOD_FOURIER_TABLE || m_idFilterMethod == FILTER_METHOD_FFT
+#if HAVE_FFTW
+      || m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW
+#endif
+      ) {
     m_nFilterPoints = m_nSignalPoints;
     if (m_zeropad > 0) {
       double logBase2 = log(m_nSignalPoints) / log(2);
@@ -321,10 +332,12 @@ SignalFilter::convertFilterMethodNameToID (const char* const filterMethodName)
     fmID = FILTER_METHOD_FOURIER_TABLE;
   else if (strcasecmp (filterMethodName, FILTER_METHOD_FFT_STR) == 0)
     fmID = FILTER_METHOD_FFT;
+#if HAVE_FFTW
   else if (strcasecmp (filterMethodName, FILTER_METHOD_FFTW_STR) == 0)
     fmID = FILTER_METHOD_FFTW;
   else if (strcasecmp (filterMethodName, FILTER_METHOD_RFFTW_STR) == 0)
     fmID = FILTER_METHOD_RFFTW;
+#endif
 
   return (fmID);
 }
@@ -342,10 +355,12 @@ SignalFilter::convertFilterMethodIDToName (const FilterMethodID fmID)
     return (FILTER_METHOD_FOURIER_TABLE_STR);
   else if (fmID == FILTER_METHOD_FFT)
     return (FILTER_METHOD_FFT_STR);
+#if HAVE_FFTW
   else if (fmID == FILTER_METHOD_FFTW)
     return (FILTER_METHOD_FFTW_STR);
   else if (fmID == FILTER_METHOD_RFFTW)
     return (FILTER_METHOD_RFFTW_STR);
+#endif
 
   return (name);
 }
@@ -391,10 +406,10 @@ SignalFilter::filterSignal (const float input[], double output[]) const
       inputSignal[i] = 0;  // zeropad
     complex<double> fftSignal[m_nFilterPoints];
     finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);
-    complex<double> filteredSignal[m_nFilterPoints];
-    dotProduct (m_vecFilter, fftSignal, filteredSignal, m_nFilterPoints);
+    for (int i = 0; i < m_nFilterPoints; i++)
+      fftSignal[i] *= m_vecFilter[i];
     double inverseFourier[m_nFilterPoints];
-    finiteFourierTransform (filteredSignal, inverseFourier, m_nFilterPoints, 1);
+    finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1);
     for (int i = 0; i < m_nSignalPoints; i++) 
       output[i] = inverseFourier[i];
   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
@@ -405,10 +420,10 @@ SignalFilter::filterSignal (const float input[], double output[]) const
       inputSignal[i] = 0;  // zeropad
     complex<double> fftSignal[m_nFilterPoints];
     finiteFourierTransform (inputSignal, fftSignal, -1);
-    complex<double> filteredSignal[m_nFilterPoints];
-    dotProduct (m_vecFilter, fftSignal, filteredSignal, m_nFilterPoints);
+    for (int i = 0; i < m_nFilterPoints; i++)
+      fftSignal[i] *= m_vecFilter[i];
     double inverseFourier[m_nFilterPoints];
-    finiteFourierTransform (filteredSignal, inverseFourier, 1);
+    finiteFourierTransform (fftSignal, inverseFourier, 1);
     for (int i = 0; i < m_nSignalPoints; i++) 
       output[i] = inverseFourier[i];
   }
@@ -930,9 +945,3 @@ SignalFilter::finiteFourierTransform (const complex<double> input[], double outp
 }
 
 
-void 
-SignalFilter::dotProduct (const double v1[], const complex<double> v2[], complex<double> output[], const int n)
-{
-    for (int i = 0; i < n; i++)
-       output[i] = v1[i] * v2[i];
-}