projects / ctsim.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
r328: *** empty log message ***
[ctsim.git] / libctsim / procsignal.cpp
diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp
index a7933b21eb2a4ffb508c88d5caa35ffc834d3723..41895ac06a37308f297e000bf8cf9b1ad5b43551 100644 (file)
--- a/libctsim/procsignal.cpp
+++ b/libctsim/procsignal.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: procsignal.cpp,v 1.11 2000/12/29 15:45:06 kevin Exp $
+**  $Id: procsignal.cpp,v 1.13 2001/01/02 05:34:57 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
@@ -189,7 +189,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
         pEZPlot->plot (pSGP);
       }
 #endif     
-      shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
+      Fourier::shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
 #ifdef HAVE_SGP
       if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
         pEZPlot->ezset ("title Filter Response: Fourier Order");
@@ -199,7 +199,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
         pEZPlot->plot (pSGP);
       }
 #endif
-      ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, -1);
+      ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, FORWARD);
       delete adFrequencyFilter;\r
 #ifdef HAVE_SGP
       if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
@@ -210,7 +210,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
         pEZPlot->plot (pSGP);
       }
 #endif
-      shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints);
+      Fourier::shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints);
 #ifdef HAVE_SGP
       if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
         pEZPlot->ezset ("title Inverse Fourier Frequency: Natural Order");
@@ -304,7 +304,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
         pEZPlot->plot (pSGP);
       }
 #endif
-      shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
+      Fourier::shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
 #ifdef HAVE_SGP
       if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
         pEZPlot->ezset ("title Filter Filter: Fourier Order");
@@ -370,10 +370,10 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
       
       m_adFilter = new double [m_nFilterPoints];
       std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];\r
-      finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, 1);
+      finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
       delete adSpatialFilter;\r
       for (i = 0; i < m_nFilterPoints; i++)
-        m_adFilter[i] = std::abs(acInverseFilter[i]) * m_dSignalInc;
+        m_adFilter[i] = std::abs (acInverseFilter[i]) * m_dSignalInc;
       delete acInverseFilter;\r
 #ifdef HAVE_SGP
       if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
@@ -553,12 +553,12 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const
     for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
       inputSignal[i] = 0;  // zeropad
     std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
-    finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);\r
+    finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);\r
     delete inputSignal;
     for (i = 0; i < m_nFilterPoints; i++)
       fftSignal[i] *= m_adFilter[i];
     double* inverseFourier = new double [m_nFilterPoints];
-    finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1);\r
+    finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);\r
     delete fftSignal;
     for (i = 0; i < m_nSignalPoints; i++) 
       output[i] = inverseFourier[i];\r
@@ -570,12 +570,12 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const
     for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
       inputSignal[i] = 0;  // zeropad
     std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
-    finiteFourierTransform (inputSignal, fftSignal, -1);\r
+    finiteFourierTransform (inputSignal, fftSignal, FORWARD);\r
     delete inputSignal;
     for (i = 0; i < m_nFilterPoints; i++)
       fftSignal[i] *= m_adFilter[i];
     double* inverseFourier = new double [m_nFilterPoints];
-    finiteFourierTransform (fftSignal, inverseFourier, 1);\r
+    finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);\r
     delete fftSignal;
     for (i = 0; i < m_nSignalPoints; i++) 
       output[i] = inverseFourier[i];\r
@@ -733,8 +733,8 @@ ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::
     std::complex<double> sum (0,0);
     for (int j = 0; j < n; j++) {
       double angle = i * j * angleIncrement;
-      std::complex<double> exponentTerm (cos(angle), sin(angle));
-      sum += input[j] * exponentTerm;
+      std::complex<double> exponentTerm (cos(angle), sin(angle));\r
+      sum += input[j] * exponentTerm;\r
     }
     if (direction < 0) {
       sum /= n;
@@ -855,119 +855,3 @@ ProcessSignal::finiteFourierTransform (const std::complex<double> input[], doubl
   }
 }
 
