Librosa like Stft in C using FFTW

Question

I am trying to write a Stft routine in C using the FFTW library that results in the same output of Librosa's Stft in Python. Unfortunately, the result seems to scaled on the bins axis. What am I do wrong here? and how can I fix the bins axis? This is my Stft function, that I use in a larger script that takes in a Wave file and saves the Stft in a CSV that I read later in python and compare with librosa.

void stft(float* signal, int signal_length, fftwf_complex** stft_result, int* num_frames, int nfft, float* wnd, int padding) {
    // Allocate input and output arrays
    fftwf_complex* in = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex) * FFT_SIZE);
    fftwf_complex* out = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex) * FFT_SIZE);

    // Create FFTW plan
    fftwf_plan plan = fftwf_plan_dft_1d(FFT_SIZE, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

    // Padding copy x to x_padded: padding by nfft /2 on both sides
    int pad_width = FFT_SIZE / 2;
    int input_size = signal_length;
    int padded_signal_length = signal_length + pad_width * 2;

    float* padded_signal = (float*)calloc(padded_signal_length, sizeof(float));
    // // Zero padding
    // memcpy(padded_signal + pad_width, signal, signal_length * sizeof(float));

    // Reflect padding
    pad_reflect(signal, signal_length, padded_signal, padded_signal_length, pad_width);

    // update lenghts
    *num_frames  = (padded_signal_length - FFT_SIZE) / HOP_SIZE + 1;
    *stft_result = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex) * (*num_frames) * nfft);

    // Process frames
    for (int i = 0; i < *num_frames; ++i) {
        int start_idx = i * HOP_SIZE;

        // Apply Hanning window and zero-padding
        for (int j = 0; j < FFT_SIZE; ++j) {
            if (start_idx + j < padded_signal_length) {
                in[j][0] = padded_signal[start_idx + j] * wnd[j];
                in[j][1] = 0.0;  // Assuming imaginary part is zero
            }
        }

        // FFT using fftw
        fftwf_execute(plan);
        int idx;  // Declare a variable to store the result
        // Normalize the output and rearrange frequencies
        for (int j = 0; j < nfft; ++j) {
            // Check if j is less than nfft - 1
            if (j < nfft - 1) {
                // If true, assign the current value of j to idx
                idx = j;
            } else {
                // If false (j is greater than or equal to nfft - 1), assign 0 to idx
                idx = 0;
            }
            
            (*stft_result)[i * nfft + j][0] = out[idx][0];
            (*stft_result)[i * nfft + j][1] = out[idx][1];
        }
    }

    fftwf_destroy_plan(plan);
    fftwf_free(in);
    fftwf_free(out);
    free(padded_signal);
}

The resulting plots are the following