Source code in the common programming language

(defun fft (a &key (inverse nil) &aux (n (length a)))
  "Perform the FFT recursively on input vector A.
   Vector A must have length N of power of 2."
  (declare (type boolean inverse)
           (type (integer 1) n))
  (if (= n 1)
      a
      (let* ((n/2 (/ n 2))
             (2iπ/n (complex 0 (/ (* 2 pi) n (if inverse -1 1))))
             (⍵_n (exp 2iπ/n))
             (⍵ #c(1.0d0 0.0d0))
             (a0 (make-array n/2))
             (a1 (make-array n/2)))
        (declare (type (integer 1) n/2)
                 (type (complex double-float) ⍵ ⍵_n))
        (symbol-macrolet ((a0[j]  (svref a0 j))
                          (a1[j]  (svref a1 j))
                          (a[i]   (svref a i))
                          (a[i+1] (svref a (1+ i))))
          (loop :for i :below (1- n) :by 2
                :for j :from 0
                :do (setf a0[j] a[i]
                          a1[j] a[i+1])))
        (let ((â0 (fft a0 :inverse inverse))
              (â1 (fft a1 :inverse inverse))
              (â (make-array n)))
          (symbol-macrolet ((â[k]     (svref â k))
                            (â[k+n/2] (svref â (+ k n/2)))
                            (â0[k]    (svref â0 k))
                            (â1[k]    (svref â1 k)))
            (loop :for k :below n/2
                  :do (setf â[k]     (+ â0[k] (* ⍵ â1[k]))
                            â[k+n/2] (- â0[k] (* ⍵ â1[k])))
                  :when inverse
                    :do (setf â[k]     (/ â[k] 2)
                              â[k+n/2] (/ â[k+n/2] 2))
                  :do (setq ⍵ (* ⍵ ⍵_n))
                  :finally (return â)))))))


;;; This is adapted from the Python sample; it uses lists for simplicity.
;;; Production code would use complex arrays (for compiler optimization).
;;; This version exhibits LOOP features, closing with compositional golf.
(defun fft (x &aux (length (length x)))
  ;; base case: return the list as-is
  (if (<= length 1) x
    ;; collect alternating elements into separate lists...
    (loop for (a b) on x by #'cddr collect a into as collect b into bs finally
          ;; ... and take the FFT of both;
          (let* ((ffta (fft as)) (fftb (fft bs))
                 ;; incrementally phase shift each element of the 2nd list
                 (aux (loop for b in fftb and k from 0 by (/ pi length -1/2)
                            collect (* b (cis k)))))
            ;; finally, concatenate the sum and difference of the lists
            (return (mapcan #'mapcar '(+ -) `(,ffta ,ffta) `(,aux ,aux)))))))


;;; Demonstrates printing an FFT in both rectangular and polar form:
CL-USER> (mapc (lambda (c) (format t "~&~6F~6@Fi = ~6Fe^~6@Fipi"
                                   (realpart c) (imagpart c) (abs c) (/ (phase c) pi)))
               (fft '(1 1 1 1 0 0 0 0)))

   4.0  +0.0i =    4.0e^  +0.0ipi
   1.0-2.414i = 2.6131e^-0.375ipi
   0.0  +0.0i =    0.0e^  +0.0ipi
   1.0-0.414i = 1.0824e^-0.125ipi
   0.0  +0.0i =    0.0e^  +0.0ipi
   1.0+0.414i = 1.0824e^+0.125ipi
   0.0  +0.0i =    0.0e^  +0.0ipi
   1.0+2.414i = 2.6131e^+0.375ipi
;;; MAPC also returns the FFT data, which looks like this:
(#C(4.0 0.0) #C(1.0D0 -2.414213562373095D0) #C(0.0D0 0.0D0)
 #C(1.0D0 -0.4142135623730949D0) #C(0.0 0.0)
 #C(0.9999999999999999D0 0.4142135623730949D0) #C(0.0D0 0.0D0)
 #C(0.9999999999999997D0 2.414213562373095D0))