vDSP FFT2d Swift 结果虚部错误答案

【问题标题】：vDSP FFT2d Swift wrong imaginary part on the resultvDSP FFT2d Swift 结果虚部错误
【发布时间】：2016-02-09 08:25:13
【问题描述】：

我正在使用 Accelerate 框架中的 vDSP 在来自网格的二维数组中执行 fft2d 操作。

问题是我在虚部得到了一个 0 数组，它与 python 中使用 pylab.fft2 的相同操作不匹配。

如果我增加数组大小，结果不会为零，但无论如何都不匹配，所以我做错了。

谁能帮我一把？这是我的第一个堆栈溢出问题，但我现在被困了两个星期。

这是网格（本例为 4x8）

[
    [1.80485138784544e-35, 2.61027906966774e-23, 1.26641655490943e-14, 2.06115362243857e-09, 1.1253517471926e-07, 2.06115362243857e-09, 1.26641655490943e-14, 2.61027906966774e-23],
    [2.93748211171084e-30, 4.24835425529162e-18, 2.06115362243857e-09, 0.000335462627902512, 0.0183156388887342, 0.000335462627902512, 2.06115362243857e-09, 4.24835425529162e-18],
    [1.60381089054866e-28, 2.31952283024359e-16, 1.1253517471926e-07, 0.0183156388887342, 1.0, 0.0183156388887342, 1.1253517471926e-07, 2.31952283024359e-16],
    [2.93748211171084e-30, 4.24835425529162e-18, 2.06115362243857e-09, 0.000335462627902512, 0.0183156388887342, 0.000335462627902512, 2.06115362243857e-09, 4.24835425529162e-18]
]

这里是fft2函数：

func fft2(arr: [[Complex<Double>]]) -> [[Complex<Double>]] {
    let nRows = arr.count
    let nCols = arr[0].count
    let N = nRows * nCols

    let radix = FFTRadix(FFT_RADIX2)
    let pass = vDSP_Length(Int(log2(Double(N))))

    // Create FFTSetup
    let setup = vDSP_create_fftsetupD(pass, radix)

    // Direction
    let dir = FFTDirection(FFT_FORWARD)

    // Get real and imag doubles from the [Complex]
    // (all imag parts are 0.0 on this example)
    let (real, imag) = complex2DArrayToDouble(arr)

    // Pack 2d arrays as 1d (function bellow)
    var realArray = pack2dArray(real, rows: nRows, cols: nCols)
    var imagArray = pack2dArray(imag, rows: nRows, cols: nCols)

    // Create the split complex with the packed arrays
    var splitComplex = DSPDoubleSplitComplex(
        realp: &realArray,
        imagp: &imagArray)

    let log2n0c = vDSP_Length(Int(log2(Double(nCols))))
    let log2n1r = vDSP_Length(Int(log2(Double(nRows))))

    let rowStride = vDSP_Stride(nRows)
    let colStride = vDSP_Stride(1) // Use all cols

    // Perform the fft2d
    vDSP_fft2d_zipD(setup, &splitComplex, rowStride, colStride, log2n0c, log2n1r, dir)

    // Destroy setup
    vDSP_destroy_fftsetupD(setup)

    // Pack the 1d arrays on 2d arrays again
    let resultReal = unpack2dArray(realArray, rows: nRows, cols: nCols)
    let resultImag = unpack2dArray(imagArray, rows: nRows, cols: nCols)

    // Ignore this...
    return complexFrom2DArray([[Double]](), imag: [[Double]]())
}

最后，这里是我用于将数组从/到 2d 到 1d 打包和解包的函数

func pack2dArray(arr: [[Double]], rows: Int, cols: Int) -> [Double] {
    var resultArray = zeros(rows * cols)
    for Iy in 0...cols-1 {
        for Ix in 0...rows-1 {
            let index = Iy * rows + Ix
            resultArray[index] = arr[Ix][Iy]
        }
    }
    return resultArray
}

func unpack2dArray(arr: [Double], rows: Int, cols: Int) -> [[Double]] {
    var resultArray = [[Double]](count: rows, repeatedValue: zeros(cols))
    for Iy in 0...cols-1 {
        for Ix in 0...rows-1 {
            let index = Iy * rows + Ix
            resultArray[Ix][Iy] = arr[index]
        }
    }
    return resultArray
}

我会很感激有关此的任何信息，如果最容易使其像在 python 中一样工作，我可以将其更改为 C 或 Objective-C。

快速结果：

