A few helper functions

src/dsp/FFT.cpp

@@ -1851,7 +1851,7 @@ public:
                 dbuf[i] = realIn[i];
             }
         fftw_execute(m_dplanf);
-        v_convert(complexOut, (fft_double_type *)m_dpacked, sz + 2);
+        v_convert(complexOut, (const fft_double_type *)m_dpacked, sz + 2);
     }
 
     void forwardPolar(const double *R__ realIn, double *R__ magOut, double *R__ phaseOut) {
@@ -1866,7 +1866,7 @@ public:
             }
         fftw_execute(m_dplanf);
         v_cartesian_interleaved_to_polar
-            (magOut, phaseOut, (fft_double_type *)m_dpacked, m_size/2+1);
+            (magOut, phaseOut, (const fft_double_type *)m_dpacked, m_size/2+1);
     }
 
     void forwardMagnitude(const double *R__ realIn, double *R__ magOut) {
@@ -1880,11 +1880,8 @@ public:
                 dbuf[i] = realIn[i];
             }
         fftw_execute(m_dplanf);
-        const int hs = m_size/2;
+        v_cartesian_interleaved_to_magnitudes
-        for (int i = 0; i <= hs; ++i) {
+            (magOut, (const fft_double_type *)m_dpacked, m_size/2+1);
-            magOut[i] = sqrt(m_dpacked[i][0] * m_dpacked[i][0] +
-                             m_dpacked[i][1] * m_dpacked[i][1]);
-        }
     }
 
     void forward(const float *R__ realIn, float *R__ realOut, float *R__ imagOut) {
@@ -1912,7 +1909,7 @@ public:
                 fbuf[i] = realIn[i];
             }
         fftwf_execute(m_fplanf);
-        v_convert(complexOut, (fft_float_type *)m_fpacked, sz + 2);
+        v_convert(complexOut, (const fft_float_type *)m_fpacked, sz + 2);
     }
 
     void forwardPolar(const float *R__ realIn, float *R__ magOut, float *R__ phaseOut) {
@@ -1941,11 +1938,8 @@ public:
                 fbuf[i] = realIn[i];
             }
         fftwf_execute(m_fplanf);
-        const int hs = m_size/2;
+        v_cartesian_interleaved_to_magnitudes
-        for (int i = 0; i <= hs; ++i) {
+            (magOut, (const fft_float_type *)m_fpacked, m_size/2+1);
-            magOut[i] = sqrtf(m_fpacked[i][0] * m_fpacked[i][0] +
-                              m_fpacked[i][1] * m_fpacked[i][1]);
-        }
     }
 
     void inverse(const double *R__ realIn, const double *R__ imagIn, double *R__ realOut) {
@@ -1978,14 +1972,8 @@ public:
 
     void inversePolar(const double *R__ magIn, const double *R__ phaseIn, double *R__ realOut) {
         if (!m_dplanf) initDouble();
-        const int hs = m_size/2;
+        v_polar_to_cartesian_interleaved
-        fftw_complex *const R__ dpacked = m_dpacked;
+            ((fft_double_type *)m_dpacked, magIn, phaseIn, m_size/2+1);
-        for (int i = 0; i <= hs; ++i) {
-            dpacked[i][0] = magIn[i] * cos(phaseIn[i]);
-        }
-        for (int i = 0; i <= hs; ++i) {
-            dpacked[i][1] = magIn[i] * sin(phaseIn[i]);
-        }
         fftw_execute(m_dplani);
         const int sz = m_size;
         fft_double_type *const R__ dbuf = m_dbuf;
@@ -2048,14 +2036,8 @@ public:
 
     void inversePolar(const float *R__ magIn, const float *R__ phaseIn, float *R__ realOut) {
         if (!m_fplanf) initFloat();
-        const int hs = m_size/2;
+        v_polar_to_cartesian_interleaved
-        fftwf_complex *const R__ fpacked = m_fpacked;
+            ((fft_float_type *)m_fpacked, magIn, phaseIn, m_size/2+1);
-        for (int i = 0; i <= hs; ++i) {
-            fpacked[i][0] = magIn[i] * cosf(phaseIn[i]);
-        }
-        for (int i = 0; i <= hs; ++i) {
-            fpacked[i][1] = magIn[i] * sinf(phaseIn[i]);
-        }
         fftwf_execute(m_fplani);
         const int sz = m_size;
         fft_float_type *const R__ fbuf = m_fbuf;

src/system/VectorOpsComplex.h

@@ -246,6 +246,63 @@ void v_cartesian_to_polar_interleaved_inplace(T *const R__ srcdst,
     }
 }
 
+template<typename S, typename T> // S source, T target
+void v_cartesian_to_magnitudes(T *const R__ mag,
+                               const S *const R__ real,
+                               const S *const R__ imag,
+                               const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        mag[i] = T(sqrt(real[i] * real[i] + imag[i] * imag[i]));
+    }
+}
+
+template<typename S, typename T> // S source, T target
+void v_cartesian_interleaved_to_magnitudes(T *const R__ mag,
+                                           const S *const R__ src,
+                                           const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        mag[i] = T(sqrt(src[i*2] * src[i*2] + src[i*2+1] * src[i*2+1]));
+    }
+}
+
+#ifdef HAVE_IPP
+template<>
+inline void v_cartesian_to_magnitudes(float *const R__ mag,
+                                      const float *const R__ real,
+                                      const float *const R__ imag,
+                                      const int count)
+{
+    ippsMagnitude_32f(real, imag, mag, count);
+}
+
+template<>
+inline void v_cartesian_to_magnitudes(double *const R__ mag,
+                                      const double *const R__ real,
+                                      const double *const R__ imag,
+                                      const int count)
+{
+    ippsMagnitude_64f(real, imag, mag, count);
+}
+
+template<>
+inline void v_cartesian_interleaved_to_magnitudes(float *const R__ mag,
+                                                  const float *const R__ src,
+                                                  const int count)
+{
+    ippsMagnitude_32fc((const Ipp32fc *)src, mag, count);
+}
+
+template<>
+inline void v_cartesian_interleaved_to_magnitudes(double *const R__ mag,
+                                                  const double *const R__ src,
+                                                  const int count)
+{
+    ippsMagnitude_64fc((const Ipp64fc *)src, mag, count);
+}
+#endif
+
 }
 
 #endif