1
$\begingroup$

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: convolution 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: convolution output 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).

$\endgroup$

1 Answer 1

1
$\begingroup$

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).

$\endgroup$
8
  • $\begingroup$ It's at step 4 that I'm stuck, believe me I've tried to fix it but the truth is I have no idea what's causing it. $\endgroup$
    – Ott
    Sep 23, 2019 at 14:39
  • $\begingroup$ I've added some code from my AudioData class to the question $\endgroup$
    – Ott
    Sep 23, 2019 at 15:17
  • $\begingroup$ What is your buffer_size and what's your fft_size ? During the convolution you should be doing a complex multiplication but I don't see any complex variables in there? Inspect your impulse response post-FFT buffer. There should all ones (real part = 1, imaginary part = 0) in there. $\endgroup$
    – Hilmar
    Sep 23, 2019 at 15:45
  • $\begingroup$ buffer_size is 448, block_size is 512, and fft_size is 1024. I'm storing imaginary numbers as simple pairs of floats; my frequency content vectors have twice as many numbers as normal. $\endgroup$
    – Ott
    Sep 23, 2019 at 16:13
  • $\begingroup$ But how do you do complex multiplication? A simple vector float multiply doesn't work here. See mathworld.wolfram.com/ComplexMultiplication.html $\endgroup$
    – Hilmar
    Sep 23, 2019 at 17:05

Your Answer

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

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