function [xa,ya,z]=kdfft2(A,k,inc,H,g);

% KDFFT2      KDE plot for 2D data using FFTs.
%
%             Example: load kddata
%                      h1=hste(mask4(:,1),'knorm');
%                      h2=hste(mask4(:,2),'knorm');
%                      kdfft2(mask4,'knorm',128,[h1,h2]);
%
%             KDFFT2(A,K,INC,H,G) plots the Kernel Density Estimate with 
%             Normal kernel.  
%             A: The data is stored in the rows of A.  For 2D data A 
%                will have two columns.
%             K: The kernel function, one of 'KNORM','KEPAN','KEPV2','KEPV3'. 
%             INC: The number of increments to use in plotting the 
%                  density function.  This should be a power of 2, for example 
%                  32, 64 and 128 give good results.  
%             H: Window width.  If H=0, an 'optimal' value is used.
%                If a single value of H is used then smoothing is the same 
%                in both directions.
%                If H=[H1, H2] smoothing is by H1 in the x and H2 in the
%                y directions.
%             G: Optional, G=0 turns graphics off.
%
%             Christian C. Beardah 1994

if k=='knorm',
  r=4;
  else,
    r=1;
end;

P=2*inc;

if nargin == 4,
  g=1;
end;

if length(H)==1,
  h1=H;
  h2=H;
  elseif length(H)>=2,
    h1=H(1);
    h2=H(2);
end;

[n d]=size(A);    % Establish dimensions of A,
                  % n=number of data points.

x=A(:,1);         % Extract the first column of A into x.
xmin=min(x);      % Find the minimum value of x.
xmax=max(x);      % Find the maximum value of x.
xrange=xmax-xmin; % Find the range of x.

% xa holds the x 'axis' vector, defining a grid of x values where 
% the k.d. function will be evaluated and plotted.

ax=xmin-xrange/4;
bx=xmax+xrange/4;

xa=linspace(ax,bx,inc); 

y=A(:,2);         % Extract the second column of A into y.
ymin=min(y);
ymax=max(y);
yrange=ymax-ymin;

% Establish an 'axis' vector for the y values.

ay=ymin-yrange/4;
by=ymax+yrange/4;

ya=linspace(ay,by,inc);

% If H=0 was specified, use the normal scale rule to find h1, h2.

if H==0,
  h1=hns(x,k);
  h2=hns(y,k);
end;

% Find the binned kernel weights, c.

c=zeros(inc,inc);

deltax=(bx-ax)/(inc-1);
binx=floor((x-ax)/deltax)+1;

deltay=(by-ay)/(inc-1);
biny=floor((y-ay)/deltay)+1;

for i=1:n, % Loop over data points.
  w=(deltax*deltay);

  weight11=(xa(binx(i)+1)-x(i))*(ya(biny(i)+1)-y(i))/w;
  weight12=(x(i)-xa(binx(i)))*(ya(biny(i)+1)-y(i))/w;  
  weight21=(xa(binx(i)+1)-x(i))*(y(i)-ya(biny(i)))/w;
  weight22=(x(i)-xa(binx(i)))*(y(i)-ya(biny(i)))/w;  

  c(binx(i),biny(i))=c(binx(i),biny(i))+weight11;
  c(binx(i)+1,biny(i))=c(binx(i)+1,biny(i))+weight12;
  c(binx(i),biny(i)+1)=c(binx(i),biny(i)+1)+weight21;
  c(binx(i)+1,biny(i)+1)=c(binx(i)+1,biny(i)+1)+weight22;
end;

c=c';

% Obtain the kernel weights.

L1=max(floor(r*h1*(inc-1)/(bx-ax)),floor(r*h2*(inc-1)/(by-ay)));

L=min(L1,inc-1);

kw=zeros(2*inc,2*inc);

[X1,Y1]=meshgrid((bx-ax)*[-L:L]/((inc-1)*h1),(by-ay)*[-L:L]/((inc-1)*h2));

% kw((inc-L+1):(inc+L+1),(inc-L+1):(inc+L+1))=feval(k,X1,Y1)/(n*h1*h2);

if strcmp(k,'knorm')==1;
   kw((inc-L+1):(inc+L+1),(inc-L+1):(inc+L+1)) = knorm(X1,Y1)/(n*h1*h2);
elseif strcmp(k,'kepan')==1;
   kw((inc-L+1):(inc+L+1),(inc-L+1):(inc+L+1)) = kepan(X1,Y1)/(n*h1*h2);
elseif strcmp(k,'krect')==1;
   kw((inc-L+1):(inc+L+1),(inc-L+1):(inc+L+1)) = krect(X1,Y1)/(n*h1*h2);
elseif strcmp(k,'ktria')==1;
   kw((inc-L+1):(inc+L+1),(inc-L+1):(inc+L+1)) = ktria(X1,Y1)/(n*h1*h2);
elseif strcmp(k,'kbiwe')==1;
   kw((inc-L+1):(inc+L+1),(inc-L+1):(inc+L+1)) = kbiwe(X1,Y1)/(n*h1*h2);
elseif strcmp(k,'ktriw')==1;
   kw((inc-L+1):(inc+L+1),(inc-L+1):(inc+L+1)) = ktriw(X1,Y1)/(n*h1*h2);
end;


% Apply 'fftshift' to kw.

kw=fftshift(kw);

% Perform the convolution.

z=real(ifft2(fft2(c,P,P).*fft2(kw)));

z=z(1:inc,1:inc).*(z(1:inc,1:inc)>0);

if g~=0,
%  mesh(xa,ya,z);
  surf(xa,ya,z);
  title(['h = [',num2str(h1),',',num2str(h2),']']);   
end;
