r7939: Initial working version of fftw3 conversion -- tested only with pjrec

[ctsim.git] / libctsim / procsignal.cpp

diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp
index bdb9fa00626813ca37509601bd7eb7698d10c470..2206ae1838cf41ebd3f2e13eb34353fdb060670e 100644 (file)
--- a/libctsim/procsignal.cpp
+++ b/libctsim/procsignal.cpp
@@ -418,21 +418,26 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
   }
   
   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
   }
   
   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {

-    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 ];
+    m_adRealFftInput = static_cast<double*>(fftw_malloc (sizeof(double) * m_nFilterPoints));
+    m_adRealFftOutput = static_cast<double*>(fftw_malloc (sizeof(double) * m_nFilterPoints));
+    m_realPlanForward = fftw_plan_r2r_1d (m_nFilterPoints, m_adRealFftInput, m_adRealFftOutput, FFTW_R2HC, FFTW_ESTIMATE);
+    m_adRealFftSignal = static_cast<double*>(fftw_malloc (sizeof(double) *  m_nOutputPoints));
+    m_adRealFftBackwardOutput = static_cast<double*>(fftw_malloc (sizeof(double) * m_nOutputPoints));
+    m_realPlanBackward = fftw_plan_r2r_1d (m_nOutputPoints, m_adRealFftSignal, m_adRealFftBackwardOutput, FFTW_HC2R, FFTW_ESTIMATE);

     for (i = 0; i < m_nFilterPoints; i++) 
       m_adRealFftInput[i] = 0;
   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
     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 | 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 ];
+    m_adComplexFftInput = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nFilterPoints));
+    m_adComplexFftOutput = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nFilterPoints));
+    m_complexPlanForward = fftw_plan_dft_1d (m_nFilterPoints, m_adComplexFftInput, m_adComplexFftOutput, FFTW_FORWARD,  FFTW_ESTIMATE);
+    m_adComplexFftSignal = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nOutputPoints));
+    m_adComplexFftBackwardOutput = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nOutputPoints));
+    m_complexPlanBackward = fftw_plan_dft_1d (m_nOutputPoints, m_adComplexFftSignal, m_adComplexFftBackwardOutput, FFTW_BACKWARD,  FFTW_ESTIMATE);
+

     for (i = 0; i < m_nFilterPoints; i++) 
     for (i = 0; i < m_nFilterPoints; i++) 

-      m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
+      m_adComplexFftInput[i][0] = m_adComplexFftInput[i][1] = 0;

     for (i = 0; i < m_nOutputPoints; i++) 
     for (i = 0; i < m_nOutputPoints; i++) 

-      m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
+      m_adComplexFftSignal[i][0] = m_adComplexFftSignal[i][1] = 0;

   }
 #endif
   
   }
 #endif
   


@@ -448,14 +453,18 @@ ProcessSignal::~ProcessSignal (void)
   if (m_idFilterMethod == FILTER_METHOD_FFTW) {
     fftw_destroy_plan(m_complexPlanForward);
     fftw_destroy_plan(m_complexPlanBackward);
   if (m_idFilterMethod == FILTER_METHOD_FFTW) {
     fftw_destroy_plan(m_complexPlanForward);
     fftw_destroy_plan(m_complexPlanBackward);

-    delete [] m_adComplexFftInput;
-    delete [] m_adComplexFftSignal;
+    fftw_free (m_adComplexFftInput);
+    fftw_free (m_adComplexFftOutput);
+    fftw_free (m_adComplexFftSignal);
+    fftw_free (m_adComplexFftBackwardOutput);

   }
   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
   }
   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {

-    rfftw_destroy_plan(m_realPlanForward);
-    rfftw_destroy_plan(m_realPlanBackward);
-    delete [] m_adRealFftInput;
-    delete [] m_adRealFftSignal;
+    fftw_destroy_plan(m_realPlanForward);
+    fftw_destroy_plan(m_realPlanBackward);
+    fftw_free (m_adRealFftInput);
+    fftw_free (m_adRealFftOutput);
+    fftw_free (m_adRealFftSignal);
+    fftw_free (m_adRealFftBackwardOutput);

   }
 #endif
 }
   }
 #endif
 }


@@ -595,35 +604,27 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const
     for (i = 0; i < m_nSignalPoints; i++)
       m_adRealFftInput[i] = input[i];
     
     for (i = 0; i < m_nSignalPoints; i++)
       m_adRealFftInput[i] = input[i];
     

-    fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ];
-    rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
+    fftw_execute (m_realPlanForward);

     for (i = 0; i < m_nFilterPoints; i++)
     for (i = 0; i < m_nFilterPoints; i++)

-      m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
-    delete [] fftOutput;
+      m_adRealFftSignal[i] = m_adFilter[i] * m_adRealFftOutput[i];

     for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
             m_adRealFftSignal[i] = 0;
     
     for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
             m_adRealFftSignal[i] = 0;
     

-    fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ];
-    rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
+    fftw_execute (m_realPlanBackward);

     for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
     for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)

-      output[i] = ifftOutput[i];
-    delete [] ifftOutput;
+      output[i] = m_adRealFftBackwardOutput[i];

   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
     for (i = 0; i < m_nSignalPoints; i++)
   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
     for (i = 0; i < m_nSignalPoints; i++)

-      m_adComplexFftInput[i].re = input[i];
+      m_adComplexFftInput[i][0] = input[i];

     
     

-    fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ];
-    fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
+    fftw_execute (m_complexPlanForward);

     for (i = 0; i < m_nFilterPoints; i++) {
     for (i = 0; i < m_nFilterPoints; i++) {

-      m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
-      m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
+      m_adComplexFftSignal[i][0] = m_adFilter[i] * m_adComplexFftOutput[i][0];
+      m_adComplexFftSignal[i][1] = m_adFilter[i] * m_adComplexFftOutput[i][1];

     }
     }

-    delete [] fftOutput;
-    fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ];
-    fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
+    fftw_execute (m_complexPlanBackward);

     for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) 
     for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) 

-      output[i] = ifftOutput[i].re;
-    delete [] ifftOutput;
+      output[i] = m_adComplexFftBackwardOutput[i][0];

   }
 #endif
   delete input;
   }
 #endif
   delete input;