Skip to content

Commit a3cac0b

Browse files
dasludaslu
committed
fft comparison wip
1 parent 22e9311 commit a3cac0b

1 file changed

Lines changed: 59 additions & 77 deletions

File tree

‎src/dsp/fft_comparison.clj‎

Lines changed: 59 additions & 77 deletions
Original file line numberDiff line numberDiff line change
@@ -41,40 +41,42 @@
4141
;; [Apache Commons Math](https://commons.apache.org/proper/commons-math/) is a comprehensive mathematical library for Java. It provides [`FastFourierTransformer`](https://commons.apache.org/proper/commons-math/javadocs/api-3.6.1/org/apache/commons/math3/transform/FastFourierTransformer.html) with various transform types and [normalization](https://en.wikipedia.org/wiki/Discrete_Fourier_transform#Normalization) options.
4242
;;
4343
;; **Key features:**
44-
;; - Part of Apache Commons (widely used, stable)
44+
;; - Part of Apache Commons (widely used, stable, maintained)
4545
;; - Supports multiple normalization conventions
46+
;; - 1D FFT, [Hadamard](https://en.wikipedia.org/wiki/Hadamard_transform), and sine/cosine transforms
4647
;; - Requires input length to be a [power of 2](https://en.wikipedia.org/wiki/Power_of_two)
47-
;; - Returns [`Complex[]`](https://commons.apache.org/proper/commons-math/javadocs/api-3.6.1/org/apache/commons/math3/complex/Complex.html) arrays (Java objects)
48+
;; - Returns [`Complex[]`](https://commons.apache.org/proper/commons-math/javadocs/api-3.6.1/org/apache/commons/math3/complex/Complex.html) arrays (allocation overhead)
4849

4950
;; ### jdsp
5051
;;
51-
;; [jdsp](https://jdsp.dev) is a Java [Digital Signal Processing](https://en.wikipedia.org/wiki/Digital_signal_processing) library providing [filters](https://en.wikipedia.org/wiki/Digital_filter), transforms, [peak detection](https://en.wikipedia.org/wiki/Peak_detection), and more. It uses Apache Commons Math internally for FFT computation.
52+
;; [jdsp](https://jdsp.dev) by [Sambit Paul](https://github.com/psambit9791) is a Java [Digital Signal Processing](https://en.wikipedia.org/wiki/Digital_signal_processing) library providing [filters](https://en.wikipedia.org/wiki/Digital_filter), transforms, [peak detection](https://en.wikipedia.org/wiki/Peak_detection), and more. It uses Apache Commons Math internally for FFT computation.
5253
;;
5354
;; **Key features:**
5455
;; - Convenient wrapper around Apache Commons Math
5556
;; - Simple API: [`new FastFourier(signal).transform()`](https://javadoc.io/doc/com.github.psambit9791/jdsp/latest/com/github/psambit9791/jdsp/transform/FastFourier.html)
56-
;; - Includes additional DSP utilities ([filters](https://en.wikipedia.org/wiki/Digital_filter), [wavelets](https://en.wikipedia.org/wiki/Wavelet), [convolution](https://en.wikipedia.org/wiki/Convolution))
57-
;; - Good for projects needing broader DSP functionality
57+
;; - Includes comprehensive DSP utilities ([filters](https://en.wikipedia.org/wiki/Digital_filter), [wavelets](https://en.wikipedia.org/wiki/Wavelet), [convolution](https://en.wikipedia.org/wiki/Convolution), [STFT](https://en.wikipedia.org/wiki/Short-time_Fourier_transform))
58+
;; - Good for projects needing broader DSP functionality beyond FFT
5859

5960
;; ### JTransforms
6061
;;
61-
;; [JTransforms](https://github.com/wendykierp/JTransforms) by Piotr Wendykier is the first open-source, multithreaded FFT library written in pure Java. It's optimized for performance with [parallel processing](https://en.wikipedia.org/wiki/Parallel_computing) support.
62+
;; [JTransforms](https://github.com/wendykierp/JTransforms) by [Piotr Wendykier](https://github.com/wendykierp) is the first open-source, multithreaded FFT library written in pure Java. It's optimized for performance with [parallel processing](https://en.wikipedia.org/wiki/Parallel_computing) support.
6263
;;
6364
;; **Key features:**
6465
;; - **Parallelized** [split-radix](https://en.wikipedia.org/wiki/Split-radix_FFT_algorithm) and [mixed-radix](https://en.wikipedia.org/wiki/Mixed-radix) algorithms
66+
;; - Supports **1D, 2D, and 3D** transforms (FFT, [DCT](https://en.wikipedia.org/wiki/Discrete_cosine_transform), [DST](https://en.wikipedia.org/wiki/Discrete_sine_transform), [DHT](https://en.wikipedia.org/wiki/Discrete_Hartley_transform))
6567
;; - [In-place mutations](https://en.wikipedia.org/wiki/In-place_algorithm) (efficient but not functional)
66-
;; - Supports FFT, [DCT](https://en.wikipedia.org/wiki/Discrete_cosine_transform), [DST](https://en.wikipedia.org/wiki/Discrete_sine_transform), [DHT](https://en.wikipedia.org/wiki/Discrete_Hartley_transform) transforms
68+
;; - Mixed-radix support: works with arbitrary sizes (not just power-of-2)
6769
;; - Used internally by fastmath and dtype-next
6870

6971
;; ### fastmath
7072
;;
71-
;; [fastmath](https://github.com/generateme/fastmath) (version 3.x) is a Clojure library for fast primitive-based mathematics. Its [`fastmath.transform`](https://generateme.github.io/fastmath/fastmath.transform.html) namespace wraps JTransforms with an idiomatic Clojure API.
73+
;; [fastmath](https://github.com/generateme/fastmath) (version 3.x) by [Tomasz Sulej](https://github.com/genmeblog) is a Clojure library for fast primitive-based mathematics. Its [`fastmath.transform`](https://generateme.github.io/fastmath/fastmath.transform.html) namespace wraps JTransforms with an idiomatic Clojure API.
7274
;;
7375
;; **Key features:**
7476
;; - Immutable, [functional](https://en.wikipedia.org/wiki/Functional_programming) API (no in-place mutations)
75-
;; - Leverages JTransforms' parallelized performance
77+
;; - Leverages JTransforms' parallelized performance and mixed-radix algorithms
7678
;; - [Protocol-based](https://clojure.org/reference/protocols) design: `transformer` → `forward-1d`/`reverse-1d`
77-
;; - Supports 1D and 2D transforms, multiple transform types
79+
;; - Supports **1D and 2D** transforms, multiple transform types (FFT, DCT, DST, DHT, wavelets)
7880

7981
;; ## Test Signal: Two-Tone Sine Wave
8082

@@ -105,7 +107,6 @@
105107

106108
;; Let's visualize the signal:
107109

108-
^:kindly/hide-code
109110
(-> (tc/dataset {:time (take 128 time-points)
110111
:amplitude (take 128 signal)})
111112
(plotly/base {:=x :time
@@ -118,6 +119,34 @@
118119
(plotly/layer-line {:=mark-color "steelblue"})
119120
plotly/plot)
120121

122+
;; ## Visualization Helper
123+
124+
;; To visualize FFT results throughout this post, we'll use a helper function:
125+
126+
(defn plot-fft-spectrum
127+
"Plot FFT magnitude spectrum.
128+
129+
Parameters:
130+
- magnitudes: sequence of magnitude values
131+
- title: plot title
132+
- color: line color (e.g., 'darkgreen', 'darkorange')
133+
134+
Returns Plotly visualization."
135+
[magnitudes title color]
136+
(let [n-bins (count magnitudes)
137+
freq-bins (dfn/* (range n-bins) (/ sample-rate n-samples))]
138+
(-> (tc/dataset {:frequency freq-bins
139+
:magnitude magnitudes})
140+
(plotly/base {:=x :frequency
141+
:=y :magnitude
142+
:=x-title "Frequency (Hz)"
143+
:=y-title "Magnitude"
144+
:=title title
145+
:=width 700
146+
:=height 300})
147+
(plotly/layer-line {:=mark-color color})
148+
plotly/plot)))
149+
121150
;; ## FFT Implementation #1: Apache Commons Math
122151

123152
;; Apache Commons Math requires creating a [`FastFourierTransformer`](https://commons.apache.org/proper/commons-math/javadocs/api-3.6.1/org/apache/commons/math3/transform/FastFourierTransformer.html) with a [normalization](https://en.wikipedia.org/wiki/Discrete_Fourier_transform#Normalization) type, then calling `transform` with the signal.
@@ -146,20 +175,10 @@
146175

147176
;; Visualize the first half (positive frequencies):
148177

149-
^:kindly/hide-code
150-
(let [freq-bins (dfn/* (range (quot (count commons-magnitudes) 2))
151-
(/ sample-rate n-samples))]
152-
(-> (tc/dataset {:frequency freq-bins
153-
:magnitude (take (quot (count commons-magnitudes) 2) commons-magnitudes)})
154-
(plotly/base {:=x :frequency
155-
:=y :magnitude
156-
:=x-title "Frequency (Hz)"
157-
:=y-title "Magnitude"
158-
:=title "FFT Spectrum (Apache Commons Math)"
159-
:=width 700
160-
:=height 300})
161-
(plotly/layer-line {:=mark-color "darkgreen"})
162-
plotly/plot))
178+
(plot-fft-spectrum
179+
(take (quot (count commons-magnitudes) 2) commons-magnitudes)
180+
"FFT Spectrum (Apache Commons Math)"
181+
"darkgreen")
163182

164183
;; ## FFT Implementation #2: jdsp
165184

@@ -176,22 +195,11 @@
176195
(def jdsp-result
177196
(time (fft-jdsp signal)))
178197

179-
;; Visualize:
180-
181-
^:kindly/hide-code
182-
(let [freq-bins (dfn/* (range (count jdsp-result))
183-
(/ sample-rate n-samples))]
184-
(-> (tc/dataset {:frequency freq-bins
185-
:magnitude jdsp-result})
186-
(plotly/base {:=x :frequency
187-
:=y :magnitude
188-
:=x-title "Frequency (Hz)"
189-
:=y-title "Magnitude"
190-
:=title "FFT Spectrum (jdsp)"
191-
:=width 700
192-
:=height 300})
193-
(plotly/layer-line {:=mark-color "darkorange"})
194-
plotly/plot))
198+
; Visualize
199+
(plot-fft-spectrum
200+
jdsp-result
201+
"FFT Spectrum (jdsp)"
202+
"darkorange")
195203

196204
;; ## FFT Implementation #3: JTransforms (Direct)
197205

@@ -224,22 +232,11 @@
224232
(def jtransforms-magnitudes
225233
(jtransforms-magnitude jtransforms-result))
226234

227-
;; Visualize:
228-
229-
^:kindly/hide-code
230-
(let [freq-bins (dfn/* (range (count jtransforms-magnitudes))
231-
(/ sample-rate n-samples))]
232-
(-> (tc/dataset {:frequency freq-bins
233-
:magnitude jtransforms-magnitudes})
234-
(plotly/base {:=x :frequency
235-
:=y :magnitude
236-
:=x-title "Frequency (Hz)"
237-
:=y-title "Magnitude"
238-
:=title "FFT Spectrum (JTransforms)"
239-
:=width 700
240-
:=height 300})
241-
(plotly/layer-line {:=mark-color "purple"})
242-
plotly/plot))
235+
; Visualize
236+
(plot-fft-spectrum
237+
jtransforms-magnitudes
238+
"FFT Spectrum (JTransforms)"
239+
"purple")
243240

244241
;; ## FFT Implementation #4: fastmath
245242

@@ -269,22 +266,11 @@
269266
(def fastmath-magnitudes
270267
(fastmath-magnitude fastmath-result))
271268

272-
;; Visualize:
273-
274-
^:kindly/hide-code
275-
(let [freq-bins (dfn/* (range (count fastmath-magnitudes))
276-
(/ sample-rate n-samples))]
277-
(-> (tc/dataset {:frequency freq-bins
278-
:magnitude fastmath-magnitudes})
279-
(plotly/base {:=x :frequency
280-
:=y :magnitude
281-
:=x-title "Frequency (Hz)"
282-
:=y-title "Magnitude"
283-
:=title "FFT Spectrum (fastmath)"
284-
:=width 700
285-
:=height 300})
286-
(plotly/layer-line {:=mark-color "crimson"})
287-
plotly/plot))
269+
; Visualize
270+
(plot-fft-spectrum
271+
fastmath-magnitudes
272+
"FFT Spectrum (fastmath)"
273+
"crimson")
288274

289275
;; ## Performance Comparison
290276

@@ -310,7 +296,6 @@
310296
:jtransforms (benchmark-fft fft-jtransforms signal n)
311297
:fastmath (benchmark-fft fft-fastmath signal n)}))
312298

313-
^:kindly/hide-code
314299
(kind/table
315300
[{:library "Apache Commons Math"
316301
:time-per-fft (format "%.3f ms" (get-in bench-small-128 [:apache-commons :per-iter-ms]))}
@@ -332,7 +317,6 @@
332317
:jtransforms (benchmark-fft fft-jtransforms signal-large n)
333318
:fastmath (benchmark-fft fft-fastmath signal-large n)}))
334319

335-
^:kindly/hide-code
336320
(kind/table
337321
[{:library "Apache Commons Math"
338322
:time-per-fft (format "%.3f ms" (get-in bench-large-131k [:apache-commons :per-iter-ms]))}
@@ -393,7 +377,6 @@
393377
(ConcurrencyUtils/setNumberOfThreads (.availableProcessors (Runtime/getRuntime)))
394378

395379
; Visualize results
396-
^:kindly/hide-code
397380
(-> (tc/dataset thread-performance)
398381
(plotly/base {:=x :threads
399382
:=y :per-iter-ms
@@ -441,7 +424,6 @@
441424
:available-processors (.availableProcessors (Runtime/getRuntime))
442425
:clojure-version (clojure-version)})
443426

444-
^:kindly/hide-code
445427
(kind/table
446428
[(assoc machine-info :property "Machine Configuration")])
447429

0 commit comments

Comments
 (0)