MATLAB Language
Trasformate di Fourier e Trasformate di Fourier inverse

Sintassi

1. Y = fft (X)% calcola la FFT di Vector o Matrix X utilizzando una lunghezza di trasformazione predefinita di 256 (da confermare per la versione)

2. Y = fft (X, n)% calcola la FFT di X usando n come Transform Length, n deve essere un numero basato su 2-power. Se la lunghezza di X è minore di n, allora Matlab calcola automaticamente X con zeri tali che la lunghezza (X) = n

3. Y = fft (X, n, dim)% calcola la FFT di X usando n come Transform Length lungo dim dimensionale (può essere 1 o 2 per orizzontale o verticale rispettivamente)

4. Y = fft2 (X)% Calcola la FFT 2D di X

5. Y = fftn (X, dim)% Calcola la FFT dim-dimensionale di X, rispetto al vettore di dimensioni dim.

6. y = ifft (X)% calcola l'inverso di FFT di X (che è una matrice / vettore di numeri) utilizzando la lunghezza di trasformazione 256 predefinita

7. y = ifft (X, n)% calcola l'IFFT di X usando n come Transform Length

8. y = ifft (X, n, dim)% calcola l'IFFT di X usando n come Transform Length over dim dim (può essere 1 o 2 per orizzontale o verticale rispettivamente)

9. y = ifft (X, n, dim, 'symmetric')% L'opzione Symmetric fa sì che ifft tratti X come coniugato simmetrico lungo la dimensione attiva. Questa opzione è utile quando X non è esattamente simmetrico coniugato, semplicemente a causa di un errore di arrotondamento.

10. y = ifft2 (X)% Calcola i metri 2D inversi di X

11. y = ifftn (X, dim)% Calcola il fft dim dimensionale inverso di X.

Parametri

Osservazioni

ifft(fft(X)) = X

Implementa una semplice trasformazione di Fourier in Matlab

Fourier Transform è probabilmente la prima lezione in Digital Signal Processing, la sua applicazione è ovunque ed è uno strumento potente quando si tratta di analizzare dati (in tutti i settori) o segnali. Matlab ha un set di potenti toolbox per la trasformata di Fourier. In questo esempio, useremo la trasformata di Fourier per analizzare un segnale sinusoidale di base e generare quello che a volte è noto come un periodogramma usando FFT:

