为啥 Swift 中的 FFT 与 Python 中的不同?
Posted
技术标签:
【中文标题】为啥 Swift 中的 FFT 与 Python 中的不同?【英文标题】:Why is FFT different in Swift than in Python?为什么 Swift 中的 FFT 与 Python 中的不同? 【发布时间】:2018-08-11 22:50:31 【问题描述】:我正在尝试使用 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]
关于我应该做些什么不同的任何想法?
【问题讨论】:
【参考方案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]
【讨论】:
行得通! vDSP_fft_zrip 是否相当于 fft.rfft? 不完全是。我们预计vDSP_fft_zrip
与 vDSP_fft_zip
的工作方式相同,所有虚部均为 0.0,但 vDSP_fft_zrip
无法按预期工作。看来vDSP_fft_zrip
和fft.rfft
使用了不同的优化方法,所以很难得到完全相同的结果。
有趣的是,看起来 fft.rfft 只保留前半部分加上一个日志二的值。因此,对于 4,它保留前 3,对于 8,它保留前 5,对于 16,它保留前 9,依此类推,因为随后的数字有点抵消。这是来自 scipy 的注释“当针对纯实输入计算 DFT 时,输出是 Hermitian 对称的,即负频率项只是相应正频率项的复共轭,因此负频率项是多余的。”
是的。 vDSP_fft_zrip
对其结果进行了类似的截断。但它们有点不同。也许需要一些缩放和重新排列,但我现在不能说一些精确的东西。一些 fft 专家可能知道如何正确使用vDSP_fft_zrip
。给定所有实数时,它应该更有效。【参考方案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 的数据打包 缩放傅里叶变换【讨论】:
非常感谢您的回答。您能否更新它以包括我们如何也可以反转该过程(FFT_INVERSE)?过去三天我一直在互联网上搜索,但没有任何运气。 感谢@Ibrahim!很高兴这对你有用。考虑将您的反向查询作为一个新问题提交,如果可以,请标记我 - 或在此处回复新问题的链接。 感谢@lentil 的好意!我终于完成了准备和发布我的问题。我决定在 Swift 中询问iSTFT
。 FFT_INVERSE
现在是问题的一部分。如果你有使用 Swift 的 iSTFT
的经验,请考虑检查一下。链接:***.com/questions/66847092/…以上是关于为啥 Swift 中的 FFT 与 Python 中的不同?的主要内容,如果未能解决你的问题,请参考以下文章
为啥我的 FFT 提供的可视化工具输出与 Windows Media Player 不同?
为啥 FFT 的计算时间比 intel MKL 中的元素到元素乘法要短?