polyinterpol.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

function pp=PolyInterPol(x,y,bs,deg,cont,fixval)
% PolyInterPol Polynomial interpolation of a function given as points. 
%              The function will be approximated by N polynomials, given in
% the pp-form (piecewise polynomial). ex: fnplt(pp); % plot function
% Note the lengths of bs, deg and cont which are (N+1), N and (N-1) respectively.
%  
% pp=PolyInterPol(x,y,bs,deg);
% pp=PolyInterPol(x,y,bs,deg,cont,fixval);
%-----------------------------------------------------------------------------------
% arguments:
%   x        - the x values of which the function is given, size Lx1
%              x should be in ascending order
%   y        - the function value for the corresponding x, size Lx1
%   bs       - break sequence, the start and end points for the N polynomials.
%              bs is a vector of length (N+1). Values of x should be within
%              range of bs: bs(1) <= x(i) <= bs(N+1)
%              bs should be in ascending order
%   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 bs.
%              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 
% 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  15.01.2002  KS: function made 
% Ver. 1.1  18.12.2002  KS: moved from ..\Frames\ to ..\FrameTools
%-----------------------------------------------------------------------------------


Mfile='PolyInterPol';
Display=1;
pp=0;

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

x=x(:);
y=y(:);
L=length(x);
l=length(y);
if (l~=L)
   error([Mfile,': size of x and y are not equal, see help.']);
end
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
K=size(fixval,1);
if (size(fixval,2)~=3)
   error([Mfile,': fixval is of wrong size, see help.']);
end

% two linear equation systems are built:
% a is the vector of unknowns for both equations, size Mx1
% 1) X*a=y      % X is LxM, a is Mx1, y is Lx1 and L>M
% this system should be approximated in the least squares sense.
% 2) Z*a=b      % Z is JxM, a is Mx1, b is Jx1 and J<M
% this system should be completly satisfied
% coefs will be an array of size Nx(maxorder) and contain
% the unknows, a, in row order possible with some zeros in the
% beginning of some of the lines.

K=size(fixval,1);
order=deg+1;       % number of unknowns in each segment
M=sum(order);      % number of unknowns (a)
N=length(bs)-1;    % number of segments
maxorder=max(order);
coefsMAPa=zeros(N,maxorder); % map from coefs to a, i.e. coefs(n,j)=a(coefsMAPa(n,j))
m=0;
for n=1:N
   for j=(maxorder-order(n)+1):maxorder
      m=m+1;
      coefsMAPa(n,j)=m;
   end
end
if m~=M
   error([Mfile,': logical error, line 106 m~=M']);
end
%
J=sum(cont)+K;     % number of conditions (equations) that should be fully satisfied
if Display
   disp([Mfile,': number of line (polynomial) segments is ',int2str(N)]);
   disp(['and total number coefficients is ',int2str(M),', but ',int2str(J),' extra']);
   disp(['conditions reduce this to ',int2str(M-J),' free coefficients.']);
