r131: *** empty log message ***

[ctsim.git] / libctsim / filter.cpp

diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp
index 27d1829380ae5224feff90fb221a801079884f20..d3a3fbdfaf01ed399ba8e2923a55ee0fd75e3443 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.8 2000/07/04 18:33:35 kevin Exp $
+**  $Id: filter.cpp,v 1.9 2000/07/04 22:21:01 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,9 +134,9 @@ 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++) 
@@ -143,7 +144,7 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
     for (int i = 0; i < halfFilter; i++)
        m_vecFilter[m_nFilterPoints - i - 1] = static_cast<double>(i) / (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,9 +157,9 @@ 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++) 
@@ -166,7 +167,7 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
     for (int i = 0; i < halfFilter; i++)
        m_vecFilter[m_nFilterPoints - i - 1] = static_cast<double>(i) / (halfFilter - 1) /  (2 * m_signalInc);
     if (halfFilter % 2) // odd
-      m_vecFilter[halfFilter] = 1;
+      m_vecFilter[halfFilter] = 1 / (2 * m_signalInc);
 
 #if HAVE_FFTW
     m_planForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
@@ -364,11 +365,24 @@ 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++) {
+    for (int i = 0; i < m_nSignalPoints; i++) 
       output[i] = abs( complexOutput[i] );
-    }
   } 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++) {
@@ -379,11 +393,32 @@ 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();
+    }
     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);
   }