MATLAB Language
Transformaty Fouriera i odwrotne transformaty Fouriera

Składnia

1. Y = fft (X)% oblicza FFT wektora lub macierzy X przy użyciu domyślnej długości transformacji 256 (do potwierdzenia dla wersji)

2. Y = fft (X, n)% oblicza FFT X przy użyciu n jako długości transformacji, n musi być liczbą opartą na 2 potęgach. Jeśli długość X jest mniejsza niż n, wówczas Matlab automatycznie wypełni X zerami, tak że długość (X) = n

3. Y = fft (X, n, dim)% oblicza FFT X za pomocą n jako Długość transformacji wzdłuż wymiaru dim (może wynosić 1 lub 2 dla odpowiednio poziomego lub pionowego)

4. Y = fft2 (X)% Oblicz 2D FFT z X

5. Y = fftn (X, dim)% Oblicz dim-wymiarową FFT X w odniesieniu do wektora wymiarów dim.

6. y = ifft (X)% oblicza odwrotność FFT X (która jest macierzą / wektorem liczb) przy użyciu domyślnej długości transformaty 256

7. y = ifft (X, n)% oblicza IFFT X, używając n jako Długość transformacji

8. y = ifft (X, n, dim)% oblicza IFFT X za pomocą n jako Długość transformacji w stosunku do dim wymiaru (może wynosić 1 lub 2 odpowiednio dla poziomu lub pionu)

9. y = ifft (X, n, dim, 'symetryczny')% Opcja Symetryczny powoduje, że ifft traktuje X jako sprzężoną symetryczną wzdłuż aktywnego wymiaru. Ta opcja jest przydatna, gdy X nie jest dokładnie sprzężony symetrycznie, tylko z powodu błędu zaokrąglenia.

10. y = ifft2 (X)% Oblicz odwrotną stopę 2D X

11. y = ifftn (X, dim)% Oblicz odwrotną fft wymiarową dla X.

Parametry

Uwagi

ifft(fft(X)) = X

Zaimplementuj prostą transformację Fouriera w Matlabie

Transformacja Fouriera to prawdopodobnie pierwsza lekcja cyfrowego przetwarzania sygnałów, jej zastosowanie jest wszędzie i jest potężnym narzędziem, jeśli chodzi o analizę danych (we wszystkich sektorach) lub sygnałów. Matlab ma zestaw potężnych zestawów narzędzi do transformacji Fouriera. W tym przykładzie użyjemy transformacji Fouriera do analizy podstawowego sygnału fali sinusoidalnej i wygenerowania tak zwanego periodiodogramu za pomocą 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

Zauważ, że Dyskretna Transformacja Fouriera jest implementowana przez Szybką Transformację Fouriera (fft) w Matlabie, obie dadzą ten sam wynik, ale FFT jest szybką implementacją 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 ]')

Odwrotne transformaty Fouriera

Jedną z głównych korzyści transformacji Fouriera jest jej zdolność do odwracania się w dziedzinie czasu bez utraty informacji. Rozważmy ten sam sygnał, którego użyliśmy w poprzednim przykładzie:

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

Jeśli otworzymy F w Matlabie, okaże się, że jest to macierz liczb zespolonych, część rzeczywista i część urojona. Z definicji, aby odzyskać oryginalny sygnał w dziedzinie czasu, potrzebujemy zarówno rzeczywistej (która reprezentuje zmianę wielkości), jak i wyobrażonej (reprezentującej zmianę fazy), więc aby powrócić do dziedziny czasu, można po prostu chcieć:

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

Zauważ, że zwrócony TD miałby długość 256, ponieważ ustawiliśmy NFFT na 256, jednak długość x wynosi tylko 64, więc Matlab będzie wstawiał zera do końca transformacji TD. Na przykład, jeśli NFFT wynosi 1024, a długość 64, wówczas zwracane TD będzie wynosić 64 + 960 zer. Zauważ też, że z powodu zaokrąglania zmiennoprzecinkowego możesz otrzymać coś takiego jak 3.1 * 10e-20, ale do ogólnego zastosowania: Dla każdego X, ifft (fft (X)) równa się X w granicach błędu zaokrąglenia.

Powiedzmy przez chwilę, że po transformacji coś zrobiliśmy i pozostała nam PRAWDZIWA część 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

