2. Costruzione di un segnale onda pressoria

Sommario

Costruzione di un segnale onda pressoria, corruzione dello stesso con rumore gaussiano, filtraggio ed analisi dei residui. Il nostro obiettivo è quello di costruire un segnale onda pressoria arteriosa tramite l’assegnazione di istanti temporali chiave e corrispondenti valori all’interno di un intervallo temporale definito e campionato in maniera regolare. Questa operazione si chiama interpolazione a tratti in italiano o “piecewise linear interpolation” in inglese.

Funzioni di base:

# FUNZIONE: Prodotto scalare
`%ps%` <- function(A, B) mean(A*B)

# FUNZIONE: Norma di un vettore
norm_ps <- function(A) sqrt(A %ps% A)


# FUNZIONE: Bodeplot
ggbodeplot <- function(tf, fmin=1, fmax=1e4, df=0.01) {
  # vector of points for each order of magnitude (OOM):
  pts <- 10^seq(0, 1, df) %>% tail(-1)
  # vector of OOMs:
  ooms <- 10^(floor(log10(fmin)):ceiling(log10(fmax)-1))
  # combine pts and ooms:
  freqs <- as.vector(pts %o% ooms)
  # warning: bode wants pulsation!
  bode(tf, freqs*2*pi) %>% {
    tibble(f=.$w/(2*pi), `magnitude (dB)`=.$mag, `phase (deg)`=.$phase)} %>%
    pivot_longer(-f) %>% 
    ggplot(aes(x=f, y=value)) +
    geom_line() +
    scale_x_log10(
      minor_breaks=scales::minor_breaks_n(10), 
      labels= ~ latex2exp::TeX(paste0("$10^{", log10(.), "}$"))) +
    facet_wrap(~name, nrow=2, scales="free") +
    labs(x="frequency (Hz)")
}

1 Sommario

Costruzione di un segnale onda pressoria, corruzione dello stesso con rumore gaussiano, filtraggio ed analisi dei residui. Il nostro obiettivo è quello di costruire un segnale onda pressoria arteriosa tramite l’assegnazione di istanti temporali chiave e corrispondenti valori all’interno di un intervallo temporale definito e campionato in maniera regolare. Questa operazione si chiama interpolazione a tratti in italiano o “piecewise linear interpolation” in inglese.

Opereremo come segue:

1. cerchiamo una figura di onda pressoria in rete per capire che istanti temporali chiave dovremmo poi assegnare
2. costruiamo una funzione in grado di eseguire l’interpolazione a tratti facendoci aiutare dall’intelligenza artificiale (AI)
3. approssimiamo la funzione lineare a tratti con N armoniche
4. aggiungiamo rumore gaussiano
5. filtriamo il segnale
6. eseguiamo l’analisi dei residui: sia statistica che in frequenza
7. in base al punto precedente cercheremo di definire una frequenza di taglio ottimale

2 Interpolazione a tratti

2.1 Risposta Chat Gpt

Domanda:

“se voglio iterare mediante purr su un vettore andando a usare degli indici di un altro vettore per fare dei conti. Ad esempio, ho un vettore”t” base dei tempi da 1 a 10. Poi un vettore tind = c(0, 3, 6, 10), un vettore tval = c(0, 1, 1, 0) ed applicare tramite programmazione funzionale una interpolazione dei valori tind e tval sulla base dei tempi?”

t     <- 1:10                      # base dei tempi
tind  <- c(0, 3, 6, 10)            # istanti "chiave"
tval  <- c(0, 1, 1, 0)             # valori ai tempi chiave

# Funzione per interpolare un singolo punto t_i
interp_fun <- function(ti, tind, tval) {
  # Trova l'intervallo [tind[k], tind[k+1]] contenente ti
  k <- max(which(tind <= ti & c(tind[-1], Inf) > ti))
  
  # Coordinate degli estremi
  t0 <- tind[k]
  t1 <- tind[k + 1]
  y0 <- tval[k]
  y1 <- tval[k + 1]
  
  # Interpolazione lineare
  y <- y0 + (ti - t0) * (y1 - y0) / (t1 - t0)
  return(y)
}

interp_fun(t[3], tind, tval)

[1] 1

Vedere anche Esercitazione scrittura Funzionale.pdf per una discussione su come interagire in maniera corretta con l’AI.

2.2 Elaborazione umana

Fare riferimento al documento Esercitazione scrittura Funzionale.pdf e alla nota soprastante.

Prove e visualizzazioni risultati parziali delle funzioni base:

t     <- 1:10                      # base dei tempi
tind  <- c(0, 3, 6, 10)            # istanti "chiave"
tval  <- c(0, 0, 3, 0)             # valori ai tempi chiave
ti <- t[4]

Dt <- ti - tind
Dt

[1]  4  1 -2 -6

Dt <- Dt[Dt>0]
Dt

[1] 4 1

ii <- length(Dt)
ii

[1] 2

