Main Content

fft

Trasformata di Fourier veloce

Descrizione

esempio

Y = fft(X) calcola la trasformata di Fourier discreta (DFT) di X utilizzando un algoritmo della trasformata di Fourier veloce (FFT). Y ha la stessa grandezza di X.

  • Se X è un vettore, fft(X) restituisce la trasformata di Fourier del vettore.

  • Se X è una matrice, fft(X) tratta le colonne di X come vettori e restituisce la trasformata di Fourier di ciascuna colonna.

  • Se X è un array multidimensionale, fft(X) tratta i valori lungo la prima dimensione dell’array, la cui grandezza non è uguale a 1, come vettori e restituisce la trasformata di Fourier di ciascun vettore.

esempio

Y = fft(X,n) restituisce la DFT a n punti.

  • Se X è un vettore e la lunghezza di X è minore di n, X viene riempita con zeri finali fino alla lunghezza di n.

  • Se X è un vettore e la lunghezza di X è maggiore di n, X viene troncata alla lunghezza di n.

  • Se X è una matrice, ciascuna colonna è trattata come nel caso vettoriale.

  • Se X è un array multidimensionale, la prima dimensione dell’array, la cui dimensione non è uguale a 1, viene trattata come nel caso vettoriale.

esempio

Y = fft(X,n,dim) restituisce la trasformata di Fourier lungo la dimensione dim. Ad esempio, se X è una matrice, fft(X,n,2) restituisce la trasformata di Fourier a n punti per ciascuna riga.

Esempi

comprimi tutto

Individuare i componenti di frequenza di un segnale immerso nel rumore e le ampiezze dei picchi di frequenza utilizzando la trasformata di Fourier.

Specificare i parametri di un segnale con una frequenza di campionamento di 1 kHz e una durata del segnale di 1,5 secondi.

Fs = 1000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector

Creare un segnale contenente un offset CC di ampiezza 0,8, una sinusoide a 50 Hz di ampiezza 0,7 e una sinusoide a 120 Hz di ampiezza 1.

S = 0.8 + 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

Corrompere il segnale con un rumore bianco casuale a media zero e varianza pari a 4.

X = S + 2*randn(size(t));

Tracciare il segnale rumoroso nel dominio del tempo. I componenti di frequenza non sono visivamente evidenti nel grafico.

plot(1000*t,X)
title("Signal Corrupted with Zero-Mean Random Noise")
xlabel("t (milliseconds)")
ylabel("X(t)")

Figure contains an axes object. The axes object with title Signal Corrupted with Zero-Mean Random Noise, xlabel t (milliseconds), ylabel X(t) contains an object of type line.

Calcolare la trasformata di Fourier del segnale.

Y = fft(X);

Poiché le trasformate di Fourier implicano numeri complessi, tracciare la magnitudine complessa dello spettro fft.

plot(Fs/L*(0:L-1),abs(Y),"LineWidth",3)
title("Complex Magnitude of fft Spectrum")
xlabel("f (Hz)")
ylabel("|fft(X)|")

Figure contains an axes object. The axes object with title Complex Magnitude of fft Spectrum, xlabel f (Hz), ylabel |fft(X)| contains an object of type line.

Il grafico mostra cinque picchi di frequenza, compreso il picco a 0 Hz per l'offset CC. In questo esempio, si presume che il segnale avrà tre picchi di frequenza a 0 Hz, 50 Hz e 120 Hz. In questo caso, la seconda metà del grafico è il riflesso speculare della prima metà, senza includere il picco a 0 Hz. Il motivo è che la trasformata discreta di Fourier di un segnale nel dominio del tempo ha una natura periodica, in cui la prima metà del relativo spettro è in frequenze positive e la seconda metà è in frequenze negative, con il primo elemento riservato alla frequenza zero.

Per i segnali reali, lo spettro fft è uno spettro a due lati, dove lo spettro nelle frequenze positive è il coniugato complesso dello spettro nelle frequenze negative. Per mostrare lo spettro fft nelle frequenze positive e negative, è possibile utilizzare fftshift. Per una lunghezza pari di L, il dominio di frequenza inizia dal negativo -Fs/2 della frequenza di Nyquist fino a Fs/2-Fs/L, con una spaziatura o risoluzione di frequenza di Fs/L.

