
function x1=pTfit(x,M,P)
N=length(x);
x1=zeros(1,N);

x0_B=x(1:2*M+1);
xM=pMfit(x0_B,M,P);
x1(1:M+1)=xM(1:M+1);

x0_E=x(end-2*M:end);
xM=pMfit(x0_E,M,P);
x1(end-M:end)=xM(end-M:end);

for J=M+2:N-M-1
  x0_M=x(J-M:J+M);
  xM=pMfit(x0_M,M,P);
  x1(J)=xM(M+1);
end
end

function xM=pMfit(x,M,P);
F=zeros(P+1,1);
S=zeros(P+1,P+1);
for r=0:P;
    I=-M:1:M;
    F(r+1)=sum(x(I+M+1).*(I.^r));
    for k=0:P;
      S(r+1,k+1)=sum(I.^(k+r));
    end    
end
a=inv(S)*F;
a1=flipud(a);
aM=a1';
xM=polyval(aM,I);
end