Trasformate di Fourier

La trasformata di Fourier è una formula matematica che trasforma un segnale campionato nel tempo o nello spazio nello stesso segnale campionato in frequenza temporale o spaziale. Nell'elaborazione dei segnali, la trasformata di Fourier può rivelare caratteristiche importanti di un segnale, ovvero i suoi componenti di frequenza.

La trasformata di Fourier è definita per un vettore $$x$$ con $$n$$ punti uniformemente campionati da

$$y_{k+1}=\sum _{j=0}^{n-1}{\omega }^{jk}{x}_{j+1}.$$

$$\omega =e^{-2\pi i/n}$$ è una delle radici complesse $$n$$ dell'unità, dove $$i$$ è l'unità immaginaria. Per $$x$$ e $$y$$, gli indici $$j$$ e $$k$$ vanno da $$0$$ a $$n-1$$.

La funzione fft in MATLAB® utilizza un algoritmo della trasformata di Fourier veloce per calcolare la trasformata di Fourier dei dati. Si consideri un segnale sinusoidale x che è una funzione del tempo t con componenti di frequenza di 15 Hz e 20 Hz. Utilizzare un vettore di tempo campionato con incrementi di 1/50 di secondo per un periodo di 10 secondi.

Ts = 1/50;
t = 0:Ts:10-Ts;
x = sin(2*pi*15*t) + sin(2*pi*20*t);
plot(t,x)
xlabel('Time (seconds)')
ylabel('Amplitude')

Calcolare la trasformata di Fourier del segnale e creare il vettore f che corrisponde al campionamento del segnale nello spazio di frequenza.

y = fft(x);
fs = 1/Ts;
f = (0:length(y)-1)*fs/length(y);

Tracciare la magnitudine del segnale trasformato in funzione della frequenza.

plot(f,abs(y))
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Magnitude')

Il grafico mostra quattro picchi di frequenza, anche se il segnale dovrebbe avere due picchi di frequenza a 15 Hz e 20 Hz. In questo caso, la seconda metà del grafico è il riflesso speculare della prima metà. La trasformata discreta di Fourier di un segnale nel dominio del tempo ha natura periodica, dove la prima metà del relativo spettro è in frequenze positive e la seconda metà in frequenze negative. I componenti di frequenza a 30 Hz e 35 Hz nel grafico corrispondono ai componenti di frequenza a –20 Hz e –15 Hz. Per visualizzare meglio questa periodicità, è possibile utilizzare la funzione fftshift, che esegue uno spostamento circolare centrato sullo zero sulla trasformata.

n = length(x);
fshift = (-n/2:n/2-1)*(fs/n);
yshift = fftshift(y);
plot(fshift,abs(yshift))
xlabel('Frequency (Hz)')
ylabel('Magnitude')

Segnali rumorosi

Nelle applicazioni scientifiche, i segnali sono spesso distorti da rumori casuali, che ne mascherano i componenti di frequenza. La trasformata di Fourier può elaborare il rumore casuale e rivelare le frequenze. Ad esempio, creare un nuovo segnale xnoise introducendo un rumore gaussiano nel segnale originale x.

rng('default')
xnoise = x + 2.5*randn(size(t));

La potenza del segnale in funzione della frequenza è una metrica comunemente utilizzata nell'elaborazione dei segnali. La potenza è la magnitudine al quadrato della trasformata di Fourier di un segnale, normalizzata per il numero di campioni di frequenza. Calcolare e tracciare lo spettro di potenza del segnale rumoroso centrato sulla frequenza zero. Nonostante il rumore, è ancora possibile distinguere le frequenze del segnale grazie ai picchi di potenza.

ynoise = fft(xnoise);
ynoiseshift = fftshift(ynoise);
power = abs(ynoiseshift).^2/n; 
plot(fshift,power)
title('Power')
xlabel('Frequency (Hz)')
ylabel('Power')

Efficienza computazionale

