For vDSP_fft_zropD the packed results seem to be wrong. Is this a bug or am I doing something wrong

Question

I have the following code which does a DFT on a ramp using Swift and Accelerate

import Foundation
import Accelerate
let N = 16
var xdtar = UnsafeMutablePointer<Double>.allocate(capacity: N)
var xdtai = UnsafeMutablePointer<Double>.allocate(capacity: N)
xdtar.initialize(to: 0.0, count: N )
xdtai.initialize(to: 0.0, count: N )
var x = DSPDoubleSplitComplex(realp: xdtar, imagp: xdtai)
var ydtar = UnsafeMutablePointer<Double>.allocate(capacity: N)
var ydtai = UnsafeMutablePointer<Double>.allocate(capacity: N)
ydtar.initialize(to: 0.0, count: N )
ydtai.initialize(to: 0.0, count: N )
var y = DSPDoubleSplitComplex(realp: ydtar, imagp: ydtai)
for i in 0..<N {
    xdtar[i] = Double(i)
}
let setup = vDSP_DFT_zop_CreateSetupD(nil, vDSP_Length(N),
       vDSP_DFT_Direction.FORWARD)
vDSP_DFT_ExecuteD(setup!, x.realp, x.imagp, y.realp, y.imagp)
vDSP_DFT_DestroySetupD(setup)
for i in 0..<N {
    print(String(format: "%2d \t-> in \t(%5.4f, %5.4fi)\t out \t(%5.4f, %5.4fi)",
                i, xdtar[i], xdtai[i], ydtar[i], ydtai[i]))
}
xdtar.deinitialize(count: N); xdtar.deallocate(capacity: N)
xdtai.deinitialize(count: N); xdtai.deallocate(capacity: N)
ydtar.deinitialize(count: N); ydtar.deallocate(capacity: N)
ydtai.deinitialize(count: N); ydtai.deallocate(capacity: N)

The output is

  0     -> in   (0.0000, 0.0000i)    out    (120.0000, 0.0000i)
  1     -> in   (1.0000, 0.0000i)    out    (-8.0000, 40.2187i)
  2     -> in   (2.0000, 0.0000i)    out    (-8.0000, 19.3137i)
  3     -> in   (3.0000, 0.0000i)    out    (-8.0000, 11.9728i)
  4     -> in   (4.0000, 0.0000i)    out    (-8.0000, 8.0000i)
  5     -> in   (5.0000, 0.0000i)    out    (-8.0000, 5.3454i)
  6     -> in   (6.0000, 0.0000i)    out    (-8.0000, 3.3137i)
  7     -> in   (7.0000, 0.0000i)    out    (-8.0000, 1.5913i)
  8     -> in   (8.0000, 0.0000i)    out    (-8.0000, 0.0000i)
  9     -> in   (9.0000, 0.0000i)    out    (-8.0000, -1.5913i)
 10     -> in   (10.0000, 0.0000i)   out    (-8.0000, -3.3137i)
 11     -> in   (11.0000, 0.0000i)   out    (-8.0000, -5.3454i)
 12     -> in   (12.0000, 0.0000i)   out    (-8.0000, -8.0000i)
 13     -> in   (13.0000, 0.0000i)   out    (-8.0000, -11.9728i)
 14     -> in   (14.0000, 0.0000i)   out    (-8.0000, -19.3137i)
 15     -> in   (15.0000, 0.0000i)   out    (-8.0000, -40.2187i)

Which I believe to be correct. The above was a complex to complex formulation. Basically checking the code below which is a real to complex FFT.