plot(Fs/L*(-L/2:L/2-1),abs(fftshift(Y)),"LineWidth",3)
title("fft Spectrum in the Positive and Negative Frequencies")
xlabel("f (Hz)")
ylabel("|fft(X)|")

Figure contains an axes object. The axes object with title fft Spectrum in the Positive and Negative Frequencies, xlabel f (Hz), ylabel |fft(X)| contains an object of type line.

Per trovare le ampiezze dei tre picchi di frequenza, convertire lo spettro fft in Y nello spettro ad ampiezza unilaterale. Poiché la funzione fft include un fattore di scala L tra il segnale originale e quello trasformato, ridimensionare Y dividendolo per L. Prendere la magnitudine complessa dello spettro fft. Lo spettro di ampiezza a due lati P2, dove lo spettro nelle frequenze positive è il coniugato complesso dello spettro nelle frequenze negative, ha la metà delle ampiezze di picco del segnale nel dominio del tempo. Per convertire lo spettro unilaterale, prendere la prima metà dello spettro a due lati P2. Moltiplicare lo spettro delle frequenze positive per 2. Non è necessario moltiplicare P1(1) e P1(end) per 2 in quanto queste ampiezze corrispondono, rispettivamente, alle frequenze zero e Nyquist e non possiedono coppie di coniugati complessi nelle frequenze negative.

P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

Definire il dominio di frequenza f per lo spetto unilaterale. Tracciare lo spettro di ampiezza unilaterale P1. Come previsto, le ampiezze sono vicine a 0,8, 0,7 e 1, ma non sono esatte a causa del rumore aggiunto. Nella maggior parte dei casi, i segnali più lunghi producono migliori approssimazioni di frequenza.

f = Fs/L*(0:(L/2));
plot(f,P1,"LineWidth",3) 
title("Single-Sided Amplitude Spectrum of X(t)")
xlabel("f (Hz)")
ylabel("|P1(f)|")

Figure contains an axes object. The axes object with title Single-Sided Amplitude Spectrum of X(t), xlabel f (Hz), ylabel |P1(f)| contains an object of type line.

Ora, prendere la trasformata di Fourier del segnale originale non corrotto e recuperare le ampiezze esatte a 0,8, 0,7 e 1,0.

Y = fft(S);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

plot(f,P1,"LineWidth",3) 
title("Single-Sided Amplitude Spectrum of S(t)")
xlabel("f (Hz)")
ylabel("|P1(f)|")

Figure contains an axes object. The axes object with title Single-Sided Amplitude Spectrum of S(t), xlabel f (Hz), ylabel |P1(f)| contains an object of type line.

Convertire un impulso gaussiano dal dominio del tempo al dominio della frequenza.

Specificare i parametri di un segnale con una frequenza di campionamento di 44,1 kHz e una durata del segnale di 1 ms. Creare un impulso gaussiano con una deviazione standard di 0,1 ms.

Fs = 44100;         % Sampling frequency
T = 1/Fs;           % Sampling period
t = -0.5:T:0.5;     % Time vector
L = length(t);      % Signal length

X = 1/(0.4*sqrt(2*pi))*(exp(-t.^2/(2*(0.1*1e-3)^2)));

Tracciate l'impulso nel dominio del tempo.

plot(t,X)
title("Gaussian Pulse in Time Domain")
xlabel("Time (t)")
ylabel("X(t)")
axis([-1e-3 1e-3 0 1.1]) 

Figure contains an axes object. The axes object with title Gaussian Pulse in Time Domain, xlabel Time (t), ylabel X(t) contains an object of type line.

Il tempo di esecuzione della fft dipende dalla lunghezza della trasformata. Le trasformate che presentano solo fattori primi piccoli comportano tempi di esecuzione significativamente più rapidi rispetto a quelle che presentano fattori primi grandi.

