Implementing convolution in C++ using fftw 3

Question

UPDATE

See my fundamental based question on DSP stackexchange here

UPDATE

I am still experiencing crackling in the output. These crackles are now less pronounced and are only audible when the volume is turned up

UPDATE

Following the advice given here has removed the crackling sound from my output. I will test with other available HRIRs to see if the convolution is indeed working properly and will answer this question once I've verified that my code now works

UPDATE

I have made some progress, but I still think there is an issue with my convolution implementation. The following is my revised program:

#define HRIR_LENGTH 512
#define WAV_SAMPLE_SIZE 256

    while (signal_input_wav.read(&signal_input_buffer[0], WAV_SAMPLE_SIZE) >= WAV_SAMPLE_SIZE)
    {
#ifdef SKIP_CONVOLUTION
        // Copy the input buffer over
        std::copy(signal_input_buffer.begin(),
                  signal_input_buffer.begin() + WAV_SAMPLE_SIZE,
                  signal_output_buffer.begin());

        signal_output_wav.write(&signal_output_buffer[0], WAV_SAMPLE_SIZE);
#else
        // Copy the first segment into the buffer
        // with zero padding
        for (int i = 0; i < HRIR_LENGTH; ++i)
        {
            if (i < WAV_SAMPLE_SIZE)
            {
                signal_buffer_fft_in[i] = signal_input_buffer[i];
            }
            else
            {
                signal_buffer_fft_in[i] = 0; // zero pad
            }
        }

        // Dft of the signal segment
        fftw_execute(signal_fft);

        // Convolve in the frequency domain by multiplying filter kernel with dft signal
        for (int i = 0; i < HRIR_LENGTH; ++i)
        {
            signal_buffer_ifft_in[i] = signal_buffer_fft_out[i] * left_hrir_fft_out[i]
                - signal_buffer_fft_out[HRIR_LENGTH - i] * left_hrir_fft_out[HRIR_LENGTH - i];

            signal_buffer_ifft_in[HRIR_LENGTH - i] = signal_buffer_fft_out[i] * left_hrir_fft_out[HRIR_LENGTH - i]
                + signal_buffer_fft_out[HRIR_LENGTH - i] * left_hrir_fft_out[i];

            //double re = signal_buffer_out[i];
            //double im = signal_buffer_out[BLOCK_OUTPUT_SIZE - i];
        }

        // inverse dft back to time domain
        fftw_execute(signal_ifft);

        // Normalize the data
        for (int i = 0; i < HRIR_LENGTH; ++i)
        {
            signal_buffer_ifft_out[i] = signal_buffer_ifft_out[i] / HRIR_LENGTH;
        }

        // Overlap-add method
        for (int i = 0; i < HRIR_LENGTH; ++i)
        {
            if (i < WAV_SAMPLE_SIZE)
            {
                signal_output_buffer[i] = signal_overlap_buffer[i] + signal_buffer_ifft_out[i];
            }
            else
            {
                signal_output_buffer[i] = signal_buffer_ifft_out[i];
                signal_overlap_buffer[i] = signal_output_buffer[i]; // record into the overlap buffer
            }
        }

        // Write the block to the output file
        signal_output_wav.write(&signal_output_buffer[0], HRIR_LENGTH);

#endif
    }

The resulting output sound file contains crackling sounds; presumably artefacts left from the buggy fftw implementation. Also, writing blocks of 512 (HRIR_LENGTH) seems to result in some aliasing, with the soundfile upon playback sounding like a vinyl record being slowed down. Writing out blocks of size WAV_SAMPLE_SIZE (256, half of the fft output) seems to playback at normal speed. However, irrespective of this the crackling sound remains.

ORIGINAL

I'm trying to implement convolution using the fftw library in C++. I can load my filter perfectly fine, and am zero padding both the filter (of length 512) and the input signal (of length 513) in order to get a signal output block of 1024 and using this as the fft size.

Here is my code:

#define BLOCK_OUTPUT_SIZE 1024
#define HRIR_LENGTH 512

#define WAV_SAMPLE_SIZE 513
#define INPUT_SHIFT 511

while (signal_input_wav.read(&signal_input_buffer[0], WAV_SAMPLE_SIZE) >= WAV_SAMPLE_SIZE)
{
#ifdef SKIP_CONVOLUTION
    // Copy the input buffer over
    std::copy(signal_input_buffer.begin(),
              signal_input_buffer.begin() + WAV_SAMPLE_SIZE,
              signal_output_buffer.begin());

    signal_output_wav.write(&signal_output_buffer[0], WAV_SAMPLE_SIZE);
#else
    // Zero pad input
    for (int i = 0; i < INPUT_SHIFT; ++i)
        signal_input_buffer[WAV_SAMPLE_SIZE + i] = 0;

    // Copy to the signal convolve buffer
    for (int i = 0; i < BLOCK_OUTPUT_SIZE; ++i)
    {
        signal_buffer_in[i] = signal_input_buffer[i];
    }

    // Dft of the signal segment
    fftw_execute(signal_fft);

    // Convolve in the frequency domain by multiplying filter kernel with dft signal
    for (int i = 1; i < BLOCK_OUTPUT_SIZE; ++i)
    {
        signal_buffer_out[i] = signal_buffer_in[i] * left_hrir_fft_in[i]
            - signal_buffer_in[BLOCK_OUTPUT_SIZE - i] * left_hrir_fft_in[BLOCK_OUTPUT_SIZE - i];

        signal_buffer_out[BLOCK_OUTPUT_SIZE - i]
            = signal_buffer_in[BLOCK_OUTPUT_SIZE - i] * left_hrir_fft_in[i]
                    + signal_buffer_in[i] * left_hrir_fft_in[BLOCK_OUTPUT_SIZE - i];

        double re = signal_buffer_out[i];
        double im = signal_buffer_out[BLOCK_OUTPUT_SIZE - i];
    }

    // inverse dft back to time domain
    fftw_execute(signal_ifft);

    // Normalize the data
    for (int i = 0; i < BLOCK_OUTPUT_SIZE; ++i)
    {
        signal_buffer_out[i] = signal_buffer_out[i] / i;
    }

    // Overlap and add with the previous block
    if (first_block)
    {
        first_block = !first_block;
        for (int i = 0; i < BLOCK_OUTPUT_SIZE; ++i)
        {
            signal_output_buffer[i] = signal_buffer_out[i];
        }
    }
    else
    {
        for (int i = WAV_SAMPLE_SIZE; i < BLOCK_OUTPUT_SIZE; ++i)
        {
            signal_output_buffer[i] = signal_output_buffer[i] + signal_buffer_out[i];
        }
    }

    // Write the block to the output file
    signal_output_wav.write(&signal_output_buffer[0], BLOCK_OUTPUT_SIZE);
#endif
}

In the end, the resulting output file contains garbage, but is not all zeros. Things I have tried:

1) Using the standard complex interface fftw_plan_dft_1d with the appropriate fftw_complex type. Same issues arise.

2) Using a smaller input sample size and iterating over the zero padded blocks (overlap-add).

I also note that its not a fault of libsndfile; toggling SKIP_CONVOLUTION does successfully result in copying the input file to the output file.