Reimplement vDSP_deq22 for Biquad IIR Filter by hand

Question

I'm porting a filterbank that currently uses the Apple-specific (Accelerate) vDSP function vDSP_deq22 to Android (where Accelerate is not available). The filterbank is a set of bandpass filters that each return the RMS magnitude for their respective band. Currently the code (ObjectiveC++, adapted from NVDSP) looks like this:

- (float) filterContiguousData: (float *)data numFrames:(UInt32)numFrames channel:(UInt32)channel {

    // Init float to store RMS volume
    float rmsVolume = 0.0f;

    // Provide buffer for processing
    float tInputBuffer[numFrames + 2];
    float tOutputBuffer[numFrames + 2];

    // Copy the two frames we stored into the start of the inputBuffer, filling the rest with the current buffer data 
    memcpy(tInputBuffer, gInputKeepBuffer[channel], 2 * sizeof(float));
    memcpy(tOutputBuffer, gOutputKeepBuffer[channel], 2 * sizeof(float));
    memcpy(&(tInputBuffer[2]), data, numFrames * sizeof(float));

    // Do the processing
    vDSP_deq22(tInputBuffer, 1, coefficients, tOutputBuffer, 1, numFrames);
    vDSP_rmsqv(tOutputBuffer, 1, &rmsVolume, numFrames);

    // Copy the last two data points of each array to be put at the start of the next buffer.
    memcpy(gInputKeepBuffer[channel], &(tInputBuffer[numFrames]), 2 * sizeof(float));
    memcpy(gOutputKeepBuffer[channel], &(tOutputBuffer[numFrames]), 2 * sizeof(float));

    return rmsVolume;
}

As seen here, deq22 implements a biquad filter on a given input vector via a recursive function. This is the description of the function from the docs:

• A =: Single-precision real input vector
  
• IA =: Stride for A.
  
• B =: 5 single-precision inputs (filter coefficients), with stride 1.
  
• C =: Single-precision real output vector.
  
• IC =: Stride for C.
  
• N =: Number of new output elements to produce.

This is what I have so far (it's in Swift, like the rest of the codebase that I already have running on Android):

// N is fixed on init to be the same size as buffer.count, below
// 'input' and 'output' are initialised with (N+2) length and filled with 0s

func getFilteredRMSMagnitudeFromBuffer(var buffer: [Float]) -> Float {
    let inputStride = 1 // hardcoded for now
    let outputStride = 1

    input[0] = input[N]
    input[1] = input[N+1]
    output[0] = output[N]
    output[1] = output[N+1]

    // copy the current buffer into input
    input[2 ... N+1] = buffer[0 ..< N]

    // Not sure if this is neccessary, just here to duplicate NVDSP behaviour:
    output[2 ... N+1] = [Float](count: N, repeatedValue: 0)[0 ..< N]

    // Again duplicating NVDSP behaviour, can probably just start at 0:
    var sumOfSquares = (input[0] * input[0]) + (input[1] * input[1])

    for n in (2 ... N+1) {
        let sumG = (0...2).reduce(Float(0)) { total, p in
            return total + input[(n - p) * inputStride] * coefficients[p]
        }

        let sumH = (3...4).reduce(Float(0)) { total, p in
            return total + output[(n - p + 2) * outputStride] * coefficients[p]
        }

        let filteredFrame = sumG - sumH
        output[n] = filteredFrame
        sumOfSquares = filteredFrame * filteredFrame
    }

    let meanSquare = sumOfSquares / Float(N + 2) // we added 2 values by hand, before the loop
    let rootMeanSquare = sqrt(meanSquare)
    return rootMeanSquare
}

The filter gives a different magnitude output to deq22 though, and seems to have a cyclic wurbling circular 'noise' in it (with a constant input tone, that frequency's magnitude pumps up and down).

I've checked to ensure the coefficients arrays are identical between each implementation. Each filter actually seems to "work" in that it picks up the correct frequency (and only that frequency), it's just this pumping, and that the RMS magnitude output is a lot quieter than vDSP's, often by orders of magnitude:

   Naive    |    vDSP
3.24305e-06   0.000108608
1.57104e-06   5.53645e-05
1.96445e-06   4.33506e-05
2.05422e-06   2.09781e-05
1.44778e-06   1.8729e-05
4.28997e-07   2.72648e-05

Can anybody see an issue with my logic?

Edit: here is a gif video of the result with a constant 440Hz tone. The various green bars are the individual filter bands. The 3rd band (shown here) is the one tuned to 440Hz.

The NVDSP version just shows a constant (non-fluctuating) magnitude reading proportional to the input volume, as expected.

Answer 1

Ok, the line sumOfSquares = filteredFrame * filteredFrame should be a +=, not an assignment. So only the last frame was being calculated, explains a lot ;)

Feel free to use this if you want to do some biquad filtering in Swift. MIT Licence like NVDSP before it.