我有以下代码,它使用Swift和Accelerate在坡道上执行DFT
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)输出为
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)我相信这是正确的。以上是一个复杂到复杂的配方。基本上检查下面的代码,这是一个实数到复数的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)它的输出是
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)所以在我看来,条目0中的压缩输出应该是120,-8,但它不是。有人对这里发生的事情有什么建议吗?我正在研究这些函数,所以我很容易出错,但输出的大部分是正确的,以至于我无法理解打包的条目。
发布于 2017-07-09 10:39:36
基于一个非常仓促的阅读,我发现了两个可能的bug,这两个都在the documentation中讨论过。
具体地说,您的输入数据布局不正确(请参阅“Real FFT的数据打包”);所有的偶数项都应该在x的“实数”部分,而奇数项应该在“虚数”部分。在修复之后,我得到了:
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)更好,但这是差了两倍;如果我们看一下“傅立叶变换的缩放”,我们会看到这样的注释:
为了提供最佳的执行速度,vDSP库的函数并不总是严格遵循教科书上的傅里叶变换公式,必须相应地进行缩放。以下各节指定由vDSP库实现的每种傅立叶变换的缩放。标度因数也在参考章节中函数定义所附的公式中明确说明。
..。
实数正向变换: RF_imp = RF_math *2
所以这就是额外的两个因素的来源。
https://stackoverflow.com/questions/44988069
复制相似问题