end
if J>0
   Z=zeros(J,M);
   b=zeros(J,1);
   % first we do the continuity conditions
   j=0;
   for n=1:(N-1)
      if cont(n)>0
         % it should at least be continuous
         j=j+1;
         % disp([Mfile,': polynomial ',int2str(n),' and polynomial ',int2str(n+1),...
         %       ' should be continuous, cont(',int2str(n),')=',int2str(cont(n))]);
         % value for polynomial n at end
         xx=bs(n+1)-bs(n);
         p=0;xp=1;   % x to the power of p
         for oi=1:order(n)
            ai=coefsMAPa(n,maxorder-oi+1);
            Z(j,ai)=xp;
            xp=xp*xx;
         end
         % subtract value for polynomial (n+1) at beginning
         ai=coefsMAPa(n+1,maxorder);
         Z(j,ai)=-1;
         % b value for this equation is 0 which is ok!
      end
      if cont(n)>1
         % and the first derivate should be continuous
         j=j+1;
         % value for derivate of polynomial n at end
         xx=bs(n+1)-bs(n);
         p=0;xp=1;   % x to the power of p
         for oi=2:order(n)
            ai=coefsMAPa(n,maxorder-oi+1);
            Z(j,ai)=xp*(oi-1);
            xp=xp*xx;
         end
         % subtract value for polynomial (n+1) at beginning
         ai=coefsMAPa(n+1,maxorder-1);
         Z(j,ai)=-1;
         % b value for this equation is 0 which is ok!
      end
      if cont(n)>2
         % and the second derivate should be continuous
         j=j+1;
         % value for derivate of polynomial n at end
         xx=bs(n+1)-bs(n);
         p=0;xp=1;   % x to the power of p
         for oi=3:order(n)
            ai=coefsMAPa(n,maxorder-oi+1);
            Z(j,ai)=xp*(oi-1)*(oi-2);
            xp=xp*xx;
         end
         % subtract value for polynomial (n+1) at beginning
         ai=coefsMAPa(n+1,maxorder-2);
         Z(j,ai)=-2;
         % b value for this equation is 0 which is ok!
      end
   end
   % and then we do the additional conditions
   for k=1:K
      % do the special conditions
      j=j+1;
      xx=fixval(k,1)-bs;
      I=find(xx>=0);
      n=length(I);
      if (n>N); n=N; end;
      xx=xx(I(n));
      p=0;xp=1;   % x to the power of p
      if fixval(k,3)==0
         for oi=1:order(n)
            ai=coefsMAPa(n,maxorder-oi+1);
            Z(j,ai)=xp;
            b(j)=fixval(k,2);
            xp=xp*xx;
         end
      elseif fixval(k,3)==1
         for oi=2:order(n)
            ai=coefsMAPa(n,maxorder-oi+1);
            Z(j,ai)=xp*(oi-1);
            b(j)=fixval(k,2);
            xp=xp*xx;
         end
      elseif fixval(k,3)==2
         for oi=3:order(n)
            ai=coefsMAPa(n,maxorder-oi+1);
            Z(j,ai)=xp*(oi-1)*(oi-2);
            b(j)=fixval(k,2);
            xp=xp*xx;
         end
      end
   end
   % now equation 2 is ok, 2) Z*a=b
   [U,S,V]=svd(Z);   % do singular value decomposition
   m=1;
   while (S(m,m)>(2*eps)); 
      m=m+1; 
      if (m>J); break; end;
   end
   nZ=V(:,m:M);    % the nullspace of Z
   a0=Z\b;         % one simple solution, all solutions a=a0+nZ*am  (am is (M-m+1)x1)
else
   % no exact conditions
   a0=zeros(M,1);   % and am is Mx1
   nZ=eye(M);
end

% now build the first equation: 1) X*a=y, X is LxM, a is Mx1, y is Lx1 and L>M
% but the equation is solved to take the exact demands into considerations
X=zeros(L,M);
l2=0;
for n=1:N
   l1=l2+1;
   while (x(l1)<bs(n)); l1=l1+1; end;   % now x(l1) >= bs(n)
   l2=l1;
   while (x(l2)<bs(n+1)); l2=l2+1; if l2>L; break; end; end; 
   l2=l2-1;  %  now x(l2) < bs(n+1)
   if ((n==N) & (x(L)==bs(n+1))); l2=l2+1; end;
   % now x(l1:l2) should be the x values within segment n
   xx=x(l1:l2)-bs(n);
   xp=ones(l2-l1+1,1);
   for oi=1:order(n)
      ai=coefsMAPa(n,maxorder-oi+1);
      X(l1:l2,ai)=xp;
      xp=xp.*xx;
   end
end

% solve 1) X*a=y in least square sense with the extra condition of 2) Z*a=b
if 0     % check size of variables
   size(a0)
   size(nZ)
   size(X)
   size(y)
end
a=a0+nZ*((X*nZ)\(y-X*a0));   % here the solution is found

% now use this to build the pp form
coefs=zeros(N,maxorder);
for n=1:N
   for j=1:maxorder
      if coefsMAPa(n,j)
         coefs(n,j)=a(coefsMAPa(n,j));
      end
   end
end

pp=ppmak(bs,coefs,1);

return;