pdfpoly.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

function pdffun=PDFpoly(data,L,bs,deg,cont,fixval);
% PDFpoly   Estimate the pdf (probability density function) by polynomials
%           i.e, a function in the pp-form (piecewise polynomial).  
% The input data should be some one-dimensional outcomes (results)
% of the random process for which the pdf is to be estimated.
% The pdf function is returned in piecewise polynomial (pp) form, 
% the pp-form and the polynomial approximation is done by the PolyInterPol function.
% The number of line (polynomial) segments is N.
% Note the lengths of bs, deg and cont which are (N+1), N and (N-1) respectively.
% A problem with this function is that the resulting function may be negative,
% this should never happen to a real probability density function.
%
% pdffun=PDFpoly(data,L,bs,deg);
% pdffun=PDFpoly(data,L,bs,deg,cont,fixval);
% -----------------------------------------------------------------------------------
% Arguments:
%   pdffun   - probability density function in pp-form
%   data     - many (real) outcomes of the random process
%              data out of range are not used, range is  bs(1) <= data(i) <= bs(N+1)
%   L        - number of points to use in the range, the pdf is first estimated
%              at these L points evenly distributed in range, this is
%              done based on 'histogram' of data.
%   bs       - break sequence, N+1 points, bs(1) < bs(2) < ... < bs(N+1)
%   deg      - degree of each polynomial, deg is a vector of length N.
%              0 <= deg(n) <= maxdeg, where maxdeg is the maximal degree of
%              the polynomials. (order is degree+1)
%   cont     - continuity properties at the break points, a vector of length (N-1)
%              if 0 - no continuity demands for this break point
%              1 - continuity, 2 also continuious first derivate, 
%              and 3 for continuious second derivate. 
%              if this variable is omitted no continuity properties is assumed
%   fixval   - a table of points where the function should have special values
%              fixval is of size Kx3 where fixval(k,1) is the x-value,
%              the x-value must be within range of breaks.
%              fixval(k,2) the y-value (or the value of the derivate), and
%              fixval(k,3) is the type or derivate degree, i.e. 0 for value
%              1 for first derivate, 2 for second derivate.
%              if this variable is omitted no fixed values is assumed
% -----------------------------------------------------------------------------------

%----------------------------------------------------------------------
% Copyright (c) 2002.  Karl Skretting.  All rights reserved.
% Hogskolen in Stavanger (Stavanger University), Signal Processing Group
% Mail:  [e-mail]   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  17.01.2002  KS, the first version made
% Ver. 1.1  05.02.2002  KS, y(end) divided by 2 in line 107
% Ver. 1.2  18.12.2002  KS: moved from ..\Frames\ to ..\FrameTools
%----------------------------------------------------------------------

Mfile='PDFpoly';
Display=1;
pdffun=0;

if nargin<4
   error([Mfile,': wrong number of arguments, see help.']);
end

data=data(:);
K=length(data);
bs=bs(:);
N=length(bs)-1;      % number of segments
deg=deg(:);
n=length(deg);
if (n~=N)
   error([Mfile,': size of bs and deg do not correspond, see help.']);
end
if nargin<5
   cont=zeros(N-1,1);
end
cont=cont(:);
if (length(cont)~=(N-1))
   error([Mfile,': cont is of wrong length, see help.']);
end
if nargin<6
   fixval=zeros(0,3);
end
if (size(fixval,2)~=3)
   error([Mfile,': fixval is of wrong size, see help.']);
end

% first find the L x values for which the pdf first should be estimated
delta=(bs(N+1)-bs(1))/(L+1);
x=linspace(bs(1)+delta/2,bs(N+1)-delta/2,L)';   % a column vector
% exclude data out of range
I1=find(data<bs(1));
I2=find(data>bs(N+1));
if length(I1)+length(I2)
   I=setdiff((1:length(data))',[I1;I2]);
   data=data(I);
end
% find y by moving each point of data to the neighboring points of x
y=zeros(L,1);
data=(data-bs(1))/delta;   % scale and translate data 0 <= data(i) <= L
I=find(data>(L-0.5001));
if length(I); data(I)=L-0.5001; end;
I=find(data<0.5);
if length(I); data(I)=0.5; end;
data=data+0.5;                     % now    1 <= data < L
for l=1:length(data)
   f=floor(data(l));   % 1 <= f <= (L-1)
   r=data(l)-f;
   y(f)=y(f)+(1-r);
   y(f+1)=y(f+1)+r;
end
if (y(end)>1.5*y(end-1))         % a practical adaption
   y(end)=y(end)/2;
end
% plot(x,y)
if length(data)<(10*L)
   % smooth data
   K=7;   % odd
   lp1=hamming(K)/sum(hamming(K));
   h1=conv(y,lp1);
   h1(((K+1)/2):(K-1))=h1(((K+1)/2):(K-1))+flipud(h1(1:((K-1)/2)));
   h1((L+1):(L+((K-1)/2)))=h1((L+1):(L+((K-1)/2)))+flipud(h1((L+((K+1)/2)):(L+K-1)));
   h1=h1(((K+1)/2):(L+((K-1)/2)));
   y=h1;
end

% now we approximate (x,y) by polynomials
pdffun=PolyInterPol(x,y,bs,deg,cont,fixval);
ppi=fnint(pdffun);
integral=fnval(ppi,bs(N+1));
if isstruct(pdffun)
   pdffun.coefs=pdffun.coefs/integral;
else
   pdffun((N+6):end)=pdffun((N+6):end)/integral;
end

if Display
   figure(1);clf;
   fnplt(pdffun);
   hold on;
   plot(x,y/integral,'g-');
   grid on;
end

return;