Real Time Partitioned Convolution Not Working

Question

EDIT: Scroll down for actual working code.

I'm working on implementing a real-time convolvution reverb JACK client on C and I've been trying to follow a number of sources (including Gardner and Wefers (pages 110-111) ) to no avail. I've tried to follow the texts to the best of my ability and honestly, I can't get rid of the horrible sounding artifacts. I looked them up on baudline (spectrogram/waveform visualizer) and you can see a combination of discontinuities and skips where the output is just zero. I tried windowing to smooth out the discontinuities generated by Wefers' description of the uniform partitioning algorithm and even if I get rid of the discontinuities I still get the weirdest aliasing. Nowhere on the internet I've found a single source outlining clearly the code needed for this task.

I'm creating the partitions like so:

// ir is being read from a .wav file using libsndfile and is zero padded to
// have a length that is a multiple of nframes (nframes*partitions, exactly).
for (int k = 0; k < partitions; k++){
            for (int i = 0; i < nframes; i++){
                    // create zero padded partitions
                    ir_time[i]           = ir[k*nframes + i];
                    ir_time[nframes + i] = 0.0;
            }

            fftw_execute(ir_forward);
            for (int i = 0; i < 2*nframes; i++){
                    // write filter partitions
                    fir[k][i] = ir_fft[i];
                    // initialize FDL to zero
                    fdl[k][i] = 0.0 + I*0.0;
            }

    }

And the jack_callback() function where the processing takes place looks like this (copied the stream processing instructions from Wefers p.111):

