How can I port code that uses numpy.fft.rfft from python to C++?

Question

I have code written in python. It computes positive part of FFT of real input using numpy. I need to port this code to C++.

import numpy as np
interp=[131.107, 133.089, 132.199, 129.905, 132.977]
res=np.fft.rfft(interp)
print res

Result of rfft is [ 659.27700000+0.j, 1.27932533-1.4548977j, -3.15032533+2.1158917j]

I tried to use OpenCV for 1D DFT:

std::vector<double> fft;
std::vector<double> interpolated = {131.107, 133.089, 132.199, 129.905, 132.977};
cv::dft( interpolated, fft );
for( auto it = fft.begin(); it != fft.end(); ++it ) {
    std::cout << *it << ' ';
}
std::cout << std::endl;

Result of cv::dft is {1.42109e-14, -127.718, -94.705, 6.26856, 23.0231}. It is much different from numpy.fft.rfft. It looks strange that DC value (zero element) is near zero on all inputs after OpenCV's dft computed.

Usage of FFTW3 library gave me the same results as OpenCV's results:

std::vector<double> interpolated = {131.107, 133.089, 132.199, 129.905, 132.977};
fftw_complex* out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * 3 );
fftw_plan plan = fftw_plan_dft_r2c_1d( interpolated.size( ), interpolated.data( ), out, FFTW_ESTIMATE );
fftw_execute(plan);
fftw_destroy_plan(plan);
for( size_t i = 0; i < interpolated.size( ); ++i ) {
    std::cout << " (" << out[ i ][ 0 ] << ", " << out[ i ][ 1 ] << ")";
}
fftw_free(out);

This code gives me the same results as OpenCV. It prints: (1.42109e-14, 0) (-127.718, -94.705) (6.26856, 23.0231).

Why do I get different results of dft in C++ and in python? What am I doing wrong?

Thanks!

Answer 1

I'm using gcc 4.6 at the moment, which doesn't have C++11, so I tried this version of your code, using OpenCV 2.4.8:

#include <iostream>
#include "opencv2/core/core.hpp"

int main(int argc, char *argv[])
{
    double data[] = {131.107, 133.089, 132.199, 129.905, 132.977};
    std::vector<double> interpolated (data, data + sizeof(data) / sizeof(double));
    std::vector<double> fft;

    cv::dft(interpolated, fft);

    for (std::vector<double>::const_iterator it = fft.begin(); it != fft.end(); ++it) {
        std::cout << *it << ' ';
    }
    std::cout << std::endl;
}

The output

659.277 1.27933 -1.4549 -3.15033 2.11589

agrees with numpy and with the cv2 python module:

In [55]: np.set_printoptions(precision=3)

In [56]: x
Out[56]: array([ 131.107,  133.089,  132.199,  129.905,  132.977])

In [57]: np.fft.rfft(x)
Out[57]: array([ 659.277+0.j   ,    1.279-1.455j,   -3.150+2.116j])

In [58]: cv2.dft(x)
Out[58]: 
array([[ 659.277],
       [   1.279],
       [  -1.455],
       [  -3.15 ],
       [   2.116]])

I don't know why your code is not working, so I guess this is more of a long comment than an answer.

Answer 2

Please check the documentation. The libfft rfft method transforms a vector of real inputs into the complex Fourier coefficients. Using the conjugacy of Fourier coefficients for real signals, the output can be given in an array of the same length as the input.

The generic fft and dft methods transform vectors of complex numbers into vectors of complex coefficients. The older codes use arrays of doubles for input and output where the real and imaginary parts of the complex numbers are given alternatingly, i.e., one array of even length. What happens to odd input lengths may be undocumented behavior.