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