vDSP_conv() 结果不正确

Incorrect results with vDSP_conv()

与 MATLAB 实现相比,尝试使用来自 Accelerate 的 vDSP_conv() 进行卷积时,我得到的结果不一致。在使用此函数计算卷积时,有几篇关于奇怪结果的 Whosebug 帖子,但据我所知,我正确使用了该框架,并采纳了其他 Stack Overflow 帖子的建议。这是我的代码:

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))
}

据我所知,我正在执行 Accelerate 框架文档中指定的正确零填充 here

我定义了两个测试数组A: [Float] = [0, 0, 1, 0, 0]B: [float] = [1, 0, 0]

在 MATLAB 中,当我 运行 conv(A, B) 时,我得到 [0, 0, 1, 0, 0, 0, 0]

然而,当我 运行 上面的 vDSP conv() 我得到,[1, 0, 0, 0, 0, 0, 0]

我的实现有什么问题?我已经多次检查过这个问题并查看了我能找到的所有 SO 帖子,但仍然无法解释这种不一致。

除此之外,是否有比我这里有的更有效的方法来对数组进行零填充?为了保持 x 不可变,我创建了新的 xPadded 数组,但毫无疑问还有更有效的方法来执行此填充。

** 编辑 ** 正如 Martin R 所建议的那样,我在数组的开头和结尾平等地填充了 k.count -1,如下所示。

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
}

使用此代码,conv(A, B) 仍然 returns [1, 0, 0, 0, 0, 0, 0]。

我调用的函数如下所示:

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

对于长度为 mn 的两个数组 AB, Accelerate 框架中的 vDSP_conv() 函数计算长度为 m - n + 1.

的新数组

这对应于形状为 conv() 的 MATLAB 函数的结果 参数设置为 "valid":

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

要获得与来自 MATLAB 的 "full" 卷积相同的结果 您必须在 A 数组的开头和结尾用 n-1 元素进行零填充,这会给出长度为 m + n - 1.

的结果数组

应用于您的函数:

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

使用 Repeat() 可能 稍微 更高效,因为它 创建一个序列而不是一个数组。但最终,一个新的阵列 必须创建为 vDSP_conv() 函数的参数, 所以没有太大的改进空间。

@MartinR 我想通了为什么我的代码不适用于数组。我在一个使用 Surge 作为链接框架的项目中编写这段代码。 Surge 为 [Float] 和 [Double] 数组重载 + 运算符,使其成为数组元素的逐元素相加。因此,当我执行 x + xPad 时,它并没有按预期扩展数组的大小,它只是返回 x,因为 xPad 仅包含零。但是,Surge 没有为序列重载 + 运算符,因此使用 Repeat() 成功扩展了数组。感谢您的帮助 - 从来没有想过要尝试序列!

对下一个遇到这个问题的可怜人的一些澄清: Apple 提供了 some sample code 如何使用 vDSP_conv,但它毫无用处。事实上,这让我感到困惑,因为该代码中的注释说需要填充输入缓冲区,而没有指定实际输入样本应该放置的位置:

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;

所以,上面的公式给你一个不同于 xPad + x + xPad 的长度(更大),其中 xPad 是 k.count - 1.

重要的是在填充缓冲区中复制输入(信号)样本的位置:它需要位于 k.count - 1

所以,上面接受的解决方案有效。但是如果你相信 Apple 示例中的评论(顺便说一句,官方文档中没有出现)那么你可以做出妥协:使用他们的公式(上面的 SignalLength)来计算和分配填充缓冲区(这会有点更大)并使用 k.count - 1(即滤波器长度 - 1)作为信号的起始偏移量(在本例中为 x)。我这样做了,结果现在匹配 ippsConvolve_32f 和 Matlab。

(抱歉,这应该是评论,但我没有足够的声誉)。