(Note that the b[1] and b[0] buffers have length 3*nsamples so that they overlap-add in the middle, I know it's weird but I know for a fact that it works with other fftw based code. I know I should change them to 2*nframes length buffers and I will get around to that sometime. They are only used as input buffers here though.)

int jack_callback (jack_nframes_t nframes, void *arg){
    jack_default_audio_sample_t *in, *out;
    int i, j;

    in  = (jack_default_audio_sample_t *)jack_port_get_buffer (input_port , nframes);
    out = (jack_default_audio_sample_t *)jack_port_get_buffer (output_port, nframes);

    for (i = 0; i < nframes; i++){
            // 1. Input buffer acts as a 2B-point sliding sliding window of
            //    the input signal. With each new input block, the right half of the input
            //    buffer is shifted to the left and the new block is stored in the right
            //    half.
            
            // tried windowing the 2*nframes length buffer but it
            // didn't work
            // shift right half of input buffer to the left
            b[1][nframes + i    ] = b[1][nframes + i];

            // store input to right hand side of buffer
            b[1][two_nframes + i] = in[i];

            // prepare buffer for FFT
            
            // tried writing in[i] to in_time[i] and zero padding the rest but
            // it also didn't work
            i_time[i]         = b[1][nframes + i];
            i_time[nframes+i] = b[1][two_nframes + i];
    }
    // 2. All conttents (DFT spectra) in the FDL are shifted up by one slot.
    for (int k = 0; k < partitions - 1; k++){
            for (int i = 0; i < two_nframes; i++){
                    fdl[k+1][i] = fdl[k][i];
            }
    }

    // 3. a 2B-point real-to-complex FFT is computed from the input buffer,
    //    resulting in 2B DFT coefficients. The result is stored in the 
    //    first FDL slot.

    // taking R2C-FFT
    fftw_execute(i_forward);

    // 4. The P sub-filter spectra are pairwisely multiplied with the input
    //    spectra in the FDL. The results are accumulated in the frequency
    //    domain.
    
    for (int i = 0; i < two_nframes; i++){
            fdl[0][i] = i_fft[i];
    }

    for (int i = 0;  i < two_nframes; i++){
            // reset o_fft[i] to erase previous callback buffer
            o_fft[i] = 0;
            for (int k = 0; k < partitions; k++){
                    // accumulation stage
                    o_fft[i] += fir[k][i] * fdl[k][i];
            }
    }

    // 5. Of the accumulated spectral convolutions, an 2B-point C2R-IFFT
    //    is computed. From the resulting 2B samples, the left half is 
    //    discarded and the right half is returned as the next output block 

    fftw_execute(o_inverse);
    for (i = 0; i < nframes; i++){
            // tried this with and without windowing but it didn't work
            //b[1][    nframes + i] = creal( o_time[i]        ) / two_nframes;
            //b[1][two_nframes + i] = creal( o_time[nframes+i]) / two_nframes;
            
            // tried to do overlap-add like this but it also didn't work
            out[i] = b[0][nframes + i]+b[1][nframes + i]; // + conv[i];

            // discarding left side of IFFT and returning rigth half as output.
            //out[i] = creal(o_time[nframes + i]) / two_nframes;
             
            // more overlap-add stuff
            // b[0][i]           = b[1][    nframes + i];
            // b[0][nframes + i] = b[1][two_nframes + i];

            //b[1][nframes+i] = in[i];
    }
    return 0;
}

If you'd like to take a look at the rest of the code I'll gladly upload it to github and add the link here. I'm at my wits end with this. Any help would be massively appreciated!

EDIT: I've since found the problem with my code. I was doing the circular shift wrong. The index for the FDL should decrease when overwriting it to avoid unwanted loops. Here is the actual working code:

int jack_callback (jack_nframes_t nframes, void *arg){
    jack_default_audio_sample_t *in, *out;
    int i, j, k;

    in = (jack_default_audio_sample_t *)jack_port_get_buffer (input_port, nframes);
    out = (jack_default_audio_sample_t *)jack_port_get_buffer (output_port, nframes);

    for (i = 0; i < nframes; i++){
            // nframes come in and are then saved in the right part of the input buffer
            b[1][nframes + i] = in[i];
            i_time[i] = b[1][i];
            i_time[nframes+i] = b[1][nframes+i];
    }
    // take the FFT of the input:
    fftw_execute(i_forward);

    // circular shift (done right!):
    for (i = 0; i < two_nframes; i++){
            for (k = partitions - 1; k > 0; k--){
                    fdl[k][i] = fdl[k-1][i];
            }
    }

    // write the most recent FFT to the first slot of the FDL:
    for (i = 0; i < two_nframes; i++){
            fdl[0][i] = i_fft[i];
            o_fft[i] = 0.0 + I*0.0;
    }
    // Processing block: multiply-add stage across the FDL
    for (i = 0; i < two_nframes; i++){
            for (k = 0; k < partitions; k++){
                    o_fft[i] += fdl[k][i] * fir[k][i];
            }
    }

    // Taking the ifft and returning to the time domain.
    fftw_execute(o_inverse);
    for (i = 0; i < nframes; i++){
            out[i] = vol*creal(o_time[nframes+i])/two_nframes;
            // shift the input buffer to the left.
            b[1][i] = in[i];
    }

    return 0;

}

Answer 1

In step 5 it looks like you are implementing "overlap-save" not "overlap-add". In this case you DON'T zero pad the input, but do an FFT over two full frames of the input.

Overlap add isn't all that complex but there is a lot of details in the buffering, indexing, handling the real values of DC & Nyquist, etc.

I suggest standard debugging procedures. Start with simple signals with known answer and than increase complexity until something breaks. Than you can single step to the point where the break occurs.

1. Add some assertions: check whether all your spectra are conjugate symmetric, whether the real part of your inverse FFT is actually 0 (or sufficiently small), etc.
2. Start with a unit impulse for both signal and impulse response (IR)
3. delay the unit impulse for each by less than a half a frame
4. delay it the unit impulses my more than half a frame but less than one frame
5. Put the unit impulses into the last sample of the first frame
6. Delay them by more than one frame
7. Add two more unit impulses to the first frame of the IR
8. Spread the IR unit impulses over multiple frames
9. add a few unit impulses to the signal as well.
10. Use a unit impulse as IR and a sine wave as signal. Make sure that the period is not an integer divider of the frame size.
11. Repeat steps 3-8 with a sine wave signal

I suggest writing automated unit tests for all these cases. In all test cases, the correct results can be easily calculated by hand. You can also use Matlab, Octave or Python as a reference .

Don't advance to the next step until you pass the test for the current step. Once you make a change, rerun all tests to make sure you didn't break anything else.

Chances are if you successfully complete step 11, your code is working fine, so something is going to break in between. Your next step depends on where exactly it breaks, but the process should narrow it down.

Answer 2

It seems you are implementing a uniform partitioned convolution, which is not very hard to realize. You should stick to the signal flow of Fig 5.2 from Wefers p110 and the stream processing on p111. Partitioned convolution doesn't need windowing. A simple matlab code could be:

NFFT = blockSize * 2;
X_fdl = zeros(NFFT / 2 + 1, hBlocks); % input blocks frequency delay line
H = zeros(NFFT / 2 + 1, hBlocks); % impulse response blocks FFT results
x_tdl = zeros(NFFT, 1); % temp buffer for FFT

%% split very long impulse response into small blocks and convert to FFT domain
startIdx = 1;
endIdx = blockSize;
for hBlock = 1:hBlocks
    H(:, hBlock) = fftr(h(startIdx:endIdx), NFFT);
    startIdx = endIdx + 1;
    endIdx = endIdx + blockSize;
end
    
%% process signal stream x frame by frame and apply UPOLS
startIdx = 1;
endIdx = blockSize;

while endIdx <= length(x)
    % 1. right half of the input buffer is shifted to the left
    x_tdl(1:blockSize) = x_tdl(blockSize + 1:end); 
    
    % and the new block is stored in the right half
    x_tdl(blockSize + 1:end) = x(startIdx:endIdx); 
    
    % 2. circular shift X_fdl, so that 1st column is current block
    X_fdl = circshift(X_fdl, [0, 1]); 

    % 3. take FFT for current block input
    X_fdl(:, 1) = fftr(x_tdl, NFFT);
    
    % 4. point-wise multiplication and accumulate the result
    Y = sum(X_fdl .* H, 2);
    
    % 5. take IFFT for current block output
    yifft = ifftr(Y, NFFT); 
    
    % left half is discarded and the right half is returned
    y(startIdx:endIdx) = yifft(blockSize + 1:end);
    
    % update loop variable
    startIdx = endIdx + 1;
    endIdx = endIdx + blockSize;
end
```