An FFT for Swift and the Xcode 6 Playground

Here's a preliminary way to view some FFT results using the new Apple Xcode6-beta Swift Playground (a recent Mac and an Apple Developer account required), perhaps useful for a bit of interactive "Matlab-like" DSP experimentation.  I can't yet figure out how to pass numeric vectors to Apple's really fast Accelerate/vDSP functions using the Swift Playground.  So, instead, I translated some ancient not-really-that-fast FFT code from Chipmunk Basic to Swift.

Enter the following Swift code in a new Playground. Wait a bit for the Playground to compute. Then click on the QuickLook bubble to the right of the m[i] variable near the end of this Swift script to see a plot of the FFT results. Enjoy.

import Foundation
// import Accelerate

var str = "Hello, playground"

var str2 = "My First Swift FFT"
println("\(str2)")

var len = 32 // radix-2 FFT length must be a power of 2

var myArray1 = [Double](count: len, repeatedValue: 0.0) // for Real Component
var myArray2 = [Double](count: len, repeatedValue: 0.0) // for Imaginary

var f0 = 3.0  // test input frequency

func myFill1 (inout a : [Double], inout b: [Double], n: Int, f0: Double) -> () {
 for i in 0 ..< n {
 // some test data
 let x = cos(2.0 * M_PI * Double(i) * f0 / Double(n))
 myArray1[i] = Double(x)  // Quicklook here to see a plot of the input waveform
 myArray2[i] = 0.0
 //
 }
 println("length = \(n)")
}

var sinTab = [Double]()

// Canonical in-place decimation-in-time radix-2 FFT
func myFFT (inout u : [Double], inout v : [Double], n: Int, dir : Int) -> () {
 
 var flag = dir // forward
 
 if sinTab.count != n {   // build twiddle factor lookup table
  while sinTab.count > n {
 sinTab.removeLast()
  }
  sinTab = [Double](count: n, repeatedValue: 0.0)
  for i in 0 ..< n {
 let x = sin(2.0 * M_PI * Double(i) / Double(n))
 sinTab[i] = x
  }
  println("sine table length = \(n)")
 }
 
 var m : Int = Int(log2(Double(n)))
 for k in 0 ..< n {
  // rem *** generate a bit reversed address vr k ***
  var ki = k
  var kr = 0
  for i in 1...m { // =1 to m
   kr = kr << 1  //  **  left shift result kr by 1 bit
   if ki % 2 == 1 { kr = kr + 1 }
   ki = ki >> 1   //  **  right shift temp ki by 1 bit
  }
  // rem *** swap data vr k to bit reversed address kr
  if (kr > k) {
   var tr = u[kr] ; u[kr] = u[k] ; u[k] = tr
   var ti = v[kr] ; v[kr] = v[k] ; v[k] = ti
  }
 }
 
 var istep = 2
 while ( istep <= n ) { //  rem  *** layers 2,4,8,16, ... ,n ***
  var is2 = istep / 2
  var astep = n / istep
  for km in 0 ..< is2 { // rem  *** outer row loop ***
   var a  = km * astep  // rem  twiddle angle index
   // var wr = cos(2.0 * M_PI * Double(km) / Double(istep))
   // var wi = sin(2.0 * M_PI * Double(km) / Double(istep))
   var wr =  sinTab[a+(n/4)] // rem  get sin from table lookup
   var wi =  sinTab[a]       // rem  pos for fft , neg for ifft
   if (flag == -1) { wi = -wi }
   for var ki = 0; ki <= (n - istep) ; ki += istep { //  rem  *** inner column loop ***
    var i = km + ki
    var j = (is2) + i
    var tr = wr * u[j] - wi * v[j]  // rem ** butterfly complex multiply **
    var ti = wr * v[j] + wi * u[j]  // rem ** using a temp variable **
    var qr = u[i]
    var qi = v[i]
    u[j] = qr - tr
    v[j] = qi - ti
    u[i] = qr + tr
    v[i] = qi + ti
   } // next ki
  } // next km
  istep = istep * 2
 }
 var a = 1.0 / Double(n)
 for i in 0 ..< n {
  u[i] = u[i] * a
  v[i] = v[i] * a
 }
}
// compute magnitude vector
func myMag (u : [Double], v : [Double], n: Int) -> [Double] {
 var m = [Double](count: n, repeatedValue: 0.0)
 for i in 0 ..< n {
 m[i] = sqrt(u[i]*u[i]+v[i]*v[i])   // Quicklook here to see a plot of the results
 
 }
 return(m)
}

myFill1(&myArray1, &myArray2, len, f0)
myFFT(  &myArray1, &myArray2, len, 1)

var mm = myMag(myArray1, myArray2, len)


// Feel free to treat the above code as if it were under an MIT style Open Source license.
// rhn 2014-June-08

// modified 2014-Jul-10 for Xcode6-beta3 newer Swift syntax

Software Defined Radio and IQ sampling

Why is quadrature or IQ sampling used for most Software Defined Radio interfaces and software algorithms? It has to do with the sampling rate, and how the sampling clock (the local oscillator or LO) relates to the signal frequency of interest. The Nyquist frequency is twice the highest frequency. But in practice, given finite length signals, and thus non-mathematically perfectly bandlimited signals, the sampling frequency for DSP has to be higher than twice the highest signal frequency. Thus doubling the number of samples by doubling the sample rate (2X LO) would still be too low. Quadrupling the sample rate (4X LO) would put you nicely above Nyquist rate, but using that much higher frequency would be more expensive in terms of circuit components, DSP data rates, megaflops required, and etc. So most IQ sampling is done with a local oscillator at (or relatively very near) the same frequency as the signal, which is obviously way too low a sampling frequency according to Nyquist. One sample per cycle of sine wave could be at the zero crossings, or at the tops, or at any point in between. You will learn almost nothing about a sinusoidal signal so sampled. But lets call this, by itself useless, set of samples the I of an IQ sample set. But how about increasing the number of samples, not by simply doubling the sample rate, but by taking an additional sample bit a little after the first one each cycle. Two samples per cycle a little bit apart would allow one to estimate the slope or derivative. If one sample was at a zero crossing the other one wouldn't be. So you would be far better off in figuring out the signal being sampled. Two points, plus knowledge that the signal is roughly periodic at the sample rate is usually enough to accurately estimate the unknowns of a canonical sinewave equation (amplitude and phase). But if you go too far apart with the second sample, to halfway between the first set of samples, you end up with the same problem as 2X sampling (one sample could be at a positive zero crossing, the other at a negative, telling you nothing). It's the same problem as 2X being too low a sample rate. But somewhere between two samples of the first set (the "I" set) there's a sweet spot. Not redundant, as with sampling at the same time, and not evenly spaced (which is equivalent to doubling the sample rate), there's an offset which gives you maximum information about the signal, with the cost being an accurate delay for the second sample instead of a much higher sample rate. Turns out that that delay is 90 degrees. That gives you a very useful "Q" set of samples, which together with the "I" set, tells you far more about a signal than either alone. Perhaps enough to demodulate AM, FM, SSB, QAM, etc., etc. (Also posted on http://electronics.stackexchange.com/ on 2014-Feb-12.)