Skip to content

Commit 26fcfd4

Browse files
dasludaslu
authored
Merge pull request #100 from ClojureCivitas/hrv-for-macroexpand-noj-2
Hrv for macroexpand noj 2
2 parents d37a1a6 + d3442cf commit 26fcfd4

File tree

2 files changed

+240
-391
lines changed

2 files changed

+240
-391
lines changed

‎site/data_analysis/heart_rate_variability/exploring_heart_rate_variability.qmd‎

Lines changed: 159 additions & 328 deletions
Large diffs are not rendered by default.

‎src/data_analysis/heart_rate_variability/exploring_heart_rate_variability.clj‎

Lines changed: 81 additions & 63 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,11 @@
11
^{:kindly/hide-code true
2-
:clay {:title "Exploring Heart Rate Variability"
3-
:external-requirements ["WESAD dataset at /workspace/datasets/WESAD/"]
4-
:quarto {:author :ludgersolbach
5-
:draft true
6-
:type :post
7-
:date "2025-10-17"
8-
:tags [:data-analysis :noj]}}}
2+
:clay {:title "Exploring Heart Rate Variability"
3+
:external-requirements ["WESAD dataset at /workspace/datasets/WESAD/"]
4+
:quarto {:author :ludgersolbach
5+
:draft true
6+
:type :post
7+
:date "2025-10-17"
8+
:tags [:data-analysis :noj]}}}
99
(ns data-analysis.heart-rate-variability.exploring-heart-rate-variability
1010
(:require [tablecloth.api :as tc]
1111
[scicloj.tableplot.v1.plotly :as plotly]
@@ -27,7 +27,6 @@
2727
[scicloj.tableplot.v1.plotly :as plotly]
2828
[java-time.api :as jt]))
2929

30-
3130
;; # Exploring HRV - DRAFT 🛠
3231

3332
(ns data-analysis.heart-rate-variability.exploring-heart-rate-variability
@@ -50,15 +49,13 @@
5049
(:import [com.github.psambit9791.jdsp.signal CrossCorrelation]
5150
[com.github.psambit9791.jdsp.signal.peaks FindPeak]
5251
[com.github.psambit9791.jdsp.filter Butterworth]
53-
[com.github.psambit9791.jdsp.transform DiscreteFourier]))
54-
52+
[com.github.psambit9791.jdsp.transform DiscreteFourier]
53+
[com.github.psambit9791.jdsp.windows Hanning]))
5554

5655
;; ## My pulse-to-pulse intervals
5756

58-
5957
;; (extracted from PPG data)
6058

61-
6259
(def my-ppi
6360
(-> (tc/dataset "src/data_analysis/heart_rate_variability/ppi-series.csv"
6461
{:key-fn keyword})
@@ -75,10 +72,9 @@
7572
:=height 300 :=width 700})
7673
(plotly/layer-bar {:=y :ppi})))
7774

78-
7975
(def compute-measures
8076
(fn [ppi-ds {:keys [sampling-rate
81-
window-size-in-sec ]}]
77+
window-size-in-sec]}]
8278
(let [spline (interp/interpolation
8379
:cubic
8480
(:t ppi-ds)
@@ -92,8 +88,8 @@
9288
bw (com.github.psambit9791.jdsp.filter.Butterworth.
9389
sampling-rate)
9490
n (tc/row-count resampled-ppi)
95-
window-size (* window-size-in-sec sampling-rate)
96-
hop-size 8
91+
window-size (* window-size-in-sec sampling-rate)
92+
hop-size (* 5 sampling-rate) ; 5-second hops for reasonable overlap
9793
n-windows (int (/ (- n window-size)
9894
hop-size))
9995
ranges (-> ppi-ds
@@ -126,60 +122,70 @@
126122
window (-> resampled-ppi
127123
:ppi
128124
(dtype/sub-buffer start-idx window-size))
129-
window-standardized (stats/standardize window)
125+
;; Filter first, then standardize for FFT normalization
130126
window-filtered (.bandPassFilter bw
131-
(double-array window-standardized)
127+
(double-array window)
132128
4
133-
0
129+
0.04 ; Lower cutoff at 0.04 Hz to remove baseline drift
134130
0.4)
135-
fft (DiscreteFourier. (double-array window-filtered))
131+
window-standardized (stats/standardize window-filtered)
132+
;; Apply Hann window to reduce spectral leakage
133+
hann-window (-> (Hanning. window-size) .getWindow)
134+
window-windowed (map * window-standardized hann-window)
135+
fft (DiscreteFourier. (double-array window-windowed))
136136
_ (.transform fft)
137137
whole-magnitude (.getMagnitude fft true)]
138138
{:t (* w hop-size (/ 1.0 sampling-rate))
139139
:whole-magnitude whole-magnitude
140140
:magnitude (take 40 whole-magnitude)})))
141141
vec)]
142142
{:sampling-rate sampling-rate
143+
:window-size window-size ;; Added: FFT window size needed for freq resolution
143144
:resampled-ppi resampled-ppi
144145
:rmssds rmssds
145146
:spectrograms spectrograms})))
146147

