Incorrect result when using int16 fixed point FFT

[Ryuk17]

I met a problem when using int16 fixed point FFT, here is my mainly code

    kiss_fft_cpx cin[NFFT];
    kiss_fft_cpx cout[NFFT];
    kiss_fft_cfg  kiss_fft_state;
    memset(cin, 0, sizeof(short) * NFFT);
    memset(cout, 0, sizeof(short) * NFFT);
    while(fread(input, sizeof(short), NFFT, in))
    {
        for(int i=0; i<NFFT;i++)
        {
            cin[i].r = input[i];
            cin[i].i = 0;
        }
        kiss_fft_state = kiss_fft_alloc(NFFT, 0, 0, 0);
        kiss_fft(kiss_fft_state, cin, cout);
        kiss_fft_state = kiss_fft_alloc(NFFT, 1, 0, 0);
        kiss_fft(kiss_fft_state, cout, cin);

        for(int i=0; i<NFFT;i++)
        {
            output[i] = cin[i].r;
        }
        ret = fwrite(output, sizeof(short), NFFT, out);
    }

The input audio spectrum is

The output audio spectrum is

I have no idea where is wrong, can you give me some advice?

[JulienMaille]

Are you plotting Y=fftinv(fft(x)) and expecting to find Y=X? I think you will obtain Y=NFFT.X Also if you work with real signals (ie. not complex) then use the dedicated kiss_fftr_**** functions.

Other problems:

• if output type is short, you might overflow SHRT_MAX (32767)
• memset(cin, 0, sizeof(short) * NFFT); cin type is kiss_fft_cpx so you would need 2*sizeof(float) but you probably don't even need that memset

[Ryuk17]

Yes, I junst want to test the process of Y=fftinv(fft(x)) and expecting to find Y=X. Does Y=NFFT.X mean NFFT multiply X, I divide the result by NFFT, it doesn't work. My input are audio samples which is int16 in Q15 format. I have changed the related macro, thus, all kiss_fft related variants are in short type.
I use kissfftr**** functions and get the same result.
Any other suggestions, bro?

[JulienMaille]

output[i] = cin[i].r/NFFT; doesn't help? No I don't have other suggestions. Have you checked the content of cin[i].r with a debugger?

[Ryuk17]

Using output[i] = cin[i].r/NFFT will make the output be samller again. Here is the complete code, you can try it yourself.

#include <string.h>
#include <stdio.h>
#include "kiss_fft.h"
#include "_kiss_fft_guts.h"

#define NFFT        (512)

int main() {

    int ret = 0;
    char inFileName[512];
    char outFileName[512];

    short input[NFFT];
    short output[NFFT];

    FILE *in;
    FILE *out;

    sprintf(inFileName, "./sample/test.wav");
    printf("Open %s\n", inFileName);
    in = fopen(inFileName, "rb");
    if(!in)
    {
        perror(inFileName);
        return -1;
    }

    fread(input, sizeof(char), 44, in);

    sprintf(outFileName, "./sample/fftout.wav");
    printf("Open %s\n", outFileName);
    out = fopen(outFileName, "wb");
    if(!out)
    {
        perror(outFileName);
        return -1;
    }
    fwrite(input, sizeof(char), 44, out);

    kiss_fft_cpx cin[NFFT];
    kiss_fft_cpx cout[NFFT];
    kiss_fft_cfg  kiss_fft_state;
    memset(cin, 0, sizeof(short) * NFFT);
    memset(cout, 0, sizeof(short) * NFFT);
    while(fread(input, sizeof(short), NFFT, in))
    {
        for(int i=0; i<NFFT;i++)
        {
            cin[i].r = input[i];
            cin[i].i = 0;
        }
        kiss_fft_state = kiss_fft_alloc(NFFT, 0, 0, 0);
        kiss_fft(kiss_fft_state, cin, cout);
        kiss_fft_state = kiss_fft_alloc(NFFT, 1, 0, 0);
        kiss_fft(kiss_fft_state, cout, cin);

        for(int i=0; i<NFFT;i++)
        {
            output[i] = cin[i].r;
        }
        ret = fwrite(output, sizeof(short), NFFT, out);
    }

    fclose(in);
    fclose(out);
    printf("Finished\n");
    return 0;
}

[mborgerding]

Looks like you're getting the peaks of the spectrogram, but a lot of the details have vanished. Is it possible to increase the input magnitude? You will lose log2(nfft) low bits of detail due to "Scaling is done both ways for the fixed-point version (for overflow prevention)." (from README.md) This is owed to lack of support for saturating math in ANSI C.

[mborgerding]

also, move the alloc functions outside your loop and free your structs.

This question might be better suited to stackoverflow.com