Signal Transforms: A Comprehensive Guide to FFT, DCT, and Wavelets
Introduction: Why Learn Signal Transforms?
About This Tutorial
This is the first in a planned series on digital signal processing in Clojure. We focus here on core transforms—the fundamental tools for moving between time and frequency domains. Future posts will explore specialized topics we deliberately omit here.
What’s covered: FFT (Fast Fourier Transform), DCT (Discrete Cosine Transform), wavelets, and their practical applications in filtering, compression, and denoising. We’ll also cover critical FFT artifacts like spectral leakage, aliasing, and windowing functions.
What’s intentionally left out for future posts:
- Filter design - Designing custom filters (lowpass, highpass, bandpass)
- Real-world data - Loading and processing audio/image files
- Advanced methods - Better spectral estimation, adaptive filtering
- Multi-rate processing - Changing sample rates, efficient resampling
The goal is to build a solid foundation in transforms before tackling these specialized topics.
Why Learn Signal Transforms?
Imagine you’re listening to a musical note. You can hear it’s an A (440 Hz), but how would you write code to discover that frequency? Or suppose you have 10 seconds of audio—how do you automatically remove a 60 Hz electrical hum without affecting the music?
These questions lead us to signal transforms—mathematical tools that reveal hidden structure in data. Classic JPEG images use the Discrete Cosine Transform for compression. MP3 audio uses a variant called the Modified DCT. Speech recognition, radar systems, medical imaging—all depend on transforming signals between different representations.
At their core, transforms answer a simple question: can we express this signal as a weighted sum of simpler basis functions? The Fourier transform uses sines and cosines. Wavelets use localized, scaled functions. Each transform gives us a different lens for understanding and manipulating data.
What Tools We’ll Use
We’ll use JTransforms—a fast, battle-tested, pure Java library—for computing the actual transforms (FFT, DCT, wavelets, etc.). It’s wrapped in fastmath, which provides a clean, idiomatic Clojure API for signal processing.
For array operations on signals, we’ll use dtype-next, which provides efficient typed arrays and functional operations like dfn/+, dfn/*, and dfn/sq. This gives us the speed of Java libraries with the expressiveness of Clojure’s functional style.
What You’ll Learn
This tutorial combines theory with practice, building from foundations to real applications. We’ll start with signal generation, creating test signals with different characteristics. Then we’ll explore the core transforms—DFT (via FFT), DCT, and wavelets—along with specialized variants. You’ll see practical applications like filtering, compression, and denoising, with efficient dtype-next array operations throughout. Each section demonstrates transforms through runnable code with visual feedback, and we’ll use test-driven validation to ensure correctness, making abstract mathematics concrete and verifiable.
Prerequisites
You should be comfortable with basic Clojure (functions, maps, sequences) and high school mathematics (sine, cosine, basic algebra). Understanding that higher pitch means higher frequency is helpful. No digital signal processing background required—we’ll build intuition through visualization and concrete examples.
Tutorial Structure
We’ll progress from simple signal generation to advanced multi-dimensional transforms. We start with signal generation, creating test signals with known properties. Then we explore the DFT/FFT for discovering frequencies, the DCT for compression (used in JPEG and MP3), and wavelets for time-frequency localization and denoising. We’ll also cover specialized transforms like DST, DHT, and 2D variants.
The later sections focus on practical applications—filtering, compression, and analysis tools—followed by a testing framework for validation and best practices for choosing the right transform. Throughout, we’ll introduce dtype-next operations as needed, learning efficient array manipulation in the context of solving real problems. Each section includes visualizations and round-trip tests to verify correctness.
Let’s begin by learning to generate the signals we’ll be transforming!
Signal Generation and Visualization
Before we can analyze frequencies, we need signals to analyze. This section introduces signal generation and visualization—the foundation for everything that follows.
We’ll also meet dtype-next, Clojure’s efficient array programming library. Rather than learning it abstractly, we’ll discover dtype-next operations as we need them for generating and manipulating signals.
Introduction to fastmath.signal and dtype-next
(require '[fastmath.signal :as sig])(require '[tech.v3.datatype :as dtype])(require '[tech.v3.datatype.functional :as dfn])(require '[tech.v3.datatype.argops :as argops])The fastmath.signal library provides signal generation through oscillators. An oscillator is simply a function mapping time to amplitude: time → amplitude. This functional approach makes it easy to compose complex signals from simple components.
Create a 440 Hz sine wave oscillator
(def osc-440 (sig/oscillator :sin 440.0 1.0 0.0))Parameters: type, frequency, amplitude, phase
Oscillators are functions - evaluate at specific times (t=0 and t=0.1s):
[(osc-440 0.0) (osc-440 0.1)][0.0 -3.9198601261660535E-14]Generate discrete samples: 44.1kHz sample rate, 10ms duration
(sig/oscillator->signal osc-440 44100.0 0.01)[0.0, 0.06264832417874365, 0.12505052369452813, 0.18696144082725336,
0.24813784794373792, 0.30833940305910035, 0.36732959406137883,
0.4248766678898383, 0.4807545410165316, 0.5347436876541294,
0.5866320022005456, 0.6362156325320929, 0.6832997808714387,
0.7276994690840094, 0.7692402653962486, 0.8077589696806925,
0.8431042546155975, 0.8751372602002244, 0.9037321392901133,
0.9287765520091285, 0.9501721070958867, 0.9678347484506665,
0.9816950853641806, 0.9916986651308531, 0.9978061869755916,
0.9999936564536084, 0.9982524797167027, 0.9925894972756649,
0.9830269571261583, 0.969602427343661, 0.9523686484908521,
0.9313933264172862, 0.9067588662653749, 0.8785620487276838,
0.8469136498274177, 0.8119380057158566, 0.7737725241965071,
0.7325671448950286, 0.6884837501957588, 0.6416955292590556,
0.5923862976180476, 0.5407497750278533, 0.486988824404367,
0.43131465484258247, 0.37394599184551275, 0.3151082180236203,
0.2550324876406644, 0.19395481848461169, 0.13211516463136738,
0.06975647374412507, 0.007123732611891674, -0.05553699532303968,
-0.11797953666743062, -0.17995857521140243, -0.24123061569515214,
-0.30155494042174963, -0.36069455495780295, -0.4184171192066551,
-0.47449586019620704, -0.5287104629953577, -0.5808479362589035,
-0.6307034490004968, -0.6780811353062296, -0.722794863827391,
-0.7646689690293306, -0.8035389413235673, -0.839252073371868,
-0.8716680600231562, -0.9006595495263089, -0.9261126438533076,
-0.9479273461671316, -0.9660179536764468, -0.980313394333687,
-0.990757506053757, -0.997309257356394, -0.999942908565361,
-0.9986481129311697, -0.9934299572800557, -0.9843089420295121,
-0.9713209006488901, -0.954516858881484, -0.933962834281163,
-0.9097395768511173, -0.8819422518036459, -0.8506800656873393,
-0.8160758373504611, -0.7782655154260821, -0.7373976442346147,
-0.6936327802020241, -0.6471428610864333, -0.5981105304912157,
-0.5467284203183442, -0.493198393980991, -0.43773075334857275,
-0.38054341253986085, -0.3218610418100651, -0.2619141848953075,
-0.20093835328207732, -0.13917310096006505, -0.076861083293317,
-0.014247103707103909, 0.04842284806422461, 0.1109025623900097,
0.1729465770199549, 0.2343111414245828, 0.2947551744109035,
0.3540412112511461, 0.41193633660356815, 0.46821309956024515,
0.5226504072268894, 0.5750343933241133, 0.6251592583977146,
0.672828078337052, 0.7178535780251143, 0.7600588670808603,
0.7992781348033673, 0.8353573015875321, 0.8681546242521786,
0.8975412529023928, 0.9234017371383941, 0.9456344796221872,
0.9641521352200855, 0.9788819541530247, 0.9897660678065152,
0.9967617160773985, 0.9998414153642309, 0.9989930665413147,
0.994220002492184, 0.9855409750158024, 0.9729900811569101,
0.9566166292499495, 0.936484945202834, 0.9126741197816123,
0.8852776978888555, 0.8544033110564996, 0.8201722545969549,
0.7827190110737126, 0.7421907219635712, 0.6987466095861785,
0.6525573515718863, 0.6038044103254782, 0.5526793201200532,
0.49938293462186145, 0.44412463780229566, 0.3871215213371295,
0.3285975317247166, 0.2687825904738485, 0.20791169081776067,
0.14622397450298433, 0.08396179228009795, 0.0213697517873014,
-0.04130624343253929, -0.10381996000604393, -0.16592580209926813,
-0.2273797762829596, -0.28794045010251973, -0.3473699005866891,
-0.4054346489686925, -0.4619065779475033, -0.5165638278857105,
-0.569191668423076, -0.6195833420814977, -0.6675408765471482,
-0.7128758624385862, -0.7554101935052618, -0.7949767663483981,
-0.831420136915289, -0.8645971311878683, -0.8943774076663372,
-0.9206439694380529, -0.9432936238199248, -0.9622373877685285,
-0.9774008374652203, -0.9887244007028445, -0.9961635909253453,
-0.9996891820008162, -0.9992873230413596, -0.9949595928186701,
-0.9867229935615607, -0.9746098841597972, -0.9586678530366604,
-0.9389595311896838, -0.9155623461340598, -0.8885682177153923,
-0.8580831969868598, -0.8242270495694983, -0.7871327851324629,
-0.7469461348417843, -0.7038249788305342, -0.6579387259397131,
-0.6094676481666195, -0.5586021724354684, -0.5055421324726097,
-0.4504959837255544, -0.3936799844100698, -0.3353173459027642,
-0.27563735581699855, -0.21487447720725855, -0.15326742744091648,
-0.09105824035620122, -0.028491315390845918, 0.0341875425829326,
0.09673208894492837, 0.158896606741102, 0.22043687202466516,
0.2811111133316544, 0.3406809615215133, 0.3989123862510416,
0.455576615402663, 0.5104510338548267, 0.5633200580636206,
0.6139759830196183, 0.6622197982525427, 0.7078619696779649,
0.750723184214399, 0.7906350542454076, 0.8274407791591957,
0.8609957613666314, 0.8911681743776243, 0.9178394807040574,
0.9409048975545885, 0.9602738084917828, 0.9758701194342947,
0.9876325576054856, 0.9955149122540149, 0.9994862162006879,
0.9995308674983134, 0.995648690726628, 0.9878549376814615,
0.9761802274554481, 0.960670426145682, 0.9413864666609029,
0.9184041093361411, 0.8918136442952654, 0.8617195367307948,
0.8282400164945015, 0.7915066136112013, 0.7516636415405575,
0.7088676302169401, 0.6632867110947972, 0.615099956615462,
0.5644966766904359, 0.5116756749649672, 0.45684446778395543,
0.4002184689284223, 0.34202014332566893, 0.2824781330576857,
0.22182635910162835, 0.16030310233122602, 0.09815006738958451,
0.035611433111129936, -0.027067106776728054, -0.08963930890343103,
-0.151859347664545, -0.21348278098965978, -0.2742675106749298,
-0.33397473350729945, -0.39236987944377644, -0.44922353315999364,
-0.5043123353475932, -0.5574198602194665, -0.6083374657755144,
-0.6568651134884359, -0.702812154189337, -0.7459980770656012,
-0.7862532188285168, -0.8234194302645385, -0.8573506975515395,
-0.8879137158991022, -0.9149884132591922, -0.9384684220497617,
-0.9582614970379841, -0.9742898777414488, -0.9864905939235211,
-0.9948157129826897, -0.9992325282639968, -0.9997236875527314,
-0.9962872612455728, -0.988936749931366, -0.9777010313517511,
-0.9626242469500127, -0.9437656284538599, -0.921199265173458,
-0.8950138129288846, -0.8653121457505724, -0.8322109517210728,
-0.7958402745459633, -0.756343002654873, -0.7138743078398655,
-0.6686010356364607, -0.6207010498424329, -0.570362533749413,
-0.5177832508326333, -0.4631697678032811, -0.4067366430758031,
-0.34870558383846073, -0.28930457503870327, -0.22876698370530607,
-0.16733064212600626, -0.10523691348271454, -0.04272974361491648,
0.019945297363298146, 0.08254197982744962, 0.1448143819978928,
0.20651785608561196, 0.26740998943292893, 0.3272515568731176,
0.3858074605673982, 0.4428476536270481, 0.4981480438919785,
0.5514913743150955, 0.6026680764938378, 0.6514770939954935,
0.6977266722417961, 0.7412351118496027, 0.7818314824680316,
0.8193562943075774, 0.853662124723059, 0.8846141973887109,
0.912090911790027, 0.9359843209522112, 0.9562005555283323,
0.9726601925811199, 0.9852985676095651, 0.9940660285944805,
0.9989281310649515, 0.9998657734193406, 0.9968752719691923,
0.9899683754112395, 0.9791722186706318, 0.964529216296742,
0.9460968958303518, 0.9239476717968579, 0.8981685612134207,
0.8688608417277132, 0.8361396537313613, 0.8001335480112053,
0.7609839807155286, 0.7188447576193768, 0.6738814298722485,
0.62627064360207, 0.5761994459306503, 0.5238645501270671,
0.4694715627858897, 0.4132341760665302, 0.35537332816698547,
0.2961163353303466, 0.23569599879401382, 0.17434969019019303,
0.1123184189908455, 0.04984588566069994, -0.012822475761727311,
-0.07544046189372612, -0.13776206726054913, -0.19954245077000266,
-0.26053889761258414, -0.3205117728081487, -0.37922546265292867,
-0.43644930036830915, -0.49195847231472983, -0.5455349012105506,
-0.5969681028859549, -0.6460560132060685, -0.6926057819144502,
-0.7364345302783594, -0.7773700695591514, -0.8152515774851836,
-0.8499302300696602, -0.8812697862911105, -0.9091471233395767,
-0.9334527203256257, -0.954091088551925, -0.9709811466569138,
-0.9840565391568189, -0.9932658971345277, -0.9985730400511654,
-0.9999571178875374, -0.9974126930569968, -0.9909497617679348,
-0.9805937147519644, -0.9663852375120869, -0.9483801504827333,
-0.9266491897296226, -0.9012777290510405, -0.8723654445722485,
-0.8400259231507723, -0.8043862161309906, -0.7655863402011965,
-0.7237787273140511, -0.6791276258316024, -0.631808455247462,
-0.5820071170213084, -0.5299192642332037, -0.4757495329269612,
-0.4197107381624439, -0.36202303793517704, -0.30291306824797226,
-0.24261305273263473, -0.1813598903196625, -0.11939422454024687,
-0.05695949811699878, 0.005699003442465559, 0.06833511549159066,
0.13070276134486217, 0.19255691903216382, 0.2536545839095089,
0.313755723344368, 0.3726242197249783, 0.43002879808874256,
0.485743934725536, 0.5395507431861802, 0.5912378342153216,
0.6406021462303005, 0.6874497430833184, 0.7315965759728014,
0.7728692065106026, 0.8111054881043979, 0.8461552029783131,
0.8778806523291781, 0.9061571972998319, 0.930873748644203,
0.9519332031604308, 0.9692528251774035, 0.9827645715960043,
0.9924153592080726, 0.9981672732428718, 0.9999977163217538,
0.9978994972358264, 0.9918808591978509, 0.9819654474573699,
0.9681922164062949, 0.9506152765399138, 0.9293036818745531,
0.9043411586570455, 0.8758257764318569, 0.8438695627580838,
0.808598063090023, 0.770149847550363, 0.7286759665337543,
0.6843393572794684, 0.6373142037445788]Available Oscillator Types
Fastmath supports: :sin, :square, :saw, :triangle, :noise, :constant
(def sample-rate 100.0)(def duration 1.0)(def oscillator-examples
{:sine (sig/oscillator :sin 3.0 1.0 0.0)
:square (sig/oscillator :square 3.0 1.0 0.0)
:saw (sig/oscillator :saw 3.0 1.0 0.0)
:triangle (sig/oscillator :triangle 3.0 1.0 0.0)})Visualize different types
(def osc-viz-data
(for [[osc-type osc] oscillator-examples
[i value] (map-indexed vector
(take 50 (sig/oscillator->signal osc sample-rate duration)))]
{:time (/ i sample-rate)
:amplitude value
:type (name osc-type)}))(-> (tc/dataset osc-viz-data)
(plotly/base {:=x :time
:=y :amplitude
:=color :type
:=x-title "Time (s)"
:=y-title "Amplitude"
:=title "Oscillator Types (3 Hz)"
:=width 700
:=height 250})
(plotly/layer-line))Combining Multiple Frequencies
Use oscillators-sum for composite signals: Three sine waves: 3 Hz, 7 Hz, and 15 Hz
(def composite-osc
(sig/oscillators-sum
(sig/oscillator :sin 3.0 1.0 0.0)
(sig/oscillator :sin 7.0 0.5 0.0)
(sig/oscillator :sin 15.0 0.3 0.0)))(def composite-signal
(sig/oscillator->signal composite-osc sample-rate duration))Teaching Signals: Designed for Transform Comparison
Now let’s create a suite of test signals, each designed to highlight different transform strengths. These aren’t arbitrary—each signal has specific characteristics that make certain transforms more effective than others.
Why multiple test signals? Different transforms excel at different tasks. The FFT easily finds exact frequencies in pure tones. The DCT concentrates energy efficiently for smooth signals (good for compression). Wavelets capture time-localized features like transient events. By testing each transform on all signals, we’ll build intuition for when to use which tool.
(defn teaching-signals
"Standard test signals with known properties.
Each signal is designed to test different transform characteristics:
- :pure-sine - Single frequency (FFT baseline)
- :two-tones - Multiple frequencies (tests superposition)
- :smooth - Low frequency content (DCT compression test)
- :chirp - Time-varying frequency (wavelet localization test)
- :impulse - Sharp transient (wavelet vs FFT comparison)"
[]
{:pure-sine
{:signal (sig/oscillator->signal
(sig/oscillator :sin 10.0 1.0 0.0)
sample-rate duration)
:description "Single frequency - FFT baseline"
:best-transform :fft}
:two-tones
{:signal (sig/oscillator->signal
(sig/oscillators-sum
(sig/oscillator :sin 10.0 1.0 0.0)
(sig/oscillator :sin 25.0 0.5 0.0))
sample-rate duration)
:description "Two frequencies - superposition principle"
:best-transform :fft}
:smooth
{:signal (sig/oscillator->signal
(sig/oscillators-sum
(sig/oscillator :sin 5.0 1.0 0.0)
(sig/oscillator :sin 10.0 0.3 0.0)
(sig/oscillator :sin 15.0 0.1 0.0))
sample-rate duration)
:description "Smooth signal - DCT compresses well"
:best-transform :dct}
:chirp
{:signal (let [n (int (* sample-rate duration))]
(double-array
(for [i (range n)]
(let [t (/ i sample-rate)
freq (+ 5.0 (* 45.0 t))] ;; 5→50 Hz sweep
(Math/sin (* 2 Math/PI freq t))))))
:description "Frequency sweep - wavelets capture time-varying content"
:best-transform :wavelet}
:impulse
{:signal (double-array
(concat [5.0] (repeat 99 0.0)))
:description "Localized event - wavelets excel here"
:best-transform :wavelet}})(def signals (teaching-signals))Signal Visualization with dtype-next
Now we’ll visualize signals using dtype-next for data preparation and Tableplot for rendering. This is where dtype-next shines: we can efficiently extract slices and perform element-wise division to compute time values.
(defn visualize-signal
"Plot signal using dtype-next for data preparation.
dtype-next operations used:
- dtype/ecount: Get array length
- dtype/sub-buffer: Zero-copy slice extraction
- dfn//: Element-wise division for time computation"
[signal title]
(let [display-n (min 200 (dtype/ecount signal))
time-vals (dfn// (range display-n) sample-rate) ; Element-wise division
amp-vals (dtype/sub-buffer signal 0 display-n) ; Zero-copy slice
dataset (tc/dataset {:time time-vals
:amplitude amp-vals})]
(-> dataset
(plotly/base {:=x :time
:=y :amplitude
:=x-title "Time (s)"
:=y-title "Amplitude"
:=title title
:=width 700
:=height 200})
(plotly/layer-line))))(visualize-signal (:signal (:two-tones signals)) "Two Tones: 10 Hz + 25 Hz")What we just learned: Signal generation with fastmath oscillators and basic dtype-next operations (element-wise math, slicing, counting). We now have test signals that will help us compare different transforms.
Next, we’ll explore the FFT—but first, we need to understand why its output is complex-valued even though our input signals are real numbers.
Understanding Complex Transform Outputs
Before we dive into the FFT, we need to address a puzzling fact: when we transform a real-valued signal (just regular numbers like 1.5, -0.3, 2.7), the FFT returns complex numbers (numbers with “real” and “imaginary” parts).
Why would a transform of real data produce complex output? This seems unnecessarily complicated. Let’s build intuition for why complex numbers are not just mathematical elegance—they’re computational necessity.
Note on dtype-next operations: In the examples below, you’ll see operations like dfn/* and dfn//. These are element-wise operations from dtype-next:
dfn/*multiplies arrays element-by-element:[1 2 3] * 2 → [2 4 6]dfn//divides arrays element-by-element:[2 4 6] / 2 → [1 2 3]
dfn/+,dfn/-work similarly for addition and subtractiondfn/cos,dfn/sinapply trigonometric functions to each element
We’ll explore these operations more deeply as we use them throughout the tutorial.
The Problem: Sines and Cosines Are Both Needed
Here’s the core issue. Suppose we want to decompose a signal into pure frequency components. Should we use sines or cosines as our basis functions?
If we use only cosines: We can represent cosine-like signals perfectly, but sine-like signals require infinite cosine terms
If we use only sines: We can represent sine-like signals perfectly, but cosine-like signals require infinite sine terms
Solution: Use BOTH sines and cosines. For each frequency, we need two numbers: 1. How much of the cosine wave at that frequency 2. How much of the sine wave at that frequency
A complex number packages these two values together: the real part stores the cosine amplitude, and the imaginary part stores the sine amplitude. This isn’t arbitrary—it’s the most compact way to store both components.
Euler’s Formula: The Mathematical Foundation
The mathematical underpinning is Euler’s formula, which connects complex exponentials to sine and cosine:
e^(iθ) = cos(θ) + i·sin(θ)
Where i = √(-1) is the imaginary unit. This beautiful equation tells us that a rotating complex number naturally encodes both sine and cosine components.
Let’s visualize this to make it concrete.
Visualizing Sine and Cosine Together
(def theta-points (dfn// (dfn/* 2.0 Math/PI (range 100)) 100.0))Real part: cos(θ)
(def cosine-vals (dfn/cos theta-points))Imaginary part: sin(θ)
(def sine-vals (dfn/sin theta-points))Create dataset for visualization
(def euler-dataset
(tc/dataset {:angle (dfn// theta-points (* 2 Math/PI))
:cosine cosine-vals
:sine sine-vals}))Using Tableplot (ggplot2-style layered grammar)
(-> euler-dataset
(plotly/base {:=x :angle
:=y :cosine
:=x-title "Angle (cycles)"
:=y-title "Value"
:=title "Euler's Formula: e^(iθ) = cos(θ) + i·sin(θ)"
:=width 700
:=height 300})
(plotly/layer-line {:=color "steelblue"
:=name "cos(θ) - Real Part"})
(plotly/layer-line {:=y :sine
:=color "orange"
:=name "sin(θ) - Imaginary Part"}))As the angle θ increases from 0 to 2π (one full rotation), cosine and sine trace out their familiar wave shapes. The key insight: we need both functions to completely describe a sinusoid at any phase.
Concrete Example: Pure Cosine vs Pure Sine Signals
Let’s demonstrate why we need complex output by transforming pure cosine and pure sine signals.
Create a simple 5 Hz cosine wave
(def pure-cosine-time (dfn// (range 100) 100.0))(def pure-cosine-signal
(dfn/cos (dfn/* 2.0 Math/PI 5.0 pure-cosine-time)))FFT gives complex output (stored as alternating real, imaginary pairs)
(def fft-basic (t/transformer :real :fft))(def cosine-spectrum (t/forward-1d fft-basic pure-cosine-signal))The spectrum has 100 values: 50 complex numbers stored as [real₀, imag₀, real₁, imag₁, …]
(dtype/ecount cosine-spectrum)100Look at frequency bin 5 (the 5 Hz component we put in) Large value in real part, ≈0 in imaginary part
(def freq-bin-5-real (dtype/get-value cosine-spectrum (* 2 5)))(def freq-bin-5-imag (dtype/get-value cosine-spectrum (+ 1 (* 2 5))))Pure cosine = nearly all energy in REAL part (cosine basis function) Imaginary part ≈ 0 (no sine component needed, apart from numerical noise)
Now try a pure sine wave
(def pure-sine-signal
(dfn/sin (dfn/* 2.0 Math/PI 5.0 pure-cosine-time)))(def sine-spectrum (t/forward-1d fft-basic pure-sine-signal))≈0 in real part, large value in imaginary part
(def sine-freq-bin-5-real (dtype/get-value sine-spectrum (* 2 5)))(def sine-freq-bin-5-imag (dtype/get-value sine-spectrum (+ 1 (* 2 5))))Pure sine = nearly all energy in IMAGINARY part (sine basis function) Real part ≈ 0 (no cosine component needed, apart from numerical noise)
Key Takeaway:
- Pure cosine signal → FFT has large real part, near-zero imaginary part
- Pure sine signal → FFT has near-zero real part, large imaginary part
- General signal → FFT has both parts, encoding the phase relationship
This is why FFT output is complex even though input is real!
Phase: Encoding the Sine/Cosine Balance
Any sinusoid can be written two ways: as amplitude + phase (A·cos(2πft + φ)) or as cosine + sine (a·cos(2πft) + b·sin(2πft)). The complex number (a + ib) from the FFT encodes both representations. The magnitude √(a² + b²) gives us the amplitude A, and the phase arctan(b/a) gives us φ.
This is why we compute magnitude as √(real² + imag²)—we’re combining the cosine and sine contributions to find the total amplitude at each frequency.
Visualization: Complex Plane Representation
(defn complex-plane-viz
"Visualize how cosine and sine components form complex numbers."
[signal-type]
(let [freq 5.0
t-vals (dfn// (range 100) 100.0)
signal (case signal-type
:cosine (dfn/cos (dfn/* 2.0 Math/PI freq t-vals))
:sine (dfn/sin (dfn/* 2.0 Math/PI freq t-vals))
:mixed (dfn/+ (dfn/cos (dfn/* 2.0 Math/PI freq t-vals))
(dfn/sin (dfn/* 2.0 Math/PI freq t-vals))))
spectrum (t/forward-1d fft-basic signal)
;; Extract complex numbers for positive frequencies
freq-bins (range 0 50)
complex-data (map (fn [k]
(let [real-part (dtype/get-value spectrum (* 2 k))
imag-part (dtype/get-value spectrum (+ 1 (* 2 k)))]
{:frequency k
:real real-part
:imag imag-part
:magnitude (Math/sqrt (+ (* real-part real-part)
(* imag-part imag-part)))}))
freq-bins)
dataset (tc/dataset complex-data)]
(-> dataset
(plotly/base {:=x :real
:=y :imag
:=x-title "Real Part (Cosine)"
:=y-title "Imaginary Part (Sine)"
:=title (str "Complex Spectrum: " (name signal-type) " signal")
:=width 400
:=height 400})
(plotly/layer-point {:=mark-size :magnitude
:=color :frequency}))))(complex-plane-viz :cosine)(complex-plane-viz :sine)(complex-plane-viz :mixed)What we just learned: Complex FFT outputs aren’t mysterious—they’re the natural way to encode both sine and cosine components at each frequency. The real part tells us “how much cosine,” the imaginary part tells us “how much sine.”
Now that we understand the complex output format, let’s put the FFT to work!
DFT and FFT - Discovering Frequencies in Signals
We’re finally ready to answer our opening question: how do we write code to discover the frequencies in a signal? The answer is the Discrete Fourier Transform (DFT), computed efficiently via the Fast Fourier Transform (FFT) algorithm.
Terminology Note:
- DFT (Discrete Fourier Transform) = the mathematical transform
- FFT (Fast Fourier Transform) = the Cooley-Tukey algorithm for computing DFT efficiently (O(n log n) instead of O(n²))
In practice, “FFT” is often used colloquially to mean both the transform and the algorithm. We’ll be technically precise here, but don’t be surprised if you see papers and codebases using “FFT” for both!
Introduction to fastmath.transform
(require '[fastmath.transform :as t])The fastmath.transform library provides a unified API for all signal transforms. The pattern is consistent across transforms, making it easy to experiment:
Step 1: Create a transformer (specifies input type and transform algorithm)
(def fft (t/transformer :real :fft)):real means input is real-valued (not complex) :fft specifies the algorithm
Step 2: Forward transform (signal → spectrum)
(t/forward-1d fft (range 16))[120.0, -8.0, -7.999999999999995, 40.218715937006785, -8.0,
19.31370849898476, -8.0, 11.972846101323913, -8.0, 8.0, -8.0,
5.345429103354391, -8.0, 3.313708498984761, -8.000000000000004,
1.591298939037264]Step 3: Reverse transform (spectrum → signal) (t/reverse-1d fft spectrum)
This same pattern works for DCT, DST, DHT, and wavelets—just change the transform type!
Performance: JTransforms and Parallelism
Under the hood, fastmath uses JTransforms, a high-performance Java library that automatically parallelizes large transforms.
For large arrays (typically >8-16K elements), JTransforms automatically splits work across available CPU cores using Java’s ForkJoinPool. It uses SIMD optimizations via Java’s auto-vectorization and can operate in-place directly on input arrays (fastmath wraps this safely). It also selects the optimal FFT algorithm (split-radix, mixed-radix, or Bluestein) based on array size.
You don’t need to think about parallelism—just pass large arrays and JTransforms automatically uses available CPU cores. For smaller signals, the overhead of parallelization exceeds the benefit, so it runs single-threaded.
Performance tip: Reuse transformer objects! The constructor pre-computes lookup tables (twiddle factors), so creating a transformer once and reusing it is much faster than creating a new one for each transform.
Understanding the DFT: Signal Decomposition
The Discrete Fourier Transform (DFT) answers the question: “What frequencies are in this signal, and how strong is each one?”
Mathematically, it expresses any signal as a weighted sum of sine and cosine waves:
Signal = a₁·cos(2π·f₁·t) + b₁·sin(2π·f₁·t) + a₂·cos(2π·f₂·t) + b₂·sin(2π·f₂·t) + …
The DFT finds the weights (amplitudes) a₁, b₁, a₂, b₂, etc. We compute it using the Fast Fourier Transform (FFT) algorithm, which reduces complexity from O(n²) to O(n log n).
Let’s see this in action by building a signal from known frequencies, then recovering them.
Concrete Example: Build and Decompose
Step 1: Manual construction - we choose the frequencies and amplitudes
(def freq-1 10.0)(def amp-1 1.0)(def freq-2 25.0)(def amp-2 0.5)Generate time samples using dtype-next
(def n-samples 1000)(def time-points (dfn// (range n-samples) sample-rate))Build each frequency component
(def component-1
(dfn/* amp-1 (dfn/sin (dfn/* 2.0 Math/PI freq-1 time-points))))(def component-2
(dfn/* amp-2 (dfn/sin (dfn/* 2.0 Math/PI freq-2 time-points))))The signal is just the sum (linear superposition!)
(def constructed-signal
(dfn/+ component-1 component-2))Step 2: FFT recovers the decomposition
(def fft-transformer (t/transformer :real :fft))(def spectrum (t/forward-1d fft-transformer constructed-signal))Step 3: Extract frequency peaks with dtype-next
(defn fft-magnitude
"Compute magnitudes from FFT spectrum using dtype-next.
FFT returns interleaved [r0, i0, r1, i1, ...] format.
Magnitude = sqrt(real² + imag²)"
[spectrum]
(let [n-bins (quot (dtype/ecount spectrum) 2)]
;; Functional approach: map over bin indices
(dtype/emap
(fn [i]
(let [real (dtype/get-value spectrum (* 2 i))
imag (dtype/get-value spectrum (inc (* 2 i)))]
(Math/sqrt (+ (* real real) (* imag imag)))))
:float64
(range n-bins))))(def magnitudes (fft-magnitude spectrum))(def freq-bins (dfn/* (range (dtype/ecount magnitudes)) (/ sample-rate n-samples)))Find peaks
(def peak-1-idx (argops/argmax magnitudes))(def peak-1-freq (dtype/get-value freq-bins peak-1-idx))(def peak-1-mag (dtype/get-value magnitudes peak-1-idx))Remove first peak, find second (functional approach)
(def mags-without-peak-1
(dtype/emap (fn [i val] (if (= i peak-1-idx) 0.0 val))
:float64
(range (dtype/ecount magnitudes))
magnitudes))(def peak-2-idx (argops/argmax mags-without-peak-1))(def peak-2-freq (dtype/get-value freq-bins peak-2-idx))(def peak-2-mag (dtype/get-value mags-without-peak-1 peak-2-idx))Verification: FFT recovered our exact construction!
{:constructed {:freq-1 freq-1 :amp-1 amp-1
:freq-2 freq-2 :amp-2 amp-2}
:recovered {:freq-1 peak-1-freq :magnitude-1 (/ peak-1-mag (/ n-samples 2))
:freq-2 peak-2-freq :magnitude-2 (/ peak-2-mag (/ n-samples 2))}}{:constructed {:freq-1 10.0, :amp-1 1.0, :freq-2 25.0, :amp-2 0.5},
:recovered
{:freq-1 10.0,
:magnitude-1 0.9999999999999997,
:freq-2 25.0,
:magnitude-2 0.49999999999999994}}Key Insight: This is linear decomposition in action! We built: signal = 1.0×sin(2π·10·t) + 0.5×sin(2π·25·t) FFT found: peaks at 10 Hz (mag ≈ 1.0) and 25 Hz (mag ≈ 0.5)
The FFT perfectly recovered the frequencies we put in. This is the power of the Fourier Transform—it reveals the hidden frequency structure of any signal.
Practical Function: Finding Dominant Frequency
Now let’s package this into a reusable function. This is the kind of utility you’d use in real applications: “What’s the main frequency in this signal?”
(defn find-dominant-frequency
"Extract strongest frequency using FFT and dtype-next."
[signal sample-rate]
(let [fft (t/transformer :real :fft)
spectrum (t/forward-1d fft signal)
n (dtype/ecount signal)
mags (fft-magnitude spectrum)
;; Skip DC component
max-idx (inc (argops/argmax (dtype/sub-buffer mags 1 (dec (dtype/ecount mags)))))
freq (* max-idx (/ sample-rate n))]
{:frequency freq
:magnitude (dtype/get-value mags max-idx)}))(find-dominant-frequency (:signal (:pure-sine signals)) sample-rate){:frequency 10.0, :magnitude 50.0}FFT Visualization
(defn visualize-fft
"Show frequency spectrum."
[signal sample-rate title]
(let [fft (t/transformer :real :fft)
spectrum (t/forward-1d fft signal)
mags (fft-magnitude spectrum)
n (dtype/ecount signal)
freqs (dfn/* (range (dtype/ecount mags)) (/ sample-rate n))
display-n (min 100 (dtype/ecount mags))
dataset (tc/dataset {:frequency (dtype/sub-buffer freqs 0 display-n)
:magnitude (dtype/sub-buffer mags 0 display-n)})]
(-> dataset
(plotly/base {:=x :frequency
:=y :magnitude
:=x-title "Frequency (Hz)"
:=y-title "Magnitude"
:=title title
:=width 700
:=height 250})
(plotly/layer-line))))(visualize-fft (:signal (:two-tones signals)) sample-rate
"FFT Spectrum: Clear Peaks at 10 Hz and 25 Hz")Round-Trip Testing - Validation
(defn test-fft-roundtrip
"Verify FFT → IFFT recovers original signal.
The 'tolerance' parameter sets how much error we accept due to
floating-point rounding - typically very small (1e-10 or less).
Returns test results with pass/fail boolean."
[signal tolerance]
(let [fft (t/transformer :real :fft)
original (vec signal)
spectrum (t/forward-1d fft signal)
recovered (t/reverse-1d fft spectrum)
;; Use dtype-next for error calculation
error (Math/sqrt (dfn/mean (dfn/sq (dfn/- signal recovered))))]
{:rmse error
:tolerance tolerance
:passed? (< error tolerance)
:interpretation (if (< error tolerance)
"✓ Perfect reconstruction"
"✗ Error exceeds tolerance")}))Test with our teaching signals
(def fft-tests
{:pure-sine (test-fft-roundtrip (:signal (:pure-sine signals)) 1e-10)
:two-tones (test-fft-roundtrip (:signal (:two-tones signals)) 1e-10)
:smooth (test-fft-roundtrip (:signal (:smooth signals)) 1e-10)})| signal | rmse | passed? | status |
|---|---|---|---|
| pure-sine | 1.69e-16 | true | ✓ PASS |
| two-tones | 1.97e-16 | true | ✓ PASS |
| smooth | 1.80e-16 | true | ✓ PASS |
Result: All signals pass - FFT perfectly recovers the original!
This round-trip property is crucial. It means we can: 1. Transform signal → frequency domain 2. Manipulate frequencies (filter, compress, analyze) 3. Transform back → time domain 4. Get our original signal back (within the limits of floating-point arithmetic)
But before we use FFT in practice, we need to understand its limitations and artifacts.
Understanding FFT Artifacts and Solutions
The FFT examples above worked perfectly because we used “nice” signals—exact integer frequencies (10 Hz, 25 Hz) that fit evenly into our sample window. Real signals aren’t so cooperative. Let’s explore what happens when signals don’t perfectly align with the FFT’s assumptions, and how to fix the resulting artifacts.
The Nyquist Frequency and Aliasing
Before we can understand FFT artifacts, we need to establish a fundamental limit of digital signal processing: the Nyquist frequency.
The Nyquist-Shannon Sampling Theorem states: to accurately represent a signal, you must sample at more than twice the highest frequency present. The Nyquist frequency is sample-rate / 2—the maximum frequency that can be unambiguously represented.
Why does this matter for FFT?
Reason 1: Spectrum Symmetry
When you FFT a real-valued signal, the output has Hermitian symmetry—the second half mirrors the first half. The meaningful frequency content is only in the first half, from 0 Hz to Nyquist (sample-rate/2).
Example: 100 Hz sample rate → FFT bins represent 0 to 50 Hz
{:sample-rate sample-rate
:nyquist-frequency (/ sample-rate 2)
:interpretation "Only analyze frequencies up to Nyquist (50 Hz)"}{:sample-rate 100.0,
:nyquist-frequency 50.0,
:interpretation "Only analyze frequencies up to Nyquist (50 Hz)"}This is why in our earlier visualize-fft function we used: (min 50 (dtype/ecount mags))
We’re limiting display to roughly the first half of the spectrum, where meaningful content lives.
Reason 2: Aliasing
What if your signal contains frequencies above Nyquist? They don’t disappear—they alias (fold back) into lower frequencies, creating false peaks.
Demonstrate aliasing:
(def nyquist-demo-rate 100.0)(def below-nyquist-signal
(sig/oscillator->signal
(sig/oscillator :sin 40.0 1.0 0.0) ; 40 Hz, below Nyquist (50 Hz)
nyquist-demo-rate 1.0))(def above-nyquist-signal
(sig/oscillator->signal
(sig/oscillator :sin 80.0 1.0 0.0) ; 80 Hz, above Nyquist (50 Hz)
nyquist-demo-rate 1.0))Find dominant frequencies
(def below-result (find-dominant-frequency below-nyquist-signal nyquist-demo-rate))(def above-result (find-dominant-frequency above-nyquist-signal nyquist-demo-rate)){:below-nyquist {:input-freq 40.0 :detected-freq (:frequency below-result) :correct? true}
:above-nyquist {:input-freq 80.0 :detected-freq (:frequency above-result)
:aliased-to (- nyquist-demo-rate 80.0)
:correct? false
:explanation "80 Hz aliases to 20 Hz (100 - 80)"}}{:below-nyquist
{:input-freq 40.0, :detected-freq 40.0, :correct? true},
:above-nyquist
{:input-freq 80.0,
:detected-freq 20.0,
:aliased-to 20.0,
:correct? false,
:explanation "80 Hz aliases to 20 Hz (100 - 80)"}}Key Takeaway: Always ensure your signal is bandlimited below Nyquist, or use an anti-aliasing filter before sampling.
Spectral Leakage - The Picket Fence Effect
Now for the most important artifact: spectral leakage. Our earlier examples used frequencies like 10 Hz and 25 Hz. What if the signal contains 10.3 Hz—a frequency that doesn’t fall exactly on an FFT bin?
The Problem: FFT assumes the signal is periodic within the window. If the frequency doesn’t complete an integer number of cycles, the signal appears discontinuous when wrapped around. This discontinuity creates artificial high-frequency content that “leaks” into adjacent bins.
Demonstrate with a frequency that doesn’t align with FFT bins
(def leakage-sample-rate 100.0)(def leakage-n-samples 100)Perfect bin alignment: 10 Hz completes exactly 10 cycles in 100 samples
(def aligned-signal
(sig/oscillator->signal
(sig/oscillator :sin 10.0 1.0 0.0)
leakage-sample-rate (/ leakage-n-samples leakage-sample-rate)))Poor alignment: 10.3 Hz doesn’t complete an integer number of cycles
(def leaky-signal
(sig/oscillator->signal
(sig/oscillator :sin 10.3 1.0 0.0)
leakage-sample-rate (/ leakage-n-samples leakage-sample-rate)))Compare spectra
(defn spectrum-comparison [signal1 signal2 label1 label2 sample-rate]
(let [fft (t/transformer :real :fft)
spec1 (fft-magnitude (t/forward-1d fft signal1))
spec2 (fft-magnitude (t/forward-1d fft signal2))
n (dtype/ecount signal1)
freqs (dfn/* (range (quot (dtype/ecount spec1) 2))
(/ sample-rate n))
dataset (tc/dataset
(concat
(map (fn [i] {:frequency (dtype/get-value freqs i)
:magnitude (dtype/get-value spec1 i)
:signal label1})
(range (dtype/ecount freqs)))
(map (fn [i] {:frequency (dtype/get-value freqs i)
:magnitude (dtype/get-value spec2 i)
:signal label2})
(range (dtype/ecount freqs)))))]
(-> dataset
(plotly/base {:=x :frequency
:=y :magnitude
:=color :signal
:=x-title "Frequency (Hz)"
:=y-title "Magnitude"
:=title "Spectral Leakage: Aligned vs Misaligned Frequency"
:=width 700
:=height 300})
(plotly/layer-line))))(spectrum-comparison aligned-signal leaky-signal
"10.0 Hz (aligned)" "10.3 Hz (leaky)"
leakage-sample-rate)Observation:
- 10.0 Hz (aligned): Clean single peak—all energy in one bin
- 10.3 Hz (leaky): Energy spreads across many bins—the “picket fence effect”
This is spectral leakage. The true frequency (10.3 Hz) falls between FFT bins, so energy leaks into adjacent bins through the sidelobes of the rectangular window.
Windowing Functions - The Solution to Leakage
The root cause of leakage is the rectangular window—our signal is multiplied by a rectangle (1 inside the window, 0 outside). This sharp cutoff creates the discontinuity.
Solution: Multiply the signal by a smoother window function that tapers to zero at the edges. This reduces the discontinuity and suppresses sidelobes.
Common window functions:
- Rectangular: No windowing (what we’ve been using)
- Hann (Hanning): Good general-purpose window, smooth taper
- Hamming: Similar to Hann, slightly different sidelobe tradeoff
- Blackman: Lower sidelobes, wider main lobe
Implement window functions
(defn hann-window
"Generate Hann window of length n."
[n]
(dtype/emap
(fn [i]
(let [angle (* 2.0 Math/PI (/ i n))]
(* 0.5 (- 1.0 (Math/cos angle)))))
:float64
(range n)))(defn hamming-window
"Generate Hamming window of length n."
[n]
(dtype/emap
(fn [i]
(let [angle (* 2.0 Math/PI (/ i n))]
(- 0.54 (* 0.46 (Math/cos angle)))))
:float64
(range n)))(defn blackman-window
"Generate Blackman window of length n."
[n]
(dtype/emap
(fn [i]
(let [angle (* 2.0 Math/PI (/ i n))]
(- 0.42
(* 0.5 (Math/cos angle))
(* 0.08 (Math/cos (* 2.0 angle))))))
:float64
(range n)))Visualize window shapes
(def window-n 100)(def window-viz-data
(concat
(map (fn [i] {:index i :amplitude 1.0 :window "Rectangular"}) (range window-n))
(map-indexed (fn [i v] {:index i :amplitude v :window "Hann"})
(hann-window window-n))
(map-indexed (fn [i v] {:index i :amplitude v :window "Hamming"})
(hamming-window window-n))
(map-indexed (fn [i v] {:index i :amplitude v :window "Blackman"})
(blackman-window window-n))))(-> (tc/dataset window-viz-data)
(plotly/base {:=x :index
:=y :amplitude
:=color :window
:=x-title "Sample Index"
:=y-title "Window Amplitude"
:=title "Window Functions - Taper to Zero at Edges"
:=width 700
:=height 300})
(plotly/layer-line))Key Observation: All non-rectangular windows taper to zero at edges, reducing discontinuity when the signal wraps around.
Applying Windows to Reduce Leakage
(defn apply-window
"Multiply signal by window function."
[signal window]
(dfn/* signal window))Apply Hann window to our leaky signal
(def windowed-leaky-signal
(apply-window leaky-signal (hann-window (dtype/ecount leaky-signal))))Compare: no window vs Hann window
(spectrum-comparison leaky-signal windowed-leaky-signal
"No window (rectangular)" "Hann window"
leakage-sample-rate)Result: The Hann window dramatically reduces sidelobe leakage! Energy is concentrated in the main lobe around 10.3 Hz, with much less spillover into distant bins.
Window Tradeoffs
Windows aren’t free—they trade frequency resolution for leakage suppression.
| window | main-lobe-width | sidelobe-level | use-case |
|---|---|---|---|
| Rectangular | Narrowest | Highest (-13 dB) | Perfect bin alignment |
| Hann | Medium | Medium (-32 dB) | General purpose |
| Hamming | Medium | Medium (-43 dB) | Better sidelobe suppression |
| Blackman | Widest | Lowest (-58 dB) | Maximum sidelobe suppression |
Choosing a window:
- Known exact frequency: Rectangular (no windowing needed)
- Unknown frequency: Hann (good default)
- Nearby strong signals: Blackman (best isolation)
- Need precise frequency: Rectangular or Hamming (narrow main lobe)
Frequency Resolution vs Window Size
The FFT’s frequency resolution depends on window size and sample rate:
Frequency resolution = sample-rate / window-size
Demonstrate resolution tradeoff
(defn resolution-demo [window-size sample-rate]
(let [;; Two close frequencies: 10 Hz and 12 Hz
signal (sig/oscillator->signal
(sig/oscillators-sum
(sig/oscillator :sin 10.0 1.0 0.0)
(sig/oscillator :sin 12.0 1.0 0.0))
sample-rate (/ window-size sample-rate))
fft (t/transformer :real :fft)
spectrum (t/forward-1d fft signal)
mags (fft-magnitude spectrum)
n (dtype/ecount signal)
freqs (dfn/* (range (min 30 (dtype/ecount mags))) (/ sample-rate n))
freq-resolution (/ sample-rate n)]
{:window-size window-size
:frequency-resolution freq-resolution
:can-resolve-2hz? (< freq-resolution 2.0)
:peaks (take 5 (sort-by :magnitude >
(map (fn [i] {:freq (dtype/get-value freqs i)
:magnitude (dtype/get-value mags i)})
(range (dtype/ecount freqs)))))}))| window-size | frequency-resolution | can-resolve-2hz? | peaks |
|---|---|---|---|
| 50 | 2.0 | false | |
| 100 | 1.0 | true | |
| 200 | 0.5 | true | |
0.5 Hz resolution - easily resolves them
Key Insight:
- Longer window → better frequency resolution, worse time resolution
- Shorter window → better time resolution, worse frequency resolution
This is the fundamental uncertainty principle in time-frequency analysis—you can’t have perfect resolution in both domains simultaneously.
Practical Guidelines for FFT Analysis
| question | answer |
|---|---|
| Do I need a window? | Always use windowing unless frequency exactly matches an FFT bin |
| Which window? | Hann for general use, Blackman for strong nearby signals |
| What window size? | Larger for better frequency resolution, smaller for time localization |
| What about DC and Nyquist? | Skip bin 0 (DC), only analyze 0 to sample-rate/2 (Nyquist) |
| How to avoid aliasing? | Ensure signal has no content above sample-rate/2 |
Improved Frequency Analysis with Windowing
(defn find-dominant-frequency-windowed
"Extract strongest frequency using FFT with Hann window."
[signal sample-rate]
(let [;; Apply Hann window to reduce leakage
window (hann-window (dtype/ecount signal))
windowed-signal (dfn/* signal window)
fft (t/transformer :real :fft)
spectrum (t/forward-1d fft windowed-signal)
n (dtype/ecount signal)
mags (fft-magnitude spectrum)
;; Only search up to Nyquist frequency
nyquist-bin (quot n 2)
;; Skip DC component (bin 0)
max-idx (inc (argops/argmax (dtype/sub-buffer mags 1 (dec nyquist-bin))))
freq (* max-idx (/ sample-rate n))]
{:frequency freq
:magnitude (dtype/get-value mags max-idx)
:frequency-resolution (/ sample-rate n)
:nyquist-frequency (/ sample-rate 2)}))Test on our leaky signal (10.3 Hz)
(find-dominant-frequency-windowed leaky-signal leakage-sample-rate){:frequency 10.0,
:magnitude 23.582009001846984,
:frequency-resolution 1.0,
:nyquist-frequency 50.0}Compare to original (no windowing)
(find-dominant-frequency leaky-signal leakage-sample-rate){:frequency 10.0, :magnitude 42.70305711359141}Improvement: Windowed version gives more accurate peak detection and magnitude!
What we just learned: Real-world FFT requires understanding Nyquist limits, spectral leakage, and windowing. Always window your data unless you know frequencies align perfectly with FFT bins. Choose window size based on the frequency resolution vs time resolution tradeoff your application needs.
Now let’s apply these insights to build robust practical tools.
Practical Application: Notch Filter
(defn notch-filter
"Remove specific frequency using FFT."
[signal sample-rate target-freq bandwidth]
(let [fft (t/transformer :real :fft)
spectrum (t/forward-1d fft signal)
n (dtype/ecount signal)
;; Functional approach: map over spectrum indices
filtered (dtype/emap
(fn [idx val]
(let [freq (* (quot idx 2) (/ sample-rate n))]
(if (< (Math/abs (- freq target-freq)) bandwidth)
0.0
val)))
:float64
(range (dtype/ecount spectrum))
spectrum)]
(t/reverse-1d fft filtered)))(def signal-with-25hz (:signal (:two-tones signals)))(def signal-25hz-removed (notch-filter signal-with-25hz sample-rate 25.0 2.0))Verify 25 Hz is gone, 10 Hz remains
(find-dominant-frequency signal-25hz-removed sample-rate){:frequency 10.0, :magnitude 45.1942088429381}Time-Frequency Analysis: Spectrograms
Limitation of FFT: It tells us WHAT frequencies are present, but not WHEN. For time-varying signals (like chirps), we need to see how frequency content changes over time.
Solution: Short-Time Fourier Transform (STFT) - apply FFT to overlapping windows.
Note on windowing: In production spectrograms, you should apply a Hann or Hamming window to each FFT window (as discussed in the FFT Artifacts section above) to reduce spectral leakage. We omit it here for simplicity, but the pattern is:
(apply-window (dtype/sub-buffer signal start-idx window-size) (hann-window window-size))
Window size choice: Remember the tradeoff from the artifacts section:
- Small windows (e.g., 32) → good time resolution, poor frequency resolution
- Large windows (e.g., 512) → poor time resolution, good frequency resolution
(defn spectrogram
"Compute spectrogram using sliding-window FFT.
Returns sequence of {:time t :frequency f :magnitude m} maps.
For production use, apply windowing to each segment:
(apply-window (dtype/sub-buffer ...) (hann-window window-size))"
[signal sample-rate window-size hop-size]
(let [fft (t/transformer :real :fft)
n-windows (quot (- (dtype/ecount signal) window-size) hop-size)]
(for [win-idx (range n-windows)
:let [start-idx (* win-idx hop-size)
window (dtype/sub-buffer signal start-idx window-size)
spectrum (t/forward-1d fft window)
mags (fft-magnitude spectrum)
time-center (/ (+ start-idx (/ window-size 2)) sample-rate)]
;; Only show up to Nyquist frequency
freq-idx (range (min 50 (quot (dtype/ecount mags) 2)))]
{:time time-center
:frequency (* freq-idx (/ sample-rate window-size))
:magnitude (dtype/get-value mags freq-idx)})))Visualize chirp signal (frequency increases over time)
(def chirp-spectrogram
(spectrogram (:signal (:chirp signals)) sample-rate 32 8))(-> (tc/dataset chirp-spectrogram)
(plotly/base {:=x :time
:=y :frequency
:=color :magnitude
:=x-title "Time (s)"
:=y-title "Frequency (Hz)"
:=title "Spectrogram: Chirp Signal (5→50 Hz sweep)"
:=width 700
:=height 350})
(plotly/layer-point {:=mark-symbol "square"}))Key Observation: The spectrogram reveals the frequency sweep! Brighter colors = higher magnitude. We can see frequency increasing linearly over time.
Limitation of FFT: It tells us WHAT frequencies exist, but loses information about WHEN they occur. For signals with time-varying content, we need different tools—wavelets provide this naturally, as we’ll see later.
What we just learned: The FFT decomposes signals into frequency components, enabling frequency analysis, filtering, and spectral visualization. Its invertibility makes it perfect for round-trip transformations.
But the FFT has a weakness: it’s not optimal for compression. Let’s explore why the Discrete Cosine Transform (DCT) outperforms FFT for smooth signals.
DCT - The Compression Transform
The FFT is excellent for finding frequencies, but it’s not the best choice for compression. Here’s why: the FFT uses both sines AND cosines (complex output), which means we’re storing twice as much information as we might need.
For smooth signals—like natural images, audio, or our “smooth” teaching signal—most energy concentrates in low frequencies. The Discrete Cosine Transform (DCT) uses only cosines, which for smooth signals concentrates energy even more efficiently than FFT.
This is why classic JPEG (baseline) compresses images with DCT, not FFT. Let’s see why.
Understanding DCT as Cosine Decomposition
Core Concept: The Discrete Cosine Transform (DCT) expresses a signal as a weighted sum of cosine waves.
Signal = a₀ + a₁·cos(π·1·t) + a₂·cos(π·2·t) + …
Unlike DFT (sines + cosines), DCT uses only cosines. This makes it perfect for smooth signals → energy concentrates in first few coefficients.
Why DCT for Compression?
Smooth signals have most energy in low frequencies. DCT with cosine-only basis captures this more efficiently than DFT. Result: JPEG (baseline), MP3 (using Modified DCT), and many modern codecs use DCT or DCT variants!
Energy Concentration Example
(def dct-n-samples 64)(def dct-time (dfn// (range dct-n-samples) dct-n-samples))Create smooth signal (mostly low frequency)
(def smooth-signal
(dfn/+ 10.0
(dfn/* 5.0 (dfn/cos (dfn/* 2.0 Math/PI 2.0 dct-time)))
(dfn/* 2.0 (dfn/cos (dfn/* 2.0 Math/PI 3.0 dct-time)))))Apply DCT
(def dct (t/transformer :real :dct))(def dct-coeffs (t/forward-1d dct smooth-signal))Measure energy concentration with dtype-next
(def energy-by-coeff (dfn/sq dct-coeffs))(def total-energy (dfn/sum energy-by-coeff))(def energy-first-10 (dfn/sum (dtype/sub-buffer energy-by-coeff 0 10))){:total-coeffs dct-n-samples
:energy-in-first-10 (format "%.1f%%" (* 100.0 (/ energy-first-10 total-energy)))
:interpretation "Smooth signal → energy concentrated in few DCT coefficients!"}{:total-coeffs 64,
:energy-in-first-10 "100.0%",
:interpretation
"Smooth signal → energy concentrated in few DCT coefficients!"}DCT vs FFT Energy Concentration
(defn measure-energy-concentration
"Compare energy concentration across transforms."
[signal transform-type keep-ratios]
(let [transformer (t/transformer :real transform-type)
coeffs (t/forward-1d transformer signal)
total-energy (dfn/sum (dfn/sq coeffs))]
(for [ratio keep-ratios]
(let [keep-n (int (* ratio (dtype/ecount coeffs)))
kept-energy (dfn/sum (dfn/sq (dtype/sub-buffer coeffs 0 keep-n)))]
{:keep-ratio (* 100 ratio)
:energy-retained (* 100 (/ kept-energy total-energy))
:transform transform-type}))))(def energy-comparison
(concat
(measure-energy-concentration (:signal (:smooth signals)) :fft [0.1 0.2 0.3 0.5])
(measure-energy-concentration (:signal (:smooth signals)) :dct [0.1 0.2 0.3 0.5])))(-> (tc/dataset energy-comparison)
(plotly/base {:=x :keep-ratio
:=y :energy-retained
:=color :transform
:=x-title "% Coefficients Kept"
:=y-title "% Energy Retained"
:=title "Energy Concentration: DCT vs FFT on Smooth Signal"
:=width 600
:=height 300})
(plotly/layer-line)
(plotly/layer-point))Observation: DCT retains ~95% energy with 10% of coefficients! FFT needs ~30% of coefficients for same energy.
Lossy Compression with Quality Metrics
(defn compress-with-dct
"Compress signal by keeping only top N% of DCT coefficients."
[signal keep-ratio]
(let [dct (t/transformer :real :dct)
coeffs (t/forward-1d dct signal)
n (dtype/ecount coeffs)
keep-n (int (* keep-ratio n))
;; Functional approach: zero out high-frequency coefficients
compressed (dtype/emap
(fn [i val] (if (< i keep-n) val 0.0))
:float64
(range n)
coeffs)]
(t/reverse-1d dct compressed)))(def original-smooth (:signal (:smooth signals)))(def compressed-50pct (compress-with-dct original-smooth 0.5))(def compressed-25pct (compress-with-dct original-smooth 0.25))(def compressed-10pct (compress-with-dct original-smooth 0.1))Calculate quality metrics with dtype-next
(defn compression-quality [original compressed keep-ratio]
(let [error (Math/sqrt (dfn/mean (dfn/sq (dfn/- original compressed))))
signal-power (Math/sqrt (dfn/mean (dfn/sq original)))
;; Quality measure: ratio of signal power to error (higher = better quality)
;; Measured in [decibels](https://en.wikipedia.org/wiki/Decibel) (dB) - a logarithmic scale where 20 dB ≈ 10x better
snr-db (* 20 (Math/log10 (/ signal-power error)))]
{:compression-ratio (format "%.0f:1" (/ 1.0 keep-ratio))
:quality-db (format "%.1f dB" snr-db)
:keep-ratio (format "%.0f%%" (* 100 keep-ratio))}))| compression-ratio | quality-db | keep-ratio |
|---|---|---|
| 2:1 | 27.1 dB | 50% |
| 4:1 | 17.9 dB | 25% |
| 10:1 | 3.8 dB | 10% |
Test DCT Energy Concentration Property
(defn test-dct-energy-concentration
"Verify DCT concentrates energy better than FFT for smooth signals."
[signal keep-ratio min-energy-threshold]
(let [dct (t/transformer :real :dct)
fft (t/transformer :real :fft)
;; DCT energy
dct-coeffs (t/forward-1d dct signal)
dct-total (dfn/sum (dfn/sq dct-coeffs))
dct-keep-n (int (* keep-ratio (dtype/ecount dct-coeffs)))
dct-kept (dfn/sum (dfn/sq (dtype/sub-buffer dct-coeffs 0 dct-keep-n)))
dct-ratio (/ dct-kept dct-total)
;; FFT energy
fft-spec (t/forward-1d fft signal)
fft-mags (fft-magnitude fft-spec)
fft-total (dfn/sum (dfn/sq fft-mags))
fft-keep-n (int (* keep-ratio (dtype/ecount fft-mags)))
fft-kept (dfn/sum (dfn/sq (dtype/sub-buffer fft-mags 0 fft-keep-n)))
fft-ratio (/ fft-kept fft-total)]
{:dct-energy-retained dct-ratio
:fft-energy-retained fft-ratio
:dct-advantage (format "%.1fx better" (/ dct-ratio fft-ratio))
:passed? (and (> dct-ratio fft-ratio)
(> dct-ratio min-energy-threshold))}))Test that DCT concentrates energy better than FFT For smooth signals, keeping 10% of DCT coefficients should retain significantly more energy than FFT
(test-dct-energy-concentration original-smooth 0.1 0.5){:dct-energy-retained 0.581972554902848,
:fft-energy-retained 3.370527878784732E-31,
:dct-advantage "1726651064262024200000000000000.0x better",
:passed? true}What we just learned: The DCT uses only cosines (instead of FFT’s sines + cosines), which concentrates energy more efficiently for smooth signals. This makes it ideal for compression—classic JPEG and MP3 (via Modified DCT) both exploit this property.
But both FFT and DCT have a fundamental limitation: they use global basis functions that span the entire signal. A sine wave at 10 Hz in the FFT basis oscillates across the whole signal—it can’t tell us that a frequency appears only briefly. Let’s explore wavelets, which solve this problem.
Wavelets - Time-Frequency Localization
Remember our chirp signal—frequency increasing from 5 to 50 Hz over time? The FFT told us “this signal contains frequencies between 5 and 50 Hz,” but it couldn’t tell us WHEN each frequency occurred. The spectrogram helped by windowing, but that was a workaround.
Wavelets provide time-frequency localization naturally by using basis functions that are inherently localized in time. Instead of global sine waves, wavelets use small, localized “wavelets” (hence the name!) that can detect features at specific times AND specific frequencies.
Understanding Wavelets as Localized Decomposition
Core Concept: Wavelets decompose a signal using scaled and shifted basis functions.
Signal = a₁·ψ(t-1) + a₂·ψ(2t-1) + a₃·ψ(2t-2) + …
⚠️ Important: Most wavelet transforms (Haar, Daubechies, Symlets) require power-of-2 array sizes (64, 128, 256, 512, etc.). This is because wavelets work by recursively splitting signals in half. If your signal isn’t a power of 2, you’ll need to pad it.
We’ll show the padding pattern in our examples below.
^kind/md “Padding to power-of-2 size:
;; Pad to next power of 2
(let [n (dtype/ecount signal)
target-n (int (Math/pow 2 (Math/ceil (/ (Math/log n) (Math/log 2)))))]
;; ... pad signal to target-n ...)
```"
Where ψ is a localized wavelet function at different scales and positions.
See also: [Wavelet Transform](https://en.wikipedia.org/wiki/Wavelet_transform)
### The Key Difference: Global vs Local
**DFT/DCT**: Basis functions span entire signal
- Answer: "What frequencies are present?"
- Cannot tell WHERE events occur
**Wavelets**: Localized basis at multiple scales
- Answer: "What frequencies, and WHEN?"
- Perfect for time-varying signals
### Visualizing Wavelet Basis Functions
To understand the difference, let's visualize the actual basis functions.
Haar wavelet basis function (simplest wavelet)
::: {.sourceClojure}
```clojure
(defn haar-wavelet
"Generate Haar wavelet at given position and scale."
[t position scale]
(let [t-scaled (/ (- t position) scale)]
(cond
(and (>= t-scaled 0) (< t-scaled 0.5)) 1.0
(and (>= t-scaled 0.5) (< t-scaled 1.0)) -1.0
:else 0.0))):::
Generate basis functions for comparison
(def basis-comparison-data
(let [n 200
t-vals (dfn// (range n) n)]
(concat
;; Sine basis (global)
(map (fn [t] {:time t :amplitude (Math/sin (* 2 Math/PI 5 t)) :type "Sine (global)"})
t-vals)
;; Cosine basis (global)
(map (fn [t] {:time t :amplitude (Math/cos (* 2 Math/PI 5 t)) :type "Cosine (global)"})
t-vals)
;; Haar wavelet at scale 0.1 (fine detail)
(map (fn [t] {:time t :amplitude (haar-wavelet t 0.3 0.1) :type "Haar scale=0.1 (fine)"})
t-vals)
;; Haar wavelet at scale 0.2 (coarse detail)
(map (fn [t] {:time t :amplitude (haar-wavelet t 0.6 0.2) :type "Haar scale=0.2 (coarse)"})
t-vals))))(-> (tc/dataset basis-comparison-data)
(plotly/base {:=x :time
:=y :amplitude
:=color :type
:=x-title "Time"
:=y-title "Amplitude"
:=title "Basis Functions: Global (Sine/Cosine) vs Local (Haar Wavelets)"
:=width 700
:=height 350})
(plotly/layer-line))Key Observation: Sine and cosine extend across the entire time axis (global). Haar wavelets are localized to specific time windows and can be scaled to capture both fine details (small scale) and coarse features (large scale).
Example: Localized Event Detection
Localized event signal: quiet, brief pulse, quiet
(def event-signal
(double-array
(concat (repeat 25 0.0)
(repeat 5 5.0)
(repeat 20 0.0)
(repeat 14 0.0))))FFT spreads event across all frequencies (global view)
(def event-fft (t/forward-1d (t/transformer :real :fft) event-signal))Wavelet localizes event in time AND frequency
(def haar (t/transformer :fast :haar))(def event-wavelet (t/forward-1d haar event-signal))Compare: where is the energy?
{:fft-view "Energy distributed across all frequencies - WHERE is unclear"
:wavelet-view "Energy concentrated at specific time-scale location"
:key-insight "Wavelets = localized basis functions at multiple scales"}{:fft-view
"Energy distributed across all frequencies - WHERE is unclear",
:wavelet-view "Energy concentrated at specific time-scale location",
:key-insight "Wavelets = localized basis functions at multiple scales"}Wavelet Families
Fastmath supports many wavelet families via JWave:
- Haar:
:haar - Daubechies:
:daubechies-4,:daubechies-8 - Symlets:
:symlet-4 - Coiflets:
:coiflet-2
(defn compare-wavelets [signal]
(let [;; Ensure signal is power of 2
n (dtype/ecount signal)
target-n (int (Math/pow 2 (Math/ceil (/ (Math/log n) (Math/log 2)))))
;; Functional padding: use original values or 0
padded (dtype/emap
(fn [i] (if (< i n)
(dtype/get-value signal i)
0.0))
:float64
(range target-n))
wavelets {:haar :haar
:daubechies-4 :daubechies-4
:daubechies-8 :daubechies-8
:symlet-4 :symlet-4}]
(for [[name wavelet-key] wavelets]
(let [trans (t/transformer :fast wavelet-key)
coeffs (t/forward-1d trans padded)
;; Sparsity = how many small coefficients
sparsity (count (filter #(< (Math/abs %) 0.01) coeffs))]
{:wavelet (clojure.core/name name)
:sparsity-pct (format "%.0f%%" (* 100.0 (/ sparsity (alength coeffs))))
:max-coeff (format "%.2f" (dfn/reduce-max (dfn/abs coeffs)))}))))| wavelet | sparsity-pct | max-coeff |
|---|---|---|
| haar | 23% | 1.51 |
| daubechies-4 | 17% | 1.73 |
| daubechies-8 | 19% | 2.39 |
| symlet-4 | 16% | 1.89 |
Observation: Different wavelets give different sparsity. Smoother wavelets (Daubechies, Symlet) → more zeros → better compression!
Wavelet Denoising
(defn add-noise [signal noise-level]
;; Functional approach: generate noise and add to signal
(let [n (dtype/ecount signal)
noise (dtype/emap (fn [_] (* noise-level (- (rand) 0.5)))
:float64
(range n))]
(dfn/+ signal noise)))(def clean-sig (sig/oscillator->signal (sig/oscillator :sin 5.0 1.0 0.0) 128.0 1.0))(def noisy-sig (add-noise clean-sig 0.3))Wavelet denoising: transform → threshold → inverse
Important: Denoising requires careful threshold tuning. The default threshold may be too aggressive (removing signal) or too weak (leaving noise). In practice, you’d analyze the wavelet coefficients and choose a threshold based on the noise characteristics. This is an active area of research with methods like:
- Universal threshold: σ√(2 log N)
- SURE threshold (Stein’s Unbiased Risk Estimate)
- Empirical Bayes methods
(def db8 (t/transformer :fast :daubechies-8))(def noisy-coeffs (t/forward-1d db8 noisy-sig))(def denoised-coeffs (t/denoise db8 noisy-coeffs :soft))(def denoised-sig (t/reverse-1d db8 denoised-coeffs))Quality comparison with dtype-next
(defn signal-quality [clean noisy denoised]
{:noisy-rmse (format "%.4f" (Math/sqrt (dfn/mean (dfn/sq (dfn/- clean noisy)))))
:denoised-rmse (format "%.4f" (Math/sqrt (dfn/mean (dfn/sq (dfn/- clean denoised)))))
:note "Denoising effectiveness depends on threshold tuning"})(signal-quality clean-sig noisy-sig denoised-sig){:noisy-rmse "0.0901",
:denoised-rmse "0.4080",
:note "Denoising effectiveness depends on threshold tuning"}Round-Trip Testing for Wavelets
(defn test-wavelet-roundtrip
"Verify wavelet transform is invertible."
[signal wavelet-type tolerance]
(let [wavelet (t/transformer :fast wavelet-type)
coeffs (t/forward-1d wavelet signal)
recovered (t/reverse-1d wavelet coeffs)
error (Math/sqrt (dfn/mean (dfn/sq (dfn/- signal recovered))))]
{:wavelet wavelet-type
:rmse error
:passed? (< error tolerance)}))(def wavelet-tests
[(test-wavelet-roundtrip (take 64 (:signal (:smooth signals))) :haar 1e-10)
(test-wavelet-roundtrip (take 64 (:signal (:smooth signals))) :daubechies-8 1e-10)
(test-wavelet-roundtrip (take 64 (:signal (:smooth signals))) :symlet-4 1e-10)])| wavelet | rmse | passed? |
|---|---|---|
| haar | 4.911240235747681E-16 | true |
| daubechies-8 | 7.675184834973428E-12 | true |
| symlet-4 | 9.536902080421727E-13 | true |
Other Transforms - DST and DHT
Discrete Sine Transform (DST)
Use case: Solving PDEs with Dirichlet boundary conditions (signal = 0 at endpoints)
The Discrete Sine Transform (DST) is like DCT but uses sine basis instead of cosine. Implicitly assumes signal goes to zero at boundaries.
(def dst (t/transformer :real :dst))Signal that naturally goes to zero at edges Signal naturally zero at boundaries (t=0 and t=1)
(def boundary-signal
(let [n 64
t (dfn// (range n) n)]
(dfn/* (dfn/sin (dfn/* Math/PI t))
(dfn/sin (dfn/* 5.0 Math/PI t)))))(def dst-coeffs (t/forward-1d dst boundary-signal))(def dst-reconstructed (t/reverse-1d dst dst-coeffs))Verify round-trip
{:rmse (format "%.2e" (Math/sqrt (dfn/mean (dfn/sq (dfn/- boundary-signal dst-reconstructed)))))
:passed? (< (Math/sqrt (dfn/mean (dfn/sq (dfn/- boundary-signal dst-reconstructed)))) 1e-10)}{:rmse "0.00e+00", :passed? true}Discrete Hartley Transform (DHT)
Use case: Real-valued alternative to DFT for convolution
The Discrete Hartley Transform (DHT) is similar to the DFT but output is real-valued (no complex numbers). Useful for applications where you only care about real operations.
(def dht (t/transformer :real :dht))(def test-sig (take 64 (:signal (:pure-sine signals))))(def dht-result (t/forward-1d dht test-sig))All real values (unlike FFT which has complex output)
{:all-real? (every? #(not (Double/isNaN %)) dht-result)
:invertible? (< (Math/sqrt (dfn/mean (dfn/sq (dfn/- test-sig (t/reverse-1d dht dht-result))))) 1e-10)}{:all-real? true, :invertible? true}Transform Selection Guide
| transform | use-case | basis | output | algorithm |
|---|---|---|---|---|
| DFT | Frequency analysis | Sine + Cosine | Complex | FFT |
| DCT | Compression (JPEG/MP3) | Cosine only | Real | Fast DCT |
| DST | PDEs, boundary problems | Sine only | Real | Fast DST |
| DHT | Real-valued convolution | Hartley | Real | Fast DHT |
| Wavelet | Time-frequency, denoising | Localized | Real | Fast Wavelet |
2D Transforms for Images
Creating a 2D Signal
(defn generate-2d-signal [rows cols]
(let [data (make-array Double/TYPE rows cols)]
(dotimes [i rows]
(dotimes [j cols]
;; Create pattern: horizontal 3Hz + vertical 5Hz
(aset data i j
(+ (Math/sin (* 2 Math/PI 3.0 (/ j cols)))
(Math/sin (* 2 Math/PI 5.0 (/ i rows)))))))
data))(def signal-2d (generate-2d-signal 32 32))Visualize 2D signal
(-> (tc/dataset (for [i (range 32) j (range 32)]
{:row i :col j :value (aget signal-2d i j)}))
(plotly/base {:=x :col
:=y :row
:=color :value
:=title "2D Signal: Horizontal 3Hz + Vertical 5Hz"
:=width 400
:=height 400})
(plotly/layer-point {:=mark-symbol "square"}))2D FFT
(def fft-2d (t/transformer :real :fft))(def spectrum-2d (t/forward-2d fft-2d signal-2d))Compute 2D magnitude
(defn magnitude-2d [spectrum]
(let [rows (alength spectrum)
cols (alength (aget spectrum 0))
result (make-array Double/TYPE rows cols)]
(dotimes [i rows]
(dotimes [j cols]
(aset result i j (Math/abs (aget spectrum i j)))))
result))(def mag-2d (magnitude-2d spectrum-2d))Visualize 2D spectrum (log scale for visibility)
(-> (tc/dataset (for [i (range 32) j (range 32)]
{:row i :col j :value (Math/log (+ 1 (aget mag-2d i j)))}))
(plotly/base {:=x :col
:=y :row
:=color :value
:=x-title "Frequency X"
:=y-title "Frequency Y"
:=title "2D FFT Spectrum - Bright Spots Show Frequency Components"
:=width 400
:=height 400})
(plotly/layer-point {:=mark-symbol "square"}))Interpretation: Just like 1D FFT, 2D FFT decomposes image into frequency components. The bright spots in the spectrum correspond to our horizontal 3Hz and vertical 5Hz patterns!
2D DCT Basis Functions
Understanding how baseline JPEG compression works: visualizing the 8×8 DCT basis images. Each 8×8 block in baseline JPEG is decomposed into a weighted sum of these 64 basis patterns.
2D DCT-II basis function
(defn dct-2d-basis
"Generate a single 2D DCT basis function for frequency (u, v)."
[size u v]
(let [data (make-array Double/TYPE size size)
pi Math/PI
cu (if (zero? u) (/ 1.0 (Math/sqrt 2.0)) 1.0)
cv (if (zero? v) (/ 1.0 (Math/sqrt 2.0)) 1.0)]
(dotimes [i size]
(dotimes [j size]
(aset data i j
(* cu cv
(Math/cos (* (/ pi size) u (+ i 0.5)))
(Math/cos (* (/ pi size) v (+ j 0.5)))))))
data))Generate a grid of low-frequency DCT basis images (like JPEG uses)
(def dct-basis-grid
(let [size 8
grid-size 4]
(for [u (range grid-size)
v (range grid-size)
i (range size)
j (range size)]
{:freq-u u
:freq-v v
:row i
:col j
:value (aget (dct-2d-basis size u v) i j)})))Visualize the first 16 DCT basis functions (0-3 in each direction)
(-> (tc/dataset dct-basis-grid)
(plotly/base {:=x :col
:=y :row
:=color :value
:=facet-x :freq-u
:=facet-y :freq-v
:=title "2D DCT Basis Functions (Baseline JPEG uses these!)"
:=width 600
:=height 600})
(plotly/layer-point {:=mark-symbol "square"}))Key Insight:
- Top-left (0,0) = DC component (constant, average brightness)
- Moving right → increasing horizontal frequency
- Moving down → increasing vertical frequency
- Baseline JPEG keeps low-frequency basis (top-left) and discards high-frequency (bottom-right)
- This is why JPEG compression works well for natural images (most energy in low frequencies)
Practical Applications
Now that we understand the theory behind FFT, DCT, and wavelets, let’s see how to apply these transforms to solve real-world problems. We’ll build practical tools for filtering, compression, and spectral analysis.
Application 1: Bandpass Filtering
(defn bandpass-filter
"Keep only frequencies in [low-freq, high-freq] range."
[signal low-freq high-freq sample-rate]
(let [fft (t/transformer :real :fft)
spectrum (t/forward-1d fft signal)
n (dtype/ecount signal)
low-bin (int (* n (/ low-freq sample-rate)))
high-bin (int (* n (/ high-freq sample-rate)))
;; Functional approach: keep frequencies in range, zero others
filtered (dtype/emap
(fn [idx val]
(let [bin (quot idx 2)]
(if (and (>= bin low-bin) (<= bin high-bin))
val
0.0)))
:float64
(range (dtype/ecount spectrum))
spectrum)]
(t/reverse-1d fft filtered)))Test: remove high frequency noise Mixed signal: 8 Hz (signal) + 25 Hz (noise) + random noise
(def mixed-signal
(add-noise
(sig/oscillator->signal
(sig/oscillators-sum
(sig/oscillator :sin 8.0 1.0 0.0)
(sig/oscillator :sin 25.0 0.3 0.0))
100.0 1.0)
0.1))(def filtered-signal (bandpass-filter mixed-signal 5.0 15.0 100.0)){:original-dominant (find-dominant-frequency mixed-signal 100.0)
:filtered-dominant (find-dominant-frequency filtered-signal 100.0)
:interpretation "High frequency (25 Hz) successfully removed!"}{:original-dominant {:frequency 8.0, :magnitude 47.9344063286437},
:filtered-dominant {:frequency 8.0, :magnitude 47.59174688032908},
:interpretation "High frequency (25 Hz) successfully removed!"}Application 2: Spectral Analysis
(defn analyze-spectrum
"Find all significant frequency components."
[signal sample-rate threshold-ratio]
(let [fft (t/transformer :real :fft)
spectrum (t/forward-1d fft signal)
mags (fft-magnitude spectrum)
n (dtype/ecount signal)
max-mag (dfn/reduce-max mags)
threshold (* threshold-ratio max-mag)
;; Find all peaks above threshold
peaks (filter #(> (dtype/get-value mags %) threshold) (range (dtype/ecount mags)))]
{:frequencies (mapv #(* % (/ sample-rate n)) peaks)
:magnitudes (mapv #(dtype/get-value mags %) peaks)
:num-components (count peaks)}))(analyze-spectrum (:signal (:two-tones signals)) sample-rate 0.1){:frequencies [10.0 25.0 40.0],
:magnitudes [45.19420884293811 19.611582314123716 4.805791157061837],
:num-components 3}Application 3: Compression Pipeline
(defn compression-pipeline
"Complete compression: signal → DCT → quantize → compress."
[signal keep-ratio]
(let [dct (t/transformer :real :dct)
;; Step 1: Transform to frequency domain
coeffs (t/forward-1d dct signal)
;; Step 2: Sort by magnitude, keep largest
indexed-coeffs (map-indexed vector coeffs)
sorted-coeffs (sort-by #(- (Math/abs (second %))) indexed-coeffs)
keep-n (int (* keep-ratio (count sorted-coeffs)))
kept-indices (set (map first (take keep-n sorted-coeffs)))
;; Step 3: Functional approach - zero out small coefficients
compressed (dtype/emap
(fn [idx val]
(if (kept-indices idx) val 0.0))
:float64
(range (dtype/ecount coeffs))
coeffs)]
;; Step 4: Inverse transform
(t/reverse-1d dct compressed)))Compare compression strategies: sequential vs magnitude-based
(def original (:signal (:smooth signals)))(def compressed-sequential (compress-with-dct original 0.25))(def compressed-magnitude (compression-pipeline original 0.25)){:sequential-rmse (format "%.4f" (Math/sqrt (dfn/mean (dfn/sq (dfn/- original compressed-sequential)))))
:magnitude-rmse (format "%.4f" (Math/sqrt (dfn/mean (dfn/sq (dfn/- original compressed-magnitude)))))
:note "Keeping largest coefficients works better than keeping first N!"}{:sequential-rmse "0.0910",
:magnitude-rmse "0.0297",
:note
"Keeping largest coefficients works better than keeping first N!"}Testing and Validation Framework
Transform operations can seem like black boxes. How do we know they’re working correctly? This section builds a comprehensive testing framework to validate transform properties: invertibility (can we get the original back?), stability across different signal sizes, and mathematical correctness.
Utility Functions
(defn compute-rmse
"Root mean square error using dtype-next."
[signal1 signal2]
(Math/sqrt (dfn/mean (dfn/sq (dfn/- signal1 signal2)))))(defn generate-random-signal
"Random signal with seeded RNG for reproducibility."
[n seed]
(let [rng (java.util.Random. seed)]
(double-array (repeatedly n #(- (* 2.0 (.nextDouble rng)) 1.0)))))Generic Round-Trip Test
(defn test-transform-roundtrip
"Generic test for any transform's invertibility."
[signal transformer tolerance]
(let [forward (t/forward-1d transformer signal)
inverse (t/reverse-1d transformer forward)
error (compute-rmse signal inverse)]
{:signal-size (dtype/ecount signal)
:rmse error
:tolerance tolerance
:passed? (< error tolerance)}))Test Multiple Signal Types
(defn test-all-signal-types
"Test transform on various signal types."
[transformer tolerance]
(let [test-signals
{:pure-sine (sig/oscillator->signal (sig/oscillator :sin 440.0 1.0 0.0) 1000.0 0.1)
:composite (sig/oscillator->signal
(sig/oscillators-sum
(sig/oscillator :sin 220.0 1.0 0.0)
(sig/oscillator :sin 440.0 0.5 0.0))
1000.0 0.1)
:square (sig/oscillator->signal (sig/oscillator :square 100.0 1.0 0.0) 1000.0 0.1)
:noise (generate-random-signal 100 42)
:chirp (:signal (:chirp signals))}]
(into {}
(for [[sig-type signal] test-signals]
[sig-type (test-transform-roundtrip signal transformer tolerance)]))))Run comprehensive FFT tests
(def fft-comprehensive (test-all-signal-types (t/transformer :real :fft) 1e-10))| signal | rmse | passed? | status |
|---|---|---|---|
| pure-sine | 2.20e-16 | true | ✓ |
| composite | 1.96e-16 | true | ✓ |
| square | 1.71e-16 | true | ✓ |
| noise | 1.17e-16 | true | ✓ |
| chirp | 1.93e-16 | true | ✓ |
Numerical Stability Test
(defn test-numerical-stability
"Test across wide range of amplitudes - verifies transform works equally
well on tiny signals (1e-6) and huge signals (1e6)."
[base-signal transformer scale-factors]
(for [scale scale-factors]
(let [scaled-signal (dfn/* scale base-signal)
;; Use relative tolerance - bigger signals need bigger error allowance
tolerance (* 1e-10 scale)
result (test-transform-roundtrip scaled-signal transformer tolerance)]
(assoc result :scale-factor scale))))(def stability-results
(test-numerical-stability
(take 100 (:signal (:pure-sine signals)))
(t/transformer :real :fft)
[1e-6 1e-3 1.0 1e3 1e6]))All should pass
(every? :passed? stability-results)trueComprehensive Test Suite
(defn comprehensive-test-suite
"Run all tests and return boolean: true if all pass."
[]
(let [fft (t/transformer :real :fft)
dct (t/transformer :real :dct)
haar (t/transformer :fast :haar)
test-signal (take 64 (:signal (:smooth signals)))
results
{:fft-roundtrip (every? :passed? (vals (test-all-signal-types fft 1e-10)))
:dct-roundtrip (every? :passed? (vals (test-all-signal-types dct 1e-10)))
:wavelet-roundtrip (:passed? (test-wavelet-roundtrip test-signal :haar 1e-10))
:numerical-stability (every? :passed? stability-results)
:dct-energy-concentration (:passed? (test-dct-energy-concentration original-smooth 0.1 0.5))}]
{:individual-results results
:all-passed? (every? identity (vals results))}))(def final-test-results (comprehensive-test-suite))(:all-passed? final-test-results)trueBest Practices and Decision Guide
dtype-next Best Practices Summary
Key Functions Actually Used:
| category | functions |
|---|---|
| Element-wise | dfn/+, dfn/-, dfn/*, dfn//, dfn/sq, dfn/sqrt, dfn/cos, dfn/sin, dfn/abs |
| Reductions | dfn/sum, dfn/mean, dfn/standard-deviation, dfn/reduce-max, dfn/reduce-min |
| Functional | dtype/emap - map function over arrays |
| Array access | dtype/ecount, dtype/get-value, dtype/sub-buffer |
| Search | argops/argmax - find index of maximum |
Common dtype-next Patterns - Quick Reference
Throughout the tutorial, you’ve seen these dtype-next patterns repeatedly. Here’s a quick reference with working examples you can use in your own code.
Measuring signal difference - compute how far apart two signals are:
(defn compute-error [signal1 signal2]
;; This computes the "root mean square error" - square the differences,
;; average them, then take the square root. Smaller = more similar.
(Math/sqrt (dfn/mean (dfn/sq (dfn/- signal1 signal2)))))Example: Compare original and noisy signal
(let [original (range 100)
noisy (dfn/+ original 0.1)]
(compute-error original noisy))0.0999999999999986Signal Energy - total power in a signal:
(defn signal-energy [signal]
(dfn/sum (dfn/sq signal)))Example: Energy before and after filtering
(signal-energy (:signal (:pure-sine signals)))50.0Normalization - scale signal to have average of zero and consistent spread:
(defn normalize [signal]
;; Subtract the average (centers signal around zero)
;; Then divide by standard deviation (makes the spread consistent)
;; Useful for comparing signals with different scales
(dfn// (dfn/- signal (dfn/mean signal))
(dfn/standard-deviation signal)))Example: Normalize before processing
(let [normalized (normalize (:signal (:two-tones signals)))]
{:mean (dfn/mean normalized)
:spread (dfn/standard-deviation normalized)}){:mean 3.552713678800501E-17, :spread 0.9999999999999997}Signal-to-Noise Ratio (SNR) in dB:
(defn snr-db [signal noise]
;; Measures signal quality: how much stronger is the signal than the noise?
;; Higher numbers = better quality. Measured in decibels (dB).
;; See: https://en.wikipedia.org/wiki/Signal-to-noise_ratio
(* 20 (Math/log10 (/ (dfn/mean (dfn/sq signal))
(dfn/mean (dfn/sq noise))))))^kind/md “Example usage - measuring compression quality:
(let [signal (:signal (:smooth signals))
noise (dfn/- signal compressed-signal)]
(snr-db signal noise))
```"
### Transform Selection Flowchart
::: {.sourceClojure}
```clojure
(defn recommend-transform
"Decision tree for transform selection."
[signal-properties]
(cond
;; Time-varying frequency
(:time-varying-frequency signal-properties)
{:transform :wavelet
:reason "Frequency changes over time → need time-frequency localization"
:recommended-type :daubechies-8}
;; Compression goal + smooth signal
(and (:smooth signal-properties) (:compress signal-properties))
{:transform :dct
:reason "Smooth signal + compression → DCT energy concentration"
:compression-ratio "10:1 to 100:1 typical"}
;; Frequency analysis with constant frequencies
(:constant-frequencies signal-properties)
{:transform :fft
:reason "Constant frequencies → DFT (via FFT algorithm) is fastest and most accurate"
:note "Use :real variant for real-valued signals (2x faster)"}
;; Denoising
(:noisy signal-properties)
{:transform :wavelet
:reason "Noise removal → wavelet thresholding"
:recommended-type :daubechies-8
:threshold-method :soft}
;; Boundary conditions (PDEs)
(:boundary-conditions signal-properties)
{:transform :dst
:reason "Signal must be zero at boundaries → use DST"
:application "Finite element methods, heat equation"}
;; Default
:else
{:transform :fft
:reason "General purpose - when in doubt, start with DFT (via FFT)"
:note "Analyze spectrum first, then choose specialized transform"})):::
Examples:
(recommend-transform #{:constant-frequencies}){:transform :fft,
:reason
"Constant frequencies → DFT (via FFT algorithm) is fastest and most accurate",
:note "Use :real variant for real-valued signals (2x faster)"}(recommend-transform #{:time-varying-frequency}){:transform :wavelet,
:reason
"Frequency changes over time → need time-frequency localization",
:recommended-type :daubechies-8}(recommend-transform #{:smooth :compress}){:transform :dct,
:reason "Smooth signal + compression → DCT energy concentration",
:compression-ratio "10:1 to 100:1 typical"}(recommend-transform #{:noisy}){:transform :wavelet,
:reason "Noise removal → wavelet thresholding",
:recommended-type :daubechies-8,
:threshold-method :soft}Performance Tips
| tip | reason | example | speedup |
|---|---|---|---|
| Use :real transforms | 2x faster than complex for real-valued signals | (t/transformer :real :fft) | |
| Reuse transformers | Pre-computes lookup tables in constructor | (def fft (t/transformer :real :fft)) | |
| Power-of-2 sizes | Enables fastest split-radix algorithm | ~2-3x vs non-power-of-2 | |
| Use dtype-next | Vectorized operations, avoids boxing | Often 10-100x vs boxed seq operations | |
| Avoid boxing | Type hints for primitive arrays | ^doubles or use dtype-next |
Common Pitfalls
| pitfall | solution | symptom |
|---|---|---|
| Forgetting to scale inverse FFT | Use scaled=true or divide by N manually | Results N× too large |
| Non-power-of-2 for wavelets | Pad signal to next power of 2 | Exception or wrong results |
| Using seqs for array math | Use dtype-next functional operations | Slow performance, high memory |
| Interpreting FFT magnitude wrong | Remember interleaved [real, imag] format | Incorrect frequency peaks |
| Not testing round-trip | Always verify transform → inverse recovers original | Lossy operations when expecting lossless |
Recipe Collection
Recipe 1: Find dominant frequency
(defn recipe-find-frequency [signal sample-rate]
(let [fft (t/transformer :real :fft)
spectrum (t/forward-1d fft signal)
mags (fft-magnitude spectrum)
max-idx (inc (argops/argmax (dtype/sub-buffer mags 1 (dec (dtype/ecount mags)))))
freq (* max-idx (/ sample-rate (dtype/ecount signal)))]
freq))Recipe 2: Remove frequency
(defn recipe-notch-filter [signal sample-rate target-freq bandwidth]
(let [fft (t/transformer :real :fft)
spectrum (t/forward-1d fft signal)
n (dtype/ecount signal)
;; Functional filtering: zero out target frequency band
filtered (dtype/emap
(fn [idx val]
(let [freq (* (quot idx 2) (/ sample-rate n))]
(if (< (Math/abs (- freq target-freq)) bandwidth)
0.0
val)))
:float64
(range (dtype/ecount spectrum))
spectrum)]
(t/reverse-1d fft filtered)))Recipe 3: Compress with quality target
(defn recipe-compress-to-quality [signal target-snr-db]
;; Compress signal until we reach a target quality level (measured in dB)
;; Higher target-snr-db = better quality = less compression
;; Example: 30 dB is good quality, 20 dB is acceptable, 10 dB is poor
(let [dct (t/transformer :real :dct)
coeffs (t/forward-1d dct signal)]
(loop [keep-ratio 0.9]
(let [compressed (compress-with-dct signal keep-ratio)
error (Math/sqrt (dfn/mean (dfn/sq (dfn/- signal compressed))))
signal-power (Math/sqrt (dfn/mean (dfn/sq signal)))
snr-db (* 20 (Math/log10 (/ signal-power error)))]
(if (or (>= snr-db target-snr-db) (< keep-ratio 0.05))
{:keep-ratio keep-ratio :snr-db snr-db :compressed compressed}
(recur (- keep-ratio 0.05)))))))Recipe 4: Denoise signal
(defn recipe-denoise [signal]
(let [;; Pad to power of 2
n (dtype/ecount signal)
padded-n (int (Math/pow 2 (Math/ceil (/ (Math/log n) (Math/log 2)))))
;; Functional padding
padded (dtype/emap
(fn [i] (if (< i n)
(dtype/get-value signal i)
0.0))
:float64
(range padded-n))
;; Wavelet denoise
wavelet (t/transformer :fast :daubechies-8)
coeffs (t/forward-1d wavelet padded)
denoised-coeffs (t/denoise wavelet coeffs :soft)
denoised-full (t/reverse-1d wavelet denoised-coeffs)]
;; Return only original length
(dtype/sub-buffer denoised-full 0 n)))Recipe 5: Spectral analysis with thresholding
(defn recipe-find-all-frequencies [signal sample-rate min-magnitude]
(let [fft (t/transformer :real :fft)
spectrum (t/forward-1d fft signal)
mags (fft-magnitude spectrum)
n (dtype/ecount signal)
peaks (filter #(> (dtype/get-value mags %) min-magnitude)
(range 1 (dtype/ecount mags)))]
(map (fn [idx]
{:frequency (* idx (/ sample-rate n))
:magnitude (dtype/get-value mags idx)})
peaks)))Summary and Next Steps
What We Learned
We used dtype-next for efficient array operations, leveraging vectorization and functional patterns for clean, performant code.
The core transforms each serve different purposes. The DFT (via FFT) provides frequency analysis using global sine/cosine basis functions. The DCT excels at compression with its cosine-only basis and energy concentration properties. Wavelets offer time-frequency localization with localized basis functions, perfect for denoising.
We also explored specialized transforms: the DST for boundary problems with sine-only basis, the DHT as a real-valued alternative for convolution, and 2D transforms for image processing and frequency analysis.
Applications included bandpass filtering, lossy compression with quality metrics, wavelet denoising, and spectral analysis. We validated everything through round-trip testing, stability testing across tiny and huge signals (1e-6 to 1e6), and property testing for characteristics like energy concentration and linearity.
Key Takeaways
| concept | insight |
|---|---|
| Linear Decomposition | All transforms express signals as weighted sums of basis functions |
| dtype-next First | Use vectorized operations for performance and clarity |
| Global vs Local | DFT/DCT see whole signal, wavelets see time AND frequency |
| Test-Driven | Always validate with round-trip tests |
| Choose Wisely | Match transform to signal properties and use case |
Next Steps
- Experiment: Try these transforms on your own signals
- Extend: Implement spectrograms (sliding window FFT)
- Optimize: Profile dtype-next vs traditional approaches
- Apply: Build audio visualizer, image compressor, or signal processor
- Deep Dive: Explore transform theory, basis functions, and mathematics
Further Resources
- dtype-next Guide: https://cnuernber.github.io/dtype-next/
- Fastmath Transform: https://generateme.github.io/fastmath/clay/transform.html
- JTransforms: https://github.com/wendykierp/JTransforms
- Clay Documentation: https://scicloj.github.io/clay/
End of Comprehensive Tutorial
You now have: theory, practice, validation, and recipes. Go forth and transform signals with confidence! 🌊📈