Oznacza to, że tracimy wyimaginowaną część naszego FFT, a zatem tracimy informacje w tym odwrotnym procesie. Aby zachować oryginał bez utraty informacji, zawsze powinieneś zachować wyimaginowaną część FFT za pomocą imag i zastosować swoje funkcje do obu lub do rzeczywistej części.

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')

Obrazy i wielowymiarowe FT

W obrazowaniu medycznym, spektroskopii, przetwarzaniu obrazów, kryptografii i innych obszarach nauki i techniki często zdarza się, że chce się obliczyć wielowymiarowe transformacje Fouriera obrazów. W Matlabie jest to całkiem proste: obrazy (wielowymiarowe) są w końcu matrycami n-wymiarowymi, a transformaty Fouriera są operatorami liniowymi: jeden tylko iteracyjnie transformaty Fouriera wzdłuż innych wymiarów. Matlab zapewnia fft2 i ifft2 aby to zrobić w 2-d lub fftn w n-wymiarach.

Jedną z potencjalnych pułapek jest to, że transformacja Fouriera obrazów jest zwykle pokazywana „uporządkowana centrycznie”, tj. Z początkiem przestrzeni K na środku obrazu. Matlab udostępnia polecenie fftshift aby odpowiednio zamienić położenie komponentów DC transformacji Fouriera. Ta notacja porządkowa znacznie ułatwia wykonywanie typowych technik przetwarzania obrazu, z których jedna jest zilustrowana poniżej.

Zero napełniania

Jednym z „szybkich i brudnych” sposobów interpolacji małego obrazu do większego rozmiaru jest transformacja Fouriera, wypełnienie transformaty Fouriera zerami, a następnie wykonanie transformacji odwrotnej. To skutecznie interpoluje każdy piksel z funkcją podstawową w kształcie cynku i jest powszechnie stosowane do skalowania danych obrazowania medycznego o niskiej rozdzielczości. Zacznijmy od załadowania przykładowego obrazu

%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

Możemy teraz uzyskać transformatę Fouriera I. Aby zilustrować działanie fftshift , porównajmy dwie metody:

% 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). 

Otrzymaliśmy teraz 2D FT przykładowego obrazu. Aby wypełnić go zerami, chcemy wziąć każdą przestrzeń k, uzupełnić krawędzie zerami, a następnie wziąć transformację wsteczną:

%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). 

W tym momencie wynik jest dość niezwykły:

Kiedy weźmiemy transformacje wsteczne, możemy zobaczyć, że (poprawnie!) Dane wypełniające zero stanowią sensowną metodę interpolacji:

% 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'); 

Pamiętaj, że rozmiar obrazu wypełnionego zerą jest dwa razy większy niż oryginału. Można wypełnić zero więcej niż dwa razy w każdym wymiarze, chociaż oczywiście nie powoduje to arbitralnego zwiększenia rozmiaru obrazu.

Wskazówki, porady, 3D i nie tylko

Powyższy przykład odnosi się do obrazów 3D (jak często są generowane za pomocą technik obrazowania medycznego lub mikroskopu konfokalnego, na przykład), ale wymagają fft2 być zastąpione przez fftn(I, 3) , na przykład. Ze względu na nieco nieporęczną naturę pisania fftshift(fft(fftshift(... kilka razy) dość często lokalnie definiuje się funkcje takie jak fft2c aby zapewnić lokalnie łatwiejszą składnię - takie jak:

function y = fft2c(x)

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

Należy pamiętać, że FFT jest szybki, ale duże, wielowymiarowe transformacje Fouriera nadal będą wymagały czasu na nowoczesnym komputerze. Jest to dodatkowo z natury skomplikowane: wielkość przestrzeni k pokazano powyżej, ale faza jest absolutnie niezbędna; tłumaczenia w dziedzinie obrazów są równoważne rampie fazowej w dziedzinie Fouriera. Istnieje szereg znacznie bardziej złożonych operacji, które można chcieć wykonywać w dziedzinie Fouriera, takich jak filtrowanie wysokich lub niskich częstotliwości przestrzennych (przez pomnożenie ich przez filtr) lub maskowanie dyskretnych punktów odpowiadających szumowi. Odpowiednio duża ilość kodu generowanego przez społeczność do obsługi typowych operacji Fouriera jest dostępna w głównej witrynie repozytorium społeczności Matlaba, File Exchange .