function [pdx,pdf,h,gx,gv,crest]=pdff(data,npts,h)
%
% Name: pdff
%
% Usage: [pdx,pdf,h,gx,gv,crest]=pdff(data,npts,h)
%
%
% COMPUTES AN ESTIMATED PROBABILITY DENSITY FUNCTION BASED ON DATA INPUT.
% BASED ON SILVERMAN, DENSITY ESTIMATION FOR STATISTICS AND DATA ANALYSIS.
% USING A GAUSSIAN KERNAL PP.65.
%
% data is the input vector.
% npts if supplied is the no. of pts to use in the computation of pdf.
% pdf is the returned pdf.
% pdx is the returned amplitude grid.
% if h is not supplied h=npts^-1/5*std*0.8
% if npts not supplied the whole of data is used.
% gx=gaussian x axis
% gv= gaussian with std equivilant to data.
% crest=data crest factor=peak/std
%
%

%
%  VERSION SWD970806
%
%  Copyright (c) Norman Hunter, Los Alamos National Laboratory
%
%  Used by permission.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This matlab source code was originally     %
% developed as part of "DIAMOND" at          %
% Los Alamos National Laboratory. It may     %
% be copied, modified, and distributed in    %
% any form, provided:                        %
%  a) This notice accompanies the files and  %
%     appears near the top of all source     %
%     code files.                            %
%  b) No payment or commercial services are  %
%     received in exchange for the code.     %
%                                            %
% Original copyright is reserved by the      %
% Regents of the University of California,   %
% in addition to Scott W. Doebling, Phillip  %
% J. Cornwell, Erik G. Straser, and Charles  %
% R. Farrar.                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%


sigdat=std(data);                       % standard deviation of data input.

if nargin <2

     npts=max(size(data));

end

if nargin <3

     mh=.80;

end

% basic formula,slight undersmooth.

hinit=2*mh*sigdat/(npts)^.2;              % initial value of h.

h=hinit/1.3;

dmin=min(data);

dmax=max(data);                         % range of data.

delta=dmax-dmin;

a=dmin-0.2*delta;                           % start of window interval.

b=dmax+0.2*delta;                           % end of window interval.

gridsize=2048;

%transfer data to grid.

hmax=2.0*npts^(-1/5)*sigdat;            % widest possible h.

ygrid=zeros(gridsize,1);

dm=(b-a)./(gridsize-1);                 % data grid increment.



for id=1:npts

     da=(data(id)-a)./dm;

     gindex=floor(da);                  % quantize data to grid.

     dl=da-gindex;                      % dist. to lower grid pt.

     dp=1.0-dl;                         % dist. to upper grid pt.

     ial=1.0/(npts*dm^2).*dp;           % lower and upper data.

     iau=1.0/(npts*dm^2).*dl;           % p. 63, Silverman.

     ygrid(gindex+1)=ygrid(gindex+1)+ial;

     ygrid(gindex+2)=ygrid(gindex+2)+iau;

     % zero index on data is 1 in array.

end

% ygrid has data on a 256 size grid. Can be any 2^n.

%disp('view data on grid, ygrid, vs raw data, data')

%keyboard

yl=1.0/(gridsize)*fft(ygrid);

l=[-gridsize/2:gridsize/2];

sl=2*pi*l/(b-a);

zeta=exp(-0.5*h^2*sl(1:gridsize).^2).*fftshift(yl)';

% Transform of Gaussian Kernal.

pdf=ifft(fftshift(zeta));

pdf=abs(pdf);

scale=1.0/(sum(pdf)*dm);

pdf=scale*pdf;

pdx=[a:dm:b];

pdf=fliplr(pdf);

sigma=std(data);

gx=[a:dm:b];

crest=max(data)/sigma;

sigma2=2*sigma^2;

gv=1.0/(sigma*sqrt(2*pi))*exp(-gx.^2/sigma2);


return





