2015-11-08 48 views
1

我正在使用Accelerate框架中的vDSP在來自網格網格的二維數組中執行fft2d操作。vDSP FFT2d Swift錯誤結果的想象部分

問題是我得到一個虛數部分的數組0與Python中使用pylab.fft2的相同操作不匹配。

如果我增加數組的大小,結果不是零,但不匹配,所以我做的不好。

有人能幫我一個忙嗎?這是我的第一個堆棧溢出問題,但現在我陷入了兩個星期。

這是網格(4×8在這個例子中)

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

我將不勝感激任何有關此信息,我可以改變這一點到C或Objective-C如果最容易讓它像python一樣工作。

斯威夫特結果:

[ 
    [(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

這裏是關於相同的代碼℃的版本: http://pastebin.com/C9RPgu68

這裏是Python代碼: http://pastebin.com/rr4e6rku

+1

您是否在討論像第二項中的'0.0.i'(Swift)與'+ 1.98409020e-17j'(Python)的區別?這似乎在64位浮點數的精度範圍內。不同的輸出甚至可能是由Swift使用定點表示和Python使用科學表示造成的。 –

+0

我想在此之前,但如果我做一個4x32數組例如我得到虛部的值,但他們是不正確的(不是0,但不正確)。無論如何,非常感謝我的信息,我會再次檢查。 –

+0

也許我需要使用vDSP_fft2d_zripD而不是vDSP_fft2d_zipD,因爲輸入數據是真實的,而不是複雜的格式,但我有問題需要爲vDSP正確打包。 –

回答

3

的不同的輸出如

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

並不表示錯誤的結果。首先,Swift代碼顯然是 使用定點表示法,因此1.98409020 -17四捨五入爲0.0。其次,即使您期望 結果完全爲零,但由於二進制浮點數的精度有限(對於64位Double,約16位十進制數),所以會出現一個非零的小數值。