lisp-random / fft.lisp

;;;; Copyright (c) 2012 Robert Smith
;;;; Simple FFT
;;;; Requires QTILITY

;;; Constants

(defconstant ii #C(0 1))

(defconstant 2*pi*i (* 2 pi ii))


;;; Utilities

(defmacro define-pointwise-operation (name function)
  (let ((v1 (gensym "V1-"))
        (v2 (gensym "V2-")))
    `(defun ,name (,v1 ,v2)
       (map 'vector ,function ,v1 ,v2))))

(define-pointwise-operation .+ #'+)
(define-pointwise-operation .- #'-)
(define-pointwise-operation .* #'*)
(define-pointwise-operation ./ #'/)
(define-pointwise-operation .^ #'expt)

(defun slice (v indexes)
  (loop
    :for i :in indexes
    :collect (aref v i) :into s
    :finally (return (coerce s 'vector))))

(defun catvec (v1 v2)
  (concatenate 'vector v1 v2))

(defun primitive-nth-root (n &optional (k 1))
  (if (= n k)
      1
      (exp (/ (* 2*pi*i k) n))))

;;; FFT
(defun fft (x)
  ;; This has a normalization factor of 1
  (let ((N (length x)))
    (if (= 1 N)
        x
        (let* ((top   (fft (slice x (qtl:range 0 N 2))))
               (bot   (fft (slice x (qtl:range 1 N 2))))
               (roots (qtl:tabulate (lambda (k) (primitive-nth-root N (- k)))
                                    (floor N 2)))
               (z     (.* roots bot)))
          (catvec (.+ top z) (.- top z))))))


;;; DFT
(defun dft (x &key (normalization-factor (/ (sqrt (length x))))
                   inversep)
  (let ((N (length x)))
    (labels ((coeff (k)
               (loop
                 :with inverter := (if inversep 1 -1)
                 :for j :below N
                 :sum (* (aref x j)
                         (primitive-nth-root N (* j k inverter))))))
      (loop :for k :below N
            :collect (* (coeff k) normalization-factor) :into dft-array
            :finally (return (coerce dft-array 'vector))))))
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.