t_p <- tind[ii]
t_p

[1] 3

val_p <- tval[ii]
val_p

[1] 0

t_d <- tind[ii+1]
t_d

[1] 6

val_d <- tval[ii+1]
val_d

[1] 3

Dt

[1] 4 1

Dt <- Dt[length(Dt)]
Dt

[1] 1

# Interpolazione lineare
val_interp <- val_p + Dt * (val_d - val_p) / (t_d - t_p)
val_interp

[1] 1

Sintesi della funzione di interpolazione:

t     <- 1:10                      # base dei tempi
tind  <- c(0, 3, 6, 10)            # istanti "chiave"
tval  <- c(0, 0, 3, 0)             # valori ai tempi chiave

interp_fun <- function(ti, tind, tval) {
  # "ti" istante da interpolare, "tind" istanti "chiave", "tval" valori ai tempi chiave
  Dt <- ti - tind
  Dt <- Dt[Dt>=0]
  ii <- length(Dt)

  Dt <- Dt[ii] # scarto con tempo chiave corrente

  t_p <- tind[ii] # valore temporale chiave precedente
  val_p <- tval[ii] # valore chiave precedente

  t_d <- tind[ii+1]  # valore temporale chiave successivo
  val_d <- tval[ii+1]  # valore chiave successivo

  # Interpolazione lineare e restituzione valore in uscita della funzione
  val_p + Dt * (val_d - val_p) / (t_d - t_p)
}

interp_fun(1, tind, tval)

[1] 0

2.3 Uso della funzione per costruire un’onda di pressione arteriosa

1. troviamo la figura di una onda pressoria in rete
2. la interpoliamo linearmente con dei valori chiave usando la programmazione funzionale tramite map_dbl()
3. ricaviamo ed usiamo le prime N armoniche che la approssimano

2.3.1 Punti 1 e 2 - interpolazione onda pressoria tipica

Troviamo la figura di una onda pressoria in rete e la interpoliamo linearmente con i valori chiave che ricaviamo dagli andamenti tipici trovati in rete:

t <- seq(from = 0, to = 1, by = 0.001)  # base dei tempi
tind  <- c(0, 0.1, 0.2, 0.3, 0.5, 1)   # istanti "chiave"
tval  <- c(80, 120, 100, 105, 90, 80)  # valori ai tempi chiave
# interp_fun(0, tind, tval)

p_art <- map_dbl(t[1:length(t)-1], \(t) interp_fun(t, tind, tval))
# NOTA: mappiamo su t[1:length(t)-1] altrimenti per l'ultimo valore non riusciremmo a trovare l'intervallo lineare corrispondente

t <- t[1:length(t)-1]


# Simuliamo 10 battiti cardiaci
p_art <- rep(p_art, 10)
t <- seq(from = 0, to = 0.001*(length(p_art)-1), by = 0.001)

pp <- plot_ly() %>%
  add_lines(t, p_art, name = "Onda pressoria a tratti", line= list(color= "red")) %>%
  layout(title = "Onda pressoria")
pp

2.3.2 Punto 3 - usiamo le prime N armoniche

sin_k <- function(t, f, k=1) sin(2*pi*f*k*t)
cos_k <- function(t, f, k=1) cos(2*pi*f*k*t)


K <- 120
f_fond <- 1 / t[length(t)-1]
comps_sin <- map_dbl(1:K, \(k) {
  s <- sin_k(t, f_fond, k)
  (p_art - mean(p_art)) %ps% s/norm_ps(s)
}) %>% zapsmall()

comps_sin

  [1] -0.000442 -0.000907 -0.001421 -0.002024 -0.002779 -0.003812 -0.005404
  [8] -0.008402 -0.017013  8.374429  0.016008  0.007371  0.004247  0.002437
 [15]  0.001061 -0.000251 -0.001830 -0.004366 -0.011036  3.115668  0.013123
 [22]  0.006535  0.004026  0.002472  0.001194 -0.000127 -0.001840 -0.004762
 [29] -0.012787  2.546676  0.017155  0.009288  0.006424  0.004769  0.003516
 [36]  0.002324  0.000876 -0.001500 -0.007951  1.535344  0.016188  0.010008
 [43]  0.007879  0.006788  0.006128  0.005699  0.005424  0.005291  0.005430
 [50] -0.055591  0.004052  0.004338  0.004484  0.004652  0.004901  0.005305
 [57]  0.006017  0.007492  0.012017 -0.739623 -0.005785 -0.001364  0.000182
 [64]  0.001040  0.001672  0.002274  0.003022  0.004289  0.007817 -0.477130
 [71] -0.005404 -0.002046 -0.000859 -0.000208  0.000248  0.000645  0.001085
 [78]  0.001756  0.003506 -0.200692 -0.002770 -0.001133 -0.000538 -0.000203
 [85]  0.000038  0.000252  0.000495  0.000873  0.001875 -0.103055 -0.001772
 [92] -0.000848 -0.000529 -0.000366 -0.000269 -0.000206 -0.000165 -0.000138
 [99] -0.000123 -0.000118 -0.000121 -0.000133 -0.000156 -0.000192 -0.000245
