Apple FFT giving inconsistent results

Question

I implemented in FFT algorithm using native Apple classes. I took the code directly off of their website here:

https://developer.apple.com/documentation/accelerate/vdsp/fast_fourier_transforms/finding_the_component_frequencies_in_a_composite_sine_wave

Nonetheless when I run the code it provides different results every time. I created a unit test which runs it repeatedly and compares if the results are the same at the unit test fails. My only guess is that this is a memory problem. That's the only way I can imagine that the results could be different each time.

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
    }
}

Answer 1

This:

var forwardInput = DSPSplitComplex(realp: &forwardInputReal,
                                   imagp: &forwardInputImag)
vDSP_ctoz(observed, 2, &forwardInput, 1, vDSP_Length(halfN))

does not do what you want it to do. The problem with it is a little bit subtle, especially if you're coming from a C or C++ background. Arrays in swift are not like arrays in C or C++; in particular, they do not have fixed addresses in memory. They are objects that Swift may choose to move. This is nice when you're working in Swift, but sometimes causes pain when you need to interact with C functions (and especially C types that want to persist pointers across function calls, as you have noticed).

When you call DSPSplitComplex(realp: &forwardInputReal, ...), the & implicitly creates an UnsafeMutablePointer<Float> to the memory of forwardInputReal, but that pointer is only valid during the duration of the call to init. When you pass forwardInput to vDSP_ctoz, the pointer has already passed out of scope and is no longer valid, so you are invoking undefined behavior. In particular, the compiler can assume that the call to vDSP_ctoz does not modify the contents of forwardInputReal or forwardInputImag, because the function doesn't receive a valid pointer to their contents.

The best way to work around this is to be more explicit:

forwardInputReal.withUnsafeMutableBufferPointer { r in
  forwardInputImag.withUnsafeMutableBufferPointer { i in
    var splitComplex = DSPSplitComplex(realp: r.baseAddress!, imagp: i.baseAddress!)
     vDSP_ctoz(observed, 2, &splitComplex, 1, vDSP_Length(halfN))
  }
}
// forwardInput[Real,Imag] now contain the de-interleaved data.
// splitComplex is out-of-scope and cannot be used, so the invalid pointers
// are discarded.

There's a few things that would make this easier.

First, there's a change coming to the Swift compiler that will diagnose this error for you.

Second, we can wrap the little dance I showed into some convenience functions:

/// De-interleave the real and imaginary parts of a complex buffer into two
/// new Float arrays.
func ctoz<T>(_ data: T) -> (real: [Float], imag: [Float])
where T: AccelerateBuffer, T.Element == DSPComplex {
  var imag = [Float]()
  let real = [Float](unsafeUninitializedCapacity: data.count) { r, n in
    imag = [Float](unsafeUninitializedCapacity: data.count) { i, n in
      ctoz(data, real: &r, imag: &i)
      n = data.count
    }
    n = data.count
  }
  return (real, imag)
}

/// De-interleave the real and imaginary parts of a complex buffer into two
/// caller-provided Float buffers.
///
/// - Precondition: data, real, and imag must all have the same length.
func ctoz<T, U, V>(_ data: T, real: inout U, imag: inout V)
where T: AccelerateBuffer, T.Element == DSPComplex,
      U: AccelerateMutableBuffer, U.Element == Float,
      V: AccelerateMutableBuffer, V.Element == Float
{
  precondition(data.count == real.count && data.count == imag.count)
  real.withUnsafeMutableBufferPointer { r in
    imag.withUnsafeMutableBufferPointer { i in
      var split = DSPSplitComplex(realp: r.baseAddress!, imagp: i.baseAddress!)
      data.withUnsafeBufferPointer { d in
        vDSP_ctoz(d.baseAddress!, 2, &split, 1, vDSP_Length(data.count))
      }
    }
  }
}

With these convenience functions defined, you can just do:

var forwardInputReal = [Float](repeating: 0, count: halfN)
var forwardInputImag = [Float](repeating: 0, count: halfN)
ctoz(observed, real: &forwardInputReal, imag: &forwardInputImag)

or even:

let (forwardInputReal, forwardInputImag) = ctoz(data)

I'll talk to the vDSP team and see if we can't get something like this added in the Framework for a future release so you don't need to write it yourself.