4. Inversione delle armoniche pure (DARE)

filtri
FFT
Autore/Autrice

Mariolino De Cecco e Paolo Bosetti

Data di Pubblicazione

15 maggio 2025

Sommario

Questo esercizio genera delle armoniche pure in ingresso ad un sistema di isolamento delle vibrazioni, simulazione dell’uscita e—mediante inversione della dinamica—stimare quella originaria in ingresso.

1 Parte fornita per l’esercitazione

1.1 Funzioni

# 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)")
}

# 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] ))))
  })
}

1.2 Segnale in ingresso

Definiamo un segnale d’ingresso (nel caso di sistemi di misura si tratta di un misurando) costituito da segnali armonici, più un eventuale disturbo normale:

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)
)

Nt <- 1000 # numero totale dei campioni
s <- tibble(
  t = 0:Nt * 0.001, # 1 kHz di frequenza di campionamento
  u = signal(t, pars), 
  un = u + rnorm(length(t), 0, pars$w[1] / 10)
)

La trasformata di Fourier mostra i picchi attesi:

s %>% 
  mutate(
    f = 0:(n()-1)/max(t),
    fft_u = fft(un),
    intensity_u = Mod(fft_u) / n()*2,
    phase_u = Arg(fft_u)/pi*180
  ) %>% 
  slice_head(n=as.integer(nrow(.)/2)) %>%
  ggplot(aes(x=f, y=intensity_u)) +
  geom_spoke(aes(y=0, radius=intensity_u, angle=pi/2)) +
  geom_point()

FFT del segnale armonico di riferimento

1.3 Impieghiamo un sistema per isolamento da vibrazioni per simulare l’uscita

Con il metodo delle impedenze generalizzate si ottiene:

\[ H(i\omega)=\frac{V_\mathrm{out}(i\omega)}{V_\mathrm{in}(i\omega)}=\frac{C i\omega + K}{M (i\omega)^2 + C i\omega + K} \tag{1}\]

La frequenza naturale del sistema è \(f_0=\frac{1}{2\pi}\sqrt{\frac{K}{M}}\), e l’attenuazione comincia a \(\sqrt{2}f_0\).

Possiamo definire la funzione di trasferimento in Equazione 1 con la funzione control::tf(), che prende come due argomenti due vettori con i coefficienti della Equazione 1, in ordine decrescente di grado della variabile \(i\omega\):

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(H, 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"))
Registered S3 methods overwritten by 'signal':
  method         from   
  print.freqs    gsignal
  print.freqz    gsignal
  print.grpdelay gsignal
  plot.grpdelay  gsignal
  print.impz     gsignal
  print.specgram gsignal
  plot.specgram  gsignal

Bode plot per il sistema di isolamento da vibrazioni
cat("Per la frequenza", f0, "Hz abbiamo: \nModulo:", Mod(freqresp(H, 2*pi*f0)), "\nFase:", Arg(freqresp(H, 2*pi*f0))*180/pi, "°\n")
Per la frequenza 10 Hz abbiamo: 
Modulo: 0.08539786 
Fase: -102.9892 °
Versione alternativa

O anche, usando glue, è più semplice comporre delle stringe inserendo tra graffe le espressioni da valutare:

glue("Per la frequenza {f0} Hz abbiamo: \nModulo: {Mod(freqresp(H, 2*pi*f0)) %>% round(3)}\nFase: {(Arg(freqresp(H, 2*pi*f0))*180/pi) %>% round(3)}°")
Per la frequenza 10 Hz abbiamo: 
Modulo: 0.085
Fase: -102.989°

Dunque, se inseriamo un’armonica a fase nulla come un coseno \(\cos(\omega_{0} t)\) di frequenza pari a 10 Hz avremo in uscita la stessa armonica attenuata di 0.085 e sfasata di 103°. Dovrebbe comportarsi effettivamente come un sistema di isolamento delle vibrazioni.

1.4 Simulazione uscita

Andiamo a verificarlo. Simuliamo cosa succede in uscita se imponiamo il segnale al sistema d’isolamento impiegando il simulatore integrato in R lsim().

output <- lsim(H, s$un, s$t)
s <- s %>% 
  mutate(
    # Simulazione dell'uscita
    y = output$y[1,]
  )

# Grafici
pp <- plot_ly() %>%
  add_lines(s$t, s$y, name = "output", line= list(color= "red")) %>%
  add_lines(s$t, s$un, name = "input", line= list(color= "blue")) %>%
  layout(title = "Input & Output")
pp
Versione alternativa

Attenzione: i grafici che usano plotly sono oggetti Javascript e in quanto tali non sono compatibili con LaTeX, quindi non è possibile generare anche la versione pdf se il documento contiene grafici plotly.

(s %>% 
  select(t, input=y, output=un) %>% 
  pivot_longer(-t, names_to = "segnale", values_to = "valore") %>% 
  ggplot(aes(x=t, y=valore, color=segnale)) +
  geom_line()) %>% 
  ggplotly()

Se avete usato come parametro frequenza della funzione creata 10 Hz potete verificare come il segnale in uscita si sia effettivamente attenuato di parecchio!

Esercitazione

Riconoscimento punti esame

Il completamento di questo esercizio comporta il riconoscmiento di ±1 punto.

Da fare

Stimare il segnale in ingresso originale a partire dal segnale in uscita e dal modulo e la fase della funzione di trasferimento

NOTA: impiegare funzioni che processano automaticamente tutte le armoniche che sono presenti nel segnale. In altre parole impiegare nelle funzioni che costruirete la struttura della funzione signal()

Suggerimenti

  • calcolare FFT
  • trovare i picchi FFT
  • applicare inversa di modulo e fase della funz trasferimento alle armoniche corrispondenti ai picchi del modulo della FFT
  • scrivere una funzione che stima l’ingresso da generiche componenti armoniche in uscita che avrete salvato in un contenitore dati ‘picchi’ ed H(w): signal_input <- function(t, picchi, H) {}
  • quando il vostro codice funziona provate ad aggiungere altre armoniche
  • quando il vostro codice funziona provate ad aggiungere rumore