Given a static time domain signal, why does the magnitude oscillate in the frequency domain after FFT?

Question

I'm trying to:

• Generate a sine wave in the time domain (locally within my program at 500 Hz).
• Analyze it with FFT
• Print the results in the frequency domain.

The expected results are that since the sine wave doesn't change its characteristics (I just keep increasing the t variable for each invocation of the function), and stays at 500 Hz, I should get a constant magnitude of the signal in the frequency domain.

Actual results, are that the magnitude of the 500 Hz signal in the frequency domain rather wildly bounces up and down.

Here I have two screenshots of the plot. Please notice the magnitude of the signal changing between each iteration. It pulsates up and down, almost like the sine wave itself.

And just to be sure, I've also plotted the sine wave in the time domain, and checked that it is a nice, continuous function between each invocation (each invocation of the sine generate function generates 512 bytes of data, the plot is for 1000 bytes, to see it is continuous between invocations):

Why does this happen? I've thought about normalization, but just multiplying or dividing everything by some constant wouldn't get rid of this oscillating behavior? Do I need to do something with the imaginary part of the FFT results?

To reproduce, put this in a file live.gnu and plot it with gnuplot live.gnu:

set xrange[0:22000]
set yrange[0:1.0]
set boxwidth 100.0
set style fill solid plot "outfile.log" with boxes
pause 2
reread

I have attached a minimally reproducible example of the code (screenshots are from exactly this code):

Compile with:

gcc fftw_example4.c -lm -lfftw3 -o fftw_example4

/**
 *
 * Minimal reproducible version of oscillating magnitude of signal in the frequency domain.
 *
 * 1: Generate a sine wave, with minor random variations.
 * 2: Do FFT.
 * 3: Write the results to a file (and also print to screen).
 *
 * Notice the 500 Hz signal in the frequency domain, its magnitude changes for each iteration of the while(1) loop.
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <sys/time.h>
#include <fcntl.h>
#include <fftw3.h>
#include <math.h>
#include <unistd.h>
#include <limits.h>

#define INPUT_BUFFER_SIZE 512
#define RESULT_BUFFER_SIZE (INPUT_BUFFER_SIZE / 2 + 1)
#define SAMPLE_RATE 44100

#define REAL 0
#define IMAG 1

#define FREQ 0
#define MAGN 1

#define OUTFILE "./outfile.log"

int generate_sine(double gain, double frequency, double *input_buffer, int input_buffer_size)
{
    static unsigned int t = 0;

    memset(input_buffer, 0, sizeof(double) * input_buffer_size);

    for(int i = 0; i < input_buffer_size; i++) {
        double sample = gain * sin(((double) t++ * M_PI * 2.0 * frequency) / (double) (SAMPLE_RATE));

        input_buffer[i] = sample;
    }

    return 0;
}

int analyze_init(fftw_plan *plan, double *input_buffer, int input_buffer_size, fftw_complex *result_buffer)
{
    unsigned flags = 0;

    *plan = fftw_plan_dft_r2c_1d(input_buffer_size, input_buffer, result_buffer, flags);

    return 0;
}

int analyze_frequencies(fftw_plan *plan, double *input_buffer, int input_buffer_size, fftw_complex (*result_buffer)[RESULT_BUFFER_SIZE], int result_buffer_size)
{
    fftw_execute(*plan);

    for(int i = 0; i < result_buffer_size; i++) {
        (*result_buffer)[i][MAGN] = (double) 2.0 / (double) (result_buffer_size) * fabs((*result_buffer)[i][REAL]);
        (*result_buffer)[i][FREQ] = i * (SAMPLE_RATE / (double) input_buffer_size);
    }

    return 0;
}

int print_frequencies(fftw_complex (*result_buffer)[RESULT_BUFFER_SIZE], int result_buffer_size)
{
    FILE *f = fopen(OUTFILE, "w");

    if(!f) {
        perror("fopen");

        exit(1);
    }

    for(int i = 1; i < result_buffer_size - 1; i++) {
        if((*result_buffer)[i][FREQ] >= 300 && (*result_buffer)[i][FREQ] <= 600)
            printf("%f %f\n", (*result_buffer)[i][FREQ], (*result_buffer)[i][MAGN]);

        fprintf(f, "%f %f\n", (*result_buffer)[i][FREQ], (*result_buffer)[i][MAGN]);
    }

    fclose(f);

    return 0;
}

int main(int argc, char **argv)
{
    double input_buffer[INPUT_BUFFER_SIZE];

    fftw_plan plan;
    int result_buffer_size = RESULT_BUFFER_SIZE;
    fftw_complex result_buffer[RESULT_BUFFER_SIZE];

    analyze_init(&plan, input_buffer, INPUT_BUFFER_SIZE, result_buffer);

    while(1) {
        generate_sine(1.0, 500, input_buffer, INPUT_BUFFER_SIZE);

        if(analyze_frequencies(&plan, input_buffer, INPUT_BUFFER_SIZE, &result_buffer, RESULT_BUFFER_SIZE) < 0) {
            exit(1);
        }

        print_frequencies(&result_buffer, RESULT_BUFFER_SIZE);
    }

    fftw_destroy_plan(plan);

    exit(0);
}

Answer 1

You should compute the magnitude of the complex output of the FFT instead of just taking the magnitude of its real part. If $X[k]$ is the (complex-valued) $k^{th}$ frequency bin, its magnitude is

$$\big|X[k]\big|=\sqrt{\left(\textrm{Re}\big\{X[k]\big\}\right)^2+\left(\textrm{Im}\big\{X[k]\big\}\right)^2}$$