function [xa,y]=kdfft1(A,k,inc,h,g,xr);

% KDFFT1      KDE plot for 1D data using FFTs.
%
%             Example: load kddata
%                      [x,y]=kdfft1(suicide,'knorm',128,40);
%
%             KDFFT1(A,K,INC,H,G,XR) plots the Kernel Density Estimate with 
%             Normal kernel.  
%             A: The data is stored in the rows of A.  For 1D data A 
%                will have one column.
%             K: The kernel function, one of 'KNORM', 'KEPAN', etc..
%             INC: The number of increments to use in plotting the 
%                  density function.  THIS MUST BE A POWER OF 2, for example 
%                  32, 64 and 128 give good results.  
%             H: Window width.  If omitted, an 'optimal' value (calculated using
%                the STE rule) is used.
%             G: Optional, G=0 turns graphics off.
%             XR: Optional, XR=[XMIN, XMAX] specifies the x range - XR must 
%                 include the range of the data.
%
%             Christian C. Beardah 1994

P=2*inc;

r=1;

if k=='knorm' | k=='cnorm' | k=='clapl' | k=='clogi',
  r=4;
end;
if k=='klogi' | k=='klapl',
  r=7;
end;
if k=='cbiwe' | k=='ctriw',
  r=3;
end;
if k=='cepan' | k=='ctria' | k=='crect',
  r=2;
end;

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

if nargin == 4 | nargin==5,
  xr=mean(A)*[1 1];
end;

if h==0,
%  h=hns(A,k);
  h=hste(A,k);
end;

n=length(A);      % Establish dimensions of A,
                  % n=number of data points,

xmin=min(A);      % Find the minimum value of A.
xmax=max(A);      % Find the maximum value of A.
xrange=xmax-xmin; % Find the range of A.

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

if xr(1)<xmin,
  ax=xr(1);
  else,
    ax=xmin-xrange/4;   
end;

if xr(2)>xmax,
  bx=xr(2);
  else,
    bx=xmax+xrange/4;
end;

xa=linspace(ax,bx,inc); 

c=zeros(inc,1);

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

% Obtain the grid counts.

for i=1:n, % Loop over data points in x direction.
  c(binx(i))=c(binx(i))+(xa(binx(i)+1)-A(i))/deltax;
  c(binx(i)+1)=c(binx(i)+1)+(A(i)-xa(binx(i)))/deltax;  
end;

% Obtain the kernel weights.

L=min(floor(r*h*(inc-1)/(bx-ax)),inc-1);

kw=zeros(1,2*inc);

temp=feval(k,(bx-ax)*[0:L]/((inc-1)*h))/(n*h);

kw(inc+1:inc+L+1)=temp;
kw(inc-L+1:inc)=temp([L+1:-1:2]);

% Apply 'fftshift' to kw.

kw=fftshift(kw)';

% Perform the convolution.

y=real(ifft(fft(c,P).*fft(kw)));

y=y(1:inc).*(y(1:inc)>0);

if g~=0,
  plot(xa,y,'w'); % Plot a 2D k.e. function on the range defined by xa. 
  title(['h = ',num2str(h)]);
end;