2
$\begingroup$

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;
}
$\endgroup$
4
  • 1
    $\begingroup$ I haven't used FFTW before so this may be a dumb question. Many times, the term "frame" is used in audio programming to declare a collection of samples (for example 1024 samples held in a buffer). Are you sure in your context it represents samples? $\endgroup$
    – ZaellixA
    Commented Jan 30, 2023 at 12:34
  • 1
    $\begingroup$ @ZaellixA right, a frame (buffer,window) and also a windowing function are required for properly working, however in my case the frame is the entire file. but that's not the actual problem, something else is wrong with my code. $\endgroup$
    – user174174
    Commented Jan 31, 2023 at 9:06
  • 1
    $\begingroup$ Thanks for the clarification. I am trying to understand the terminology and the way FFTW works in order to try and provide some help. I believe that there is some truth in Hilmar's answer. I am not sure why scaling the spectrum results in added noise but for the other cases it seems like time aliasing is the issue here. Imagine the simple case of zeroing out all bins above a specific index. This would result in a brick-wall filter which has an infinite impulse response, thus it would alias in time (not realisable). Something similar must be happening in your case too. (cont.) $\endgroup$
    – ZaellixA
    Commented Jan 31, 2023 at 10:05
  • $\begingroup$ What I suggest you try is to either use the overlap-add or overlap-save methods, using FIR/IIR filtering, or try and pad A LOT of zeros to your input_samples array (and similarly change your sfinfo.frames to match the value) in order to try and minimise the effects of aliasing. If by padding a lot of zeros you see the situation becoming better then this could be indeed an indication of time aliasing. Hope this helps somehow. $\endgroup$
    – ZaellixA
    Commented Jan 31, 2023 at 10:09

1 Answer 1

2
$\begingroup$

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.

$\endgroup$
1
  • 2
    $\begingroup$ TYSM for the good interpretation, however it seems I have to use FT bc I need some complex filterings (vocal boost/cut, etc). Actually I was able to debug my code and main problem was incorrect handling of complex numbers that FT generates. Now it's working fine. Being a total newbie in SP, I had no idea of firs/iirs, but again thanks to you I implemented a simple fir. The problem with those are that generating coefficients for a fir are sometimes not easy, I was able to use some tools to get a simple bandpass/stop co-efs but not further than that. $\endgroup$
    – user174174
    Commented Feb 2, 2023 at 13:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.