-// Odd Number of Points
-//   Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2
-//   Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1
-// Even Number of Points
-//   Natural Frequency Order: -n/2...-1,0,1...((n/2)-1)
-//   Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1
-
-void\r
-ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n)\r
-{\r
-  double* pdTemp = new double [n];\r
-  int i;\r
-  if (n % 2) { // Odd\r
-    int iHalfN = (n - 1) / 2;\r
-    \r
-    pdTemp[0] = pdVector[iHalfN];\r
-    for (i = 0; i < iHalfN; i++)\r
-      pdTemp[i + 1] = pdVector[i + 1 + iHalfN];\r
-    for (i = 0; i < iHalfN; i++)\r
-      pdTemp[i + iHalfN + 1] = pdVector[i];\r
-  } else {     // Even\r
-    int iHalfN = n / 2;\r
-    pdTemp[0] = pdVector[iHalfN];\r
-    for (i = 0; i < iHalfN; i++)\r
-      pdTemp[i + 1] = pdVector[i + iHalfN];\r
-    for (i = 0; i < iHalfN - 1; i++)\r
-      pdTemp[i + iHalfN + 1] = pdVector[i];\r
-  }\r
-  \r
-  for (i = 0; i < n; i++)\r
-    pdVector[i] = pdTemp[i];\r
-  delete pdTemp;\r
-}\r
-\r
-void\r
-ProcessSignal::shuffleNaturalToFourierOrder (std::complex<double>* pdVector, const int n)\r
-{\r
-  std::complex<double>* pdTemp = new std::complex<double> [n];\r
-  int i;\r
-  if (n % 2) { // Odd\r
-    int iHalfN = (n - 1) / 2;\r
-    \r
-    pdTemp[0] = pdVector[iHalfN];\r
-    for (i = 0; i < iHalfN; i++)\r
-      pdTemp[i + 1] = pdVector[i + 1 + iHalfN];\r
-    for (i = 0; i < iHalfN; i++)\r
-      pdTemp[i + iHalfN + 1] = pdVector[i];\r
-  } else {     // Even\r
-    int iHalfN = n / 2;\r
-    pdTemp[0] = pdVector[iHalfN];\r
-    for (i = 0; i < iHalfN; i++)\r
-      pdTemp[i + 1] = pdVector[i + iHalfN];\r
-    for (i = 0; i < iHalfN - 1; i++)\r
-      pdTemp[i + iHalfN + 1] = pdVector[i];\r
-  }\r
-  \r
-  for (i = 0; i < n; i++)\r
-    pdVector[i] = pdTemp[i];\r
-  delete [] pdTemp;\r
-}\r
-\r
-
-void\r
-ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n)\r
-{\r
-  double* pdTemp = new double [n];\r
-  int i;\r
-  if (n % 2) { // Odd\r
-    int iHalfN = (n - 1) / 2;\r
-    \r
-    pdTemp[iHalfN] = pdVector[0];\r
-    for (i = 0; i < iHalfN; i++)\r
-      pdTemp[i + 1 + iHalfN] = pdVector[i + 1];\r
-    for (i = 0; i < iHalfN; i++)\r
-      pdTemp[i] = pdVector[i + iHalfN + 1];\r
-  } else {     // Even\r
-    int iHalfN = n / 2;\r
-    pdTemp[iHalfN] = pdVector[0];\r
-    for (i = 0; i < iHalfN; i++)\r
-      pdTemp[i] = pdVector[i + iHalfN];\r
-    for (i = 0; i < iHalfN - 1; i++)\r
-      pdTemp[i + iHalfN + 1] = pdVector[i+1];\r
-  }\r
-  \r
-  for (i = 0; i < n; i++)\r
-    pdVector[i] = pdTemp[i];\r
-  delete pdTemp;\r
-}\r
-\r
-void\r
-ProcessSignal::shuffleFourierToNaturalOrder (std::complex<double>* pdVector, const int n)\r
-{\r
-  std::complex<double>* pdTemp = new std::complex<double> [n];\r
-  int i;\r
-  if (n % 2) { // Odd\r
-    int iHalfN = (n - 1) / 2;\r
-    \r
-    pdTemp[iHalfN] = pdVector[0];\r
-    for (i = 0; i < iHalfN; i++)\r
-      pdTemp[i + 1 + iHalfN] = pdVector[i + 1];\r
-    for (i = 0; i < iHalfN; i++)\r
-      pdTemp[i] = pdVector[i + iHalfN + 1];\r
-  } else {     // Even\r
-    int iHalfN = n / 2;\r
-    pdTemp[iHalfN] = pdVector[0];\r
-    for (i = 0; i < iHalfN; i++)\r
-      pdTemp[i] = pdVector[i + iHalfN];\r
-    for (i = 0; i < iHalfN - 1; i++)\r
-      pdTemp[i + iHalfN + 1] = pdVector[i+1];\r
-  }\r
-  \r
-  for (i = 0; i < n; i++)\r
-    pdVector[i] = pdTemp[i];\r
-  delete [] pdTemp;\r
-}\r
-\r
Computed Tomography simulator and reconstructor