[106] -0.000328 -0.000466 -0.000737 -0.001541  0.069099  0.001451  0.000680
[113]  0.000385  0.000201  0.000044 -0.000124 -0.000349 -0.000740 -0.001827
[120]  0.083465

# Con rowSums():
p_art_K = rowSums(
  sapply(1:K, \(k) {
    s_sin <- sin_k(t, f_fond, k)
    comps_sin[k] * s_sin / norm_ps(s_sin)
  }))

# Con reduce():
# da completare


pp <- pp %>%
  add_lines(t, p_art_K + mean(p_art), name = "Onda pressoria approssimata", line= list(color= "magenta"))
pp

2.4 Aggiungiamo rumore e filtriamo per ridurne l’effetto

2.4.1 Aggiungiamo rumore

Abbiamo la base dei tempi t, e l’approssimazione dell’onda pressoria p_art_K. Aggiungiamo rumore bianco con distribuzione normale:

# Aggiunta di rumore
p_art_K_noise <- p_art_K + rnorm(length(t), 0, 3)


pp <- pp %>%
  add_lines(t, p_art_K_noise + mean(p_art), name = "Onda pressoria \n approssimata con \n rumore", line= list(color= "black"))
pp

2.4.2 Filtriamo il rumore e valutiamone l’effetto

Costruzione ed applicazione del filtro:

fc <- 1/(t[2]-t[1]) # frequenza di campionamento
fn <- fc/2 # frequenza di Nyquist
ft <- 50

# Filtro Butterworth di ordine 'ord' nel dominio `z` discreto
# Il parametro w è normalizzato tra 0 e 1, con 1 = fn
ord <- 3
flt <- butter(ord, W=ft/fn, plane="z")
flt

$b
[1] 0.002898195 0.008694584 0.008694584 0.002898195

$a
[1]  1.0000000 -2.3740947  1.9293557 -0.5320754

attr(,"class")
[1] "Arma"

# Applicazione del filtro
p_art_K_noise_filt <- signal::filter(flt, p_art_K_noise)

pp %>%
  add_lines(t, p_art_K_noise_filt + mean(p_art), name = "Onda pressoria \n approssimata con \n rumore e filtrata", line= list(color= "green"))

Come si vede, l’oggetto filtro contiene i due vettori con i coefficienti polinomiali di un modello ARMA.

Plot della Funzione di trasferimento del filtro (nel dominio continuo s):

filtro_tf <- butter(ord, W=ft/fn, plane="s") %>% {tf(.$b, .$a)}
ggbodeplot(filtro_tf, fmin=0.01, fmax=10)

2.5 Analisi dei residui

Inserire qui la parte richiesta.

Esercitazione

Riconoscimento punti esame

Il completamento di questo esercizio comporta il riconoscimento di -2…+3 punti in base alla valutazione del docente (3 = eccellente, 2 = buono, 1 = mediocre o svolta parzialmente, …, -2 = non conoscenza di quanto riportato).

Assegnazioni

Scrivete in modo diverso la funzione di interpolazione [NOTA: scrivere in modo diverso la funzione può voler dire qualsiasi cosa. Se userete un approccio completamente diverso questo punto varrà parecchio, in caso modifichiate solo qualche assegnazione varrà meno]

Riscrivete con reduce() la routine seguente:

p_art_K <- rowSums(
  sapply(1:K, \(k) {
    s_sin <- sin_k(t, f_fond, k)
    comps_sin[k] * s_sin / norm_ps(s_sin)
  }))

Usare anche i coseni per ricostruire il segnale.

Trovate una frequenza di taglio opportuna.

A valle del filtraggio il segnale risulta ritardato. Fate in modo di eliminare tale ritardo.

Effettuare l’analisi statistica dei residui (residui := p_art_K - p_art_K_noise_filt).

Effettuare l’analisi in frequenza dei residui ed estrapolare le considerazioni del caso.

Scegliere una frequenza di taglio a valle dell’analisi dei residui [NOTA: al variare della frequenza di taglio vedrete che i residui cambiano. Come indicatori di tale cambiamento potrete usare: valor medio, deviazione standard, analisi della normalità, banda in frequenza, etc. Quello che si chiede è di dimostrare di aver compreso la relazione tra ft e tali indicatori e di saper indicare anche in modo non quantitativo con che criterio si può effettuare la scelta]

Suggerimenti

• Per la frequenza di taglio finale notare il fatto che il rumore aggiunto (e poi da ridurre tramite filtro) è gaussiano. Ripensate anche all’approssimazione tramite Fourier eseguita.
• Per eliminare il ritardo potete usare funzioni diverse da filter (o anche, meglio, filter stessa con una serie di comandi)