% numerical verification of the Wigner-Ville transform of
% rectangular, triangular, and Gaussian functions
% given in:
% "Influence of the analysis window on the metrological 
% performance of the grid method"
% Frederic Sur and Michel Grediac
% Journal of Mathematical Imaging and Vision
% vol. 56, no. 3, p. 472-498
% 2016.

clear all, close all;

fprintf('\n');
fprintf('numerical verification of the Wigner-Ville transform\n')
fprintf('of rectangular, triangular, and Gaussian functions.\n')
fprintf('cf "Influence of the analysis window on the metrological performance of the grid method"\n')
fprintf('Frederic Sur and Michel Grediac\n')
fprintf('Journal of Mathematical Imaging and Vision\n')
fprintf('vol. 56, no. 3, p. 472-498, 2016\n')
fprintf('\n')

% window half-width, randomly chosen (test also with an integer value)
% (this script is to be run several times)
b=10*rand(1);
fprintf(['half-width: b=',num2str(b),'\n']);

% functions:
rec = @(x) (abs(x)<=b)/2/b;
tri = @(x) (b-abs(x))/b^2.*(abs(x)<=b);
gau = @(x) 1/(sqrt(2*pi)*b)*exp(-x.^2/(2*b^2));

% variables x and lambda denote the space and frequency variables
% in the Wigner-Ville transform
% xstep and lambdastep can be set by the user
xstep=4*b/100;
lambdastep=4*b/100;
xrange=-2*b:xstep:2*b;
lambdarange=-1.5*b:lambdastep:1.5*b;


% rectangular function
%%%%%%%%%%%%%%%%%%%%%%

fprintf('\n');
fprintf('numerical estimation of the Wigner-Ville transform of the rectangular function\n');
WVnum=zeros(numel(xrange),numel(lambdarange));
WVth=zeros(numel(xrange),numel(lambdarange));
reverseStr = '';

% subintervals for numerical integration:
% numerical accuracy depends on istep
istep=.001;
irange=-2*b:istep:2*b;


for jx=1:numel(xrange),
    
    percent = 100*jx/numel(xrange);
    msg = sprintf('Percent done: %3.1f', percent); 
    fprintf([reverseStr, msg]);
    reverseStr=repmat(sprintf('\b'), 1, length(msg));
   
    for jlambda=1:numel(lambdarange),
        x=xrange(jx);
        lambda=lambdarange(jlambda);
       
        WVnum(jx,jlambda) = istep* ...
            sum(rec(x+irange/2).*rec(x-irange/2).*exp(-1i*irange*lambda)) ;
        if (lambda~=0) 
          WVth(jx,jlambda) = 1/(2*b^2*lambda) * sin(2*lambda*(b-abs(x))) * (abs(x)<=b);
        else
          WVth(jx,jlambda) = (b-abs(x))/b^2 * (abs(x)<=b);
        end;
        
    end;
    
end;

fprintf('\n');
fprintf(['relative difference between numerical and theoretical value: ',num2str(max(abs(WVnum(:)-WVth(:)))/max(WVth(:)-min(WVth(:)))),'\n']);
fprintf(['maximum imaginary part of the numerical value (theoretically real): ',num2str(max(abs(imag(WVnum(:))))),'\n']);
WVnum=real(WVnum);
figure, imagesc(lambdarange,xrange,WVnum), colorbar; title('WF(rectangle): numerical estimation')
xlabel('\lambda'), ylabel('x')
figure, imagesc(lambdarange,xrange,WVth), colorbar; title('WF(rectangle): theoretical expression')
xlabel('\lambda'), ylabel('x')
figure, imagesc(lambdarange,xrange,WVnum-WVth), colorbar; title('WF(rectangle): difference between th. and num.');   
xlabel('\lambda'), ylabel('x')

