我使用本地Apple类实现了FFT算法。 我直接从他们的网站上获取了代码:
然而,每次运行代码时,它都提供不同的结果。我创建了一个单元测试,重复运行代码并比较结果,如果单元测试失败,则说明结果不同。我唯一的猜测是这是一个内存问题。这是我想象中唯一可能导致结果每次不同的方式。
import Foundation
import Accelerate
class AppleFFT{
var windowSize = 512
var n = vDSP_Length(512)
var halfN = Int(512 / 2)
var fftSetUp : FFTSetup?
var log2n : vDSP_Length?
init(windowSize: Int){
self.windowSize = windowSize
n = vDSP_Length(windowSize)
halfN = Int(n / 2)
initialize()
}
private init(){
initialize()
}
func initialize(){
log2n = vDSP_Length(log2(Float(n)))
if log2n == nil { return }
fftSetUp = vDSP_create_fftsetup(log2n!, FFTRadix(kFFTRadix2))
}
func process(signal : [Float], n: vDSP_Length) ->DSPSplitComplex{
let window = vDSP.window(ofType: Float.self,
usingSequence: .hanningDenormalized,
count: Int(n),
isHalfWindow: false)
let signal2 = vDSP.multiply(signal, window)
let observed: [DSPComplex] = stride(from: 0, to: Int(n), by: 2).map {
return DSPComplex(real: signal[$0],
imag: signal[$0.advanced(by: 1)])
}
var forwardInputReal = [Float](repeating: 0, count: halfN)
var forwardInputImag = [Float](repeating: 0, count: halfN)
var forwardInput = DSPSplitComplex(realp: &forwardInputReal,
imagp: &forwardInputImag)
vDSP_ctoz(observed, 2,
&forwardInput, 1,
vDSP_Length(halfN))
//Create some empty arrays we can put data into
var forwardOutputReal = [Float](repeating: 0, count: halfN)
var forwardOutputImag = [Float](repeating: 0, count: halfN)
var forwardOutput = DSPSplitComplex(realp: &forwardOutputReal,
imagp: &forwardOutputImag)
//Perform actual fft, placing results in forwardOutput
vDSP_fft_zrop(fftSetUp!,
&forwardInput, 1,
&forwardOutput, 1,
log2n!,
FFTDirection(kFFTDirection_Forward))
//Do cheap analysis to figure out original frequencies
let componentFrequencies = forwardOutputImag.enumerated().filter {
$0.element < -1
}.map {
return $0.offset
}
return forwardOutput
}
}
import XCTest
import Accelerate
class testAppleFFT: XCTestCase {
func testFFTConsistency(){
let signal = genSignalWith(frequencies:[100, 500], numSamples: 512, sampleRate: 44100)
let fft = AppleFFT(windowSize: 512)
let complex1 = fft.process(signal: signal , n: 512)
for i in 0..<10{
print("i = \(i)")
let complex2 = fft.process(signal: signal, n: 512)
var complex1realp = complex1.realp
var complex1imagp = complex1.imagp
var complex2realp = complex2.realp
var complex2imagp = complex2.imagp
for j in 0..<512 {
let r1 = complex1realp.pointee
let i1 = complex1imagp.pointee
let r2 = complex2realp.pointee
let i2 = complex2imagp.pointee
XCTAssert(abs(r1 - r2) < 0.00001)
XCTAssert(abs(i1 - i2) < 0.00001)
if !(abs(r1 - r2) < 0.00001){
print(" error: i: \(i) j: \(j) r1: \(r1) r2: \(r2)")
}
if !(abs(i1 - i2) < 0.00001){
print(" error: index: \(i) i1: \(i1) i2: \(i2)")
}
complex1realp = complex1realp.advanced(by: 1)
complex1imagp = complex1imagp.advanced(by: 1)
complex2realp = complex2realp.advanced(by: 1)
complex2imagp = complex2imagp.advanced(by: 1)
}
}
}
func genSignalWith(frequencies: [Float], numSamples: Int, sampleRate: Float, amplitudes: [Float] = []) -> [Float]{
var sig : [Float] = []
for t in 0..<numSamples{
var sum : Float = 0.0
for i in 0..<frequencies.count{
let f = frequencies[i]
var a : Float = 1.0
if(amplitudes.count > i){
a = amplitudes[i]
}
let thisValue = sin(Float(t) / sampleRate * 2 * .pi * f)
sum += thisValue
}
sig.append(sum)
}
return sig
}
}