Search code examples
csignal-processingfftfftw

FFTW interlaced zero's in output (not end padding)


I can't find where this is documented or similar questions but when I perform an FFT with more than 2k elements my output interlaces zeros, for example if I double N to 4k, my output is 4k elements, with 2k data points alternating with 2k zeros i.e. {...9413.5, 0.0, 9266.2, 0.0, ...}. Can someone explain what I overlooked, thanks!

//testing the fftw compile with gcc fftw_test.c -lfftw3 -lm -o fftw_test
#include <stdlib.h>
#include <fftw3.h>
#include <math.h>

#define N 1024*4

int main(void)
{
    double input[N] = 
        {/*This is filled with 4k elements representing a basic sine wave*/};

    //declare
    double *in;
    fftw_complex *out;
    fftw_plan my_plan;

    //allocate and assign
    in = (double*) fftw_malloc(sizeof(double) * N);
    for (int i = 0; i < N; i++)
    {
        in[i] = input[i];
    }
    out = (fftw_complex*) fftw_malloc((sizeof(fftw_complex) * N) + 1);
    my_plan = fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);

    //execute
    fftw_execute(my_plan);

    //print magnitude
    FILE *log = fopen("log.txt", "w");
    for (int i = 0; i < N; i++)
    {
        input[i] = sqrt(out[i][0]*out[i][0] + out[i][1]*out[i][1]);
        fprintf(log, " %.01lf,", input[i]);
    }

    //exit
    fclose(log);
    fftw_destroy_plan(my_plan);
    fftw_free(in);
    fftw_free(out);
    return 0;
}

The code is generated using this prewritten python script (I realize it only generates 2k points, I just copied it twice):

#this program generates a sine wave and prints it to sine.txt

import numpy as np
import matplotlib.pylab as plt
file = open("sine.txt","w")

x = np.linspace(0, 2048, 2048)
y = [2048]

plt.plot(np.int16(np.sin(x/16)*2048 + 2048))

for i in x:
    file.write(str(np.int16(np.sin(i/16)*2048+2048)))
    file.write(", ")
file.close()
plt.show()

Solution

  • This is key:

    I realize it only generates 2k points, I just copied it twice

    You created a sine wave signal whose period does not evenly divide your signal length. Its Discrete Fourier Transform has values for all frequencies because of this (if it did evenly divide the signal length, you'd see only two non-zero elements in the FFT). But then you duplicated the signal, in effect creating a signal where N is exactly twice the period. Consequently, the frequency content for all odd k is zero.

    If you make four copies of your signal, you'll find three zeros between each non-zero component. If you make 8 copies, you'll find 7 zeros. The non-zero elements are the same in all these cases, but scaled by the number of copies.