147-
148148
(comment
149149
(compute-measures my-ppi
150150
{:sampling-rate 10
151151
:window-size-in-sec 60}))
152152

153-
154153
;; [An Overview of Heart Rate Variability Metrics and Norms](https://pmc.ncbi.nlm.nih.gov/articles/PMC5624990/)
155154
;; by Fred Shaffer, & J P Ginsberg.
156155
;; doi: [10.3389/fpubh.2017.00258](https://www.frontiersin.org/journals/public-health/articles/10.3389/fpubh.2017.00258/full)
157156

158157
(defn LF-to-HF [freqs spectrogram]
159158
(let [ds (tc/dataset {:f freqs
160159
:s (:magnitude spectrogram)}
161-
tc/dataset)]
162-
(/ (-> ds
163-
(tc/select-rows #(<= 0.04 (% :f) 0.15))
164-
:s
165-
tcc/sum)
166-
(-> ds
167-
(tc/select-rows #(<= 0.15 (% :f) 0.4))
168-
:s
169-
tcc/sum))))
170-
160+
tc/dataset)
161+
lf-power (-> ds
162+
(tc/select-rows #(<= 0.04 (% :f) 0.15))
163+
:s
164+
tcc/sum)
165+
hf-power (-> ds
166+
(tc/select-rows #(<= 0.15 (% :f) 0.4))
167+
:s
168+
tcc/sum)]
169+
(if (pos? hf-power)
170+
(/ lf-power hf-power)
171+
Double/NaN))) ; Return NaN when HF power is zero or negative
171172

172173
(defn plot-with-measures [{:keys [sampling-rate
174+
window-size
173175
resampled-ppi
174176
rmssds
175177
spectrograms]}]
176178
(when spectrograms
177-
(let [n (-> spectrograms first :magnitude count)
179+
(let [;; Number of frequency bins we're displaying (truncated spectrum)
180+
n (-> spectrograms first :magnitude count)
181+
;; Correct frequency resolution based on FFT window size
182+
;; freq_resolution = sampling_rate / FFT_size
183+
freq-resolution (/ sampling-rate window-size)
184+
;; Nyquist frequency (maximum representable frequency)
178185
Nyquist-freq (/ sampling-rate 2.0)
179-
freq-resolution (/ Nyquist-freq n)
180186
times (map (comp str :t) spectrograms)
181-
freqs (tcc/* (range n)
182-
freq-resolution)]
187+
;; Generate frequency values for the n bins we're displaying
188+
freqs (tcc/* (range n) freq-resolution)]
183189
{:resampled-ppi (-> resampled-ppi
184190
(plotly/base {:=height 300 :=width 700})
185191
(plotly/layer-bar (merge {:=x :t
@@ -206,7 +212,7 @@
206212
:width 700
207213
:margin {:t 25}
208214
:xaxis {:title {:text "t"}}
209-
:yaxis {:title {:text "freq"}}}})
215+
:yaxis {:title {:text "freq (Hz)"}}}})
210216
:mean-power-spectrum (-> {:freq freqs
211217
:mean-power (-> spectrograms
212218
(->> (map :magnitude))
@@ -227,14 +233,12 @@
227233
(assoc-in [:layout :yaxis :range] [0 4])
228234
(assoc-in [:layout :yaxis :title] {:text "LF/HF"}))})))
229235

230-
231236
(delay
232237
(-> my-ppi
233238
(compute-measures {:sampling-rate 10
234239
:window-size-in-sec 60})
235240
plot-with-measures))
236241

237-
238242
;; ## Analysing ECG data
239243

240244
;; ### The [WESAD](https://dl.acm.org/doi/10.1145/3242969.3242985) dataset
@@ -265,8 +269,8 @@
265269
:ECG (-> ld
266270
(get-in [:signal :chest :ECG])
267271
(py. flatten))
268-
:label (-> ld
269-
(get :label))})))))
272+
:label (-> ld
273+
(get :label))})))))
270274

271275
(delay
272276
(labelled-dataset 5))
@@ -275,7 +279,7 @@
275279

276280
;; [Pan-Tompkins Algorithm](https://en.wikipedia.org/wiki/Pan%E2%80%93Tompkins_algorithm)
277281

278-
;; [Unupervised ECG QRS Detection](https://hooman650.github.io/ECG-QRS.html) by Hooman SedghamizHooman Sedghamiz
282+
;; [Unupervised ECG QRS Detection](https://hooman650.github.io/ECG-QRS.html) by Hooman Sedghamiz
279283

280284
;; scipy: `peaks = signal.find_peaks(signal, height=mean, distance=200)`
281285
;; JDSP equivalent:
@@ -290,7 +294,6 @@
290294
final-peaks (.filterByPeakDistance peak-obj height-filtered distance)]
291295
final-peaks)) ; Returns int[] of peak row-numbers
292296

293-
294297
(delay
295298
(let [bw (com.github.psambit9791.jdsp.filter.Butterworth.
296299
WESAD-sampling-rate)
@@ -328,7 +331,6 @@
328331
(plotly/layer-point {:=y :peak
329332
:=name "peak"}))))
330333

331-
332334
(defn extract-ppi
333335
"Extract peak-to-peak intervals from ECG signal.
334336
Returns dataset with columns: :t (time in seconds), :ppi (interval in seconds)"
@@ -345,8 +347,9 @@
345347
4 5 15))
346348
;; Differentiate and square
347349
(tc/add-column :sqdiff
348-
#(tcc/- (:filtered %)
349-
(tcc/shift (:filtered %) 1))))
350+
#(tcc/sq
351+
(tcc/- (:filtered %)
352+
(tcc/shift (:filtered %) 1)))))
350353
;; Find peaks with distance constraint (200 samples = ~0.29s)
351354
peak-indices (find-peaks (:sqdiff pipeline)
352355
{:distance 200})
@@ -381,7 +384,6 @@
381384
extract-ppi
382385
(compute-measures measures-params)))))
383386

384-
385387
(delay
386388
(-> {:ppi-params {:subject-id 5
387389
:row-interval [0 1000000]}
@@ -413,24 +415,40 @@
413415
(tc/add-column :offset #(cons 0 (reductions + (:n %))))
414416
(tc/select-columns [:offset :n :label])))))
415417

416-
417418
(delay
418419
(label-intervals 5))
419420

420-
421421
(delay
422422
(let [subject 5]
423-
(kind/fragment
424-
(-> (label-intervals subject)
425-
(tc/select-rows #(not= (:label %) :ignore))
426-
#_(tc/select-rows #(= (:label %) :meditation))
427-
(tc/rows :as-maps)
428-
(->> (mapcat (fn [{:keys [offset n label]}]
429-
[label
430-
(try (-> {:ppi-params {:subject-id subject
431-
:row-interval [offset (+ offset n)]}
432-
:measures-params {:sampling-rate 10
433-
:window-size-in-sec 120}}
434-
WESAD-spectrograms
435-
plot-with-measures)
436-
(catch Exception e 'unavailable))])))))))
423+
(-> (label-intervals subject)
424+
(tc/select-rows #(not= (:label %) :ignore))
425+
#_(tc/select-rows #(= (:label %) :meditation))
426+
(tc/rows :as-maps)
427+
(->> (mapcat (fn [{:keys [offset n label]}]
428+
[label
429+
(try (-> {:ppi-params {:subject-id subject
430+
:row-interval [offset (+ offset n)]}
431+
:measures-params {:sampling-rate 10
432+
:window-size-in-sec 120}}
433+
WESAD-spectrograms
434+
plot-with-measures)
435+
(catch Exception e 'unavailable))]))))))
436+
437+
;; ## Conclusion
438+
439+
;; - measures are tricky
440+
;; - domain knowledge matters
441+
;; - workflow matters
442+
;; - visualization matters
443+
444+
445+
446+
447+
448+
449+
450+
451+
452+
453+
454+

0 commit comments

Comments
 (0)