From fcc952e81a7b607c5e949f0d74498fc11a6f8a9a Mon Sep 17 00:00:00 2001 From: Chris Cannam Date: Tue, 9 Mar 2021 17:29:21 +0000 Subject: [PATCH] A few helper functions --- src/dsp/FFT.cpp | 40 +++++++----------------- src/system/VectorOpsComplex.h | 57 +++++++++++++++++++++++++++++++++++ 2 files changed, 68 insertions(+), 29 deletions(-) diff --git a/src/dsp/FFT.cpp b/src/dsp/FFT.cpp index f6f0d97..eb5c8f6 100644 --- a/src/dsp/FFT.cpp +++ b/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; - for (int i = 0; i <= hs; ++i) { - magOut[i] = sqrt(m_dpacked[i][0] * m_dpacked[i][0] + - m_dpacked[i][1] * m_dpacked[i][1]); - } + v_cartesian_interleaved_to_magnitudes + (magOut, (const fft_double_type *)m_dpacked, m_size/2+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; - for (int i = 0; i <= hs; ++i) { - magOut[i] = sqrtf(m_fpacked[i][0] * m_fpacked[i][0] + - m_fpacked[i][1] * m_fpacked[i][1]); - } + v_cartesian_interleaved_to_magnitudes + (magOut, (const fft_float_type *)m_fpacked, m_size/2+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; - fftw_complex *const R__ dpacked = m_dpacked; - 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]); - } + v_polar_to_cartesian_interleaved + ((fft_double_type *)m_dpacked, magIn, phaseIn, m_size/2+1); 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; - fftwf_complex *const R__ fpacked = m_fpacked; - 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]); - } + v_polar_to_cartesian_interleaved + ((fft_float_type *)m_fpacked, magIn, phaseIn, m_size/2+1); fftwf_execute(m_fplani); const int sz = m_size; fft_float_type *const R__ fbuf = m_fbuf; diff --git a/src/system/VectorOpsComplex.h b/src/system/VectorOpsComplex.h index 8d4161b..aca2a37 100644 --- a/src/system/VectorOpsComplex.h +++ b/src/system/VectorOpsComplex.h @@ -246,6 +246,63 @@ void v_cartesian_to_polar_interleaved_inplace(T *const R__ srcdst, } } +template // 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 // 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