DFT (Discrete Fourier Transform) Algorithm in Swift

Question

I am looking to replicate in Swift what the FFT function does in Matlab. Essentially, it takes an arbitrary length signal (not necessarily a multiple of \$2^n\$) and gives real and complex DFT coefficients.

Since the FFT described in Accelerate can only handle sample sizes that are multiples of \$2^n\$, I wrote a brute force algorithm in Swift that produces exactly the same results as the Matlab FFT function for arbitrary sample size.

The problem: When my sample size > 15,000 sample (say), this algorithm takes about 20 s to complete. Could this be sped up?

import Foundation

public func fft(x: [Double]) -> ([Double],[Double]) {

    let N = x.count
    var Xre: [Double] = Array(repeating:0, count:N)
    var Xim: [Double] = Array(repeating:0, count:N)

    for k in 0..<N
    {
        Xre[k] = 0
        Xim[k] = 0
        for n in 0..<N {
            let q = (Double(n)*Double(k)*2.0*M_PI)/Double(N)
            Xre[k] += x[n]*cos(q) // Real part of X[k]
            Xim[k] -= x[n]*sin(q) // Imag part of X[k]
        }
    }
    return (Xre, Xim)
}

// Call FFT
let x: [Double] = [1, 2, 3, 4, 5, 6] // works rapidly
// let x = Array(stride(from: 0, through: 15000, by: (1.0))) // Will   choke it
let (fr, fi) = fft (x: x)
print("Real:", fr)
print(" ")
print("Imag:", fi)
// Call FFT

Answer 1

The fastest algorithm would be a native implementation of the "chirp z-transform" as described here with example code in Python.

Answer 2

I think your code is doing DFT (Discrete Fourier Transform) and not FFT. You are doing \$O(n^{2})\$ operations in the 2 for loops. The FFT is supposed to be \$n*\log(n)\$.

First thing to do is remove repeated multiplications. The terms Double(n)*2*M_PI/Double(N) can be calculated (as initial step) for every \$n\$ in \$0:(N-1)\$. Make a map for each \$n\$ to this value and use it to calculate \$q\$.

Then, see if using properties of \$\sin(2x)\$, \$\sin(\frac{x}{2})\$, \$\cos(2x)\$, \$\cos(\frac{x}{2})\$, etc. are faster than actually calling those functions.

For example:

$$\begin{align} \sin(2x) &= 2\sin(x)\cos(x) \\ \cos(2x) &= \cos^{2}(x) - \sin^{2}(x) \end{align} $$

or

$$\begin{align} \sin\left(\frac{x}{2}\right) &= \sqrt{\frac{1-\cos(x)}{2}} \\ \cos\left(\frac{x}{2}\right) &= \sqrt{\frac{1+\cos(x)}{2}} \end{align} $$

If these are faster, you will reduce the complexity significantly.

Answer 3

Profile

The first step of improving the performance is to profile it. I recommend running this using Xcode's profile option and see where the time is spent. I suspect (but don't know for sure) that it will be in the calls to sin() and cos().

Avoid Casts

One thing that can slow down calculations is lots of casts between types. You're using k, n, and N as integers in most of the code, but need to cast them to Double to calculate q. You could keep a parallel dk, dn, and dN that are floating point copies of k, n, and N to avoid the casts. You'll need to manually increment them in the loops, though.

Do more at once

If you look in <Accelerate/vfp.h>, you'll find vsinf() and vcosf() and more importantly, vsincosf() which calculate sine, cosine, and both at once for a whole vector of Floats. The precision is less than Double, so I don't know if it meets your precision needs, but I'd look into it. This should allow you to work on 16 elements at a time instead of only 1.

Answer 4

Did you turn on the Swift compiler's optimiser? On my early 2016 Macbook Pro (3.1 GHz dual core i7) your code runs in 15.5 seconds with it off, 4.1 seconds with on, and 3.3 seconds if I move any loop invariants outside of the relevant loops:

public func fft2(x: [Double]) -> ([Double],[Double]) {

    let N = x.count
    var Xre: [Double] = Array(repeating:0, count:N)
    var Xim: [Double] = Array(repeating:0, count:N)

    let f = 2.0 * M_PI / Double(N)    // <---   here

    for k in 0..<N
    {
        Xre[k] = 0
        Xim[k] = 0
        let kf = Double(k) * f        // <---   and here
        for n in 0..<N {
            let q = Double(n) * kf
            let (cq, sq) = (cos(q), sin(q))
            Xre[k] += x[n] * cq // Real part of X[k]
            Xim[k] -= x[n] * sq // Imag part of X[k]
        }
    }
    return (Xre, Xim)
}

I tried using the identity sq = sqrt(1 - cq * cq) to remove one trig call (by replacing it with a sqrt call) but it was slightly slower.

Using this loop instead, based on an idea from @alpehzero using the two following trigonometric identities:

cos(q + a) = cos q cos a - sin q sin a
sin(q + a) = sin q cos a + cos q sin a

gets the run time down to 0.65 seconds, albeit with a very small possible loss of precision:

for k in 0 ..< N {

    let kf = Double(k) * f
    let (cosa, sina) = (cos(kf), sin(kf))
    var (cosq, sinq) = (1.0, 0.0)

    for n in 0 ..< N {
        Xre[k] += x[n] * cosq
        Xim[k] -= x[n] * sinq
        (cosq, sinq) = (cosq * cosa - sinq * sina, sinq * cosa + cosq * sina)
    }
}

Answer 5

In many practical cases you can do a lot better than your basic approach, if $N$ has factors which are small integers. Since 15,000 = 2.2.2.3.5.5.5.5 your specific benchmark would run almost as fast as when N is a power of 2.

This approach doesn't work when N is a prime number, but there is an alternative idea based on factorizing N-1. See C. M. Rader, "Discrete Fourier transforms when the number of data samples is prime," Proc. IEEE 56, 1107–1108 (1968).

Good implementations of DFT algorithms in C (using the algorithms mentioned above) are at http://www.fftw.org/. If you can't link C code with Swift, it should be possible to translate what you need.

Answer 6

This is not a language specific problem.

The Discrete FT running time is theta(N^2), and as you note, it is not fast. However the Fast FT is named as such because it has theta(N lg N) performance. Just adapt a C FFT to swift.