import Foundation
import Accelerate
let N = 16
let log2N = vDSP_Length(4)
var xdtar = UnsafeMutablePointer<Double>.allocate(capacity: N)
xdtar.initialize(to: 0.0, count: N )
var x = DSPDoubleSplitComplex(realp: xdtar, imagp: xdtar + N/2)
var ydtar = UnsafeMutablePointer<Double>.allocate(capacity: N/2)
var ydtai = UnsafeMutablePointer<Double>.allocate(capacity: N/2)
ydtar.initialize(to: 0.0, count: N )
ydtai.initialize(to: 0.0, count: N )
var y = DSPDoubleSplitComplex(realp: ydtar, imagp: ydtai)
for i in 0..<N {
    xdtar[i] = Double(i)
}
let setup = vDSP_create_fftsetupD(log2N, Int32(2))
vDSP_fft_zropD(setup!, &x, vDSP_Stride(1), &y, vDSP_Stride(1), log2N, Int32(1))
vDSP_destroy_fftsetupD(setup)
for i in 0..<N/2 {
    print(String(format: "%2d \t-> in \t(%5.4f, %5.4fi)\t out \t(%5.4f, %5.4fi)",
                 i, x.realp[i], x.imagp[i], ydtar[i], ydtai[i]))
}
xdtar.deinitialize(count: N); xdtar.deallocate(capacity: N)
ydtar.deinitialize(count: N/2); ydtar.deallocate(capacity: N/2)
ydtai.deinitialize(count: N/2); ydtai.deallocate(capacity: N/2)

It's output is

 0  -> in   (0.0000, 8.0000i)    out    (240.0000, -128.0000i)
 1  -> in   (1.0000, 9.0000i)    out    (-8.0000, 40.2187i)
 2  -> in   (2.0000, 10.0000i)   out    (-8.0000, 19.3137i)
 3  -> in   (3.0000, 11.0000i)   out    (-8.0000, 11.9728i)
 4  -> in   (4.0000, 12.0000i)   out    (-8.0000, 8.0000i)
 5  -> in   (5.0000, 13.0000i)   out    (-8.0000, 5.3454i)
 6  -> in   (6.0000, 14.0000i)   out    (-8.0000, 3.3137i)
 7  -> in   (7.0000, 15.0000i)   out    (-8.0000, 1.5913i)

So it seems to me the packed output in entry zero should be 120, -8 but it is not. Anybody have any suggestions as to what is going on here? I am exploring these functions so I could easily have an error but enough of the output is correct that I don't understand the packed entries.

score 1 · Answer 1 · answered Jul 09 '17 at 02:39

Based on a very hasty reading, I see two probable bugs, both of which are discussed in the documentation.

Specifically, your input data layout is incorrect (see "Data packing for Real FFTs"); all the even terms should be in the "real" part of x and the odd terms in the "imaginary" part. After fixing that, I get:

0 -> in (0.0000, 1.0000i)    out (240.0000, -16.0000i)
1 -> in (2.0000, 3.0000i)    out (-16.0000, 80.4374i)
2 -> in (4.0000, 5.0000i)    out (-16.0000, 38.6274i)
3 -> in (6.0000, 7.0000i)    out (-16.0000, 23.9457i)
4 -> in (8.0000, 9.0000i)    out (-16.0000, 16.0000i)
5 -> in (10.0000, 11.0000i)  out (-16.0000, 10.6909i)
6 -> in (12.0000, 13.0000i)  out (-16.0000, 6.6274i)
7 -> in (14.0000, 15.0000i)  out (-16.0000, 3.1826i)

Better, but that's off by a factor of two; if we look at "Scaling for Fourier Transforms", we see this note:

To provide the best possible execution speeds, the vDSP library's functions don't always adhere strictly to textbook formulas for Fourier transforms, and must be scaled accordingly. The following sections specify the scaling for each type of Fourier transform implemented by the vDSP Library. The scaling factors are also stated explicitly in the formulas that accompany the function definitions in the reference chapter.

...

Real forward transforms: RF_imp = RF_math * 2

So that's where the extra factor of two is coming from.

For vDSP_fft_zropD the packed results seem to be wrong. Is this a bug or am I doing something wrong

1 Answers1