Unexpected results with non-power of 2 real DFT

[clydeclements]

I am getting unexpected results for non-power of 2 DFT. I am using KFR version 3.0.9 with Clang compiler version 6.0.0 and testing on Ubuntu 18.04.

My test code is as follows:

#include <kfr/base.hpp>
#include <kfr/dft.hpp>
#include <kfr/dsp.hpp>

using namespace kfr;

int main() {
    println(library_version());
    constexpr size_t N = 63;
    const double fs = 1000.;
    const double ts = 1 / fs;
    univector<double, N> t = linspace(0., N * ts, N);
    univector<double, N> s = 0.7 * sin(2 * constants<double>::pi * 50 * t) + sin(2 * constants<double>::pi * 120 * t);
    univector<complex<double>> S(N);
    const dft_plan_real<double> dft(N);
    univector<u8> temp(dft.temp_size);
    dft.execute(S, s, temp);
    println("Sum of signal values: ", sum(s));
    println("DC component of DFT : ", S[0].real());
    return 0;
}

When I compile and run the above code, I get the following output:

Sum of signal values: 3.24885
DC component of DFT : 2.46928

I expected the first component of the DFT to be the same as the sum of the signal values.

For comparison purposes, I created a corresponding Python test script (Python 3.7.3 with NumPy 1.16.2):

import numpy as np
fs = 1000.
ts = 1 / fs
N = 63
t = np.linspace(0., N * ts, N, endpoint=False)
s = 0.7 * np.sin(2 * np.pi * 50. * t) + np.sin(2 * np.pi * 120. * t)
S = np.fft.rfft(s)
print("Sum of signal values: ", np.sum(s))
print("DC component of DFT : ", S[0].real)

Output from the Python script is:

Sum of signal values:  3.248853363223806
DC component of DFT :  3.248853363223807

This is what I expected. Note also that the sum of the signal values in the Python case agrees with the KFR case. However, the first component of the DFT differs (and in fact, all values of the DFT differ).

If I change the value of N from 63 to 64 and rerun, I get this result from KFR:

Sum of signal values: 3.44704
DC component of DFT : 3.44704

This agrees with the result from Python:

Sum of signal values:  3.447040706601597
DC component of DFT :  3.447040706601598

So power of 2 size DFT works as expected, but non-power of 2 size DFT does not.

Am I calling the DFT functions in KFR incorrectly?

[dancasarin]

Hi,

This is caused by the way real DFT is calculated. Pair of real values are interpreted as complex for high performance, so there is limitation for real DFT size, it must be even.

The workaround is to use complex DFT instead of real one.

#include <kfr/base.hpp>
#include <kfr/dft.hpp>
#include <kfr/dsp.hpp>

using namespace kfr;

int main() {
    println(library_version());
    constexpr size_t N = 63;
    const double fs = 1000.;
    const double ts = 1 / fs;
    univector<double, N> t = linspace(0., N * ts, N);
    univector<double, N> s = 0.7 * sin(2 * constants<double>::pi * 50 * t) + sin(2 * constants<double>::pi * 120 * t);
    univector<complex<double>> S(N);
    const dft_plan<double> dft(N); // dft_plan instead of dft_plan_real
    univector<u8> temp(dft.temp_size);
    dft.execute(S, univector<complex<double>>(s), temp); // Explicit cast to complex vector
    println("Sum of signal values: ", sum(s));
    println("DC component of DFT : ", S[0].real());
    return 0;
}

Tested on Ubuntu 18.04, Clang 6.0.0, KFR 3.0.9.

KFR 3.0.9 optimized avx2 64-bit (clang-6.0.0/linux) +in +ve
Sum of signal values: 3.24885
DC component of DFT : 3.24885

KFR definitely should show warning for odd sizes for real DFT or use the above workaround internally. This will be fixed in one of future versions.

[clydeclements]

Great, that works for me. And thanks for the prompt reply!

Just a suggestion: maybe add a note to the documentation until a warning or fix is implemented.