X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprocsignal.cpp;fp=libctsim%2Fprocsignal.cpp;h=2206ae1838cf41ebd3f2e13eb34353fdb060670e;hp=bdb9fa00626813ca37509601bd7eb7698d10c470;hb=d42d3d062dd1aca92b5a2552a1f474aab0bee610;hpb=c64771ab78d25f0cbfbaac2456c8790c786a7a68

diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp
index bdb9fa0..2206ae1 100644
--- 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) {
-    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) {
-    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++) 
-      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++) 
-      m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
+      m_adComplexFftSignal[i][0] = m_adComplexFftSignal[i][1] = 0;
   }
 #endif
   
@@ -448,14 +453,18 @@ ProcessSignal::~ProcessSignal (void)
   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) {
-    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
 }
@@ -595,35 +604,27 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const
     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++)
-      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;
     
-    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++)
-      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++)
-      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++) {
-      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++) 
-      output[i] = ifftOutput[i].re;
-    delete [] ifftOutput;
+      output[i] = m_adComplexFftBackwardOutput[i][0];
   }
 #endif
   delete input;