function [rc,ic,Fx,Fy] = fouriertransform2(z,varargin)
%FOURIERTRANSFORM2 Performs the complex Fourier transform of a matrix.
% [RC,IC,Fx,Fy] = FOURIERTRANSFORM2(Z,dX,dY) Gives the components of the
% complex Fourier transform
% Cz(Fx,Fy) = RC(Fx,Fy) + i*IC(Fx,Fy) = complex(RC,IC)
% of a 2D field Z(X,Y) taken at dX,dY sampling interval, via
% FFT2. The spatial 2D frequencies Fx,Fy have the origin at the center.
% Fy (and Y) increase upwards contrary to the row index.
%
% [RC,IC,Fx,Fy] = FOURIERTRANSFORM2(Z,dX,dY,'fft') does the same thing
% but via fft. Only for academic purpose.
%
% [RC,IC,Fx,Fy] = FOURIERTRANSFORM2(Z,dX,dY,'exp') does the same thing
% but via complex exponentials. Only for academic purpose. Not so fast.
%
% Inputs: units: Size
% Z - 2D real field, for example units [u] Ny times Nx
% dX - Sampling spatial interval at x [v] Scalar
% dY - Sampling spatial interval at y [w] Scalar
%
% Outputs:
% RC - Real components [u*v*w] Ny times Nx
% IC - Imaginary components [u*v*w] Ny times Nx
% Fx - Centered x-frequency grid [cycle/v] Ny times Nx
% Fy - Centered y-frequency grid [cycle/w] Ny times Nx
%
% Note: the origin of (X,Y) is at the (Ny,1) element (left-bottom
% corner), and of (Fx,Fy) at the (Ny-floor(Ny/2),floor(Nx/2)+1) element
% (center).
%
% Note: the components at null frequencies are:
% RC(Ny-floor(Ny/2),floor(Nx/2)+1) = mean(Z(:))*(Tx*Ty);
% IC(Ny-floor(Ny/2),floor(Nx/2)+1) = 0;
% where Tx=Nx*dX and Ty=Ny*dY are the total sampling interval at x and y.
%
% Note: the Fourier spatial frequencies (cycle/v, cycle/w) are
% nx = -floor(Nx/2):(Nx-1)/2; % centralized x index
% ny = -floor(Ny/2):(Ny-1)/2; % centralized y index
% ny = fliplr(ny); % fy increase upwards
% Fx = nx/(Nx*dX); Fy = ny/(Ny*dY);
% [Fx,Fy] = meshgrid(Fx,Fy);
%
% Note: if dX = 1/Nx => normalized space at x (Bloomfield)
% if dY = 1/Ny => normalized space at y (Bloomfield)
% if dX = 1 or empty => normalized frequency at x (MATLAB)
% if dY = 1 or empty => normalized frequency at y (MATLAB)
%
% Example:
% Nx = 100; Ny = 200; dX = 120; dY = dX;
% fx1 = 0.0005; fx2 = -0.0010;
% fy1 = 0.0012; fy2 = 0.0017;
% x = (0:Nx-1)*dX; y = fliplr(0:Ny-1)*dY; [x,y] = meshgrid(x,y);
% z = 20*sin(2*pi*(fx1*x + fy1*y)) + 30*cos(2*pi*(fx2*x + fy2*y));
% z = z + rand(size(x));
% z = z-mean(z(:));
% [rc,ic,Fx,Fy] = fouriertransform2(z,dX,dY);
% subplot(311), surf(x,y,z), view(2), shading interp
% xlabel('x, m'), ylabel('y, m'), axis([0 (Nx-1)*dX 0 (Ny-1)*dY])
% title('Complex Fourier transform 2D example')
% subplot(312), surf(Fx,Fy,rc), view(2), shading interp
% hold on, axis([Fx(1,1) Fx(1,end) Fy(end,1) Fy(1,1)])
% plot3(fx1*[1 -1],fy1*[1 -1],max(rc(:))*[1 1],'ro',...
% fx2*[1 -1],fy2*[1 -1],max(rc(:))*[1 1],'bo',...
% [Fx(1,1) Fx(1,end)],[0 0],max(rc(:))*[1 1],'-k',...
% [0 0],[Fy(end,1) Fy(1,1)],max(rc(:))*[1 1],'-k')
% hold off
% xlabel('frecuency, cycle/m'), ylabel('Frequency, cycle/m')
% title('Real part')
% subplot(313), surf(Fx,Fy,ic), view(2), shading interp
% hold on, axis([Fx(1,1) Fx(1,end) Fy(end,1) Fy(1,1)])
% plot3(fx1*[1 -1],fy1*[1 -1],max(ic(:))*[1 1],'ro',...
% fx2*[1 -1],fy2*[1 -1],max(ic(:))*[1 1],'bo',...
% [Fx(1,1) Fx(1,end)],[0 0],max(ic(:))*[1 1],'-k',...
% [0 0],[Fy(end,1) Fy(1,1)],max(ic(:))*[1 1],'-k')
% hold off
% xlabel('frecuency, cycle/m'), ylabel('Frequency, cycle/m')
% title('Imaginary part')
%
% See also INVERSEFOURIERTRANSFORM2, FOURIERTRANSFORM, FFT2
% Written by
% Lic. on Physics Carlos Adrin Vargas Aguilera
% Physical Oceanography MS candidate
% UNIVERSIDAD DE GUADALAJARA
% Mexico, 2004
%
% nubeobscura@hotmail.com
% Time series information:
[dX,dY,Nx,Ny,Method] = check_arguments(z,varargin,nargin);
% Complex Fourier transform:
switch lower(Method)
case 'fft2'
[rc,ic] = fouriertransform_fft2(z,dX*dY);
case 'fft'
[rc,ic] = fouriertransform_fft(z,dX,dY);
case 'exp'
[rc,ic] = fouriertransform_exponential(z,dX*dY,Nx,Ny);
otherwise
error('Method unknown. Must be one of ''fft2'', ''fft'' or ''exp''.')
end
% Fourier frequency fields:
mx = -floor(Nx/2):(Nx-1)/2; % centered index
my = -floor(Ny/2):(Ny-1)/2; % centered index
my = fliplr(my); % increase upwards
Fx = mx/(Nx*dX); % Fourier frequency at x, centered
Fy = my/(Ny*dY); % Fourier frequency at y, centered
[Fx,Fy] = meshgrid(Fx,Fy);
function [rc,ic] = fouriertransform_fft2(z,dS)
% 2D complex Fourier Transform via FFT2.
Cz = fft2(flipud(z))*dS; Cz = fftshift(Cz); Cz = flipud(Cz);
rc = real(Cz); ic = imag(Cz);
function [rc,ic] = fouriertransform_fft(z,dX,dY)
% 2D complex Fourier Transform via FFT.
Cy = fft(flipud(z),[],1)*dY;
Cz = fft(Cy.',[],1).'*dX;
Cz = fftshift(Cz); Cz = flipud(Cz);
rc = real(Cz); ic = imag(Cz);
function [rc,ic] = fouriertransform_exponential(z,dS,Nx,Ny)
% 2D complex Fourier Transform via complex exponentials.
nx = 0:Nx-1; mx = -floor(Nx/2):(Nx-1)/2;
ny = 0:Ny-1; my = -floor(Ny/2):(Ny-1)/2;
ny = fliplr(ny); my = fliplr(my);
xn = nx; yn = ny;
[xn,yn] = meshgrid(xn,yn); % (xn*dX,yn*dY): space origin left-lower corner
fxm = mx/Nx; % fxm/dX: Fourier frequency at x centered
fym = my/Ny; % fym/dY: Fourier frequency at y centered
kxm = 2*pi*fxm; % Fourier wavenumber x-component
kym = 2*pi*fym; % Fourier wavenumber y-component
Cz = zeros(Ny,Nx);
for my = 1:length(kym)
for mx = 1:length(kxm)
Cz(my,mx) = sum( sum( z .* exp(-i*(kxm(mx)*xn + kym(my)*yn)) )) * dS;
end
end
rc = real(Cz); ic = imag(Cz);
function [dX,dY,Nx,Ny,Method] = check_arguments(z,Ventries,Nentries)
% Is matrix?
[Ny,Nx] = size(z);
% Sampling interval?
dX = [];
dY = 1; % Default (MATLAB way)
Method = 'fft2'; % Default (MATLAB way)
if Nentries > 1
if ~ischar(Ventries{1})
dX = Ventries{1}; dY = dX;
else
Method = Ventries{1};
end
end
if Nentries > 2
if ~ischar(Ventries{2})
if isempty(dX)
dX = Ventries{2}; dY = dX;
else
dY = Ventries{2};
end
else
Method = Ventries{2};
end
end
if Nentries > 3
if ~ischar(Ventries{3})
dY = Ventries{3};
else
Method = Ventries{3};
end
end
if isempty(dX)
dX = dY;
end
% Carlos Adrin. nubeobscura@hotmail.com