In questo esempio, la lunghezza del segnale L è 44,101, che è un numero primo piuttosto grande. Per migliorare la prestazione della fft, individuare una lunghezza di input che sia la potenza prossima di 2 dalla lunghezza del segnale originale. Chiamando fft con questa lunghezza di input, l’impulso X viene riempito di zeri finali fino alla lunghezza specificata della trasformata.

n = 2^nextpow2(L);

Convertire l'impulso gaussiano nel dominio della frequenza.

Y = fft(X,n);

Definire il dominio della frequenza e tracciare le frequenze uniche.

f = Fs*(0:(n/2))/n;
P = abs(Y/sqrt(n)).^2;

plot(f,P(1:n/2+1)) 
title("Gaussian Pulse in Frequency Domain")
xlabel("f (Hz)")
ylabel("|P(f)|")

Figure contains an axes object. The axes object with title Gaussian Pulse in Frequency Domain, xlabel f (Hz), ylabel |P(f)| contains an object of type line.

Confrontare le onde cosiniche nel dominio del tempo e nel dominio della frequenza.

Specificare i parametri di un segnale con una frequenza di campionamento di 1 kHz e una durata del segnale di 1 secondo.

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sampling period
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector

Creare una matrice in cui ciascuna riga rappresenta un'onda cosinica con frequenza scalata. Il risultato X è una matrice 3x1000. La prima riga ha una frequenza d'onda pari a 50, la seconda riga ha una frequenza d'onda pari a 150 e la terza riga ha una frequenza d'onda pari a 300.

x1 = cos(2*pi*50*t);          % First row wave
x2 = cos(2*pi*150*t);         % Second row wave
x3 = cos(2*pi*300*t);         % Third row wave

X = [x1; x2; x3];

Tracciare le prime 100 voci di ciascuna riga di X in una singola figura in ordine e confrontare le frequenze.

for i = 1:3
    subplot(3,1,i)
    plot(t(1:100),X(i,1:100))
    title("Row " + num2str(i) + " in the Time Domain")
end

Figure contains 3 axes objects. Axes object 1 with title Row 1 in the Time Domain contains an object of type line. Axes object 2 with title Row 2 in the Time Domain contains an object of type line. Axes object 3 with title Row 3 in the Time Domain contains an object of type line.

Specificare l’argomento dim per utilizzare fft lungo le righe di X, ossia per ciascun segnale.

dim = 2;

Calcolare la trasformata di Fourier dei segnali.

Y = fft(X,L,dim);

Calcolare lo spettro bilaterale e quello unilaterale di ciascun segnale.

P2 = abs(Y/L);
P1 = P2(:,1:L/2+1);
P1(:,2:end-1) = 2*P1(:,2:end-1);

Nel dominio della frequenza, tracciare lo spettro di ampiezza unilaterale per ciascuna riga in un'unica figura.

for i=1:3
    subplot(3,1,i)
    plot(0:(Fs/L):(Fs/2-Fs/L),P1(i,1:L/2))
    title("Row " + num2str(i) + " in the Frequency Domain")
end

Figure contains 3 axes objects. Axes object 1 with title Row 1 in the Frequency Domain contains an object of type line. Axes object 2 with title Row 2 in the Frequency Domain contains an object of type line. Axes object 3 with title Row 3 in the Frequency Domain contains an object of type line.

Creare un segnale composto da due sinusoidi di frequenza 15 Hz e 40 Hz. La prima sinusoide è un'onda cosinica con fase -π/4, mentre la seconda è un'onda cosinica con fase π/2. Campionare il segnale a 100 Hz per un 1 s.

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

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))
title("Double-Sided Amplitude Spectrum of x(t)")
xlabel("Frequency (Hz)")
ylabel("|y|")
grid

Figure contains an axes object. The axes object with title Double-Sided Amplitude Spectrum of x(t), xlabel Frequency (Hz), ylabel |y| contains an object of type stem.

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)
title("Phase Spectrum of x(t)")
xlabel("Frequency (Hz)")
ylabel("Phase/\pi")
grid

