Possible errors in fft implementation

[inkydragon]

The output of fft_forward and fft_inverse seems inconsistent with other reference implementations.
I'm not sure if this is a bug.

But it might have an effect on the precision of the final results.

CPP test code

Details

Add those code to the end of https://github.com/Mysticial/Mini-Pi/blob/master/decomposed/src/FFT.cpp

Compile and Run: g++ -o fft FFT.cpp FFT.h && ./fft

#include <iostream>
#include <vector>
#include <complex>

std::vector<std::complex<double>> generateComplexVector(int n) {
    std::vector<std::complex<double>> result;
    for (int i = 1; i <= n; ++i) {
        result.push_back(std::complex<double>(i, 0));
    }
    return result;
}
void print_vec_cf64(std::vector<std::complex<double>> vec) {
    for (const auto& num : vec) {
        // printf("%.16g+%.16gim, ", num.real(), num.imag());
        printf("%.4g+%.4gim, ", num.real(), num.imag());
    }
    puts("");
}

int main() {
    auto k = 2;
    
    for (int k=1; k <= 3; k++) {
            auto len = 0x1 << k;
        auto Ta = generateComplexVector(len);

        Mini_Pi::ensure_FFT_tables(k);

        printf("k=%d;  len=%d;\n", k, len);

        puts("fft_forward(1..len):");
        printf("    in:   "); print_vec_cf64(Ta);
            Mini_Pi::fft_forward(Ta.data(), k);
        printf("    out:  "); print_vec_cf64(Ta);

        puts("fft_inverse(fft_forward(1..len))");
        printf("    in:   "); print_vec_cf64(Ta);
            Mini_Pi::fft_inverse(Ta.data(), k);
        printf("    out:  "); print_vec_cf64(Ta);

        puts("fft_inverse(1..len):");
        Ta = generateComplexVector(len);
        printf("    in:   "); print_vec_cf64(Ta);
            Mini_Pi::fft_inverse(Ta.data(), k);
        printf("    out:  "); print_vec_cf64(Ta);
        
        puts("");
    }

    return 0;
}

Test code output:

$ g++ -o fft FFT.cpp FFT.h && ./fft
k=1;  len=2;
fft_forward(1..len):
    in:   1+0im, 2+0im,
    out:  3+0im, -1+0im,
fft_inverse(fft_forward(1..len))
    in:   3+0im, -1+0im,
    out:  2+0im, 4+0im,
fft_inverse(1..len):
    in:   1+0im, 2+0im,
    out:  3+0im, -1+0im,

k=2;  len=4;
fft_forward(1..len):
    in:   1+0im, 2+0im, 3+0im, 4+0im,
    out:  10+0im, -2+0im, -2+-2im, -2+2im,
fft_inverse(fft_forward(1..len))
    in:   10+0im, -2+0im, -2+-2im, -2+2im,
    out:  4+0im, 8+-2.288e-17im, 12+0im, 16+2.288e-17im,
fft_inverse(1..len):
    in:   1+0im, 2+0im, 3+0im, 4+0im,
    out:  10+0im, -1+1im, -4+0im, -1+-1im,

k=3;  len=8;
fft_forward(1..len):
    in:   1+0im, 2+0im, 3+0im, 4+0im, 5+0im, 6+0im, 7+0im, 8+0im,
    out:  36+0im, -4+0im, -4+-4im, -4+4im, -4+-9.657im, -4+1.657im, -4+-1.657im, -4+9.657im,
fft_inverse(fft_forward(1..len))
    in:   36+0im, -4+0im, -4+-4im, -4+4im, -4+-9.657im, -4+1.657im, -4+-1.657im, -4+9.657im,
    out:  8+0im, 16+-1.822e-15im, 24+7.966e-16im, 32+-1.731e-15im, 40+0im, 48+1.731e-15im, 56+-7.966e-16im, 64+1.822e-15im,
fft_inverse(1..len):
    in:   1+0im, 2+0im, 3+0im, 4+0im, 5+0im, 6+0im, 7+0im, 8+0im,
    out:  36+0im, -1+2.414im, -4+4im, -1+0.4142im, -16+0im, -1+-0.4142im, -4+-4im, -1+-2.414im,

issues

1. fft_inverse(fft_forward(1..len)) != 1..len
   Simply take the out * 1/len to get the correct result.
   This may be due to the lack of the necessary 1/n factor in ifft (fft_inverse).

2. fft_inverse(1..2) != [1.5 + 0.0im, -0.5 + 0.0im]
   Same as 1.

3. fft_forward give correct answer only when len=2
   The output of fft_forward is in bit-reversed order.
   But the positions of the first and last elements of the array should remain unchanged.
   So the result is still wrong.

• numpy ref output

import numpy
for i in range(1,4):
    arr_len = 1<<i
    print(numpy.fft.fft(range(1,arr_len+1)))

[ 3.+0.j -1.+0.j]
[10.+0.j -2.+2.j -2.+0.j -2.-2.j]
[36.+0.j         -4.+9.65685425j -4.+4.j         -4.+1.65685425j
 -4.+0.j         -4.-1.65685425j -4.-4.j         -4.-9.65685425j]

• matlab bitrevorder position map

>> bitrevorder((1:2))
ans =
     1     2

>> bitrevorder((1:4))
ans =
     1     3     2     4

>> bitrevorder((1:8))
ans =
     1     5     3     7     2     6     4     8

[Mysticial]

The behavior here is expected. The order of the frequency domain doesn't matter when performing convolution. And where the scaling is done doesn't matter as long as it's done somewhere. Reference implementations either don't scale at all or they scale by sqrt(length) in both the forward and inverse transforms. (which is not what we want here anyway)

The reason why it appears worse than bit-reversed is because the twiddle factors are not inverted like in reference implementations. If you take the complex conjugate of the frequency domain, it should be bit-reversed from a reference implementation. (aside from scaling factors)