Incorrect results with vDSP_conv()

Question

I am getting inconsitent results when attempting to do convolution using vDSP_conv() from Accelerate when compared to the MATLAB implementation. There have been a couple of StackOverflow posts about weird results when using this function to calculate convolution, however as far as I can tell, I am using the framework correctly and have incorporated the suggestions from the other Stack Overflow posts. Here is my code:

public func conv(x: [Float], k: [Float]) -> [Float] {  
    let resultSize = x.count + k.count - 1
    var result = [Float](count: resultSize, repeatedValue: 0)
    let kEnd = UnsafePointer<Float>(k).advancedBy(k.count - 1)
    let xPad: [Float] = [Float](count: (2*k.count)+1, repeatedValue: 0.0)
    let xPadded = x + xPad
    vDSP_conv(xPadded, 1, kEnd, -1, &result, 1, vDSP_Length(resultSize), vDSP_Length(k.count))
}

As far as I can tell, I am doing the correct zero padding as specified in the Accelerate framework documentation here

I defined two test arrays A: [Float] = [0, 0, 1, 0, 0] and B: [float] = [1, 0, 0].

In MATLAB, when I run conv(A, B), I get [0, 0, 1, 0, 0, 0, 0].

However, when I run the above vDSP conv() I get, [1, 0, 0, 0, 0, 0, 0].

What is wrong with my implementation? I have gone over this a number of times and looked through all the SO posts that I could find, and still haven't been able to account for this inconsistency.

Beyond that, is there a more efficient method to zero-pad the array then what I have here? In order to keep x immutable, I created the new xPadded array but there is undoubtedly a more efficient method of performing this padding.

** EDIT ** As suggested by Martin R, I padded k.count -1 equally at the beginning and end of the array as shown below.

public func conv(x: [Float], k: [Float]) -> [Float] {
    let resultSize = x.count + k.count - 1
    var result = [Float](count: resultSize, repeatedValue: 0)
    let kEnd = UnsafePointer<Float>(k).advancedBy(k.count - 1)
    let xPad: [Float] = [Float](count: k.count-1, repeatedValue: 0.0)
    let xPadded = xPad + x + xPad
    vDSP_conv(xPadded, 1, kEnd, -1, &result, 1, vDSP_Length(resultSize), vDSP_Length(k.count))

    return result
}

Using this code, conv(A, B) still returns [1, 0, 0, 0, 0, 0, 0].

I am calling the function as shown below:

let A: [Float] = [0, 0, 1, 0, 0]
let B: [Float] = [1, 0, 0]
let C: [Float] = conv(A, k: B)

Answer 1

For two arrays A and B of length m and n, the vDSP_conv() function from the Accelerate framework computes a new array of length m - n + 1.

This corresponds to the result of the MATLAB function conv() with the shape parameter set to "valid":

Only those parts of the convolution that are computed without the zero-padded edges. ...

To get the same result as the with "full" convolution from MATLAB you have to zero-pad the A array with n-1 elements at the beginning and the end, this gives a result array of length m + n - 1.

Applied to your function:

let xPad = Repeat(count: k.count - 1, repeatedValue: Float(0.0))
let xPadded = xPad + x + xPad 

Using Repeat() might be slightly more performant because it creates a sequence and not an array. But ultimately, a new array has to be created as an argument to thevDSP_conv() function, so there is not much room for improvement.

Answer 2

Some clarifications for the next poor soul who stumbles into this: Apple provides some sample code on how to use vDSP_conv but it's pretty useless. In fact, it was confusing me because a comment in that code says that the input buffer needs to be padded without specifying where the actual input samples should be placed:

The SignalLength defined below is used to allocate space, and it is the filter length rounded up to a multiple of four elements and added to the result length.

SignalLength = (FilterLength+3 & -4u) + ResultLength;

So, the above formula gives you a different length (bigger) than the xPad + x + xPad where xPad is the k.count - 1.

The important thing is where in that padded buffer you copy your input (signal) samples: it needs to be at k.count - 1.

So, the above accepted solution works. But if you trust that comment in Apple's example (which BTW doesn't show up in the official docs) then you can do a compromise: use their formula (the SignalLength above) to calculate and allocate the padded buffer (it will be a bit larger) and use the k.count - 1 (i.e. filter length - 1) as the starting offset for your signal (x in this case). I did this and the results now match ippsConvolve_32f and Matlab.

(Sorry, this should have been a comment but I don't have enough reputation for that).

Answer 3

@MartinR I figured out why my code doesn't work with Arrays. I was writing this code in a project that was using Surge as a linked framework. Surge overloads the + operator for [Float] and [Double] arrays so that it becomes element-wise addition of array elements. So when I was doing x + xPad it wasn't extending the size of the array as expected, it was simply returning x as xPad only contained zeros. However, Surge had not overloaded the +operator for sequences, so using Repeat() successfully extended the array. Thanks for your help - never would have thought to try sequences!