Figure contains an axes object. The axes object with title Phase Spectrum of x(t), xlabel Frequency (Hz), ylabel Phase/ pi contains an object of type stem.

Interpolare la trasformata di Fourier di un segnale tramite il riempimento con zeri.

Specificare i parametri di un segnale con una frequenza di campionamento di 80 kHz e una durata del segnale di 0,8 s.

Fs = 80;
T = 1/Fs;
L = 65;
t = (0:L-1)*T;

Creare una sovrapposizione di un segnale sinusoidale di 2 Hz e delle sue armoniche superiori. Il segnale contiene un'onda cosinica di 2 Hz, un'onda cosinica di 4 Hz e un'onda sinusoidale di 6 Hz.

X = 3*cos(2*pi*2*t) + 2*cos(2*pi*4*t) + sin(2*pi*6*t);

Tracciare il segnale nel dominio del tempo.

plot(t,X)
title("Signal superposition in time domain")
xlabel("t (ms)")
ylabel("X(t)")

Figure contains an axes object. The axes object with title Signal superposition in time domain, xlabel t (ms), ylabel X(t) contains an object of type line.

Calcolare la trasformata di Fourier del segnale.

Y = fft(X);

Calcolare lo spettro di ampiezza unilaterale del segnale.

f = Fs*(0:(L-1)/2)/L;
P2 = abs(Y/L);
P1 = P2(1:(L+1)/2);
P1(2:end) = 2*P1(2:end);

Nel dominio della frequenza, tracciare lo spettro unilaterale. Poiché il campionamento temporale del segnale è piuttosto breve, la risoluzione di frequenza della trasformata di Fourier non è abbastanza precisa per evidenziare la frequenza di picco vicino a 4 Hz.

plot(f,P1,"-o") 
title("Single-Sided Spectrum of Original Signal")
xlabel("f (Hz)")
ylabel("|P1(f)|")

Figure contains an axes object. The axes object with title Single-Sided Spectrum of Original Signal, xlabel f (Hz), ylabel |P1(f)| contains an object of type line.

Per valutare meglio le frequenze di picco, è possibile aumentare la lunghezza della finestra di analisi riempiendo di zeri il segnale originale. Questo metodo interpola automaticamente la trasformata di Fourier del segnale con una risoluzione di frequenza più precisa.

Individuare una nuova lunghezza di input che sia la potenza prossima di 2 dalla lunghezza del segnale originale. Riempire il segnale X con zeri finali per estenderne la lunghezza. Calcolare la trasformata di Fourier del segnale riempito con zeri.

n = 2^nextpow2(L);
Y = fft(X,n);

Calcolare lo spettro di ampiezza unilaterale del segnale riempito. Poiché la lunghezza del segnale n è aumentata da 65 a 128, la risoluzione di frequenza diventa Fs/n, ossia 0,625 Hz.

f = Fs*(0:(n/2))/n;
P2 = abs(Y/L);
P1 = P2(1:n/2+1);
P1(2:end-1) = 2*P1(2:end-1);

Tracciare lo spetto unilaterale del segnale riempito. Questo nuovo spettro mostra le frequenze di picco vicino a 2 Hz, 4 Hz e 6 Hz con una risoluzione di frequenza di 0,625 Hz.

plot(f,P1,"-o") 
title("Single-Sided Spectrum of Padded Signal")
xlabel("f (Hz)")
ylabel("|P1(f)|")

Figure contains an axes object. The axes object with title Single-Sided Spectrum of Padded Signal, xlabel f (Hz), ylabel |P1(f)| contains an object of type line.

Argomenti di input

comprimi tutto

Array di input, specificato come vettore, matrice o array multidimensionale.

Se X è una matrice 0x0 vuota, fft(X) restituisce una matrice 0x0 vuota.

Tipi di dati: double | single | int8 | int16 | int32 | uint8 | uint16 | uint32 | logical
Supporto numeri complessi:

