r135: *** empty log message ***

[ctsim.git] / libctsim / filter.cpp

diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp
index d3a3fbdfaf01ed399ba8e2923a55ee0fd75e3443..978a5025777a4ed4166e223242a01b853a31b0ef 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.9 2000/07/04 22:21:01 kevin Exp $
+**  $Id: filter.cpp,v 1.12 2000/07/05 17:59:26 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
@@ -142,7 +142,7 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
     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 / (2 * m_signalInc);
   } else if (m_idFilterMethod == FILTER_METHOD_FFT || m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_2 || m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_4) {
@@ -163,11 +163,11 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
     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 / (2 * m_signalInc);
+      m_vecFilter[halfFilter] = 1 / (2 * m_signalInc) / m_nSignalPoints;
 
 #if HAVE_FFTW
     m_planForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
@@ -365,25 +365,12 @@ SignalFilter::filterSignal (const float input[], double output[]) const
     complex<double> complexOutput[m_nSignalPoints];
     complex<double> filteredSignal[m_nSignalPoints];
     finiteFourierTransform (input, fftSignal, m_nSignalPoints, -1);
-    if (m_traceLevel >= TRACE_PLOT) {
-       double test[m_nSignalPoints];
-       for (int i = 0; i < m_nSignalPoints; i++) 
-           test[i] = abs(fftSignal[i]);
-       ezplot_1d(test, m_nSignalPoints);
-       cio_kb_getc();
-    }
     dotProduct (m_vecFilter, fftSignal, filteredSignal, m_nSignalPoints);
-    if (m_traceLevel >= TRACE_PLOT) {
-       double test[m_nSignalPoints];
-       for (int i = 0; i < m_nSignalPoints; i++) 
-           test[i] = abs(filteredSignal[i]);
-       ezplot_1d(test, m_nSignalPoints);
-       cio_kb_getc();
-    }
     finiteFourierTransform (filteredSignal, complexOutput, m_nSignalPoints, 1);
     for (int i = 0; i < m_nSignalPoints; i++) 
-      output[i] = abs( complexOutput[i] );
+      output[i] = complexOutput[i].real();
   } else if (m_idFilterMethod == FILTER_METHOD_FFT || FILTER_METHOD_FFT_ZEROPAD_2 || FILTER_METHOD_FFT_ZEROPAD_4) {
+#if HAVE_FFTW
     fftw_complex in[m_nFilterPoints], out[m_nFilterPoints];
     for (int i = 0; i < m_nSignalPoints; i++) {
       in[i].re = input[i];
@@ -393,35 +380,15 @@ SignalFilter::filterSignal (const float input[], double output[]) const
       in[i].re = in[i].im = 0;      // ZeroPad
     }
     fftw_one(m_planForward, in, out);
-    if (m_traceLevel >= TRACE_PLOT) {
-       double test[m_nFilterPoints];
-       for (int i = 0; i < m_nFilterPoints; i++) 
-           test[i] = sqrt(out[i].re * out[i].re + out[i].im * out[i].im);
-       ezplot_1d(test, m_nFilterPoints);
-       cio_kb_getc();
-    }
     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;
-    }
-    if (m_traceLevel >= TRACE_PLOT) {
-       double test[m_nFilterPoints];
-       for (int i = 0; i < m_nFilterPoints; i++) 
-           test[i] = sqrt(out[i].re * out[i].re + out[i].im * out[i].im);
-       ezplot_1d(test, m_nFilterPoints);
-       cio_kb_getc();
+      out[i].re = m_vecFilter[i] * out[i].re;
+      out[i].im = m_vecFilter[i] * out[i].im;
     }
     fftw_one(m_planBackward, out, in);
-    if (m_traceLevel >= TRACE_PLOT) {
-       double test[m_nFilterPoints];
-       for (int i = 0; i < m_nFilterPoints; i++) 
-           test[i] = sqrt(in[i].re * in[i].re + in[i].im * in[i].im);
-       ezplot_1d(test, m_nFilterPoints);
-       cio_kb_getc();
-    }
     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;
   }
+#endif
 }
 
 double