function W = wigner2D(Ex)
% W(y,x,v,u) = wigner2D(Ein(y,x))
% 
% W(x,y,u,v): output WDF of 4 variables
% (x,y): spatial coordinates     x-> row direction, y-> column direction
% (u,v): spatial frequency coordinates
% Ein: input signal of 2 variables
%
%
%	Notes of wigner.m:
%		W = Int(-inf..inf){E(x+y)E(x-y)exp[2ixy]}
%
%		E(x+y) & E(x-y) are calculated via a FFT (fast Fourier transform) using the
%		shift theorem. The integration is performed via a FFT. Thus it is important
%		for the data to satisfy the sampling theorem:
%		dy = 2*pi/X			X = span of all x-values	dy = y resolution
%		dx = 2*pi/Y			Y = span of all y-values	dx = x resolution
%		The data must be completely contained within the range x(0)..x(N-1) &
%		y(0)..y(N-1) (i.e. the function must fall to zero within this range).
%
%	v1.0
%
%	Currently waiting for update:
%		Remove the fft/ifft by performing this inside the last function calls
%		Allow an arbitrary output resolution
%		Allow an input vector for x (and possibly y).


% generate W template
W = repmat(0, [size(Ex,1) size(Ex,2) size(Ex,1) size(Ex,2)]);

Ny = size(Ex, 1);
Nx = size(Ex, 2);

%   Generate linear vector
xvec = [-ceil((Nx-1)/2):Nx-1-ceil((Nx-1)/2)]*2*pi/(Nx-1);
yvec = [-ceil((Ny-1)/2):Ny-1-ceil((Ny-1)/2)]*2*pi/(Ny-1);

Xvec = [-ceil((Nx-1)/2):Nx-1-ceil((Nx-1)/2)];
Yvec = [-ceil((Ny-1)/2):Ny-1-ceil((Ny-1)/2)];

%[y,x,Y,X] = ndgrid(yvec, xvec, Yvec, Xvec);
[x,y,X,Y] = ndgrid(xvec, yvec, Xvec, Yvec);

clear xvec yvec Xvec Yvec;

% repliate the fft of input signal
EX = zeros(size(W));
tmp = fftshift(fft2(ifftshift(Ex)));
for n=1:size(EX,1),
    for m=1:size(EX,2),
        EX(n,m,:,:) = tmp(n,m);
    end
end


% computing the shift correlation like term
EX1 = ifftshift(ifftshift(...
    ifft(ifft(...
    fftshift(fftshift(...
    EX.*exp(i*((x.*X/2)+(y.*Y/2))),...
        1), 2),...
    [], 1), [], 2),...
    1), 2);

EX2 = ifftshift(ifftshift(...
    ifft(ifft(...
    fftshift(fftshift(...
    EX.*exp(-i*((x.*X/2)+(y.*Y/2))),...
        1), 2),...
    [], 1), [], 2),...
    1), 2);

% computing the WDF
W = fftshift(fftshift(...
    fft(fft(...
    ifftshift(ifftshift(EX1.*conj(EX2), 3), 4),...
    [], 3), [], 4),...
    3 ), 4);


% for numerical reason, just take real values
W = real(W);