[
    [(1.07460475603902+0.0.i), (-1.06348244974363+0.0.i), (1.03663115699765+0.0.i), (-1.00978033088166+0.0.i), (0.998658491216246+0.0.i), (-1.00978033088166+0.0.i), (1.03663115699765+0.0.i), (-1.06348244974363+0.0.i)],
    [(-1.03663138619031+0.0.i), (1.02590210946989+0.0.i), (-0.999999662394501+0.0.i), (0.974097665459761+0.0.i), (-0.963368838879988+0.0.i), (0.974097665459761+0.0.i), (-0.999999662394501+0.0.i), (1.02590210946989+0.0.i)],
    [(0.998658482971633+0.0.i), (-0.988322230996495+0.0.i), (0.963368617931946+0.0.i), (-0.938415438518917+0.0.i), (0.928079620195301+0.0.i), (-0.938415438518917+0.0.i), (0.963368617931946+0.0.i), (-0.988322230996495+0.0.i)],
    [(-1.03663138619031+0.0.i), (1.02590210946989+0.0.i), (-0.999999662394501+0.0.i), (0.974097665459761+0.0.i), (-0.963368838879988+0.0.i), (0.974097665459761+0.0.i), (-0.999999662394501+0.0.i), (1.02590210946989+0.0.i)]
]

Python 结果：

[
    [ 1.07460476 +0.00000000e+00j, -1.06348245 +1.98409020e-17j, 1.03663116 +0.00000000e+00j -1.00978033 -1.97866921e-17j, 0.99865849 +0.00000000e+00j -1.00978033 -1.98409020e-17j, 1.03663116 +0.00000000e+00j -1.06348245 +1.97866921e-17j]
    [-1.03663139 +0.00000000e+00j, 1.02590211 -1.90819560e-17j, -0.99999966 +0.00000000e+00j, 0.97409767 +1.90819558e-17j, -0.96336884 +0.00000000e+00j, 0.97409767 +1.90819560e-17j, -0.99999966 +0.00000000e+00j, 1.02590211 -1.90819558e-17j]
    [ 0.99865848 +0.00000000e+00j, 0.98832223 +1.83230190e-17j, 0.96336862 +0.00000000e+00j, 0.93841544 -1.83772293e-17j, 0.92807962 +0.00000000e+00j, 0.93841544 -1.83230190e-17j, 0.96336862 +0.00000000e+00j, 0.98832223 +1.83772293e-17j]
    [-1.03663139 +0.00000000e+00j, 1.02590211 -1.90819560e-17j, -0.99999966 +0.00000000e+00j, 0.97409767 +1.90819558e-17j, -0.96336884 +0.00000000e+00j, 0.97409767 +1.90819560e-17j, -0.99999966 +0.00000000e+00j, 1.02590211 -1.90819558e-17j]
]

提前致以诚挚的谢意！

编辑 1

这是相同代码的 C 版本： http://pastebin.com/C9RPgu68

这里是python代码： http://pastebin.com/rr4e6rku

【问题讨论】：

您是在谈论第二个条目中的 0.0.i (Swift) 与 +1.98409020e-17j (Python) 之类的区别吗？这似乎在 64 位浮点数的精度范围内。不同的输出甚至可能是由于 Swift 使用定点表示而 Python 使用科学表示。
我之前想过这个问题，但是如果我制作一个 4x32 的数组，例如，我会得到虚部的值，但它们不正确（不是 0，而是不正确）。无论如何，非常感谢您提供的信息，我会再次检查。
也许我需要使用 vDSP_fft2d_zripD 而不是 vDSP_fft2d_zipD，因为输入数据是真实的，而不是复杂格式，但我在为 vDSP 正确打包时遇到问题。
嗯，也许我只是愚蠢而你是对的@Martin R ...哈哈。将其写为答案，我会将其标记为有效。也许我的错是真的痴迷于获得与 python 中完全相同的结果，并且不认为 1.984...-17j 是一个非常薄的数字。谢谢

标签： c swift fft accelerate-framework vdsp

【解决方案1】：

不同的输出如

Swift:  (-1.06348244974363+0.0.i)
Python: -1.06348245 +1.98409020e-17j

不表示错误的结果。首先，Swift 代码显然使用定点表示，因此 1.98409020 ⋅ 10^-17 舍入为0.0。其次，即使你期望结果正好为零，由于精度有限，预计会有一个小的非零值二进制浮点数（对于 64 位 Double，大约有 16 个十进制数字）。

【讨论】：