为什么 Swift 中的 FFT 与 Python 中的不同？

Question

我正在尝试使用 Accelerate 框架将一些 python numpy 代码移植到 Swift。

我用python写

import numpy as np
frames = np.array([1.0, 2.0, 3.0, 4.0])
fftArray = np.fft.fft(frames, len(frames))
print(fftArray)

输出是：

[10.+0.j -2.+2.j -2.+0.j -2.-2.j]

所以在 Swift 中，我尝试像这样计算 FFT：

import Foundation
import Accelerate

    func fftAnalyzer(frameOfSamples: [Float]) {
        // As above, frameOfSamples = [1.0, 2.0, 3.0, 4.0]            

        let analysisBuffer = frameOfSamples
        let frameCount = frameOfSamples.count

        var reals = [Float]()
        var imags = [Float]()
        for (idx, element) in analysisBuffer.enumerated() {
            if idx % 2 == 0 {
                reals.append(element)
            } else {
                imags.append(element)
            }
        }
        var complexBuffer = DSPSplitComplex(realp: UnsafeMutablePointer(mutating: reals), imagp: UnsafeMutablePointer(mutating: imags))

        let log2Size = Int(log2f(Float(frameCount)))

        guard let fftSetup = vDSP_create_fftsetup(vDSP_Length(log2Size), Int32(kFFTRadix2)) else {
            return []
        }

        // Perform a forward FFT
        vDSP_fft_zrip(fftSetup, &(complexBuffer), 1, UInt(log2Size), Int32(FFT_FORWARD))

        let realFloats = Array(UnsafeBufferPointer(start: complexBuffer.realp, count: Int(frameCount)))
        let imaginaryFloats = Array(UnsafeBufferPointer(start: complexBuffer.imagp, count: Int(frameCount)))

        print(realFloats)
        print(imaginaryFloats)

        // Release the setup
        vDSP_destroy_fftsetup(fftSetup)

        return realFloats
    }

realFloats 和 imaginaryFloats 打印如下：

[20.0, -4.0, 0.0, 0.0]
[-4.0, 4.0, 0.0, 0.0]

关于我应该做些什么不同的任何想法？

Answer 1

我不擅长 numpy，但根据the doc，fft 需要复杂输入。那么它的等价物将是vDSP_fft_zip，而不是vDSP_fft_zrip。

您的代码会导致缓冲区溢出或可能导致悬空指针，所有这些问题都已修复，我得到这个：

func fftAnalyzer(frameOfSamples: [Float]) -> [Float] {
    // As above, frameOfSamples = [1.0, 2.0, 3.0, 4.0]

    let frameCount = frameOfSamples.count

    let reals = UnsafeMutableBufferPointer<Float>.allocate(capacity: frameCount)
    defer {reals.deallocate()}
    let imags =  UnsafeMutableBufferPointer<Float>.allocate(capacity: frameCount)
    defer {imags.deallocate()}
    _ = reals.initialize(from: frameOfSamples)
    imags.initialize(repeating: 0.0)
    var complexBuffer = DSPSplitComplex(realp: reals.baseAddress!, imagp: imags.baseAddress!)

    let log2Size = Int(log2(Float(frameCount)))
    print(log2Size)

    guard let fftSetup = vDSP_create_fftsetup(vDSP_Length(log2Size), FFTRadix(kFFTRadix2)) else {
        return []
    }
    defer {vDSP_destroy_fftsetup(fftSetup)}

    // Perform a forward FFT
    vDSP_fft_zip(fftSetup, &complexBuffer, 1, vDSP_Length(log2Size), FFTDirection(FFT_FORWARD))

    let realFloats = Array(reals)
    let imaginaryFloats = Array(imags)

    print(realFloats)
    print(imaginaryFloats)

    return realFloats
}

印刷

[10.0, -2.0, -2.0, -2.0]
[0.0, 2.0, 0.0, -2.0]

Answer 2

