arm_cfft_f64 output not matching Matlab

[adwodon]

I'm relatively new to DSP so forgive me if I've missed something obvious. Currently use pocketfft and I am looking into CMSIS-DSP for a new platform and I can't seem to get the same results, this is with v1.15.0 and v1.16.0

I've created a small C++ project which demonstrates the issue, using -DHOST=ON on a windows machine, it prints out the results from arm_cfft_f64 (as well as pocketfft), I've attached in output.txt, which can be fed into Matlab.

extern "C" {
#include "pocketfft.h"
}
#include "arm_math.h"
#include "arm_const_structs.h"
#include <complex>
#include <cmath>
#include <vector>
#include <fmt/core.h>

int main(int argc, char **) {
    constexpr int kDataSize{256};
    constexpr float64_t kTs = 1.0/kDataSize;
    std::vector<float64_t> times{};

    std::vector<std::complex<float64_t>> pocketfft_data{};
    std::vector<std::complex<float64_t>> arm_data{};
    for(int i = 0; i < kDataSize; i++)
    {
        float64_t t = i * kTs;
        times.push_back(t);
        float64_t real = std::cos(2*M_PI*15*t);
        float64_t imag = 0.0;

        pocketfft_data.emplace_back(real, imag);
        arm_data.emplace_back(real, imag);
    }

    // PocketFFT
    cfft_plan plan = make_cfft_plan(kDataSize);
    cfft_forward(plan, reinterpret_cast<float64_t*>(pocketfft_data.data()), 1.0);
    destroy_cfft_plan(plan);

    fmt::print("\n\npocketFFTResults = [\n");
    for(const auto item : pocketfft_data)
    {
        fmt::print("complex({}, {})\n", item.real(), item.imag());
    }
    fmt::print("];\n");

    // CMSIS-DSP
    arm_cfft_instance_f64 inst{};
    arm_cfft_init_f64(&inst, kDataSize);
    arm_cfft_f64(&inst, reinterpret_cast<float64_t*>(arm_data.data()), 0, 0);

    fmt::print("\n\narmResults = [\n");
    for(const auto item : arm_data)
    {
        fmt::print("complex({}, {})\n", item.real(), item.imag());
    }
    fmt::print("];\n");
}

I've then tried to compare these results in Matlab:

Ts = 1/256;
t = 0:Ts:1-Ts;
S = complex(cos(2*pi*15*t), 0);
Fs = 1/Ts;
n = length(S);
f = (0:n-1)*(Fs/n);
fshift = (-n/2:n/2-1)*(Fs/n);

X = fft(S);
power = abs(X).^2/n;
Y = fftshift(X);
powershift = abs(Y).^2/n;
stem(fshift, powershift, 'DisplayName', 'Matlab');
hold on

% X = pocketFFTResults;
% power = abs(X).^2/n;
% Y = fftshift(X);
% powershift = abs(Y).^2/n;
% stem(fshift, powershift, 'DisplayName', 'PocketFFT');
% hold on

X = armResults;
power = abs(X).^2/n;
Y = fftshift(X);
powershift = abs(Y).^2/n;
stem(fshift, powershift, 'DisplayName', 'ARM DSP');
hold off
title("Double-Sided Amplitude Spectrum of x(t)")
xlabel("Frequency (Hz)")
ylabel("|y|")
grid
legend('-DynamicLegend');

I get the following figure, I've disabled the pocketfft output, but it matches the Matlab:

If I increase kDataSize to 1024 (with arm_cfft_sR_f64_len1024) I get:

Any help would be greatly appreciated.

[christophe0606]

@adwodon You need to set the bit reverse argument to 1 when calling the FFT. It is set to 0 in your example:

arm_cfft_f64(&inst, reinterpret_cast<float64_t*>(arm_data.data()), 0, 1);

This argument is a bit legacy. There are indeed a few cases when you do FFT -> linear processing -> IFFT where you can avoid doing the bit reversing if the linear processing step is done differently.

But it is such an unfrequent case, that most FFT libraries are always doing the bit reversal.

[adwodon]

When I was only comparing the raw output I did try bit reversal, but didn't look at it again when doing the analysis. It's all working as intended now, thankyou for the explanation, very helpful.