r7061: initial property settings

[ctsim.git] / libctsim / procsignal.cpp

diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp
index c503096945d4dc48ca997eb32a1cd9fc818cb530..bdb9fa00626813ca37509601bd7eb7698d10c470 100644 (file)
--- a/libctsim/procsignal.cpp
+++ b/libctsim/procsignal.cpp
@@ -1,15 +1,15 @@
 /*****************************************************************************
 ** File IDENTIFICATION
 ** 
-**     Name:                   filter.cpp
-**     Purpose:                Routines for signal-procesing filters
-**     Progammer:             Kevin Rosenberg
-**     Date Started:           Aug 1984
+**     Name:            procsignal.cpp
+**     Purpose:         Routines for processing signals and projections
+**     Progammer:          Kevin Rosenberg
+**     Date Started:    Aug 1984
 **
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: procsignal.cpp,v 1.27 2001/03/01 07:30:49 kevin Exp $
+**  $Id$
 **
 **  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
@@ -28,7 +28,7 @@
 #include "ct.h"
 
 #ifdef HAVE_WXWINDOWS
-#include "dlgezplot.h"
+#include "nographics.h"
 #endif
 
 // FilterMethod ID/Names
@@ -42,23 +42,23 @@ const int ProcessSignal::FILTER_METHOD_FFTW = 4;
 const int ProcessSignal::FILTER_METHOD_RFFTW =5 ;
 #endif
 const char* const ProcessSignal::s_aszFilterMethodName[] = {
-  {"convolution"},
-  {"fourier"},
-  {"fouier-table"},
-  {"fft"},
+  "convolution",
+  "fourier",
+  "fouier-table",
+  "fft",
 #if HAVE_FFTW
-  {"fftw"},
-  {"rfftw"},
+  "fftw",
+  "rfftw",
 #endif
 };
 const char* const ProcessSignal::s_aszFilterMethodTitle[] = {
-  {"Convolution"},
-  {"Fourier"},
-  {"Fouier Trigometric Table"},
-  {"FFT"},
+  "Convolution",
+  "Fourier",
+  "Fouier Trigometric Table",
+  "FFT",
 #if HAVE_FFTW
-  {"FFTW"},
-  {"Real/Half-Complex FFTW"},
+  "FFTW",
+  "Real/Half-Complex FFTW",
 #endif
 };
 const int ProcessSignal::s_iFilterMethodCount = sizeof(s_aszFilterMethodName) / sizeof(const char*);
@@ -68,12 +68,12 @@ const int ProcessSignal::FILTER_GENERATION_INVALID = -1;
 const int ProcessSignal::FILTER_GENERATION_DIRECT = 0;
 const int ProcessSignal::FILTER_GENERATION_INVERSE_FOURIER = 1;
 const char* const ProcessSignal::s_aszFilterGenerationName[] = {
-  {"direct"},
-  {"inverse-fourier"},
+  "direct",
+  "inverse-fourier",
 };
 const char* const ProcessSignal::s_aszFilterGenerationTitle[] = {
-  {"Direct"},
-  {"Inverse Fourier"},
+  "Direct",
+  "Inverse Fourier",
 };
 const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGenerationName) / sizeof(const char*);
 
@@ -258,22 +258,10 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
     
     if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
       // calculate number of filter points with zeropadding
-      m_nFilterPoints = m_nSignalPoints;
-      if (m_iZeropad > 0) {
-        double logBase2 = log(m_nFilterPoints) / log(2);
-        int nextPowerOf2 = static_cast<int>(floor(logBase2));
-        if (logBase2 != floor(logBase2))
-          nextPowerOf2++;
-        nextPowerOf2 += (m_iZeropad - 1);
-        m_nFilterPoints = 1 << nextPowerOf2;
-#if defined(DEBUG) || defined(_DEBUG)
-        if (m_traceLevel >= Trace::TRACE_CONSOLE)
-          sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints);
-#endif
-      }
+      m_nFilterPoints = addZeropadFactor (m_nSignalPoints, m_iZeropad);
       m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
       
-      if (m_nFilterPoints % 2) { // Odd
+      if (isOdd (m_nFilterPoints)) { // Odd
         m_dFilterMin = -1. / (2 * m_dSignalInc);
         m_dFilterMax = 1. / (2 * m_dSignalInc);
         m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
@@ -297,7 +285,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
         dlgEZPlot.ShowModal();
       }
 #endif
-
+      
       // This works fairly well. I'm not sure why since scaling for geometries is done on
       // frequency filter rather than spatial filter as it should be.
       // It gives values slightly off than freq/inverse filtering
@@ -329,9 +317,9 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
         dlgEZPlot.ShowModal();
       }
 #endif
-
+      
       // FILTERING:  FREQUENCY - INVERSE FOURIER
-
+      
     } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
       // calculate number of filter points with zeropadding
       int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
@@ -430,15 +418,15 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
   }
   
   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
-    m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
-    m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
+    m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE | FFTW_USE_WISDOM);
+    m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE | FFTW_USE_WISDOM);
     m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
     m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
     for (i = 0; i < m_nFilterPoints; i++) 
       m_adRealFftInput[i] = 0;
   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
-    m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
-    m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
+    m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD,  FFTW_ESTIMATE | FFTW_USE_WISDOM);
+    m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD,  FFTW_ESTIMATE | FFTW_USE_WISDOM);
     m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
     m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
     for (i = 0; i < m_nFilterPoints; i++) 
@@ -876,3 +864,19 @@ ProcessSignal::finiteFourierTransform (const std::complex<double> input[], doubl
   }
 }
 
+
+int
+ProcessSignal::addZeropadFactor (int n, int iZeropad)
+{
+  if (iZeropad > 0) {
+    double dLogBase2 = log(n) / log(2);
+    int iLogBase2 = static_cast<int>(floor (dLogBase2));
+    int iPaddedN = 1 << (iLogBase2 + iZeropad);
+#ifdef DEBUG
+    sys_error (ERR_TRACE, "Zeropadding %d to %d", n, iPaddedN);
+#endif
+    return iPaddedN;
+  }
+
+  return n;
+}