hfft.h

#ifndef __HFFT_H
#define __HFFT_H
 
#include <complex>
#include <valarray>
#include <complex.h>
 
#include "hwindow.h"
 
template <class T>
class HFft {
 
    private:
 
        HWindow<T>* _window;
        int _size;
 
        T* _fftBuffer;
        std::complex<double>* _hilbertBuffer;
        std::complex<T>* _iqComplex;
        std::complex<double>* _iqSpectrum;
 
        void FFT(std::valarray<std::complex<double>>& x)
        {
            const size_t N = x.size();
            if (N <= 1) return;
 
            // divide
            std::valarray<std::complex<double>> even = x[std::slice(0, N/2, 2)];
            std::valarray<std::complex<double>>  odd = x[std::slice(1, N/2, 2)];
 
            // conquer
            FFT(even);
            FFT(odd);
 
            // combine
            for (size_t k = 0; k < N/2; ++k)
            {
                std::complex<double> t = std::polar(1.0, -2 * M_PI * k / N) * odd[k];
                x[k    ] = even[k] + t;
                x[k+N/2] = even[k] - t;
            }
        }
 
    public:
 
        HFft(int size, HWindow<T>* window = nullptr):
            _window(window),
            _size(size) {
 
            // Set window size
            if( _window != nullptr ) {
                _window->SetSize(_size);
            }
 
            // Allocate a buffer for intermediate results in the fft functions
            _fftBuffer = new T[_size];
 
            // Allocate a buffer for intermediate results in the hilbert function
            _hilbertBuffer = new std::complex<double>[_size];
 
            // Allocate buffers for IQ2REAl operations
            _iqComplex = new std::complex<T>[_size];
            _iqSpectrum = new std::complex<double>[_size];
        }
 
        ~HFft() {
            delete _fftBuffer;
            delete _hilbertBuffer;
            delete _iqComplex;
            delete _iqSpectrum;
        }
 
        void FFT(std::complex<T>* src, std::complex<double>* result) {
 
            // Prepare a set of complex values
            std::valarray<std::complex<double>> x(_size);
            for( int i = 0; i < _size ; i++ )
            {
                x[i] = std::complex<double>(src[i].real(), src[i].imag());
            }
 
            // Calculate the FFT
            FFT(x);
 
            // Convert to simple array with bins for positive AND negative frequencies
            for( int i = 0; i < _size; i++ )
            {
                result[i] = x[i];
            }
        }
 
        void FFT(T* src, std::complex<double>* result, bool window = true) {
 
            // Run the input buffer through a window function, if any is given
            if( window && _window != nullptr ) {
                _window->Apply(src, _fftBuffer, _size);
            } else {
                memcpy((void*) _fftBuffer, (void*) src, _size * sizeof(T));
            }
 
            // Prepare a set of complex values
            std::valarray<std::complex<double>> x(_size);
            for( int i = 0; i < _size ; i++ )
            {
                x[i] = std::complex<double>(_fftBuffer[i], 0);
            }
 
            // Calculate the FFT
            FFT(x);
 
            // Convert to simple array with bins for positive AND negative frequencies
            for( int i = 0; i < _size; i++ )
            {
                result[i] = x[i];
            }
        }
 
        void FFT(T* src, double* spectrum, double* phase = nullptr) {
 
            // Allocate a buffer for the complex spectrum results
            std::complex<double>* result = new std::complex<double>[_size];
 
            // Run fft on the input
            FFT(src, result);
 
            // convert to real valued spectrum powers for each bin
            for( int i = 0; i < _size / 2; i++ )
            {
                // Get absolute value at point n
                spectrum[i] = std::abs(result[i]);
 
                if( phase != nullptr )
                {
                    phase[i] = std::arg(result[i]);
                }
            }
 
            // Cleanup
            delete[] result;
        }
 
        void IFFT(std::complex<double>* src, T* result) {
 
            // Prepare the input array to the fft function, filled with the complex conjugate of the input values
            std::valarray<std::complex<double>> x(_size);
            for( int i = 0; i < _size ; i++ )
            {
                x[i] = src[i];
            }
            x.apply(std::conj);
 
            // Calculate the FFT
            FFT(x);
 
            // Convert to simple array with the realvalued signal samples
            double factor = 1.0 / _size;
            x.apply(std::conj);
            for( int i = 0; i < _size; i++ )
            {
                result[i] = x[i].real() * factor;
            }
        }
 
        void HILBERT(T* src, T* dest) {
            FFT(src, _hilbertBuffer, false);
 
            for(int i = _size / 2; i < _size; i++ ) {
                _hilbertBuffer[i] = 0;
            }
 
            IFFT(_hilbertBuffer, dest);
        }
 
        void SPECTRUM(T* response, int length, std::complex<double>* spectrum) {
 
            // Sanity check
            if( length > _size / 2 ) {
                throw new HInitializationException("Length must not exceed fft-size/2 in call to HFft::SPECTRUM");
            }
 
            // Initialize buffer
            T* buffer = new T[_size];
            memset((void*) buffer, 0, sizeof(T) * _size);
 
            // Copy frequency response into buffer
            memcpy((void*) buffer, (void*) response, sizeof(T) * length);
 
            // Take FFT of the frequency response
            FFT(buffer, spectrum);
        }
 
        void IQ2REAL(T *multiplexed, T* real) {
 
            // Demultiplex into complex signal
            for( int i = 0; i < _size * 2; i += 2 ) {
                _iqComplex[i / 2] = {multiplexed[i], multiplexed[i + 1]};
            }
 
            // Convert the complex IQ buffer to real valued samples
            IQ2REAL(_iqComplex, real);
        }
 
        void IQ2REAL(std::complex<T>* iq, T* real) {
 
            // FFT -> symmetrical spectrum with positive and negative frequencies
            FFT(iq, _iqSpectrum);
 
            // Zero negative frequencies
            for( int i = _size / 2; i < _size; i++ ) {
                _iqSpectrum[i] = 0;
            }
 
            // IFFT -> real signal with negative frequencies removed
            IFFT(_iqSpectrum, real);
        }
};
 
#endif