|
| 1 | +(ns fft.core |
| 2 | + (:require [complex.core :as c])) |
| 3 | +;; complex is a jar for complex numbers |
| 4 | +;; https://github.com/alanforr/complex |
| 5 | +;; add [complex "0.1.11"] to :dependencies in your project.clj |
| 6 | +;; and run lein repl or lein deps in the terminal |
| 7 | +(defn matrix-mult |
| 8 | + "take a matrix m and a vector v which length is number of columns |
| 9 | + ,return a vector of applying dot-product between v and each row of |
| 10 | + m. the returned vector's length is the number of rows of m" |
| 11 | + [m v] |
| 12 | + (mapv (comp (partial apply c/+) |
| 13 | + (partial map c/* v)) |
| 14 | + m)) |
| 15 | +(defn dft |
| 16 | + "take a vector of real numbers and return a vector of frequency |
| 17 | + space" |
| 18 | + [vx] |
| 19 | + (let [len (count vx)] |
| 20 | + (matrix-mult |
| 21 | + (partition len |
| 22 | + (for [n (range len) |
| 23 | + k (range len)] |
| 24 | + ;; expresion below is |
| 25 | + ;; e^(n*k*2*pi*(1/len)*(-i)) |
| 26 | + (c/exp (c/* n k |
| 27 | + 2 Math/PI |
| 28 | + (/ len) |
| 29 | + (c/complex 0 -1))))) |
| 30 | + vx))) |
| 31 | +(defn fft [vx] |
| 32 | + (let [len (count vx)] |
| 33 | + (if (= len 1) |
| 34 | + vx |
| 35 | + ;;else |
| 36 | + (let [;; take values of vx in the even indices |
| 37 | + even-indices (keep-indexed #(if (even? %1) %2) vx) |
| 38 | + ;; take values in the odd indices |
| 39 | + odd-indices (keep-indexed #(if (odd? %1) %2) vx) |
| 40 | + ;; recursion |
| 41 | + even-fft (fft even-indices) |
| 42 | + odd-fft (fft odd-indices) |
| 43 | + ;; make a sequence of e^(-2pi*i*k/N) where N is the length |
| 44 | + ;; vx and k range from 0 to N/2 |
| 45 | + omegas-half (map |
| 46 | + (comp c/exp |
| 47 | + (partial c/* |
| 48 | + (/ len) |
| 49 | + 2 Math/PI |
| 50 | + (c/complex 0 -1))) |
| 51 | + (range 0 (quot len 2))) |
| 52 | + ;; take the negative of the first sequence because |
| 53 | + ;; e^(-2pi*i*(k+N/2)/N=-e^(-2pi*i*k/N) where k ranges from |
| 54 | + ;; 0 to N/2 |
| 55 | + omegas-2half (map c/- omegas-half) |
| 56 | + mult-add (partial map #(c/+ %3 (c/* %1 %2)))] |
| 57 | + (concat (mult-add omegas-half odd-fft even-fft) |
| 58 | + (mult-add omegas-2half odd-fft even-fft)))))) |
| 59 | +(defn -main [& args] |
| 60 | + (let [vx [0 1 2 3] |
| 61 | + len (count vx) |
| 62 | + ;; calculate the next power of 2 after len |
| 63 | + ;; the reason behind this is to fill them with zeros for fft |
| 64 | + next-len (->> |
| 65 | + [len 2] |
| 66 | + (map #(Math/log %)) |
| 67 | + (apply /) |
| 68 | + Math/ceil |
| 69 | + (Math/pow 2) |
| 70 | + int) |
| 71 | + ;; add zeros at the end of vx |
| 72 | + complete-vx (into vx (repeat (- next-len len) 0)) |
| 73 | + fft-cvx (fft complete-vx) |
| 74 | + dft-cvx (dft complete-vx) |
| 75 | + diffv (mapv c/- fft-cvx dft-cvx)] |
| 76 | + (println "vx:" vx) |
| 77 | + (println "complete-vx:" complete-vx) |
| 78 | + (println "result from fft:" (map c/stringify fft-cvx)) |
| 79 | + (println "result from dft:" (map c/stringify dft-cvx)) |
| 80 | + (println "difference: " (map c/stringify diffv)))) |
0 commit comments