Для vDSP_fft_zropD упакованные результаты кажутся неверными. Это ошибка или я делаю что-то не так

У меня есть следующий код, который делает DFT на рампе, используя Swift и 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)

Выход

  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)

Что я считаю правильным. Выше было сложным в сложной формулировке. В основном проверка кода ниже, который является реальным для сложных БПФ.

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)

Так что мне кажется, что упакованный вывод в нулевой записи должен быть 120, -8, но это не так. У кого-нибудь есть предложения относительно того, что здесь происходит? Я изучаю эти функции, чтобы у меня могла легко возникнуть ошибка, но достаточно вывода правильного, что я не понимаю упакованные записи.

Источник

1 ответ

Исходя из очень поспешного прочтения, я вижу две вероятные ошибки, обе из которых обсуждаются в документации.

В частности, неправильное расположение входных данных (см. "Упаковка данных для реальных БПФ"); все четные члены должны быть в "реальной" части 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

Так вот откуда берется дополнительный фактор в два.

Источник
Другие вопросы по тегам