r133: *** empty log message ***

[ctsim.git] / libctsim / filter.cpp

diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp
index 27d1829380ae5224feff90fb221a801079884f20..2d315252f43012a0f8d3ea2468ba53b0b19c6ba7 100644 (file)
--- a/libctsim/filter.cpp
+++ b/libctsim/filter.cpp
@@ -1,5 +1,5 @@
 /*****************************************************************************
-** FILE IDENTIFICATION
+** File IDENTIFICATION
 ** 
 **     Name:                   filter.cpp
 **     Purpose:                Routines for signal-procesing filters
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: filter.cpp,v 1.8 2000/07/04 18:33:35 kevin Exp $
+**  $Id: filter.cpp,v 1.11 2000/07/05 17:24:33 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
@@ -115,6 +115,7 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
     m_fail = true;
     return;
   }
+  m_traceLevel = TRACE_NONE;
   m_nameFilter = convertFilterIDToName (m_idFilter);
   m_nameDomain = convertDomainIDToName (m_idDomain);
   m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod);
@@ -133,17 +134,17 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
       m_vecFourierSinTable[i] = sin (angleIncrement * i);
     }
     m_nFilterPoints = m_nSignalPoints;
-    m_filterMin = 0;
-    m_filterMax = m_nSignalPoints * m_signalInc;
-    m_filterInc = (m_filterMax - m_filterMin) / (m_nFilterPoints - 1);
+    m_filterMin = -1. / (2 * m_signalInc);
+    m_filterMax = 1. / (2 * m_signalInc);
+    m_filterInc = (m_filterMax - m_filterMin) / m_nFilterPoints;
     m_vecFilter = new double [m_nFilterPoints];
     int halfFilter = m_nFilterPoints / 2;
     for (int i = 0; i < halfFilter; i++) 
        m_vecFilter[i] = static_cast<double>(i) / (halfFilter - 1) / (2 * m_signalInc);
     for (int i = 0; i < halfFilter; i++)
-       m_vecFilter[m_nFilterPoints - i - 1] = static_cast<double>(i) / (halfFilter - 1) / (2 * m_signalInc);
+       m_vecFilter[m_nFilterPoints - i - 1] = static_cast<double>(i+1) / (halfFilter - 1) / (2 * m_signalInc);
     if (halfFilter % 2) // odd
-      m_vecFilter[halfFilter] = 1;
+      m_vecFilter[halfFilter] = 1 / (2 * m_signalInc);
   } else if (m_idFilterMethod == FILTER_METHOD_FFT || m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_2 || m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_4) {
     m_nFilterPoints = m_nSignalPoints;
     if (m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_2 || m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_4) {
@@ -156,17 +157,17 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
       m_nFilterPoints = 1 << nextPowerOf2;
       cout << "nFilterPoints = " << m_nFilterPoints << endl;
     }
-    m_filterMin = 0;
-    m_filterMax = m_nSignalPoints * m_signalInc;
-    m_filterInc = (m_filterMax - m_filterMin) / (m_nFilterPoints - 1);
+    m_filterMin = -1. / (2 * m_signalInc);
+    m_filterMax = 1. / (2 * m_signalInc);
+    m_filterInc = (m_filterMax - m_filterMin) / m_nFilterPoints;
     m_vecFilter = new double [m_nFilterPoints];
     int halfFilter = m_nFilterPoints / 2;
     for (int i = 0; i < halfFilter; i++) 
-       m_vecFilter[i] = static_cast<double>(i) / (halfFilter - 1) /  (2 * m_signalInc);
+       m_vecFilter[i] = static_cast<double>(i) / (halfFilter - 1) /  (2 * m_signalInc) / m_nSignalPoints;
     for (int i = 0; i < halfFilter; i++)
-       m_vecFilter[m_nFilterPoints - i - 1] = static_cast<double>(i) / (halfFilter - 1) /  (2 * m_signalInc);
+       m_vecFilter[m_nFilterPoints - i - 1] = static_cast<double>(i+1) / (halfFilter - 1) /  (2 * m_signalInc) / m_nSignalPoints;
     if (halfFilter % 2) // odd
-      m_vecFilter[halfFilter] = 1;
+      m_vecFilter[halfFilter] = 1 / (2 * m_signalInc) / m_nSignalPoints;
 
 #if HAVE_FFTW
     m_planForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
@@ -366,9 +367,8 @@ SignalFilter::filterSignal (const float input[], double output[]) const
     finiteFourierTransform (input, fftSignal, m_nSignalPoints, -1);
     dotProduct (m_vecFilter, fftSignal, filteredSignal, m_nSignalPoints);
     finiteFourierTransform (filteredSignal, complexOutput, m_nSignalPoints, 1);
-    for (int i = 0; i < m_nSignalPoints; i++) {
-      output[i] = abs( complexOutput[i] );
-    }
+    for (int i = 0; i < m_nSignalPoints; i++) 
+      output[i] = complexOutput[i].real();
   } else if (m_idFilterMethod == FILTER_METHOD_FFT || FILTER_METHOD_FFT_ZEROPAD_2 || FILTER_METHOD_FFT_ZEROPAD_4) {
     fftw_complex in[m_nFilterPoints], out[m_nFilterPoints];
     for (int i = 0; i < m_nSignalPoints; i++) {
@@ -380,12 +380,12 @@ SignalFilter::filterSignal (const float input[], double output[]) const
     }
     fftw_one(m_planForward, in, out);
     for (int i = 0; i < m_nFilterPoints; i++) {
-      out[i].re = m_vecFilter[i] * out[i].re / m_nSignalPoints;
-      out[i].im = m_vecFilter[i] * out[i].im / m_nSignalPoints;
+      out[i].re = m_vecFilter[i] * out[i].re;
+      out[i].im = m_vecFilter[i] * out[i].im;
     }
     fftw_one(m_planBackward, out, in);
     for (int i = 0; i < m_nSignalPoints; i++) 
-      output[i] = sqrt (in[i].re * in[i].re + in[i].im * in[i].im);
+      output[i] = in[i].re;
   }
 }