figure, surf(ones(length(xrange),1)*lambdarange, xrange'*ones(1,length(lambdarange)), WVth)
xlabel('\lambda'), ylabel('x')
axis([lambdarange(1),lambdarange(end),xrange(1),xrange(end),])
view(-37.5,35)

drawnow;


% triangular function
%%%%%%%%%%%%%%%%%%%%%

fprintf('\n');
fprintf('numerical estimation of the Wigner-Ville transform of the triangular function\n');
WVnum=zeros(numel(xrange),numel(lambdarange));
WVth=zeros(numel(xrange),numel(lambdarange));
reverseStr = '';

% subintervals for numerical integration:
% numerical accuracy depends on istep
istep=.001;
irange=-2*b:istep:2*b;

for jx=1:numel(xrange),
    
    percent = 100*jx/numel(xrange);
    msg = sprintf('Percent done: %3.1f', percent); 
    fprintf([reverseStr, msg]);
    reverseStr=repmat(sprintf('\b'), 1, length(msg));
   
    for jlambda=1:numel(lambdarange),
        x=xrange(jx);
        lambda=lambdarange(jlambda);
        WVnum(jx,jlambda) = istep* ...
            sum(tri(x+irange/2).*tri(x-irange/2).*exp(-1i*irange*lambda)) ;
        
        if (lambda~=0)
            WVth(jx,jlambda) = 2/b^4/lambda^2 * ...
                  ( -(b-abs(x))*cos(2*lambda*(b-abs(x))) + ...
                  1/(2*lambda)*sin(2*lambda*(b-abs(x))) ) * (abs(x)<=b);
        else
            WVth(jx,jlambda) = 8*(b-abs(x))^3/(3*b^4)*(abs(x)<=b);
        end;
        
        if (abs(x)<=b/2)
            if (lambda~=0)
                 WVth(jx,jlambda) = WVth(jx,jlambda) + 2/b^4/lambda^2 * 2*cos(lambda*b) * ...
                     ( (b-2*abs(x))*cos(lambda*(b-2*abs(x))) - 1/lambda*sin(lambda*(b-2*abs(x))) );
            else
                WVth(jx,jlambda) = WVth(jx,jlambda) - 4*(b-2*abs(x))^3/(3*b^4); 
            end;     
        end;
        
    end;
    
end;

fprintf('\n');
fprintf(['relative difference between numerical and theoretical value: ',num2str(max(abs(WVnum(:)-WVth(:)))/max(WVth(:)-min(WVth(:)))),'\n']);
fprintf(['maximum imaginary part of the numerical value (theoretically real): ',num2str(max(abs(imag(WVnum(:))))),'\n']);
WVnum=real(WVnum);
figure, imagesc(lambdarange,xrange,WVnum), colorbar; title('WF(triangle): numerical estimation')
xlabel('\lambda'), ylabel('x')
figure, imagesc(lambdarange,xrange,WVth), colorbar; title('WF(triangle): theoretical expression')
xlabel('\lambda'), ylabel('x')
figure, imagesc(lambdarange,xrange,WVnum-WVth), colorbar; title('WF(triangle): difference between th. and num.');   
xlabel('\lambda'), ylabel('x')

figure, surf(ones(length(xrange),1)*lambdarange, xrange'*ones(1,length(lambdarange)), WVth)
xlabel('\lambda'), ylabel('x')
axis([lambdarange(1),lambdarange(end),xrange(1),xrange(end),])
view(-37.5,35)

drawnow;

% Gaussian function
%%%%%%%%%%%%%%%%%%%%%%

fprintf('\n');
fprintf('numerical estimation of the Wigner-Ville transform of the Gaussian function\n');
WVnum=zeros(numel(xrange),numel(lambdarange));
WVth=zeros(numel(xrange),numel(lambdarange));
reverseStr = '';

% subintervals for numerical integration:
% numerical accuracy depends on istep
istep=.001;
irange=-7*b:istep:7*b;

for jx=1:numel(xrange),
    
    percent = 100*jx/numel(xrange);
    msg = sprintf('Percent done: %3.1f', percent); 
    fprintf([reverseStr, msg]);
    reverseStr=repmat(sprintf('\b'), 1, length(msg));
   
    for jlambda=1:numel(lambdarange),
        x=xrange(jx);
        lambda=lambdarange(jlambda);
       
        WVnum(jx,jlambda) = istep* ...
            sum(gau(x+irange/2).*gau(x-irange/2).*exp(-1i*irange*lambda)) ;
          WVth(jx,jlambda) = 1/(sqrt(pi)*b)*exp(-x^2/b^2-b^2*lambda^2);
        
    end;
    
end;

fprintf('\n');
fprintf(['relative difference between numerical and theoretical value: ',num2str(max(abs(WVnum(:)-WVth(:)))/max(WVth(:)-min(WVth(:)))),'\n']);
fprintf(['maximum imaginary part of the numerical value (theoretically real): ',num2str(max(abs(imag(WVnum(:))))),'\n']);
WVnum=real(WVnum);
figure, imagesc(lambdarange,xrange,WVnum), colorbar; title('WF(Gaussian): numerical estimation')
xlabel('\lambda'), ylabel('x')
figure, imagesc(lambdarange,xrange,WVth), colorbar; title('WF(Gaussian): theoretical expression')
xlabel('\lambda'), ylabel('x')
figure, imagesc(lambdarange,xrange,WVnum-WVth), colorbar; title('WF(Gaussian): difference between th. and num.');   
xlabel('\lambda'), ylabel('x')

figure, surf(ones(length(xrange),1)*lambdarange, xrange'*ones(1,length(lambdarange)), WVth)
xlabel('\lambda'), ylabel('x')
axis([lambdarange(1),lambdarange(end),xrange(1),xrange(end),])
view(-37.5,35)

drawnow;
