Skip to main content
Tweeted twitter.com/StackSignals/status/1449253822874931215
Added the working code for reference.
Source Link

EDIT: Scroll down for actual working code.

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;

}

EDIT: Scroll down for actual working code.

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;

}

Source Link

Real Time Partitioned Convolution Not Working

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!