fft.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_FFT_H
#define THEORETICA_FFT_H
 
#include "../core/bits.h"
#include "../algebra/algebra_types.h"
#include "../algebra/algebra.h"
#include "../complex/complex.h"
 
 
namespace theoretica {
 
 
    namespace signal {
 
 
        template<typename ReturnVector = cvec, typename InputVector = cvec>
        inline ReturnVector fft(const InputVector& x, bool inverse = false) {
 
            if (x.size() == 0) {
                TH_MATH_ERROR("fft", x.size(), MathError::InvalidArgument);
                return make_error<ReturnVector>(1);
            }
 
            // Resulting vector in the frequency domain
            ReturnVector k = x;
            const unsigned int N = x.size();
            const real sign = (inverse ? 1.0 : -1.0);
 
            // Compute the logarithm of the size
            const unsigned int log2N = ilog2(N);
 
            // Enforce power of 2 vector size
            if (N != (unsigned int) (1 << log2N)) {
                return algebra::vec_error(k);
            }
 
            // Bit reverse
            bits::swap_bit_reverse(k, log2N);
 
            for (unsigned int p = 1; p <= log2N; p++) {
 
                const unsigned int m = 1 << p;
                const unsigned int offset = m / 2;
                
                complex<real> w (1.0, 0.0);
 
                // Phase shift between iterations
                const complex<real> phase = complex<real>(
                    cos(sign * 2 * PI / m),
                    sin(sign * 2 * PI / m)
                );
 
                for (unsigned int j = 0; j < offset; j++) {
 
                    for (unsigned int i = j; i < N; i += m) {
 
                        const complex<real> t = w * k[i + offset];
                        k[i + offset] = k[i] - t;
                        k[i] += t;
                    }
                    
                    w *= phase;
                }
            }
 
            // The normalization constant is 1/N
            if (inverse)
                k /= N;
 
            return k;
        }
 
 
        template<typename ReturnVector = cvec, typename InputVector = cvec>
        inline ReturnVector ifft(const InputVector& k) {
            return fft(k, true);
        }
    }
}
 
 
#endif