3

How am I supposed to use the Accelerate framework to compute the FFT of a real signal in Swift on iOS?

Available example on the web

Apple’s Accelerate framework seems to provide functions to compute the FFT of a signal efficiently.

Unfortunately, most of the examples available on the Internet, like Swift-FFT-Example and TempiFFT, crash if tested extensively and call the Objective C API.

The Apple documentation answers many questions, but also leads to some others (Is this piece mandatory? Why do I need this call to convert?).

Threads on Stack Overflow

There are few threads addressing various aspects of the FFT with concrete examples. Notably FFT Using Accelerate In Swift, DFT result in Swift is different than that of MATLAB and FFT Calculating incorrectly - Swift. None of them address directly the question “What is the proper way to do it, starting from 0”?

It took me one day to figure out how to properly do it, so I hope this thread can give a clear explanation of how you are supposed to use Apple's FFT, show what are the pitfalls to avoid, and help developers save precious hours of their time.

Eric Postpischil
  • 195,579
  • 13
  • 168
  • 312
Jeremy Cochoy
  • 2,480
  • 2
  • 24
  • 39
  • What crashes? I've been using vDSP for FFTs in dozens of iOS apps for over a decade, and have never seen an FFT related crash. For real signals you don't need to convert as your strictly real data is already half of a split complex (the other half starts out as a zero'd buffer). – hotpaw2 Feb 08 '20 at 18:41
  • See post below. Using the same DSPSplitComplex as input and output doesn't work correctly if you do hundreads of FFT per seconds. Not correctly assuring that the pointer stay valid during the computation will also gives you sporadic crashes. It is not the vDSP.FFT object in itself that has any problem, it is more related to lack of documentation. – Jeremy Cochoy Feb 10 '20 at 13:43

1 Answers1

12

TL ; DR : If you need a working implementation to copy past here is a gist.

What is FFT?

The Fast Fourier transform is an algorithm that take a signal in the time domain -- a collection of measurements took at a regular, usual small, interval of time -- and turn it into a signal expressed into the phase domain (a collection of frequency). The ability to express the signal along time lost by the transformation (the transformation is invertible, which means no information is lost by computing the FFT and you can apply a IFFT to get the original signal back), but we get the ability to distinguish between frequencies that the signal contained. This is typically used to display the spectrograms of the music you are listening to on various hardware and youtube videos.

The FFT works with complexe numbers. If you don't know what they are, lets just pretend it is a combination of a radius and an angle. There is one complex number per point on a 2D plane. Real numbers (your usual floats) can be saw as a position on a line (negative on the left, positive on the right).

Nb: FFT(FFT(FFT(FFT(X))) = X (up to a constant depending on your FFT implementation).

How to compute the FFT of a real signal.

Usual you want to compute the FFT of a small window of an audio signal. For sake of the example, we will take a small 1024 samples window. You would also prefer to use a power of two, otherwise things gets a little bit more difficult.

var signal: [Float] // Array of length 1024

First, you need to initialize some constants for the computation.

// The length of the input
length = vDSP_Length(signal.count)
// The power of two of two times the length of the input.
// Do not forget this factor 2.
log2n = vDSP_Length(ceil(log2(Float(length * 2))))
// Create the instance of the FFT class which allow computing FFT of complex vector with length
// up to `length`.
fftSetup = vDSP.FFT(log2n: log2n, radix: .radix2, ofType: DSPSplitComplex.self)!

Following apple's documentation, we first need to create a complex array that will be our input. Dont get mislead by the tutorial. What you usual want is to copy your signal as the real part of the input, and keep the complex part null.

// Input / Output arrays

var forwardInputReal = [Float](signal) // Copy the signal here
var forwardInputImag = [Float](repeating: 0, count: Int(length))
var forwardOutputReal = [Float](repeating: 0, count: Int(length))
var forwardOutputImag = [Float](repeating: 0, count: Int(length))

Be careful, the FFT function do not allow to use the same splitComplex as input and output at the same time. If you experience crashs, this may be the cause. This is why we define both an input and an output.

Now, we have to be careful and "lock" the pointer to this four arrays, as showed in the documentation example. If you simply use &forwardInputReal as argument of your DSPSplitComplex, the pointer may become invalidated at the following line and you will likely experience sporadic crash of your app.

    forwardInputReal.withUnsafeMutableBufferPointer { forwardInputRealPtr in
      forwardInputImag.withUnsafeMutableBufferPointer { forwardInputImagPtr in
        forwardOutputReal.withUnsafeMutableBufferPointer { forwardOutputRealPtr in
          forwardOutputImag.withUnsafeMutableBufferPointer { forwardOutputImagPtr in
            // Input
            let forwardInput = DSPSplitComplex(realp: forwardInputRealPtr.baseAddress!, imagp: forwardInputImagPtr.baseAddress!)
            // Output
            var forwardOutput = DSPSplitComplex(realp: forwardOutputRealPtr.baseAddress!, imagp: forwardOutputImagPtr.baseAddress!)

            // FFT call goes here
          }
        }
      }
    }

Now, the finale line: the call to your fft:

fftSetup.forward(input: forwardInput, output: &forwardOutput)

The result of your FFT is now available in forwardOutputReal and forwardOutputImag.

If you want only the amplitude of each frequency, and you don't care about the real and imaginary part, you can declare alongside the input and output an additional array:

var magnitudes = [Float](repeating: 0, count: Int(length))

add right after your fft compute the amplitude of each "bin" with:

vDSP.absolute(forwardOutput, result: &magnitudes)
Anthony Dito
  • 3,610
  • 3
  • 29
  • 56
Jeremy Cochoy
  • 2,480
  • 2
  • 24
  • 39
  • Depends on whether you need the FFT to go "fast" or not. If you need the performance (max flops or min latency), then the way to go is to manually manage memory (allocate and deallocate unsafe mutable pointers) and use the C API. This will allow using the same raw memory for FFT input and output, which will result in less 1st and 2nd level data cache spills and misses on the CPU and memory hierarchy. Possibly less memory copies as well. – hotpaw2 Feb 08 '20 at 18:27
  • Could you elaborate on how to do it ? I am interested by this result, because using in this code the same DSPSplitComplex as input and output do not work on iPhone X and emulators. :) Feel free to edit or post a new answer. – Jeremy Cochoy Feb 10 '20 at 13:44
  • I am also interested by real fft (cf numpy rfft) implementation (i.e. reduce the number of computation using symmetries) if you have an implementation of this. – Jeremy Cochoy Feb 10 '20 at 13:46
  • You can use the same buffers for input and output of a vDSP FFT. If your program crashed while doing this, you had something else wrong in it. – Eric Postpischil Feb 12 '20 at 00:07
  • 2
    You should use vDSP’s real-to-complex FFT, not a complex-to-complex FFT with explicit zeros for the imaginary parts. Although the interface for this is not great (and requires using ztoc to prepare the data), it will be faster and use less memory. – Eric Postpischil Feb 12 '20 at 00:09
  • 1
    Performing an FFT does not lose the notion of time. Time information is present in the phases of the complex numbers. (A simple proof that no information is lost is that the FFT is fully reversible, aside from arithmetic rounding errors.) – Eric Postpischil Feb 12 '20 at 00:12
  • @EricPostpischil I didn't said any information was lost. FFT is just expressing a signal in a different basis (complex exponential). Cosine transform do it using cosines as base. Would you write an answer using the vDSP real only interface? I'd be happy to accept a better answer. :) Numpy and Scipy provide fft.rfft and fft.irfft which are optimised real->complex FFT but when I looked at the apple documentation I didn't found anything about real->complex FFT. – Jeremy Cochoy Feb 13 '20 at 09:03
  • @EricPostpischil If you don't have time to write an answer, just few references to the API functions to call for real->Complex fft and it's inverse is already a very good information. :) Edit: do you mean this https://developer.apple.com/documentation/accelerate/1449739-vdsp_dft_zrop_createsetup?language=objc ? Does the swift binding https://developer.apple.com/documentation/accelerate/vdsp/dft correspond to the same thing? :) – Jeremy Cochoy Feb 13 '20 at 09:08
  • I also wonder how costly it is to re-order the output of vDSP's iRFFT into a single array (as the output is interleaved both into the real and imaginary part according to documentation). – Jeremy Cochoy Feb 13 '20 at 09:19
  • @JeremyCochoy: I might consider showing how to use the real-to-complex FFT, but it may be quite a while. I have little free time these days, and I do not use Swift, so I will have to look up some things. It would be easier in C. – Eric Postpischil Feb 16 '20 at 23:52
  • @JeremyCochoy: That [DFT class](https://developer.apple.com/documentation/accelerate/vdsp/dft) is new to me; I am not sure it had been published yet when I left Apple. But I would expect it to bind to the family of routines that includes [`vDSP_DFT_zrop_CreateSetup`](https://developer.apple.com/documentation/accelerate/1449739-vdsp_dft_zrop_createsetup?language=objc). – Eric Postpischil Feb 16 '20 at 23:55
  • For additional documentation, find the file `vDSP.h` inside the Xcode SDKs (e.g., `/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/System/Library/Frameworks/Accelerate.framework/Frameworks/vecLib.framework/Headers/vDSP.h`) and open it in a text (or source code) editor. Read the “Documentation conventions” in the beginning. After that, you will find C declarations of the vDSP routines with documentation, including some pseudocode. – Eric Postpischil Feb 16 '20 at 23:58