3
$\begingroup$

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;

}

$\endgroup$
4
  • 1
    $\begingroup$ Instead of zero-padding you need to have some overlap between the blocks, which you of course discard from the output. Overlap size is given by the filter size. This ensures that the filter always samples input data for each relevant output pixel. $\endgroup$ Sep 26 at 1:32
  • 1
    $\begingroup$ Please take a look at this. I talked with Bill about this. The issue I brought to Bill is that the inherent minimum delay of efficiently processing a buffer of samples of length $B$ must be $2B$. As best as I can tell, you cannot have a buffer length of half of the FFT length. $\endgroup$ Sep 26 at 18:25
  • $\begingroup$ Chris: Overlap at the output, right? I was using a 50% overlap with my Hann window. Should I do it differently? Robert: From what I gathered, you're right in that the FFT block convolution induces a latency of at least 2B blocks plus the time it takes to perform the calculations. Bill's paper as I understood it suggests that we calculate the convolution of the first blocks in the time domain to reduce this latency to just B samples. Honestly it's been quite tricky to get it right! $\endgroup$
    – Juan F.
    Sep 27 at 23:34
  • $\begingroup$ // Robert: From what I gathered, you're right in that the FFT block convolution induces a latency of at least 2B blocks plus the time it takes to perform the calculations. // --- NO, NO, NO! Not "plus the time it takes to perform the calculations". The latency is because of double buffering; buffering both the input and the output. The calculations of a block of samples is done at the time the samples for the following block are input into a buffer and during the time that the samples of the previous block are being output from a buffer. $\endgroup$ Oct 1 at 19:22
3
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ I have to say that even though the other answer is more clear on the code needed, I found your unit test suggestion exceedingly useful. I've been implementing these unit tests and I've managed to go from step 5 to step 8. Fix a man's code and it will compile once; teach him how to fix bugs and it will compile for life! (The IR seems to be stuck in the first columns of the FDL though!) $\endgroup$
    – Juan F.
    Sep 27 at 23:29
5
$\begingroup$

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
```
$\endgroup$
2
  • $\begingroup$ Thank you for uploading your code, I'm going to try to implement it in C. I think the problem with mine might be the circular shifting of the FDL or the initial partitioning. I'll cross reference it with mine. Do you hear any artifacts with this implementation? $\endgroup$
    – Juan F.
    Sep 27 at 23:43
  • 2
    $\begingroup$ @JuanF. No, as long as your impulse response is time-invariant, the results is identical with time-domain convolution and is artifact-free. $\endgroup$
    – ZR Han
    Sep 28 at 1:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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