Fast Fourier transform: Difference between revisions
→{{header|Common Lisp}}
(Added graphic FreeBasic demo) |
|||
Line 661:
=={{header|Common Lisp}}==
As the longer standing solution below didn't work out for me
and I don't find it very nice, I want to give another one,
that's not just a plain translation. Of course it could be
optimized in several ways.
The function uses some non ASCII symbols for better readability
and condenses also the inverse part, by a keyword.
<lang lisp>
(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*pi/n (complex 0 (/ (* 2 pi) n (if inverse -1 1))))
(⍵_n (exp 2i*pi/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 â)))))))
</lang>
From here on the old solution.
<lang lisp>
;;; This is adapted from the Python sample; it uses lists for simplicity.
| |||