1 Analisi di segnali
In questo tutorial vedremo come implementare la routine per il calcolo della trasformata di Fourier discreta. La useremo poi per riprodurre l'analisi della funzione 161, fatta nelle dispense.
Come vediamo è costituita da una sommatoria di termini esponenziali complessi dipendenti da due indici moltiplicati per i valori della funzione di cui si vuole fare la trasformata . L'indice è anche quello su cui si esegue la sommatoria, per ogni valore dell'indice fissato, rappresenterà i valori discreti della trasformata . Il tempo di campionamento invece è rappresentato da . Vediamo di seguito le righe di codice che permettono di imprementare questa formula:
1 2 3 4 5 6 7 8 9 10 11 12 13 | #definisco una funzione h di prova
tau=0.05
t=arange(-10,10,tau)
T=2.
h=sin(2*pi*t/T)
N=len(h)
H=zeros(N,dtype='complex')
for n in range(N):
for k in range(N):
H[n] += h[k]*exp(2j*pi*k*n/N)
H*=tau
|
- nelle righe 2-5 definiamo una funzione di prova molto semplice con una frequenza nota. Solitamente la funzione h di cui si vuol fare la trasformata è data da un segnale proveniente da una qualche sorgente che quindi verrà inizializzata partendo da dati registrati su un file (vedremo meglio in seguito un caso di questo tipo).
- nelle due righe successive, definiamo il numero N di componenti della funzione e poi il vettore H che conterrà la trasformata. Poichè la trasformata in genere restituisce dei valori complessi, definiamo come tali le componenti di H, usando dtype='complex'.
- nelle righe successive troviamo i due cicli for che permettono il calcolo della traasformata: fissato l'indice n dal relativo ciclo, per ogni valore dell'indice k calcoliamo l'esponenziale complesso, lo moltiplichiamo per la componente k-sima della funzione h e aggiungiamo il risultato al valore calcolato per l'indice k precedente. Una volta eseguite tutte le iterazioni, moltiplichiamo tutte le componenti del vettore H per il tempo tau.
- soffermiamoci sul termine 2j usato dentro l'esponenziale: è la forma rapida per rappresentare il numero immaginario ed è del tutto uguale a scrivere 2*1j.
Poichè la routine che calcola la trasformata è indipendente dal tipo di funzione di partenza possiamo racchiuderla all'interno di una funzione python che potremo chiamare ogni volta che vogliamo dandogli in ingresso la funzione di cui vogliamo calcolare la trasformata. Vediamo dunque come definire una funzione in python:
1 2 3 4 5 6 7 8 9 10 11 12 | def my_dft(h,tau):
N=len(h)
H=zeros(N,dtype='complex')
for n in range(N):
for k in range(N):
H[n] += h[k]*exp(2j*pi*k*n/N)
H *= tau
return H
#uso la funzione appena definita
Transf_f1 = my_dft(f1,tau1)
Transf_f2 = my_dft(f2,tau2)
|
- con def nome_funzione(arg1,arg2,...): si definisce una funzione in python, avente come argomenti in ingresso le variabili (di qualunque tipo) arg1, arg2,... etc.
- Le variabili in ingresso saranno usate all'interno della funzione, cioè in tutte le righe che dopo il def sono indentate sino alla riga return.
- Come si nota nelle ultime due righe, il vantaggio di definire una funzione come my_dft è quello che, una volta fatto, potremo richiamarla tutte le volte che vorremo dandogli in ingresso le funzioni che vogliamo trasformare senza dover cambiare ogni volta nella routine il nome delle delle variabili che contengono la funzione stessa.
- La routine della trasformata viene quindi scritta una sola volta, all'interno della funzione my_dft e posta solitamente all'inizio del file (o anche in un file a parte contenente tutte le funzioni che definiamo).
- la riga return H permette di restituire la variabile H calcolata all'interno della funzione, quando questa finisce, e quindi di poterla salvare in un'altra variabile, come si fa nelle ultime due righe di codice.
In maniera del tutto analoga a quanto fatto per la trasformata, possiamo definire una funzione che calcoli l'anti-trasformata usando la formula 163 delle dispense.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | def my_antidft(H,tau):
N=len(H)
h=zeros(N,dtype='complex')
for n in range(N):
for k in range(N):
h[n] += H[k]*exp(-2j*pi*k*n/N)
H *= 1/tau/N
return h
#uso la funzione appena definita
Transf_f1 = my_dft(f1,tau1)
f2 = my_antidft(Transf_f1,tau1)
plot(f1)
plot(f2)
|
Come vediamo la routine eà del tutto simile a quella scritta per la trasformata. Nelle ultime quattro righe facciamo la trasformata di una funzione f1 e la salviamo nel vettore Transf_f1. Poi facciamo l'anti-trasformata di quest'ultima e salviamo il risultato nel vettore f2. Cosa vi aspettate dal plot delle funzioni f1 e f2?
Ora che abbiamo gli strumenti necessari possiamo studiare la funzione 166 delle dispense nella configurazione di partenza, cioe con tempo di campionamento e numero totale di punti , quindi periodo , e frequenze multipli di 2. Vediamo il codice:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 | #supponiamo che siano state definite
#sopra le funzioni my_dft e my_antidft
tau=1/32.
T=1.
A=1.
f0=2.
f1=2*f0
f2=4*f0
t=arange(0,T,tau)
N=len(t) #o anche N=T/tau
funct = 1 + A*(sin(2*pi*f0*t)+2*sin(2*pi*f1*t)+0.5*sin(2*pi*f2*t))
Transf = my_dft(funct,tau)
freq=arange(N)/tau/N
bar(freq,abs(Transf)**2,width=0.5,align='center',color='b')
figure(2)
bar(freq,Transf.real,width=0.5,align='center',color='b',label='Reale')
bar(freq,Transf.imag,width=0.5,align='center',color='r',label=Immaginaria')
|
Come vediamo, una volta definite le variabili utili al calcolo della funzione, definiamo un vettore con i tempi, la funzione funct che vogliamo analizzare e ne facciamo la trasformata chiamando la nostra funzione definita all'inizio del file e salviamo il risultato dentro la variabile Transf.
Definiamo inoltre un vettore contenente i valori delle frequenze mediante la formula 164 delle dispense. Da notare che in questo modo, però, saranno corrette solo le prime , che sono comunque quelle di interesse. Le successive frequenze sono delle repliche traslate di (vedasi dispense e grafici seguenti).
I plot che riportiamo di seguito sono stati eseguiti mediante la funzione bar che disegna delle barre verticali in corrispondenza di ogni frequenza con altezza pari al valore della trasformata. I primi due argomenti di ingresso sono gli array delle frequenze e del potere spettrale, o della parte reale e immaginaria. Notiamo come per fare il modulo di un vettore complesso si possa ancora usare la funzione abs, già vista nei precedenti tutorials. Inoltre per estrarre dall'array complesso Transf le due componenti usiamo la forma Transf.real/.imag che di fatto genera due array.
Gli altri parametri di ingresso sono abbastanza espliciti: larghezza, allineamento centrato rispetto ai valori dell'asse x e colore delle barre.
Il grafico a sinistra rappresenta il potere spettrale, mentre in quello a destra troviamo la parte reale e immaginaria della trasformata. Osservando i grafici notiamo come le tutte le frequenze che compongono il segnale vengano identificate con precisione e si trovano nella parte immaginaria della trasformata. Nella parte reale invece troviamo la sola componente a frequenza nulla che non è altro che la costante che abbiamo nel segnale. Svolgendo analiticamente la trasformata si ha una conferma dei diversi contributi che si osservano.
Vediamo di passaggio un interessante esempio, sicuramente il più semplice, di come si possa modificare il segnale nello spazio delle frequenze. Supponiamo di voler togliere la componente a frequenza nulla dalla parte reale e tornando indietro allo spazio dei tempi, mediante l'anti-trasformata, vediamo come il segnale viene modificato. Di seguito il codice necessario:
1 2 3 4 5 6 7 | #eliminiamo la componente a frequenza zero
Transf.real[0] = 0
antiTransf = my_antidft(Transf,tau)
plot(t,antiTransf,'-s',label='funct modificata')
plot(t,funct,'-d',label='funct originale')
|
Come osserviamo la funzione modificata è stata soltanto traslata in basso di una quantità pari a . Questo è in linea con quanto abbiamo accennato in precedenza: la componente reale della trasformata contiene la parte costante della funzione, perciò eliminare quella componente nello spazio delle frequenze, significa eliminare la costante nello spazio dei tempi. Cosa ci aspettiamo dunque se dovessimo togliere dalla trasformata le qualcuna delle altre componenti?
- cambiamo il periodo di campionamento da a ;
- aumentiamo ancora il periodo portandolo a ;
- modifichiamo ora le frequenze in e proviamo diversi periodi di campionamento.
dove abbiamo fatto andare il programma per due volte, modificando il solo periodo, usando la funzione bar(freq,abs(Transf)**2) per il primo grafico e la funzione plot(freq,abs(Transf)**2) per aggiungere le linee nel secondo. Quello che osserviamo nel primo grafico è che la frequenza 2 non è più identificata e il potere spettrale si è spostato su frequenze spurie ad essa affiancate. Questo spostamento del peso spettrale si nota anche per la frequenza 4, anche se questa resta identificata insieme alla 8. Questo è dovuto al semplice fatto che il periodo non è più un multiplo della frequenza 2, mentre lo è ancora per le frequenze 4 e 8. Si provi a modificare in maniera ancora diversa il periodo e si osservi come cambia lo spettro delle frequenze.
Nel secondo grafico abbiamo fatto il grafico del potere spettrale per il caso di un periodo pari a 19.6. In questo caso osserviamo che aumentando il periodo di quasi venti volte si ottengono dei picchi in corrispondenza delle frequenze corrette, ma come effetto negativo il peso spettrale si sposta anche su frequenze adiacenti a quelle reali, in maniera differente come si nota dalla diversa larghezza dei picchi. Questo è il fenomeno del leakage. Si provi a modificare il periodo di campionamento con valori vicini a 20 (multiplo del periodo) e si osservi come varia lo spettro delle frequenze.
Passiamo ora al terzo caso, in cui abbiamo intenzionalmente posto delle frequenze che difficilmente saranno contemporaneamente suttomultipli di un periodo di campionamento.
Come osserviamo se il periodo di campionamento è corto ( per esempio) le frequenze che compongono il segnale non vengono identificate correttamente, mentre al loro posto ne compaiono delle altre. Aumentando il periodo di campionamento ( per esempio) riusciamo a identificare le corrette frequenze anche se un pò di leakage è presente. Questo mostra dunque l'importanza di avere un tempo di campionamento il più lungo possibile per assicurarsi di poter trovare le corrette frequenze di un certo segnale.
Una tecnica usata per ridurre il leakage è quella del windowing: moltiplicare il segnale per una funzione finestra che vada a zero lentamente. La funzione finestra che usiamo di seguito, tra le tante che si possono usare, è la 167 delle dispense. Di seguito riportiamo il codice per implementare il windowing:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | #definiamo la funzione finestra
r=2*pi*t/T
han=0.5-0.5*cos(r)
funct_windowed=funct*han
Transf_window=my_dft(funct_windowed,tau)
#facciamo i grafici
plot(t,han)
figure(2)
plot(t,funct)
plot(t,funct_windowed)
figure(3)
plot(freq,abs(Transf)**2,label='Senza window')
plot(freq,abs(Transf_window)**2,label='Con window')
|