关于 vDSP_fft_zip 和 vDSP_fft_zrip 之间的区别存在一些混淆。在您的情况下，最直接的影响是 zrip 输出已打包，需要按 1/2 缩放以等于标准 FFT 的正常数学输出。

这在 Apple 文档中有所隐藏：而不是出现在 vDSP_fft_zrip 的页面上，而是出现在 vDSP Programming Guide 中。但是，从指南中仍然不清楚在每种情况下如何准备输入数据 - OOPer 的回答对此有很大帮助。

import Foundation
import Accelerate

var input: [Float] = [1.0, 2.0, 3.0, 4.0]

let length = input.count
let log2n = log2(Float(length))
let fftSetup = vDSP_create_fftsetup(vDSP_Length(log2n), FFTRadix(kFFTRadix2))!

print("--vDSP_fft_zip---")
var real = input
var imag = [Float](repeating: 0.0, count: length)
var splitComplexBuffer = DSPSplitComplex(realp: &real, imagp: &imag)

vDSP_fft_zip(fftSetup, &splitComplexBuffer, 1, vDSP_Length(log2n), FFTDirection(FFT_FORWARD))

print("real: \(real)")
print("imag: \(imag)")
print("-----------------")
print("DC      : real[0]")
print("Nyquist : real[2]")
print()
print()

print("--vDSP_fft_zrip--")
// only need half as many elements because output will be packed
let halfLength = length/2
real = [Float](repeating: 0.0, count: halfLength)
imag = [Float](repeating: 0.0, count: halfLength)

// input is alternated across the real and imaginary arrays of the DSPSplitComplex structure
splitComplexBuffer = DSPSplitComplex(fromInputArray: input, realParts: &real, imaginaryParts: &imag)

// even though there are 2 real and 2 imaginary output elements, we still need to ask the fft to process 4 input samples
vDSP_fft_zrip(fftSetup, &splitComplexBuffer, 1, vDSP_Length(log2n), FFTDirection(FFT_FORWARD))

// zrip results are 2x the standard FFT and need to be scaled
var scaleFactor = Float(1.0/2.0)
vDSP_vsmul(splitComplexBuffer.realp, 1, &scaleFactor, splitComplexBuffer.realp, 1, vDSP_Length(halfLength))
vDSP_vsmul(splitComplexBuffer.imagp, 1, &scaleFactor, splitComplexBuffer.imagp, 1, vDSP_Length(halfLength))

print("real: \(real)")
print("imag: \(imag)")
print("-----------------")
print("DC      : real[0]")
print("Nyquist : imag[0]")

vDSP_destroy_fftsetup(fftSetup)

印刷：

--vDSP_fft_zip---
real: [10.0, -2.0, -2.0, -2.0]
imag: [0.0, 2.0, 0.0, -2.0]
-----------------
DC      : real[0]
Nyquist : real[2]


--vDSP_fft_zrip--
real: [10.0, -2.0]
imag: [-2.0, 2.0]
-----------------
DC      : real[0]
Nyquist : imag[0]

vDSP_fft_zip

输入放在DSPSplitComplex结构的实数数组中，虚数数组清零。

输出很复杂，需要两倍于输入的内存。然而，这个输出的大部分是对称的——这就是 zrip 的打包输出能够在一半内存中表示相同信息的方式。

vDSP_fft_zrip

输入使用 DSPSplitComplex.init(fromInputArray: ) 或使用 vDSP_ctoz 分布在 DSPSplitComplex。

fromInputArray: 方法的作用与 ctoz 相同，但它是一种更简单/更安全的方法 - 不必使用 UnsafeRawPointer 和 bindMemory。

输出很复杂。打包后，输出需要与输入相同的内存量。

缩放：结果是标准数学 FFT 的 2 倍，因此需要进行缩放。

---

见： vDSP Programming Guide

• 真实 FFT 的数据打包
• 缩放傅里叶变换