function [z,x,y] = inversefouriertransform2(Cz,varargin)
%INVERSEFOURIERTRANSFORM2 Performs the Real 2D Inverse Fourier Transform.
% [Z,X,Y] = INVERSEFOURIERTRANSFORM2(RC,IC,dFx,dFy) Gives the real 2D
% field Z(X,Y) from the components of the complex 2D Fourier transform,
% Cz(Fx,Fy) = RC(Fx,Fy) + i*IC(Fx,Fy) = complex(RC,IC)
% taken at dFx,dFy sampling frequency interval, via IFFT2. The spatial 2D
% frequencies (Fx,Fy) must have the origin at the center and Fy must
% increase upwards contrary to the row index.
%
% [Z,X,Y] = INVERSEFOURIERTRANSFORM2(Cz,dFx,dFy) does the same thing.
%
% [Z,X,Y] = INVERSEFOURIERTRANSFORM2(Cz,dFx,dFy,'ifft') does the same
% thing but via IFFT. Only for academic purpose.
%
% [Z,X,Y] = INVERSEFOURIERTRANSFORM2(Cz,dFx,dFy,'exp') does the same
% thing but via complex exponentials. Only for academic purpose. Not so
% fast as IFFT or IFFT2.
%
% Inputs: units: Size:
% RC - Real components [u*v*w] Ny times Nx
% IC - Imaginary components [u*v*w] Ny times Nx
% (Cz - Complex transform [u*v*w] Ny times Nx)
% dFx - Sampling frequency interval at Fx [cycle/v] Scalar
% dFy - Sampling frequency interval at Fy [cycle/w] Scalar
%
% Outputs:
% Z - 2D real field [u] Ny times Nx
% X - x grid origin left-lower corner [v] Ny times Nx
% Y - y grid origin left-lower corner [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) must be
% 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 dFx = 1 => normalized space at x (Bloomfield)
% if dFy = 1 => normalized space at y (Bloomfield)
% if dFx = 1/Nx or empty => normalized frequency at x (MATLAB)
% if dFy = 1/Ny 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);
% [zi,xi,yi] = inversefouriertransform2(rc,ic,1/(Nx*dX),1/(Ny*dY));
% subplot(221), surf(x,y,z), view(2), shading interp
% xlabel('x, m'), ylabel('y, m'), axis([0 (Nx-1)*dX 0 (Ny-1)*dY])
% title('Original 2D field')
% subplot(223), surf(xi,yi,zi), view(2), shading interp
% xlabel('x, m'), ylabel('y, m'), axis([0 (Nx-1)*dX 0 (Ny-1)*dY])
% title('Result 2D field')
% subplot(222), 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('frequency, cycle/m'), ylabel('Frequency, cycle/m')
% title('Real part')
% subplot(224), 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('frequency, cycle/m'), ylabel('Frequency, cycle/m')
% title('Imaginary part')
%
% See also FOURIERTRANSFORM2, INVERSEFOURIERTRANSFORM, IFFT2
% Written by
% Lic. on Physics Carlos Adrin Vargas Aguilera
% Physical Oceanography MS candidate
% UNIVERSIDAD DE GUADALAJARA
% Mexico, 2004
%
% nubeobscura@hotmail.com
% Time series information:
[Cz,dFx,dFy,Nx,Ny,Method] = check_arguments(Cz,varargin,nargin);
% Complex Fourier transform:
switch lower(Method)
case 'ifft2'
z = inversefouriertransform_ifft2(Cz,dFx*dFy,Nx*Ny);
case 'ifft'
z = inversefouriertransform_ifft(Cz,dFx,dFy,Nx,Ny);
case 'exp'
z = inversefouriertransform_exponential(Cz,dFx*dFy,Nx,Ny);
otherwise
error('Method unknown. Must be one of ''ifft2'', ''ifft'' or ''exp''.')
end
% Space grid (y increase upwards and origin at the lower-left corner):
x = (0:Nx-1)/(Nx*dFx); y = fliplr(0:Ny-1)/(Ny*dFy);
[x,y] = meshgrid(x,y);
function z = inversefouriertransform_ifft2(Cz,dFs,Ns)
% 2D real inverse Fourier transform via IFFT2.
Cz = flipud(Cz); Cz = ifftshift(Cz);
z = ifft2(Cz)*Ns*dFs; z = flipud(z); z = real(z);
function z = inversefouriertransform_ifft(Cz,dFx,dFy,Nx,Ny)
% 2D real inverse Fourier transform via IFFT.
Cz = flipud(Cz); Cz = ifftshift(Cz);
Cy = ifft(Cz,[],1)*Ny*dFy;
z = ifft(Cy.',[],1).'*Nx*dFx; z = flipud(z); z = real(z);
function z = inversefouriertransform_exponential(Cz,dFs,Nx,Ny)
% 2D real inverse 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*dX,yn*dY): space origin left-lower corner
fxm = mx/Nx; % fxm/dX: Fourier frequency at x
fym = my/Ny; % fym/dY: Fourier frequency at y
[fxm,fym] = meshgrid(fxm,fym); % Centered frequency grid
kxm = 2*pi*fxm; % Fourier wavenumber x-component
kym = 2*pi*fym; % Fourier wavenumber y-component
z = zeros(Ny,Nx);
for ny = 1:length(yn)
for nx = 1:length(xn)
z(ny,nx) = sum( sum( Cz .* exp(i*(kxm*xn(nx) + kym*yn(ny))) )) * dFs;
end
end
z = real(z);
function [Cz,dFx,dFy,Nx,Ny,Method] = check_arguments(Cz,Ventries,Nentries)
% Is matrix?
[Ny,Nx] = size(Cz);
% Sampling interval?
ic = []; % Default (complex entry)
dFx = []; % 1/Nx: Default (MATLAB way)
dFy = []; % 1/Ny: Default (MATLAB way)
Method = 'ifft2'; % Default (MATLAB way)
if Nentries > 1
if ischar(Ventries{1})
Method = Ventries{1};
elseif numel(Ventries{1}) == Nx*Ny
ic = Ventries{1};
else
dFx = Ventries{1}; dFy = dFx;
end
end
if Nentries > 2
if ischar(Ventries{2})
Method = Ventries{2};
elseif numel(Ventries{2}) == Nx*Ny
ic = Ventries{2};
elseif isempty(dFx)
dFx = Ventries{2}; dFy = dFx;
else
dFy = Ventries{2};
end
end
if Nentries > 3
if ischar(Ventries{3})
Method = Ventries{3};
elseif numel(Ventries{3}) == Nx*Ny
ic = Ventries{3};
elseif isempty(dFx)
dFx = Ventries{3}; dFy = dFx;
else
dFy = Ventries{3};
end
end
if Nentries > 4
if ischar(Ventries{4})
Method = Ventries{4};
elseif numel(Ventries{4}) == Nx*Ny
ic = Ventries{4};
else
dFy = Ventries{4};
end
end
if ~isempty(ic)
Cz = complex(Cz,ic);
end
% Default sampling frequency interval?
if isempty(dFx), dFx = 1/Nx; dFy = 1/Ny; end
% Carlos Adrin. nubeobscura@hotmail.com