L'utilizzo diretto della formula della trasformata di Fourier per calcolare ciascuno degli $$n$$ elementi di $$y$$ richiede un ordine di $$n^2$$ operazioni in virgola mobile. L'algoritmo della trasformata di Fourier veloce richiede solo un ordine di $\mathit{n}\log \mathit{n}$ operazioni per il calcolo. Questa efficienza computazionale è un grande vantaggio quando si elaborano dati con milioni di punti di dati. Molte implementazioni specializzate dell'algoritmo della trasformata di Fourier veloce sono ancora più efficienti quando $$n$$ ha fattori primi piccoli, come nel caso in cui $$n$$ è una potenza di 2.

Si considerino i dati audio acquisiti dai microfoni subacquei al largo della costa della California. Questi dati sono disponibili in una libreria gestita dal Cornell University Bioacoustics Research Program. Caricare e formattare un sottoinsieme dei dati in bluewhale.au, che contiene una vocalizzazione di una balenottera azzurra del Pacifico. Poiché i richiami emessi dalla balenottera azzurra sono suoni a bassa frequenza, sono appena udibili dall'orecchio umano. La scala temporale nei dati è compressa di un fattore 10 per aumentare il tono e rendere più chiaramente udibile il richiamo. È possibile utilizzare il comando sound(x,fs) per ascoltare l'intero file audio.

whaleFile = 'bluewhale.au';
[x,fs] = audioread(whaleFile);
whaleMoan = x(2.45e4:3.10e4);
t = 10*(0:1/fs:(length(whaleMoan)-1)/fs);

plot(t,whaleMoan)
xlabel('Time (seconds)')
ylabel('Amplitude')
xlim([0 t(end)])

Specificare una nuova lunghezza del segnale che corrisponda alla potenza di 2 immediatamente superiore alla lunghezza originale. Quindi, utilizzare fft per calcolare la trasformata di Fourier utilizzando la nuova lunghezza del segnale. fft riempie automaticamente i dati con zeri per aumentare la dimensione del campione. Questo riempimento può rendere il calcolo della trasformata significativamente più veloce, in particolare per campioni con fattori primi di grandi dimensioni.

m = length(whaleMoan);
n = pow2(nextpow2(m));
y = fft(whaleMoan,n);

Tracciare lo spettro di potenza del segnale. Il grafico indica che la vocalizzazione è composta da una frequenza fondamentale intorno ai 17 Hz e da una sequenza di armoniche, dove la seconda armonica è enfatizzata.

f = (0:n-1)*(fs/n)/10; % frequency vector
power = abs(y).^2/n;   % power spectrum
plot(f(1:floor(n/2)),power(1:floor(n/2)))
xlabel('Frequency (Hz)')
ylabel('Power')

Fase delle sinusoidi

Utilizzando la trasformata di Fourier, è inoltre possibile estrarre lo spettro di fase del segnale originale. Ad esempio, creare un segnale composto da due sinusoidi di frequenza 15 Hz e 40 Hz. La prima sinusoide è un'onda cosinica con fase $-\pi /4$, mentre la seconda è un'onda cosinica con fase $\pi /2$. Campionare il segnale a 100 Hz per un 1 secondo.

fs = 100;
t = 0:1/fs:1-1/fs;
x = cos(2*pi*15*t - pi/4) - sin(2*pi*40*t);

Calcolare la trasformata di Fourier del segnale. Tracciare la magnitudine della trasformata come una funzione di frequenza.

y = fft(x);
z = fftshift(y);

ly = length(y);
f = (-ly/2:ly/2-1)/ly*fs;

stem(f,abs(z))
xlabel("Frequency (Hz)")
ylabel("|y|")
grid

Calcolare la fase della trasformata, eliminando i valori di piccola entità della trasformata. Tracciare la fase come una funzione di frequenza.

tol = 1e-6;
z(abs(z) < tol) = 0;

theta = angle(z);

stem(f,theta/pi)
xlabel("Frequency (Hz)")
ylabel("Phase / \pi")
grid

Vedi anche

fft | fftshift | nextpow2 | ifft | fft2 | fftn | fftw

Argomenti

• 2-D Fourier Transforms

Siti web esterni

• Fourier Analysis (Analisi di Fourier) (Risorse didattiche MathWorks)