Lunghezza della trasformata, specificata come [] o come scalare intero non negativo. La definizione di uno scalare intero positivo per la lunghezza della trasformata può migliorare la prestazione della fft. La lunghezza è di norma specificata come una potenza di 2 o come un valore che può essere scomposto in un prodotto di numeri primi piccoli (con fattori primi non maggiori di 7). Se n è inferiore alla lunghezza del segnale, fft ignora i valori rimanenti del segnale dopo la n-esima voce e restituisce il risultato troncato. Se n è 0, fft restituisce una matrice vuota.

Esempio n = 2^nextpow2(size(X,1))

Tipi di dati: double | single | int8 | int16 | int32 | uint8 | uint16 | uint32 | logical

Dimensione lungo la quale operare, specificata come scalare intero positivo. Se non si specifica la dimensione, per impostazione predefinita viene considerata la prima dimensione dell'array la cui grandezza è maggiore di 1.

  • fft(X,[],1) opera lungo le colonne di X e restituisce la trasformata di Fourier di ciascuna colonna.

    fft(X,[],1) column-wise operation

  • fft(X,[],2) opera lungo le righe di X e restituisce la trasformata di Fourier di ciascuna riga.

    fft(X,[],2) row-wise operation

Se dim è maggiore di ndims(X), fft(X,[],dim) restituisce X. Quando n è specificato, fft(X,n,dim) esegue il riempimento o il troncamento di X alla lunghezza di n lungo la dimensione dim.

Tipi di dati: double | single | int8 | int16 | int32 | uint8 | uint16 | uint32 | logical

Argomenti di output

comprimi tutto

Rappresentazione nel dominio della frequenza, restituita come vettore, matrice o array multidimensionale.

Se X è di tipo single, fft calcola in modo nativo a precisione singola e anche Y è di tipo single. Diversamente, Y viene restituito come tipo double.

La dimensione di Y è la seguente:

  • Per Y = fft(X) o Y = fft(X,[],dim), la dimensione di Y è uguale alla dimensione di X.

  • Per Y = fft(X,n,dim), il valore di size(Y,dim) è uguale a n, mentre la grandezza di tutte le altre dimensioni rimane come in X.

Se X è reale, Y è simmetrico coniugato e il numero di punti unici in Y è ceil((n+1)/2).

Tipi di dati: double | single

Ulteriori informazioni

comprimi tutto

Trasformata di Fourier discreta di vettore

Y = fft(X) e X = ifft(Y) implementano rispettivamente la trasformata di Fourier e la trasformata di Fourier inversa. Per X e Y di lunghezza n, queste trasformate sono definite come segue:

Y(k)=j=1nX(j)Wn(j1)(k1)X(j)=1nk=1nY(k)Wn(j1)(k1),

dove

Wn=e(2πi)/n

è una delle radici n dell’unità.

Suggerimenti

  • Il tempo di esecuzione della fft dipende dalla lunghezza della trasformata. Le trasformate che presentano solo fattori primi piccoli (non superiori a 7) comportano tempi di esecuzione significativamente più rapidi rispetto a quelle in numero primo o che presentano fattori primi grandi.

  • Per la maggior parte dei valori di n, le DFT con input reale richiedono approssimativamente la metà del tempo di calcolo delle DFT con input complesso. Tuttavia, quando n presenta fattori primi grandi, la differenza di velocità è minima o nulla.

  • È eventualmente possibile aumentare la velocità di fft utilizzando la funzione di utilità fftw. Questa funzione controlla l'ottimizzazione dell'algoritmo utilizzato per calcolare una FFT di una particolare grandezza e dimensione.

Algoritmi

Le funzioni FFT (fft, fft2, fftn, ifft, ifft2, ifftn) sono basate su una libreria chiamata FFTW [1] [2].

Riferimenti

[2] Frigo, M., and S. G. Johnson. “FFTW: An Adaptive Software Architecture for the FFT.” Proceedings of the International Conference on Acoustics, Speech, and Signal Processing. Vol. 3, 1998, pp. 1381-1384.

Funzionalità estese

Generazione di codice GPU
Genera codice CUDA® per GPU NVIDIA® con GPU Coder™.

Cronologia versioni

Introduzione prima di R2006a