Question Regarding Twiddle Factor in Catch22 Function

[ceap29]

I reached out to Ben via email, and he suggested that I open an issue regarding the following.

My question is: shouldn’t the tmp section inside the twiddle function use the formula 0.0 - 2.0 * PI * i / size to calculate a twiddle factor, instead of 0.0 - PI * i / size?

Here’s the function I’m referring to:
https://github.com/DynamicsAndNeuralSystems/catch22/blob/main/C/fft.c

void twiddles(cplx a[], int size)
{
    double PI = 3.14159265359;

    for (int i = 0; i < size; i++) {
        // cplx tmp = { 0, -PI * i / size };
        #if defined(__GNUC__) || defined(__GNUG__)
   cplx tmp = 0.0  - PI * i / size * I;
        #elif defined(_MSC_VER)
   cplx tmp = {0.0, -PI * i / size };
        #endif
        a[i] = cexp(tmp);
        //a[i] = cexp(-I * M_PI * i / size);
    }
}

Thank you in advance for your clarification. I really appreciate it. Also, thanks for creating such an amazing package; it has been incredibly helpful.

[brendanjohnharris]

It looks like the indexing here only uses every second twiddle, so that seems to recover the conventional 2pi formula:

catch22/C/fft.c

Lines 46 to 51 in d39af48

	for (int i = 0; i < size; i += 2 * step) {
	//cplx t = tw[i] * out[i + step];
	cplx t = _Cmulcc(tw[i], out[i + step]);
	a[i / 2] = _Caddcc(out[i], t);
	a[(i + size) / 2] = _Cminuscc(out[i], t);
	}

I don't use C so maybe @hendersontrent can double check my reading of that section?

In any case, the output fft coefficients are correct; see the quick test below:

using Catch22
using FFTW
using LinearAlgebra
using CairoMakie
using Statistics
using Libdl

# Create test signal
N = 64
t = 0:(N - 1)
f1 = 5
signal = randn(64)

begin # * Now compare the raw ffts
    y = Vector{ComplexF64}(signal)
    twiddles = similar(y)
    ccall(dlsym(lib, :twiddles), Ptr{ComplexF64},
          (Ptr{ComplexF64}, Cint),
          twiddles, length(twiddles)) # In-place

    ccall(dlsym(lib, :fft), Cvoid,
          (Ptr{ComplexF64}, Cint, Ptr{ComplexF64}),
          y, length(y), twiddles)
    catch22_fft = deepcopy(y) # fft modifies signal in-place

    standard_fft = fft(signal)

    fig = Figure()
    ax = PolarAxis(fig[1, 1],
                   title = "FFT Comparison: Standard vs Catch22")
    p1 = scatter!(ax, angle.(standard_fft), (log ∘ norm).(standard_fft), marker = :x,
                  color = :red, markersize = 15)
    p2 = scatter!(ax, angle.(catch22_fft), (log ∘ norm).(catch22_fft),
                  color = :green, marker = :+, markersize = 15)
    Legend(fig[1, 2],
           [MarkerElement(color = :red, marker = :x, markersize = 15),
               MarkerElement(color = :green, marker = :+, markersize = 15)], ["Standard FFT", "Catch22 FFT"])
    display(fig)
end

begin # * Then double-check autocorrelations
    lib = dlopen(Catch22.libcatch22)
    ac_ptr = ccall(dlsym(lib, :co_autocorrs), Ptr{Cdouble},
                   (Ptr{Cdouble}, Cint), signal, length(signal))
    catch22_acf = unsafe_wrap(Array, ac_ptr, length(signal); own = false)

    # Standard autocorrelation using FFT
    function autocorr_fft(x)
        N = length(x)
        x_centered = x .- mean(x)
        nfft = nextpow(2, 2 * N)
        X = fft(vcat(x_centered, zeros(nfft - N)))
        acf = real(ifft(abs2.(X)))
        return acf[1:N] / acf[1]  # Normalize
    end

    acf_standard = autocorr_fft(signal)

    begin # Create plot
        fig = Figure()
        ax = Axis(fig[1, 1],
                  title = "Autocorrelation Comparison: Standard vs π-twiddle vs Catch22",
                  xlabel = "Lag",
                  ylabel = "Autocorrelation")

        lags = 0:20
        scatterlines!(ax, lags, acf_standard[1:21],
                      label = "Standard FFT",
                      linewidth = 2, color = :red, markersize = 20, marker = :x)

        scatter!(ax, lags, catch22_acf[1:21],
                 label = "Catch22 (actual)",
                 marker = :+,
                 markersize = 20,
                 color = :green)

        axislegend(ax, position = :rt)

        display(fig)
    end
end