FFT, why my program works but doesn't?

Question

I'm working on a simple 1D FFT and iFFT for two days now but no matter what I do, the code fails. I'm using FFTW library and libsndfile on linux and my goal is to use FFT for filtering audio.

The code reads a wave file (input.wav) performs a global FT , modify frequencies, inverse FT and writes to file (output.wav). I know it works because if I set all frequency coefficients to 0 then I get a silent file (same length, no contents or DC offset) and copying buffers gives me a duplicate but any process except these causes to change the amplitude and introducing white noise and some weird things.

First I thought it's because the wave file length isn't a power of 2, so I padded some silence at the end, now it's $2^{23}$ samples (16bit PCM 44100Hz mono). But that didn't help.

Then I thought it might be that input and output buffers are the same (although the FFTW docs says it's ok) so I used another buffer, but it also didn't change anything. I'd be so grateful for any clues or answers that why it behaves strangely.

my algorithm is:

• open file
• read all the file into input buffer
• real to real fft on input buffer results in output buffer
• modifying output buffer, then inverse fft and write:
  • all elements = 0
    • $\rightarrow$ expected result (same-length total silence wave file)
  • out = in
    • $\rightarrow$ expected result (exact copy)
  • all elements *= $x$
    • $\rightarrow$ unpredictable results, amplitude change according to $x$ (less for $0<x<1$, more for $x > 1$) + some weird harmonics + white noise + clippings + for some values a reversed time-domain signal mixed with original)
  • $x >1,$ $x <1,$ $x=0$ for 3 bands
    • $\rightarrow$ again amp change, white noise, clips, no difference in frequency harmonics.

and here is the full code if more probbing is required:

#include <stdio.h>
#include <fftw3.h>
#include <sndfile.h>
#include <stdlib.h>
#include <math.h>

int main(int argc, char *argv[])
{
    SF_INFO sfinfo, out_sfinfo;
    SNDFILE *input_file,*output_file;
    sfinfo.format = 0; // prepare for reading
    unsigned i, fft_size;
    fftw_plan forward,reverse;

    input_file = sf_open("input.wav", SFM_READ, &sfinfo);
    if (!input_file)
    {
        printf("Error openning input file...\n");
        return EXIT_FAILURE;
    }

    // Allocate memory
    double *input_samples = (double *)fftw_malloc(sizeof(double) * sfinfo.frames);
    double *out_samples = (double *)fftw_malloc(sizeof(double) * sfinfo.frames);

    // Read the input samples
    sf_readf_double(input_file, input_samples, sfinfo.frames);
    fft_size = sfinfo.frames;  // 8388608 samples
    sf_close(input_file);

    printf("%d frames read.\n", sfinfo.frames);

    forward = fftw_plan_r2r_1d(sfinfo.frames, input_samples, 
                            out_samples, FFTW_R2HC, 
                            FFTW_ESTIMATE);
    fftw_execute(forward);

    // This gives a total silence
    for (i = 0; i < sfinfo.frames; i++)
        input_samples[i]  = out_samples[i];
/*
    // This rises the amplitude & adds white noise
    for (i = 0; i < sfinfo.frames; i++)
        input_samples[i]  = 2.1*out_samples[i];
    
    // Trying bass,mid,tre has no apparent effect other than amp.
    for (i = 0; i < sfinfo.frames ; i++) {
        unsigned frq = i * sfinfo.frames / 44100;
        // frq calculated by: index * (buffersize / sample rate)
        if (frq > 8000) 
            input_samples[i] = 2.1 * out_samples[i];
        else if (frq < 100)
            input_samples[i] = 2.1 * out_samples[i];
        else
            input_samples[i] = 0.01 * out_samples[i];
    }

*/
    // Doing the inverse FFT 
    reverse = fftw_plan_r2r_1d(sfinfo.frames, input_samples, 
                            out_samples, FFTW_HC2R, 
                            FFTW_ESTIMATE);
    fftw_execute(reverse);

    // Normalize
    for (i = 0; i < sfinfo.frames; i++)
        out_samples[i] /= fft_size 
        ;
    

    // Opening and writing to the ouput file
    out_sfinfo = sfinfo;
    output_file = sf_open("output.wav", SFM_WRITE, &out_sfinfo);
    if (!output_file)
    {
        printf("Error openning output file...\n");
        return EXIT_FAILURE;
    }
    sf_write_double(output_file, out_samples, sfinfo.frames);
    sf_close(output_file);

    // give back taken resources
    fftw_destroy_plan(forward);
    fftw_destroy_plan(reverse);
    fftw_free(input_samples);
    fftw_free(out_samples);

    return 0;
}

Answer 1

FFT for filtering audio.

Frequency domain filtering is mathematically quite complicated. Manipulating the spectrum directly typically results in lots of time-domain aliasing, which I'm guessing is what you are seeing here.

If your desired filter is LTI (linear time-invariant) you can use an overlap-add or overlap-save algorithm to implement frequency domain filtering (starting with a time domain impulse response of the filter).

I know it works

No, it doesn't. It ONLY works if you scale the entire spectrum with a constant factor. You can't implement frequency dependent audio filtering this way.

In general, it's much easier to implement audio filtering directly in the time domain using suitable IIR of FIR filters.