Swift vDSP rFFT 与 Python np.fft.rfft() 不同答案

【问题标题】：Swift vDSP rFFT not same as Python np.fft.rfft()Swift vDSP rFFT 与 Python np.fft.rfft() 不同
【发布时间】：2021-03-21 20:21:23
【问题描述】：

我正在尝试在 iOS 中实现真正的 FFT，因为我正在使用 Accelerate Framework。这是我的 Swift 代码。

class FFT {

private var fftSetup: FFTSetup?
private var log2n: Float?
private var length: Int?

func initialize(count: Int){
    length = count
    log2n = log2(Float(length!))
    self.fftSetup = vDSP_create_fftsetup(vDSP_Length(log2n!), FFTRadix(kFFTRadix2))!
}

func computeFFT(input: [Float]) -> ([Float], [Float]) {
    
    var real = input
    var imag = [Float](repeating: 0.0, count: input.count)
    var splitComplexBuffer = DSPSplitComplex(realp: &real, imagp: &imag)

    let halfLength = (input.count/2) + 1
    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))
    
    return (real, imag)
    
}

func computeIFFT(real: [Float], imag: [Float]) -> [Float]{
    
    var real = [Float](real)
    var imag = [Float](imag)
    
    var result : [Float] = [Float](repeating: 0.0, count: length!)
    var resultAsComplex : UnsafeMutablePointer<DSPComplex>? = nil

    result.withUnsafeMutableBytes {
        resultAsComplex = $0.baseAddress?.bindMemory(to: DSPComplex.self, capacity: 512)
    }
    
    var splitComplexBuffer = DSPSplitComplex(realp: &real, imagp: &imag)
    
    vDSP_fft_zrip(fftSetup!, &splitComplexBuffer, 1, vDSP_Length(log2n!), FFTDirection(FFT_INVERSE));
    
    vDSP_ztoc(&splitComplexBuffer, 1, resultAsComplex!, 2, vDSP_Length(length! / 2));
    //
    //// Neither the forward nor inverse FFT does any scaling. Here we compensate for that.
    var scale : Float = 1.0/Float(length!);
    var copyOfResult = result;
    vDSP_vsmul(&result, 1, &scale, &copyOfResult, 1, vDSP_Length(length!));
    result = copyOfResult
    
    return result
}

func deinitialize(){
    vDSP_destroy_fftsetup(fftSetup)
}

}

这是我用于计算 rFFT 和 irFFT 的 Python 代码

# calculate fft of input block
    in_block_fft = np.fft.rfft(np.squeeze(in_buffer)).astype("complex64")

# apply mask and calculate the ifft
    estimated_block = np.fft.irfft(in_block_fft * out_mask)

问题：

斯威夫特 如果我计算 512 帧的 rFFT 并将 irFFT 应用于 rFFT 的结果，我会得到相同的原始数组。

Python python 也是如此，如果我采用 rFFT 和 irFFT，我会得到原始数组作为回报。

问题如果我比较 Swift rFFT 和 Python rFFT 的结果，就会出现问题。他们的结果在十进制值上是不同的。有时实部是一样的，但虚部是完全不同的。

我在 Python 中尝试了不同的框架，如 Numpy、SciPy 和 TensorFlow，它们的结果完全相同（小数部分略有不同）。但是当我使用上面的 Swift 代码在 iOS 的相同输入上计算 rfft 时，结果是不同的。

如果任何有 Accelerate Framework 经验并掌握 FFT 知识的人帮助我解决这个问题，将会非常有帮助。我对 FFT 的了解有限。

【问题讨论】：

标签： python swift numpy fft vdsp

【解决方案1】：

是的，它有区别。

我使用相同的数据将此功能从 python librosa 部分移植到 iOS。

https://github.com/dhrebeniuk/RosaKit

【讨论】：

【解决方案2】：

两个观察结果：首先，您对 DSPSplitComplex(realp: &real, imagp: &imag) 的调用会创建临时指针，而 fromInputArray 初始化程序已被弃用。查看Data Packing for Fourier Transforms 以获得最佳实践。

其次，bindMemory 容量是元素的数量。以下几行创建了 512 个复杂元素，它们应该（如果我理解正确的话）创建 256 个复杂元素来表示 512 个真实元素。

        result.withUnsafeMutableBytes {
            resultAsComplex = $0.baseAddress?.bindMemory(to: DSPComplex.self, capacity: 512)
        }

西蒙

【讨论】：