GNU Octave: esempi di Digital Signal Processing

GNU Octave: esempi di Digital Signal Processing

Ing. Andrea [email protected]

Dipartimento di Ingegneria Informatica e delle TelecomunicazioniUniversita di Catania

20 Maggio 2009

Sommario

Elaborazioni nel dominio del tempo

Esempio di analisi in frequenza

Ing. Andrea Spadaccini DIIT - Universita degli Studi di Catania


Suddivisione in trame Funzioni finestra Ricampionamento Cross-correlazione

Parte I

Elaborazioni nel dominio del tempo




La funzione buffer

Suddivisione in trame

Spesso, per elaborare un segnale puo essere utile suddividerlo in trame,eventualmente sovrapponendo parte del segnale.In GNU Octave, questo si ottiene con la funzione buffer, presente nelpackage Signal Processing.

Sintassi di buffer

buffer(signal, frame-size, overlap)

signal il segnale da suddividere in trame

frame-size dimensione della trama

overlap numero di campioni della trama precedente da riportare (oda saltare, se negativo)




buffer: Esempio

>> s = 1:20;>> b = buffer(s, 4)b =

1 5 9 13 172 6 10 14 183 7 11 15 194 8 12 16 20

>> b = buffer(s, 4, 1)b =

0 3 6 9 12 15 181 4 7 10 13 16 192 5 8 11 14 17 203 6 9 12 15 18 0




buffer: Esempio con overlap negativo

>> s = 1:20;>> b = buffer(s, 4, -1)b =

1 6 11 162 7 12 173 8 13 184 9 14 19




Funzioni finestra

Scopo

Per evitare problemi come la perdita spettrale (spectral leakage) a voltesi ricorre all’applicazioni di finestre differenti dalla semplice finestrarettangolare.

Esempi di finestra

Entrambe le funzioni richiedono come parametro solo il numero di punti.

Hamming

Hann (o Hanning)




Plot delle due finestre

>> figure;>> subplot(311);>> plot(hamming(100));>> subplot(312);>> plot(hanning(100));>> subplot(313);>> plot(hamming(100) - hanning(100));




Ricampionamento (resampling)

Scopo

L’idea alla base del ricampionamento (resampling) e quella di variare lafrequenza di campionamento di un segnale, generalmentesottocampionandolo per elaborazioni successive.

La funzione resample

Sintassi: y = resample(x, p, q)

x e il segnale da ricampionare

p e q costituiscono il numero p/q, ovvero il fattore da utilizzare pervariare la frequenza di campionamento.




Esempio di ricampionamento

>> Fs = 100;>> Fo = 4;>> Ts = 1 / Fs;>> t = 0 : Ts : 1 - Ts;>> y = sin(2 * pi * Fo * t);>> plot(t, y)>> rs_y = resample(y, 1, 2);>> new_t = 0 : 2 * Ts : 1 - Ts>> figure>> subplot(211)>> plot(t, y)>> subplot(212)>> plot(new_t, rs_y);




Cross-correlazione (I)

Un pizzico di teoria

Autocorrelazione di due segnali discreti x ed y al variare dell’intervallo(lag) k :

Rxy (k) =

1

N − k

N−k∑i=1

xiyi−k k ≥ 0

Ryx(−k) k ≤ 0

(1)




Cross-correlazione (II)

La funzione xcorr

La funzione xcorr consente di computare la cross-correlazione di duesegnali. Sintassi:[R, lag] = xcorr(X[, Y] [, maxlag] [, scale]).

X e il primo segnale

Y, se presente, e il secondo segnale

maxlag, se presente, imposta l’intervallo di variazione di lag a[−maxlag:maxlag] (altrimenti maxlag = max(length(X),length(Y)) - 1)

scale, se presente, definisce la tipologia di normalizzazione dellasomma. Utilizzare ’coeff’ per normalizzare con Rxy (0) (ovvero ilmassimo).




Esempio: periodo di un segnale periodico (I)

Autocorrelazione e segnali periodici

Nel caso di segnali reali periodici (infiniti!), il lag del primo massimopositivo dell’autocorrelazione corrisponde al periodo del segnale.

Implementeremo qualcosa di simile, anche se non possiamo generaresegnali periodici.

>> T = 10; % periodo desiderato: 10 secondi>> Fs = 1000;>> Ts = 1 / Fs;>> t = 1 : Ts : 100 - Ts;>> y = sin(2 * pi * (1/T) * t);>> [R, lag] = xcorr(y);




Esempio: periodo di un segnale periodico (II)

>> mid = ceil(length(R) / 2);>> [min_val, min_idx] = min(R(mid:end));>> [max_val, max_idx] = max(R(mid + min_idx:end));>> idx = mid + min_idx + max_idx;>> lag(idx) * Tsans = 9.9910



FFT di una cosinusoide

Parte II

Esempio di analisi in frequenza




Trasformata di Fourier di un coseno

Un pizzico di teoria (di nuovo?)

Trasformata di Fourier di un coseno:

F [cos(2πf0t)] =1

2δ(f + f0) +

1

2δ(f − f0)

FT, DFT, FFT

La normale Trasformata di Fourier (FT) si applica a segnali continui.

La DFT (Discrete Fourier Transform) e la FT per segnali discreti.

La FFT (Fast Fourier Transform) e un’implementazione (veloce!) diDFT. Variante piu utilizzata: FFTW (Fastest Fourier Transform inthe West).




Primi passi

La funzione fft

Prende in input un vettore (o una matrice), e restituisce un vettorecontenente il risultato dell’esecuzione dell’algoritmo FFT.

>> Fs = 100;>> Ts = 1 / Fs;>> F0 = 4;>> t = 0 : Ts : 1 - Ts;>> y = sin(2*pi*F0*t);>> plot(t, y);>> Fy = fft(f);>> plot(abs(Fy)); % difficile da comprendere




Il giusto ordine

La funzione fftswitch

Manipola il risultato di fft in modo da avvicinare le frequenze basse alcentro del vettore (in sostanza effettua il flipping delle due meta delrisultato di fft)

>> Fy = abs(fftswitch(Fy));>> plot(Fy)>> stem(Fy) % tipo di grafico piu indicativo




Le giuste proporzioni

L’asse delle frequenze

Ciascun elemento del vettore FFT non rappresenta una frequenza, ma uninsieme (bin) di frequenze. Dobbiamo generare un asse delle frequenze infunzione della frequenza di campionamento.

>> n = length(Fy);>> f = -n/2 : n/2 - 1;>> f = f * (Fs / n);>> Fy = Fy / n; % normalizzazione valori FFT>> stem(f, Fy) % tipo di grafico piu indicativo>> xlabel(’Frequency bins (Hz)’)



Riferimenti

Sito ufficiale di GNU Octave

http://www.octave.org

Documentazione

John W. Eaton et al. GNU Octave Manual -http://www.gnu.org/software/octave/doc/interpreter/http://en.wikipedia.org/wiki/Window functionhttp://blinkdagger.com/matlab/matlab-introductory-fft-tutorialhttp://fftw.org



Licenza

Licenza Creative Commons

Questo lavoro e rilasciato sotto licenza Creative CommonsAttribution-Share Alike 2.5 Italy(http://creativecommons.org/licenses/by-sa/2.5/it/)

Significato della licenza

La presentazione e liberamente scaricabile e redistribuibile perqualsiasi fine

E necessario citare l’autore

E necessario redistribuire secondo gli stessi termini eventuali lavoriderivati da questa presentazione



GNU Octave: esempi di Digital Signal Processing

Documents

Transcript of GNU Octave: esempi di Digital Signal Processing