summaryrefslogtreecommitdiff
path: root/libs/qm-dsp/dsp/transforms/FFT.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'libs/qm-dsp/dsp/transforms/FFT.cpp')
-rw-r--r--libs/qm-dsp/dsp/transforms/FFT.cpp261
1 files changed, 141 insertions, 120 deletions
diff --git a/libs/qm-dsp/dsp/transforms/FFT.cpp b/libs/qm-dsp/dsp/transforms/FFT.cpp
index 0ca618f40b..da476b8a8b 100644
--- a/libs/qm-dsp/dsp/transforms/FFT.cpp
+++ b/libs/qm-dsp/dsp/transforms/FFT.cpp
@@ -4,178 +4,199 @@
QM DSP Library
Centre for Digital Music, Queen Mary, University of London.
- This file is based on Don Cross's public domain FFT implementation.
*/
#include "FFT.h"
#include "maths/MathUtilities.h"
+#include "kiss_fft.h"
+#include "kiss_fftr.h"
+
#include <cmath>
#include <iostream>
-FFT::FFT(unsigned int n) :
- m_n(n),
- m_private(0)
+#include <stdexcept>
+
+class FFT::D
{
- if( !MathUtilities::isPowerOfTwo(m_n) )
- {
- std::cerr << "ERROR: FFT: Non-power-of-two FFT size "
- << m_n << " not supported in this implementation"
- << std::endl;
- return;
+public:
+ D(int n) : m_n(n) {
+ m_planf = kiss_fft_alloc(m_n, 0, NULL, NULL);
+ m_plani = kiss_fft_alloc(m_n, 1, NULL, NULL);
+ m_kin = new kiss_fft_cpx[m_n];
+ m_kout = new kiss_fft_cpx[m_n];
}
-}
-FFT::~FFT()
-{
+ ~D() {
+ kiss_fft_free(m_planf);
+ kiss_fft_free(m_plani);
+ delete[] m_kin;
+ delete[] m_kout;
+ }
-}
+ void process(bool inverse,
+ const double *ri,
+ const double *ii,
+ double *ro,
+ double *io) {
-FFTReal::FFTReal(unsigned int n) :
- m_n(n),
- m_private_real(0)
-{
- m_private_real = new FFT(m_n);
-}
+ for (int i = 0; i < m_n; ++i) {
+ m_kin[i].r = ri[i];
+ m_kin[i].i = (ii ? ii[i] : 0.0);
+ }
-FFTReal::~FFTReal()
-{
- delete (FFT *)m_private_real;
-}
+ if (!inverse) {
-void
-FFTReal::process(bool inverse,
- const double *realIn,
- double *realOut, double *imagOut)
-{
- ((FFT *)m_private_real)->process(inverse, realIn, 0, realOut, imagOut);
-}
+ kiss_fft(m_planf, m_kin, m_kout);
-static unsigned int numberOfBitsNeeded(unsigned int p_nSamples)
-{
- int i;
+ for (int i = 0; i < m_n; ++i) {
+ ro[i] = m_kout[i].r;
+ io[i] = m_kout[i].i;
+ }
- if( p_nSamples < 2 )
- {
- return 0;
- }
+ } else {
- for ( i=0; ; i++ )
- {
- if( p_nSamples & (1 << i) ) return i;
- }
-}
+ kiss_fft(m_plani, m_kin, m_kout);
-static unsigned int reverseBits(unsigned int p_nIndex, unsigned int p_nBits)
-{
- unsigned int i, rev;
+ double scale = 1.0 / m_n;
- for(i=rev=0; i < p_nBits; i++)
- {
- rev = (rev << 1) | (p_nIndex & 1);
- p_nIndex >>= 1;
+ for (int i = 0; i < m_n; ++i) {
+ ro[i] = m_kout[i].r * scale;
+ io[i] = m_kout[i].i * scale;
+ }
+ }
}
+
+private:
+ int m_n;
+ kiss_fft_cfg m_planf;
+ kiss_fft_cfg m_plani;
+ kiss_fft_cpx *m_kin;
+ kiss_fft_cpx *m_kout;
+};
+
+FFT::FFT(int n) :
+ m_d(new D(n))
+{
+}
- return rev;
+FFT::~FFT()
+{
+ delete m_d;
}
void
-FFT::process(bool p_bInverseTransform,
+FFT::process(bool inverse,
const double *p_lpRealIn, const double *p_lpImagIn,
double *p_lpRealOut, double *p_lpImagOut)
{
- if (!p_lpRealIn || !p_lpRealOut || !p_lpImagOut) return;
-
-// std::cerr << "FFT::process(" << m_n << "," << p_bInverseTransform << ")" << std::endl;
-
- unsigned int NumBits;
- unsigned int i, j, k, n;
- unsigned int BlockSize, BlockEnd;
-
- double angle_numerator = 2.0 * M_PI;
- double tr, ti;
+ m_d->process(inverse,
+ p_lpRealIn, p_lpImagIn,
+ p_lpRealOut, p_lpImagOut);
+}
+
+class FFTReal::D
+{
+public:
+ D(int n) : m_n(n) {
+ if (n % 2) {
+ throw std::invalid_argument
+ ("nsamples must be even in FFTReal constructor");
+ }
+ m_planf = kiss_fftr_alloc(m_n, 0, NULL, NULL);
+ m_plani = kiss_fftr_alloc(m_n, 1, NULL, NULL);
+ m_c = new kiss_fft_cpx[m_n];
+ }
- if( !MathUtilities::isPowerOfTwo(m_n) )
- {
- std::cerr << "ERROR: FFT::process: Non-power-of-two FFT size "
- << m_n << " not supported in this implementation"
- << std::endl;
- return;
+ ~D() {
+ kiss_fftr_free(m_planf);
+ kiss_fftr_free(m_plani);
+ delete[] m_c;
}
- if( p_bInverseTransform ) angle_numerator = -angle_numerator;
+ void forward(const double *ri, double *ro, double *io) {
- NumBits = numberOfBitsNeeded ( m_n );
+ kiss_fftr(m_planf, ri, m_c);
+ for (int i = 0; i <= m_n/2; ++i) {
+ ro[i] = m_c[i].r;
+ io[i] = m_c[i].i;
+ }
- for( i=0; i < m_n; i++ )
- {
- j = reverseBits ( i, NumBits );
- p_lpRealOut[j] = p_lpRealIn[i];
- p_lpImagOut[j] = (p_lpImagIn == 0) ? 0.0 : p_lpImagIn[i];
+ for (int i = 0; i + 1 < m_n/2; ++i) {
+ ro[m_n - i - 1] = ro[i + 1];
+ io[m_n - i - 1] = -io[i + 1];
+ }
}
+ void forwardMagnitude(const double *ri, double *mo) {
- BlockEnd = 1;
- for( BlockSize = 2; BlockSize <= m_n; BlockSize <<= 1 )
- {
- double delta_angle = angle_numerator / (double)BlockSize;
- double sm2 = -sin ( -2 * delta_angle );
- double sm1 = -sin ( -delta_angle );
- double cm2 = cos ( -2 * delta_angle );
- double cm1 = cos ( -delta_angle );
- double w = 2 * cm1;
- double ar[3], ai[3];
-
- for( i=0; i < m_n; i += BlockSize )
- {
+ double *io = new double[m_n];
- ar[2] = cm2;
- ar[1] = cm1;
+ forward(ri, mo, io);
- ai[2] = sm2;
- ai[1] = sm1;
+ for (int i = 0; i < m_n; ++i) {
+ mo[i] = sqrt(mo[i] * mo[i] + io[i] * io[i]);
+ }
- for ( j=i, n=0; n < BlockEnd; j++, n++ )
- {
+ delete[] io;
+ }
- ar[0] = w*ar[1] - ar[2];
- ar[2] = ar[1];
- ar[1] = ar[0];
+ void inverse(const double *ri, const double *ii, double *ro) {
- ai[0] = w*ai[1] - ai[2];
- ai[2] = ai[1];
- ai[1] = ai[0];
+ // kiss_fftr.h says
+ // "input freqdata has nfft/2+1 complex points"
- k = j + BlockEnd;
- tr = ar[0]*p_lpRealOut[k] - ai[0]*p_lpImagOut[k];
- ti = ar[0]*p_lpImagOut[k] + ai[0]*p_lpRealOut[k];
+ for (int i = 0; i < m_n/2 + 1; ++i) {
+ m_c[i].r = ri[i];
+ m_c[i].i = ii[i];
+ }
+
+ kiss_fftri(m_plani, m_c, ro);
- p_lpRealOut[k] = p_lpRealOut[j] - tr;
- p_lpImagOut[k] = p_lpImagOut[j] - ti;
+ double scale = 1.0 / m_n;
- p_lpRealOut[j] += tr;
- p_lpImagOut[j] += ti;
+ for (int i = 0; i < m_n; ++i) {
+ ro[i] *= scale;
+ }
+ }
- }
- }
+private:
+ int m_n;
+ kiss_fftr_cfg m_planf;
+ kiss_fftr_cfg m_plani;
+ kiss_fft_cpx *m_c;
+};
- BlockEnd = BlockSize;
+FFTReal::FFTReal(int n) :
+ m_d(new D(n))
+{
+}
- }
+FFTReal::~FFTReal()
+{
+ delete m_d;
+}
+void
+FFTReal::forward(const double *ri, double *ro, double *io)
+{
+ m_d->forward(ri, ro, io);
+}
- if( p_bInverseTransform )
- {
- double denom = (double)m_n;
+void
+FFTReal::forwardMagnitude(const double *ri, double *mo)
+{
+ m_d->forwardMagnitude(ri, mo);
+}
- for ( i=0; i < m_n; i++ )
- {
- p_lpRealOut[i] /= denom;
- p_lpImagOut[i] /= denom;
- }
- }
+void
+FFTReal::inverse(const double *ri, const double *ii, double *ro)
+{
+ m_d->inverse(ri, ii, ro);
}
+
+
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-22 16:57:38 +0000