Skip to content

Commit 5cd5329

Browse files
dasludaslu
committed
mths - wip
1 parent 9eb7518 commit 5cd5329

1 file changed

Lines changed: 127 additions & 118 deletions

File tree

‎src/dsp/mths.clj‎

Lines changed: 127 additions & 118 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,8 @@
77
:date "2025-11-21"
88
:tags [:dsp :signal-processing :biometrics :camera-ppg :heart-rate]
99
:draft true}}}
10-
(ns dsp.mths)
10+
(ns dsp.mths
11+
(:require [scicloj.kindly.v4.kind :as kind]))
1112

1213
;; ## Introduction
1314
;;
@@ -308,123 +309,131 @@
308309
;;
309310
;; Let's walk through the complete pipeline for one subject:
310311

311-
;; (let [;; Frequency range for heart rate (0.65-4 Hz = 39-240 bpm)
312-
;; [low-cutoff high-cutoff] [0.65 4]
313-
314-
;; ;; Window parameters for spectral analysis
315-
;; window-size 10 ; seconds
316-
;; window-samples (* sampling-rate window-size) ; 300 samples at 30 Hz
317-
318-
;; ;; Create Hanning window to reduce spectral leakage
319-
;; hanning (-> window-samples
320-
;; Hanning.
321-
;; .getWindow
322-
;; dtype/as-reader)
323-
324-
;; ;; Overlap parameters
325-
;; overlap-fraction 0.05 ; 5% overlap between windows
326-
;; hop (int (* overlap-fraction window-samples))
327-
328-
;; ;; Select subject to analyze
329-
;; subject 51
330-
;; sig (signal subject)
331-
;; n-samples (tc/row-count sig)
332-
333-
;; ;; Apply preprocessing pipeline: DC removal → bandpass → standardization
334-
;; std (-> sig
335-
;; (tc/update-columns [:R :G :B]
336-
;; #(-> %
337-
;; remove-dc
338-
;; (bandpass-filter {:fs sampling-rate
339-
;; :order 4
340-
;; :low-cutoff low-cutoff
341-
;; :high-cutoff high-cutoff})
342-
;; stats/standardize))
343-
;; (tc/select-columns [:R :G :B])
344-
;; ds-tensor/dataset->tensor)
345-
346-
;; ;; Calculate how many windows we can extract
347-
;; n-windows (-> n-samples
348-
;; (- window-samples)
349-
;; (quot hop))
350-
351-
;; ;; Create sliding windows: [time × windows × channels]
352-
;; windows-shape [window-samples
353-
;; n-windows
354-
;; 3]
355-
;; windows (tensor/compute-tensor
356-
;; windows-shape
357-
;; (fn [i j k]
358-
;; ;; Extract sample at time i, window j, channel k
359-
;; (std (+ (* j hop) i)
360-
;; k))
361-
;; :float32)
362-
363-
;; ;; Apply Hanning window to each segment
364-
;; hanninged-windows (tensor/compute-tensor
365-
;; windows-shape
366-
;; (fn [i j k]
367-
;; (* (windows i j k)
368-
;; (hanning i)))
369-
;; :float32)
370-
371-
;; ;; Compute power spectrum for each window
372-
;; power-spectrum (-> hanninged-windows
373-
;; ;; Rearrange to [windows × channels × time]
374-
;; (tensor/transpose [1 2 0])
375-
;; (tensor/slice 2)
376-
;; ;; Apply FFT to each window
377-
;; (->> (mapv (fn [window]
378-
;; (let [fft (-> window
379-
;; double-array
380-
;; DiscreteFourier.)]
381-
;; (.transform fft)
382-
;; ;; Get magnitude spectrum (power)
383-
;; (.getMagnitude fft true)))))
384-
;; tensor/->tensor
385-
;; ;; Reshape to [windows × channels × frequencies]
386-
;; (#(tensor/reshape %
387-
;; [n-windows
388-
;; 3
389-
;; (-> % first count)]))
390-
;; ;; Transpose to [frequencies × windows × channels] for plotting
391-
;; (tensor/transpose [2 0 1]))]
392-
;; ;; ### Visualizations
393-
;; ;;
394-
;; ;; We'll create several visualizations to understand the signal processing pipeline:
395-
396-
;; [(-> windows
397-
;; (dfn/* 100)
398-
;; plotly/imshow
399-
;; (kind/hiccup
400-
;; [:div
401-
;; [:h4 "1. Raw Windowed Signal (after preprocessing)"]
402-
;; [:p "Each row is a time sample, showing how the signal varies across windows (columns) and channels."]]))
403-
404-
;; (-> hanninged-windows
405-
;; (dfn/* 100)
406-
;; plotly/imshow
407-
;; (kind/hiccup
408-
;; [:div
409-
;; [:h4 "2. After Hanning Window"]
410-
;; [:p "Notice how the edges of each window are tapered to zero. This reduces spectral leakage in the FFT."]]))
411-
412-
;; (-> power-spectrum
413-
;; plotly/imshow
414-
;; (kind/hiccup
415-
;; [:div
416-
;; [:h4 "3. Power Spectrum"]
417-
;; [:p "Each row is a frequency bin. Bright bands indicate strong frequency components. "
418-
;; "The dominant frequency (brightest band in the 0.65-4 Hz range) corresponds to the heart rate."]]))
419-
420-
;; (-> power-spectrum
421-
;; (tensor/slice 1)
422-
;; (->> (map (fn [s]
423-
;; (-> s
424-
;; ds-tensor/tensor->dataset
425-
;; (tc/rename-columns [:R :G :B])
426-
;; plot-signal))))
427-
;; (cons (kind/md "### 4. Power Spectrum Time Series\n\nEach plot shows how the power spectrum evolves over time for one channel.")))])
312+
(let [;; Frequency range for heart rate (0.65-4 Hz = 39-240 bpm)
313+
[low-cutoff high-cutoff] [0.65 4]
314+
315+
;; Window parameters for spectral analysis
316+
window-size 10 ; seconds
317+
window-samples (* sampling-rate window-size) ; 300 samples at 30 Hz
318+
319+
;; Create Hanning window to reduce spectral leakage
320+
hanning (-> window-samples
321+
Hanning.
322+
.getWindow
323+
dtype/as-reader)
324+
325+
;; Overlap parameters
326+
overlap-fraction 0.05 ; 5% overlap between windows
327+
hop (int (* overlap-fraction window-samples))
328+
329+
;; Select subject to analyze
330+
subject 51
331+
sig (signal subject)
332+
n-samples (tc/row-count sig)
333+
334+
;; Apply preprocessing pipeline: DC removal → bandpass → standardization
335+
std (-> sig
336+
(tc/update-columns [:R :G :B]
337+
#(-> %
338+
remove-dc
339+
(bandpass-filter {:fs sampling-rate
340+
:order 4
341+
:low-cutoff low-cutoff
342+
:high-cutoff high-cutoff})
343+
stats/standardize))
344+
(tc/select-columns [:R :G :B])
345+
ds-tensor/dataset->tensor)
346+
347+
;; Calculate how many windows we can extract
348+
n-windows (-> n-samples
349+
(- window-samples)
350+
(quot hop))
351+
352+
;; Create sliding windows: [time × windows × channels]
353+
windows-shape [window-samples
354+
n-windows
355+
3]
356+
windows (tensor/compute-tensor
357+
windows-shape
358+
(fn [i j k]
359+
;; Extract sample at time i, window j, channel k
360+
(std (+ (* j hop) i)
361+
k))
362+
:float32)
363+
364+
;; Apply Hanning window to each segment
365+
hanninged-windows (tensor/compute-tensor
366+
windows-shape
367+
(fn [i j k]
368+
(* (windows i j k)
369+
(hanning i)))
370+
:float32)
371+
372+
;; Compute power spectrum for each window
373+
power-spectrum (-> hanninged-windows
374+
;; Rearrange to [windows × channels × time]
375+
(tensor/transpose [1 2 0])
376+
(tensor/slice 2)
377+
;; Apply FFT to each window
378+
(->> (mapv (fn [window]
379+
(let [fft (-> window
380+
double-array
381+
DiscreteFourier.)]
382+
(.transform fft)
383+
;; Get magnitude spectrum (power)
384+
(.getMagnitude fft true)))))
385+
tensor/->tensor
386+
;; Reshape to [windows × channels × frequencies]
387+
(#(tensor/reshape %
388+
[n-windows
389+
3
390+
(-> % first count)]))
391+
;; Transpose to [frequencies × windows × channels] for plotting
392+
(tensor/transpose [2 0 1]))]
393+
;; ### Visualizations
394+
;;
395+
;; We'll create several visualizations to understand the signal processing pipeline:
396+
397+
(kind/fragment
398+
[(kind/hiccup
399+
[:div
400+
[:h4 "1. Raw Windowed Signal (after preprocessing)"]
401+
[:p "Each row is a time sample, showing how the signal varies across windows (columns) and channels."]])
402+
403+
(-> windows
404+
(dfn/* 100)
405+
plotly/imshow)
406+
407+
(kind/hiccup
408+
[:div
409+
[:h4 "2. After Hanning Window"]
410+
[:p "Notice how the edges of each window are tapered to zero. This reduces spectral leakage in the FFT."]])
411+
412+
(-> hanninged-windows
413+
(dfn/* 100)
414+
plotly/imshow)
415+
416+
(kind/hiccup
417+
[:div
418+
[:h4 "3. Power Spectrum"]
419+
[:p "Each row is a frequency bin. Bright bands indicate strong frequency components. "
420+
"The dominant frequency (brightest band in the 0.65-4 Hz range) corresponds to the heart rate."]])
421+
422+
(-> power-spectrum
423+
plotly/imshow)
424+
425+
(kind/md "### 4. Power Spectrum Time Series\n\nEach plot shows how the power spectrum evolves over time.")
426+
427+
(-> power-spectrum
428+
(tensor/slice 1)
429+
(->> (map (fn [s]
430+
(-> s
431+
ds-tensor/tensor->dataset
432+
(tc/rename-columns [:R :G :B])
433+
plot-signal)))
434+
(into [:div {:style {:max-height "600px"
435+
:overflow-y "auto"}}])
436+
kind/hiccup))]))
428437

429438

430439

0 commit comments

Comments
 (0)