# FUNZIONE: Generazione segnale sinusoidale con diverse armoniche
signal <- function(t, pars, rad = FALSE) {
stopifnot(is.data.frame(pars))
with(pars, {
if (!rad) {
phi <- phi/180*pi
f <- 2*pi*f
}
map_dbl(t, \(t) sum( map_vec(seq_along(w) , \(i) w[i]*cos(t*f[i] + phi[i] ))))
})
}Esercitazione 5
Inversione delle armoniche pure
Per la consegna dell’esercitazione si faccia riferimento a questo link. Mentre le librerie sono rimandate nell’icona Librerie.
1 Funzioni
2 Segnale d’ingresso
f0 <- 10
# Provare con un'armonica pura
# pars <- tibble(
# w = c(1),
# f = c(f0), # 1 funziona
# phi = c(0)
# )
# Provare con tante armoniche
pars <- tibble(
w = c(1, 0.1, 0.3),
f = c(15, 20, 45),
phi = c(0, 0, 0)
)
N <- 1000 # dimensione del campione
fs <- 1000 # freq. di campionamento
s <- tibble(
t = 0:N / fs, # 1 kHz di frequenza di campionamento
u = signal(t, pars),
u_n = u + rnorm(length(t), 0, pars$w[1] / 10)
)
3 FdT del sistema di isolamento
M <- 10
K <- 1000
C <- 50
# Frequenza naturale:
fn <- 1/(2*pi) * sqrt(K/M)
num <- c(C, K)
den <- c(M, C, K)
H <- tf(num, den)
ggbodeplot_continous(H, fs, fmin=0.1, fmax=100) +
geom_vline(xintercept=c(1, sqrt(2)) * fn, color="red", linetype=2) +
labs(title=paste(
"Natural frequency:", round(fn, 2), "Hz",
" - Isolation: >", round(sqrt(2)*fn, 2), "Hz"))
4 Simulazione dell’uscita
output <- lsim(H, s$u_n, s$t)
s <- s %>%
mutate(
# Simulazione dell'uscita
y = output$y[1,]
)
# Grafici
plot_ly() %>%
add_lines(x = s$t, y = s$u_n, name = "Input") %>%
add_lines(x = s$t, y = s$y, name = "Output") %>%
layout(xaxis = list(title = "Tempo (s)"), yaxis = list(title = ""), title = "Input & Output")5 Ricostruzione dell’ingresso
Stima dell’ingresso a partire da uscita del sistema e da modulo e fase della FdT.
R è in grado di operare nel domino complesso, perciò viene effettuata la compensazione dell’ingresso nel dominio della frequenza, come in Esercitazione 4. Nello specifico, mediante la funzione di libreria dyn_compensation vengono essenzialmente svolti i segueni passaggi, presupponendo un’opportuno segnale d’uscita passatogli come parametro:
eventuale padding, impostato dall’utente, per garantire la completa periodicità del segnale;
calcolo dalla FFT del segnale d’uscita e valutazione numerica della risposta in frequenza della FdT;
inversione della caratteristica dinamica;
FFT inversa che calcola l’ingresso nel dominio del tempo.
Per effettuare la compensazione, è bene prendere atto che l’FdT inversa è un passa-alto:
[1] "TFCHK: Transfer function may not be proper and may lead to errors. Num > Den"
Sulla base dei diagrammi di bode ed osservando lo spettro dell’uscita che segue, è evidente che le componenti oltre i \(2.25\,\mathrm{Hz}\) verranno amplificate.
A questo punto viene effettuata la compensazione:
# Compensazione dell'uscita
u <- s %>% select(t, u, u_n) %>%
mutate(
u_comp = dyn_compensation(y, H, name_sig = "y")$u
)Come anticipato l’inversione della FdT comporta la ricostruzione di un segnale abbastanza rumoroso; perdipiù la coda destra presenta delle forti amplificazioni - derivanti presumibilmente dal taglio delle code durante il padding -. Per comprendere il tipo di rumore e quali componenti lo costituiscano si osserva lo spettro del segnale.
Il segnale compensato presenta un fondo spettrale diffuso e privo di struttura armonica evidente, compatibile con un rumore di tipo bianco. Pertanto, ai fini di eliminare quest’ultimo e ricostruire il segnale utile come somma delle sole componenti armoniche dominanti, viene realizzata una maschera in frequenza che annulla le componenti associate al fondo spettrale.
N <- nrow(u)
k <- 0:(N-1)
# La soglia per selezionare le componenti predominanti
# viene individuata visivamente dallo spettro, ricordando
# che lo spettro completo (bilatero), possiede metà dell'energia
# spettrale di quello monolarero
threashold <- 0.04
u_spectrum <- tibble(
f = ifelse(k <= N/2, k, k-N) * fs/N,
U = fft(u$u_comp),
mask = Mod(U)/N > threashold, # vettore di booleani
U_filtered = U*mask
)
u_spectrum %>% dplyr::filter(f>=0) %>%
mutate(mod = 2*Mod(U_filtered)/N) %>%
ggplot(aes(f, mod)) +
geom_spoke(aes(y = 0, radius = mod, angle = pi/2), color = "#F8766D") +
geom_point(color = "#F8766D") +
labs(x = "Frequenza (Hz)", y = "Intensità", title = "Spettro del input mascherato")
Nota: la componente DC viene processata dalla maschera perché visivamente si vede (e si sa a priori in questo caso) che il segnale di ingresso è a media nulla. In generale sarebbe buona pratica sottrarre tale componente.
# Segnale nel tempo mediante FFT^-1
u <- u %>% mutate(
u_new = Re(fft(u_spectrum$U_filtered, inverse = T))/N
)Si conclude visivamente come adesso il segnale compensato tenda a sovrapporsi su quello ideale, peraltro con nessuna particolare discrepanza sulle code, a differenza del caso precedente. Cionostante, si osserva un lieve ritardo che va accumulandosi nel tempo. Come prova di questo fenomeno, si può osservare la differenza di fase tra i due segnali.
phase <- tibble(
f = ifelse(k <= N/2, k, k-N) * fs/N,
U_id = fft(u$u),
U_mask = fft(u$u_new),
diff = (Arg(U_id) - Arg(U_mask))*180/pi
) %>% dplyr::filter(f>0) # fase antisimmetrica rispetto alla frequenza
# Regressione qualitativa per evidenziare l'andamento della differenza di fase
phase.fit <- loess(diff ~ f, data = phase, span = 0.5)
phase$fit <- predict(phase.fit)
phase %>% select(-U_id, -U_mask) %>% pivot_longer(-f) %>%
ggplot(aes(x = f, y = value, color = name)) +
geom_line() +
scale_color_discrete(labels = c(diff = latex2exp::TeX("$\\Delta \\varphi"), fit = "Fit")) +
labs(x = "Frequenza (Hz)", y = "Fase (°)", colour = "Legenda")
È qualitativamente evidente come la differenza di fase non sia costante, che spiega il ritardo accumulato nel tempo. La ragione di questo fatto potrebbe risiedere in errori numerici, come per esempio nella computazione di fft() e della sua inversa o di freqresp relativa alla FdT. D’altro canto, si può supporre di attribuire questo ritardo non costante ad una fase non lineare della FdT inversa in corrispondenza delle armoniche del segnale (cfr. i diagrammi di Bode della FdT inversa). Ad ogni modo, è molto contenuto tale ritardo, percui ai fini dell’esercitazione si può accettare il risultato.