Partitioned overlap-add convolution - strange behavior at buffer boundaries

Question

I've implemented a convolution reverb that operates in real-time, one audio buffer at a time (using FFTS for the fft bits). However, there's some strange behavior at the start of every buffer. Convolving a sinusoid with an impulse (a 1 followed by many zeroes), I don't get a sinusoid as the output: Instead, I get peaks that are exactly twice the amplitude they should be at the start of every buffer. In fact, even if I don't use the spectra from the actual impulse file and instead multiply the complex portion of the input by 0, I get the same result. Conversely, if I multiply the real portion of the input by 0, I get 0 at the start of every buffer: It seems like an off-by-one error or something, but then again the input is left intact if I don't modify the frequency components. I've been reading papers and going over my code for the past week, and I'd really appreciate it if someone could verify that it's correct.

class AudioEffectConvolver : public IAudioEffect
{
public:
    AudioEffectConvolver(const char *impulse_name);
    void process(AudioData *);
    void calculateImpulse(unsigned buffer_size);
    ~AudioEffectConvolver();
private:
    std::shared_ptr<AudioData> impulse;
    std::vector<float> impulse_bins;
    std::vector<std::vector<float>> partitions;
    std::vector<std::vector<float>> bin_ring;
    std::vector<float> overlap;
    ffts_plan_t *forward;
    ffts_plan_t *backward;
    unsigned ring_index = 0;
    bool impulse_calculated = false;
    unsigned block_size = 0;
};

static unsigned npo2(unsigned size)
{
    size--;
    size |= size >> 1;
    size |= size >> 2;
    size |= size >> 4;
    size |= size >> 8;
    size |= size >> 16;
    size++;

    return size;
}

void AudioEffectConvolver::calculateImpulse(unsigned buffer_size)
{
    unsigned impulse_size = impulse->frames();
    block_size = npo2(buffer_size);
    unsigned fft_size = block_size * 2;
    forward = ffts_init_1d_real(fft_size, FFTS_FORWARD);
    backward = ffts_init_1d_real(fft_size, FFTS_BACKWARD);

    overlap.resize(block_size);
    std::vector<float> window(fft_size);

    for (unsigned i = 0; i * block_size < impulse->frames(); ++i)
    {
        unsigned offset = i * block_size;
        if (impulse->frames() >= offset + block_size)
        {
            memcpy(window.data(), impulse->split(0) + offset, sizeof(float) * block_size);
            memset(window.data() + block_size, 0, sizeof(float) * (fft_size - block_size));
        }
        else
        {
            memcpy(window.data(), impulse->split(0) + offset, sizeof(float) * (impulse->frames() - offset));
            memset(window.data() + impulse->frames() - offset, 0, sizeof(float) * (fft_size - (impulse->frames() - offset)));
        }
        partitions.emplace_back(fft_size + 2); // (n / 2 + 1) * 2
        ffts_execute(forward, window.data(), partitions[i].data());

        bin_ring.resize(i + 1);
        bin_ring[i].resize(fft_size + 2);
    }
    impulse_calculated = true;
}


void AudioEffectConvolver::process(AudioData *buffer)
{
    unsigned buffer_size = buffer->frames();
    if (!impulse_calculated) calculateImpulse(buffer->frames());
    unsigned fft_size = block_size * 2;

    buffer->resize(fft_size);

    memset(bin_ring[ring_index].data(), 0, sizeof(float) * (fft_size + 2));
    ffts_execute(forward, buffer->split(0), bin_ring[ring_index].data());

    std::vector<float> convolution;
    convolution.resize(fft_size + 2);

    for (unsigned k = 0; k < partitions.size(); ++k)
    {
        int index = ring_index - k;
        while (index < 0) index += (int)bin_ring.size();
        for (unsigned i = 0; i < fft_size + 2; ++i)
        {
            convolution[i] += bin_ring[index][i] * partitions[k][i];
        }
    }

    std::vector<float> output;
    output.resize(fft_size);
    ffts_execute(backward, convolution.data(), output.data());
    //output.resize(block_size); // circular convolution; chop off the second half

    for (unsigned i = 0; i < buffer_size; ++i)
    {
        float outsample = (output[i] + overlap[i]) / (fft_size);
        buffer->sample(i) = outsample;
    }

    memcpy(overlap.data(), output.data() + buffer_size, sizeof(float) * (block_size - buffer_size));

    ring_index++;
    if (ring_index >= bin_ring.size()) ring_index = 0;
    buffer->resize(buffer_size);
}

class AudioData
{
public:
    AudioData(const char *filename);
    AudioData(unsigned frames, unsigned channels = 1, unsigned rate = 44100);
    ~AudioData();
    float sample(unsigned frame, unsigned channel = 0) const;
    float& sample(unsigned frame, unsigned channel = 0);
    float seek(float seconds, unsigned channel = 0) const;
    void resize(unsigned frames);

    float *data();
    const float *data() const;
    const float *split(unsigned channel);
    unsigned frames() const;
    unsigned rate() const;
    unsigned channels() const;

private:
    std::vector<float> audio_data;
    std::vector<float *> splits;
    unsigned frame_count;
    unsigned sampling_rate;
    unsigned channel_count;
    void clearSplits();
};

float& AudioData::sample(unsigned frame, unsigned channel)
{
    clearSplits();
    return audio_data[frame * channel_count + channel];
}

void AudioData::resize(unsigned frames)
{
    clearSplits();
    audio_data.resize(frames * channel_count);
    frame_count = frames;
}

float *AudioData::data()
{
    return audio_data.data();
}

const float *AudioData::split(unsigned channel)
{
    if (splits[channel] != nullptr) return splits[channel];
    if (channel_count == 1) return data();

    float *split = new float[frame_count];

    for (unsigned i = 0; i < frame_count; ++i)
    {
        split[i] = sample(i, channel);
    }

    splits[channel] = split;
    return split;
}

unsigned AudioData::frames() const
{
    return frame_count;
}

void AudioData::clearSplits()
{
    for (unsigned i = 0; i < splits.size(); ++i)
    {
        if (splits[i] != nullptr) delete[] splits[i];
        splits[i] = nullptr;
    }
}

In my specific case, the buffers are 448 samples (something to do with WASAPI shared mode).

Answer 1

You call a lot of methods and functions that are not included so it's hard to read. Here is how I debug this step by step.

1. Verify your audio framework. Do NOTHING in the process() function other than copying the input to the output
2. Verify simple processing. Now add multiplication with 0.5 or something simple like this.
3. Verify the FFT based processing. Do just the zero padding, forward FFT, inverse FFT, and output calculation
4. Add a "pass through" impulse response. Just a single sample at $n=0$
5. Verify your overlap handling: use an impulse response with a single tap at $n=256$
6. Verify the framing of the impulse response: use an impulse response with a single tap at $n=2000$

At each step calculated the expected results and calculate the RMS error to the expected result. For single precision floating point, that should be in the order of -130dB or so. Use both a sine wave and a unit impulse as input signal.

If you get a large error, stop and fix this step.

If that all checks out, chances are you code is good but you should still test it with a random input signal and a real room impulse response and calculate the RMS error to a known-good reference model (Matlab, Octave, Python).