Fun with Fourier series, Clojure lazy sequences and Incanter
Jun 6, 2010
The fourier series expansion of a square wave is given by:
sin(x) + (1/3)sin(3x) + (1/5)sin(5x) + ...
Let's try to generate (and plot) this series in Clojure.
We need to install a program called `Incanter' to generate plots; Incanter is an R-like statistical computing environment for Clojure; you can download it from here. The downloaded zip archive is self-contained - you need not install anything else to run Incanter. Just unzip the package, change into the `incanter' directory and start the Clojure REPL by running a script called `script/repl'. You can test the installation by executing this code fragment:
user=> (use '(incanter core charts stats)) user=> (def x (range 0 6.28 0.01)) user=> (def y (map sin x)) user=> (view (xy-plot x y))
We have plotted the first term of the fourier expansion - a simple sin curve!
What does the following code fragment do?
(defn seq-mul [n s] (map (fn [x] (* n x)) s))
It takes a number `n' and a sequence `s' and returns another sequence each element of which is multiplied by `n'.
Here is a function which will add two sequences:
(defn seq-add [a b] (map + a b))
We need a list of odd numbers:
(def odds (iterate #(+ % 2) 1))
Now we can generate the sequence: (x, 3x, 5x ... )
(def odd-multiples (map seq-mul odds (repeat x)))
`odd-multiples' is a sequence each element of which is again a sequence. We use `odd-multiples' to generate the sequence (sin(x), sin(3x), sin(5x) ...):
(def odd-harmonics (map (fn [x] (map sin x)) odd-multiples))
Now we need to `scale' each value in the sequence - we first generate the scale values 1, 1/3, 1/5, 1/7 ...
(def scale-values (map (fn [x] (/ 1.0 x)) odds))
We are now ready to generate the sequence (sin(x), (1/3)sin(3x), (1/5)sin(5x) ...):
(def fourier-series (map seq-mul scale-values odd-harmonics))
`fourier-series' represents an infinite series. We can take say the first 10 terms and add them up by:
(reduce seq-add (take 10 fourier))
Let's wrap this up in a convenient function:
(defn plot-fourier [n] (let [y (reduce seq-add (take n fourier))] (view (xy-plot x y))))
`plot-fourier' plots the sum of the first `n' terms of the series. Have fun by trying out different values for `n'!