Non-Uniformly Partitioned Convolution Implementation

Question

I've succeded in implementing the uniformly partitioned convolution algorithm and now I'm looking to implement the non-uniformly partitioned version. I've had no luck with running parallel threads on JACK audio connection kit, so I'm now shooting for a single thread implementation. The somewhat obvious solution would be to create two input buffers to work at different frequencies (meaning that one runs every 128 samples and the other one runs every 1024 samples, for example) but this will mean that the processor load will spike when both buffers are filled up. This problem is described here (page 162).

How would we go about scheduling the sub-convolutions in a single thread process so that the DSP load remains constant? Here is the starting code I'm using for the UPOLS convolution. It gives me about 10% DSP load with a 2048 samples long IR at a 48000hz sampling rate, so there's no way I can run a 2 second long IR with it. Any help (even if it's only optimizing the UPOLS code) will be greatly appreciated.

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
        buffer[nframes + i] = in[i];
        i_time[i] = buffer[i];
        i_time[nframes+i] = buffer[nframes+i];
}

// take the FFT of the input:
fftw_execute(i_forward);

// circular shift of the frequency delay line:
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,
// reset o_fft to zero to erase the previous calculations
for (i = 0; i < two_nframes; i++){
        fdl[0][i] = i_fft[i];
        o_fft[i] = 0.0 + I*0.0;
}

// multiply-add the frequency domain line (fdl) with
// the frequency domain ir partitions (fir)
for (i = 0; i < two_nframes; i++){
        for (k = 0; k < partitions; k++){
                o_fft[i] += fdl[k][i] * fir[k][i];
        }
}

// take the ifft.
fftw_execute(o_inverse);

// output the right half of the ifft, discard the rest.
for (i = 0; i < nframes; i++){
        out[i] = vol*creal(o_time[nframes+i])/two_nframes;

        // shift the input buffer to the left.
        buffer[i] = in[i];
}
return 0;

}

Answer 1

I don't think it is possible to make the CPU load remain constant for non-uniformly partitioned convolution, but a more favorable timing-dependency can be achieved.

Section 6.5 from Wefers demonstrates an optimization method to find the optimal filter partition. It is based on the minimum-load partitions and introduces some restrictions to smooth the CPU load distribution. Essentially, it is to avoid the partitions being time aligned with each other.

The optimal partition gives a slightly higher cost per sample compared with the minimum-load partition but a much more practical load distribution.

The optimal filter partition is related to the length of impulse response, the block size, and the computational costs of each operation. Thus you need to test the computational costs per sample of FFT, IFFT, complex and real multiplication and add on your target platform and you may get a different result from the author.

Answer 2

This is really a tough nut to crack. If it's a single thread, for every output sample, you are adding the results of several FIR filters running in parallel. Some FIR filters are for short impulse responses and shorter buffers. When a sample is output, you must check each filter to see if it has work to do. You do the shorter FIRs first. When a shorter FIR is done, then you have to check if the FIR next in length has work to do.