%Signal Generation
A1=10;                % Amplitude 1
A2=10;                % Amplitude 2
w1=2*pi*0.2;          % Angular frequency 1
w2=2*pi*0.225;        % Angular frequency 2
Ts=1;                 % Sampling time
N=64;                 % Number of process samples to be generated
K=5;                  % Number of independent process realizations
sgm=1;                % Standard deviation of the noise
n=repmat([0:N-1].',1,K);             % Generate resolution
phi1=repmat(rand(1,K)*2*pi,N,1);     % Random phase matrix 1
phi2=repmat(rand(1,K)*2*pi,N,1);     % Random phase matrix 2
x=A1*sin(w1*n*Ts+phi1)+A2*sin(w2*n*Ts+phi2)+sgm*randn(N,K);   % Resulting Signal

NFFT=256;            % FFT length
F=fft(x,NFFT);       % Fast Fourier Transform Result
Z=1/N*abs(F).^2;     % Convert FFT result into a Periodogram

Nota che la Trasformata di Fourier Discreta è implementata da Fast Fourier Transform (fft) in Matlab, entrambi produrranno lo stesso risultato, ma FFT è un'implementazione veloce di DFT.

figure
w=linspace(0,2,NFFT);
plot(w,10*log10(Z)),grid;
xlabel('w [\pi rad/s]')
ylabel('Z(f) [dB]')
title('Frequency Range: [ 0 ,  \omega_s ]')

Trasformate di Fourier inverse

Uno dei principali vantaggi di Fourier Transform è la sua capacità di tornare al dominio del tempo senza perdere informazioni. Consideriamo lo stesso segnale che abbiamo usato nell'esempio precedente:

A1=10;                % Amplitude 1
A2=10;                % Amplitude 2
w1=2*pi*0.2;          % Angular frequency 1
w2=2*pi*0.225;        % Angular frequency 2
Ts=1;                 % Sampling time
N=64;                 % Number of process samples to be generated
K=1;                  % Number of independent process realizations
sgm=1;                % Standard deviation of the noise
n=repmat([0:N-1].',1,K);             % Generate resolution
phi1=repmat(rand(1,K)*2*pi,N,1);     % Random phase matrix 1
phi2=repmat(rand(1,K)*2*pi,N,1);     % Random phase matrix 2
x=A1*sin(w1*n*Ts+phi1)+A2*sin(w2*n*Ts+phi2)+sgm*randn(N,K);   % Resulting Signal

NFFT=256;            % FFT length
F=fft(x,NFFT);       % FFT result of time domain signal

Se apriamo F in Matlab, scopriremo che è una matrice di numeri complessi, una parte reale e una parte immaginaria. Per definizione, al fine di recuperare il segnale del dominio del tempo originale, abbiamo bisogno sia del Real (che rappresenta la variazione di Magnitudo) sia dell'Imaginario (che rappresenta la variazione di Fase), quindi per tornare al dominio del tempo, si può semplicemente voler:

TD = ifft(F,NFFT);   %Returns the Inverse of F in Time Domain

Nota che TD restituito sarebbe lungo 256 perché impostiamo NFFT su 256, tuttavia, la lunghezza di x è solo 64, quindi Matlab eseguirà il pad degli zeri fino alla fine della trasformazione TD. Ad esempio, se NFFT era 1024 e la lunghezza era 64, TD restituito sarà 64 + 960 zeri. Si noti inoltre che a causa dell'arrotondamento in virgola mobile, è possibile ottenere qualcosa come 3.1 * 10e-20 ma per scopi generali: per ogni X, ifft (fft (X)) equivale a X fino all'errore roundoff.

Diciamo per un momento che dopo la trasformazione, abbiamo fatto qualcosa e siamo rimasti solo con la parte REALE della FFT:

R = real(F);         %Give the Real Part of the FFT
TDR = ifft(R,NFFT);  %Give the Time Domain of the Real Part of the FFT

Ciò significa che stiamo perdendo la parte immaginaria della nostra FFT e, pertanto, stiamo perdendo informazioni in questo processo inverso. Per preservare l'originale senza perdere informazioni, dovresti sempre mantenere la parte immaginaria della FFT usando imag e applicare le tue funzioni a entrambi o alla parte reale.

figure
subplot(3,1,1)
plot(x);xlabel('time samples');ylabel('magnitude');title('Original Time Domain Signal')
subplot(3,1,2)
plot(TD(1:64));xlabel('time samples');ylabel('magnitude');title('Inverse Fourier Transformed - Time Domain Signal')
subplot(3,1,3)
plot(TDR(1:64));xlabel('time samples');ylabel('magnitude');title('Real part of IFFT transformed Time Domain Signal')

Immagini e FT multidimensionali

Nell'imaging medico, nella spettroscopia, nell'elaborazione delle immagini, nella crittografia e in altre aree della scienza e dell'ingegneria, è spesso il caso che si desideri calcolare trasformazioni multidimensionali di Fourier delle immagini. Questo è abbastanza semplice in Matlab: le immagini (multidimensionali) sono solo matrici n-dimensionali, dopotutto, e le trasformate di Fourier sono operatori lineari: uno solo iterativamente trasforma Fourier lungo altre dimensioni. Matlab fornisce fft2 e ifft2 per farlo in 2-d o fftn in n-dimensioni.

Una potenziale trappola è che la trasformata di Fourier delle immagini viene solitamente mostrata "ordinata centrica", cioè con l'origine dello spazio-k nel mezzo dell'immagine. Matlab fornisce il comando fftshift per scambiare fftshift la posizione dei componenti DC della trasformata di Fourier. Questa notazione di ordinamento rende sostanzialmente più semplice l'esecuzione di tecniche di elaborazione di immagini comuni, una delle quali è illustrata di seguito.

Riempimento a zero

Un modo "veloce e sporco" di interpolare una piccola immagine ad una dimensione più grande è trasformarlo in Fourier, riempire la trasformata di Fourier con zeri e poi prendere la trasformazione inversa. Questo interpola in modo efficace tra ciascun pixel con una funzione di base a forma di sincro, ed è comunemente usato per i dati di imaging medico in bassa risoluzione su scala. Iniziamo caricando un esempio di immagine incorporata

%Load example image 
I=imread('coins.png'); %Load example data -- coins.png is builtin to Matlab 
I=double(I); %Convert to double precision -- imread returns integers 
imageSize = size(I); % I is a 246 x 300 2D image

%Display it
imagesc(I); colormap gray; axis equal; 
%imagesc displays images scaled to maximum intensity

Ora possiamo ottenere la trasformata di Fourier di I. Per illustrare cosa fa il fftshift , confrontiamo i due metodi:

% Fourier transform 
%Obtain the centric- and non-centric ordered Fourier transform of I
k=fftshift(fft2(fftshift(I))); 
kwrong=fft2(I); 

%Just for the sake of comparison, show the magnitude of both transforms: 
figure; subplot(2,1,1); 
imagesc(abs(k),[0 1e4]); colormap gray; axis equal; 
subplot(2,1,2); 
imagesc(abs(kwrong),[0 1e4]); colormap gray; axis equal;
%(The second argument to imagesc sets the colour axis to make the difference clear). 

Ora abbiamo ottenuto il FT 2D di un'immagine di esempio. Per riempirlo a zero, vogliamo prendere ogni k-spazio, riempire i bordi con zeri e quindi prendere la trasformazione posteriore:

%Zero fill
kzf = zeros(imageSize .* 2); %Generate a  492x600 empty array to put the result in 
kzf(end/4:3*end/4-1,end/4:3*end/4-1) = k; %Put k in the middle 
kzfwrong = zeros(imageSize .* 2); %Generate a  492x600 empty array to put the result in 
kzfwrong(end/4:3*end/4-1,end/4:3*end/4-1) = kwrong; %Put k in the middle 

%Show the differences again 
%Just for the sake of comparison, show the magnitude of both transforms: 
figure; subplot(2,1,1); 
imagesc(abs(kzf),[0 1e4]); colormap gray; axis equal; 
subplot(2,1,2); 
imagesc(abs(kzfwrong),[0 1e4]); colormap gray; axis equal;
%(The second argument to imagesc sets the colour axis to make the difference clear). 

A questo punto, il risultato è abbastanza insignificante:

Una volta che prendiamo le back-transforms, possiamo vedere che (correttamente!) I dati a riempimento zero forniscono un metodo ragionevole per l'interpolazione:

% Take the back transform and view 
Izf = fftshift(ifft2(ifftshift(kzf))); 
Izfwrong = ifft2(kzfwrong);

figure; subplot(1,3,1); 
imagesc(abs(Izf)); colormap gray; axis equal; 
title('Zero-filled image'); 
subplot(1,3,2); 
imagesc(abs(Izfwrong)); colormap gray; axis equal;
title('Incorrectly zero-filled image'); 
subplot(1,3,3); 
imagesc(I); colormap gray; axis equal; 
title('Original image'); 
set(gcf,'color','w'); 

Si noti che la dimensione dell'immagine riempita a zero è il doppio dell'originale. Si può riempire a zero di più di un fattore due in ciascuna dimensione, anche se ovviamente farlo non aumenta arbitrariamente la dimensione di un'immagine.

Suggerimenti, consigli, 3D e oltre

L'esempio sopra vale per le immagini 3D (come spesso sono generate dalle tecniche di imaging medico o dalla microscopia confocale, per esempio), ma richiede che fft2 sia sostituito da fftn(I, 3) , per esempio. A causa della natura alquanto scomoda della scrittura di fftshift(fft(fftshift(... più volte, è abbastanza comune definire funzioni come fft2c localmente per fornire una sintassi più semplice a livello locale, come ad esempio:

function y = fft2c(x)

y = fftshift(fft2(fftshift(x)));

Nota che la FFT è veloce, ma le trasformazioni di Fourier multidimensionali di grandi dimensioni richiederanno ancora del tempo su un computer moderno. È inoltre intrinsecamente complesso: la grandezza dello spazio-k è stata mostrata sopra, ma la fase è assolutamente vitale; le traduzioni nel dominio dell'immagine sono equivalenti a una rampa di fase nel dominio di Fourier. Ci sono molte operazioni molto più complesse che si potrebbero desiderare di fare nel dominio di Fourier, come filtrare frequenze spaziali alte o basse (moltiplicandolo con un filtro), o mascherare i punti discreti corrispondenti al rumore. Esiste corrispondentemente una grande quantità di codice generato dalla comunità per la gestione delle operazioni di Fourier comuni disponibili sul sito di repository della comunità principale di Matlab, il File Exchange .