X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Ffilter.cpp;h=978a5025777a4ed4166e223242a01b853a31b0ef;hb=cabc5359e111e9d8a9554263ece8e51d8de1ef86;hp=d3a3fbdfaf01ed399ba8e2923a55ee0fd75e3443;hpb=03c4bc505de131323cbc70a70bdceda9229269bf;p=ctsim.git

diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp
index d3a3fbd..978a502 100644
--- 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