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;
}