# DFT (Discrete Fourier Transform) Algorithm in Swift

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

• This is not an FFT! This is a discrete Fourier transform and has none of the Cooley-Tukey “fast” algorithmic features. FFTs are well-understood and you should be able to port KissFFT or Nayuki’s FFT to Swift for arbitrary-length FFTs. Like I said in your StackOverflow post, the CZT will likely be much slower than your own port of KissFFT (probably even if CZT uses Accelerate FFT under the hood). – Ahmed Fasih Feb 3 '17 at 15:50
• The FFT is a divide-and-conquer recursive algorithm, which is how it achieves O(N log N) performance. – Ahmed Fasih Feb 3 '17 at 15:56

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

• Why the fastest? The article that you link to explains that the chirp transform can be used together with 3 FFTs, but states that this is "not as efficient as using methods specifically tailored to the factors of the sequence length" – Martin R Jan 31 '17 at 12:06
• @ Martin R Thanks for pointing that out. Edited it now. – UpSampler Jan 31 '17 at 13:58

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.

• The multiplications aren't a problem, although the division might be. The idea of using a map makes no sense because a single multiplication is very cheap compared to a map look up. Those trig identities make no sense either - the only logical one would be sin(x) = sqrt(1 - cos(x)^2) but it's unlikely to be any faster. – Alnitak Jan 31 '17 at 14:07
• It should be mentioned that expressions with square roots likesin(x) = sqrt(1 - cos(x)^2) is only true when 0<x<pi. If pi < x < 2pi then there is a minus sign in front of the square-root so this also needs to be taken into account if one is to use this. The expressions withsin(x/2)holds for 0 < x < 2pi, but not for 2pi < x < 4pi and -2pi < x < 0. – Winther Jan 31 '17 at 15:24
• @Winther good point, although in this case trig identities based around 2x would seem to be of very limited utility. I did try the more common identity and it was actually slower than just calling sin and cos directly. – Alnitak Jan 31 '17 at 15:52
• You picked the wrong trig identities to use. What you want are cos(q + a) = cos q cos a - sin q sin a and sin(q + a) = sin q cos a + cos q sin a. You only need to use the trig functions once in each loop, with a = 2 pi k/N. You already know sin q and cos q at the start of each loop, because q = 0. This needs only 4 multiplications and 2 adds to calculate each pair of trig functions, and is vectorizable. – alephzero Jan 31 '17 at 17:03
• @alephzero very nice! Would you expect that to be stable within the constraints of floating point arithmetic given that a is going to be a reasonably small number (4.189e-4) – Alnitak Jan 31 '17 at 21:48

# 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.

• FWIW, with the -O flag the swift compiler appears to be able to aggregate the two trig calls into a single call to ___sincos_stret, albeit without vectorisation. – Alnitak Jan 31 '17 at 14:15
• Good find! Very cool! – user1118321 Jan 31 '17 at 17:42
• If the first step is profiling, the zeroth step is finding the best algorithm. Profiling on its own will never turn O(n^2) into O(n log n). – alephzero Jan 31 '17 at 23:23

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


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.

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.

• Thank you very much friends! I have been able to accomplish the CZT based FFT -- works well. My next challenge is to accomplish inverse FFT (IFFT) in the same manner using CZT or rather ICZT. I have opened a new Q regarding it here: [link] (stackoverflow.com/questions/42017971/…) – Pat Feb 3 '17 at 6:23