X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprocsignal.cpp;h=8c45da42d91300ef3c04902dad49ba5255a58548;hp=a7933b21eb2a4ffb508c88d5caa35ffc834d3723;hb=7ec2cd66921180a624813dff9f8bac76c6b268cc;hpb=bc5a9ca28bc4da3691fb859945d2862f6155835b

diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp
index a7933b2..8c45da4 100644
--- 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.12 2001/01/01 10:14:34 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
@@ -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;
 #ifdef HAVE_SGP
       if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
@@ -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];
-      finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, 1);
+      finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
       delete adSpatialFilter;
       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;
 #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);
+    finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
     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);
+    finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);
     delete fftSignal;
     for (i = 0; i < m_nSignalPoints; i++) 
       output[i] = inverseFourier[i];
@@ -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);
+    finiteFourierTransform (inputSignal, fftSignal, FORWARD);
     delete inputSignal;
     for (i = 0; i < m_nFilterPoints; i++)
       fftSignal[i] *= m_adFilter[i];
     double* inverseFourier = new double [m_nFilterPoints];
-    finiteFourierTransform (fftSignal, inverseFourier, 1);
+    finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);
     delete fftSignal;
     for (i = 0; i < m_nSignalPoints; i++) 
       output[i] = inverseFourier[i];
@@ -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));
+      sum += input[j] * exponentTerm;
     }
     if (direction < 0) {
       sum /= n;
@@ -878,10 +878,17 @@ ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n)
   } else {     // Even
     int iHalfN = n / 2;
     pdTemp[0] = pdVector[iHalfN];
+#if USE_BROKEN_SHUFFLE
     for (i = 0; i < iHalfN; i++)
       pdTemp[i + 1] = pdVector[i + iHalfN];
     for (i = 0; i < iHalfN - 1; i++)
       pdTemp[i + iHalfN + 1] = pdVector[i];
+#else
+    for (i = 0; i < iHalfN - 1; i++)
+      pdTemp[i + 1] = pdVector[i + iHalfN + 1];
+    for (i = 0; i < iHalfN; i++)
+      pdTemp[i + iHalfN] = pdVector[i];
+#endif
   }
   
   for (i = 0; i < n; i++)
@@ -905,10 +912,17 @@ ProcessSignal::shuffleNaturalToFourierOrder (std::complex<double>* pdVector, con
   } else {     // Even
     int iHalfN = n / 2;
     pdTemp[0] = pdVector[iHalfN];
+#if USE_BROKEN_SHUFFLE
     for (i = 0; i < iHalfN; i++)
       pdTemp[i + 1] = pdVector[i + iHalfN];
     for (i = 0; i < iHalfN - 1; i++)
       pdTemp[i + iHalfN + 1] = pdVector[i];
+#else
+    for (i = 0; i < iHalfN - 1; i++)
+      pdTemp[i + 1] = pdVector[i + iHalfN + 1];
+    for (i = 0; i < iHalfN; i++)
+      pdTemp[i + iHalfN] = pdVector[i];
+#endif
   }
   
   for (i = 0; i < n; i++)
@@ -916,7 +930,43 @@ ProcessSignal::shuffleNaturalToFourierOrder (std::complex<double>* pdVector, con
   delete [] pdTemp;
 }
 
-
+
+void
+ProcessSignal::shuffleNaturalToFourierOrder (float* pdVector, const int n)
+{
+  float* pdTemp = new float [n];
+  int i;
+  if (n % 2) { // Odd
+    int iHalfN = (n - 1) / 2;
+    
+    pdTemp[0] = pdVector[iHalfN];
+    for (i = 0; i < iHalfN; i++)
+      pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
+    for (i = 0; i < iHalfN; i++)
+      pdTemp[i + iHalfN + 1] = pdVector[i];
+  } else {     // Even
+    int iHalfN = n / 2;
+    pdTemp[0] = pdVector[iHalfN];
+#if USE_BROKEN_SHUFFLE
+    for (i = 0; i < iHalfN; i++)
+      pdTemp[i + 1] = pdVector[i + iHalfN];
+    for (i = 0; i < iHalfN - 1; i++)
+      pdTemp[i + iHalfN + 1] = pdVector[i];
+#else
+    for (i = 0; i < iHalfN - 1; i++)
+      pdTemp[i + 1] = pdVector[i + iHalfN + 1];
+    for (i = 0; i < iHalfN; i++)
+      pdTemp[i + iHalfN] = pdVector[i];
+#endif
+  }
+  
+  for (i = 0; i < n; i++)
+    pdVector[i] = pdTemp[i];
+  delete pdTemp;
+}
+
+
+
 void
 ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n)
 {
@@ -944,6 +994,7 @@ ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n)
   delete pdTemp;
 }
 
+
 void
 ProcessSignal::shuffleFourierToNaturalOrder (std::complex<double>* pdVector, const int n)
 {
@@ -971,3 +1022,33 @@ ProcessSignal::shuffleFourierToNaturalOrder (std::complex<double>* pdVector, con
   delete [] pdTemp;
 }
 
+
+
+
+void
+ProcessSignal::shuffleFourierToNaturalOrder (float* pVector, const int n)
+{
+  float* pTemp = new float [n];
+  int i;
+  if (n % 2) { // Odd
+    int iHalfN = (n - 1) / 2;
+    
+    pTemp[iHalfN] = pVector[0];
+    for (i = 0; i < iHalfN; i++)
+      pTemp[i + 1 + iHalfN] = pVector[i + 1];
+    for (i = 0; i < iHalfN; i++)
+      pTemp[i] = pVector[i + iHalfN + 1];
+  } else {     // Even
+    int iHalfN = n / 2;
+    pTemp[iHalfN] = pVector[0];
+    for (i = 0; i < iHalfN; i++)
+      pTemp[i] = pVector[i + iHalfN];
+    for (i = 0; i < iHalfN - 1; i++)
+      pTemp[i + iHalfN + 1] = pVector[i+1];
+  }
+  
+  for (i = 0; i < n; i++)
+    pVector[i] = pTemp[i];
+  delete [] pTemp;
+}
+