Spline interpolation and Savitzki-Golay smoothing

转自：http://octave.1599824.n4.nabble.com/Spline-interpolation-and-Savitzki-Golay-smoothing-td1675136.html

## natural-cubic-spline interpolation
## usage: yspline = spline(x,y,xspline)
## example:
## x = 0:10; y = sin(x);
## xspline = 0:0.1:10; yspline = spline(x,y,xspline);
## plot(x,y,"+",xspline,yspline);
## Given the vectors x and y, which tabulate a function, with
## x(1) < x(2) < x(3) <... or x(1) > x(2) > x(3) >..., and given
## the vector xspline, this function returns a natural-cubic-spline
## interpolated vector yspline.
## author: Zdenek Remes, May 22, 1999


function ynew = spline(x,y,xnew)
[x,index]=sort(x);
y=y(index);
n=length(y);
y2(1)=0.0;
y2(n)=0.0;
u(1)=0.0;
for i=2:n-1
sig=(x(i)-x(i-1))/(x(i+1)-x(i-1));
p=sig*y2(i-1)+2.0;
y2(i)=(sig-1.0)/p;
u(i)=(y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1));
u(i)=(6.0*u(i)/(x(i+1)-x(i-1))-sig*u(i-1))/p;
endfor;
k=n-1;
while (k >= 1)
y2(k)=y2(k)*y2(k+1)+u(k);
k--;
endwhile;

i1=1; in=length(xnew);

#if (xnew(1) < x(1))
# error("spline: bad xspline");
#endif;
#if (xnew(in) > x(n))
# error("spline: bad xspline");
#endif;

if (xnew(1) == x(1))
ynew(1)=y(1);
i1=2;
endif;
if (xnew(in) == x(n))
ynew(in)=y(n);
in=in-1;
endif;


for i=i1:in
khi=n;
klo=1;
while ((khi-klo) > 1)
k=floor((khi+klo)/2);
if (x(k) > xnew(i))
khi=k;
else
klo=k;
endif;
endwhile;
h=x(khi)-x(klo);
a=(x(khi)-xnew(i))/h;
b=(xnew(i)-x(klo))/h;
ynew(i)=a*y(klo)+b*y(khi)+((a^3-a)*y2(klo)+(b^3-b)*y2(khi))*(h*h)/6.0;
endfor;
endfunction;

## Savitzky-Golay smoothing filter
## usage: [xsavgol,ysavgol]=savgol(x,y,nl,nr,m)
## example: x=0:0.01:3;y1=sin(x.^3);y=y1+(rand(1,301)-0.5)/3;
## [xsavgol,ysavgol]=savgol(x,y,10,10,2);
## plot(x,y,"+",xsavgol,ysavgol,x,y1)
## Given vectors x, y containing a tabulated data y=f(x) with
## equally spaced x's this function calculates smoothed data
## ysavgol=g(xsavgol) by Savitzky-Golay smoothing filter.
## nl is the number of leftward (past) data points used, while
## nr is the number of rightward (future) data points, making
## the total number of data points used nl+nr+1. m is the order
## of the smoothing polynomial, also equal to the highest
## conserved moment; usual values are m=2 or m=4.
## The idea of Savitzky-Golay filtering is to smooth the
## underlying data y=f(x) within the moving window not by a
## constant (whose estimate is the average), but by a poly-
## nomial of higher order. Thus for a point y(i) the function
## savgol fits by a least-squares method a polynomial to
## points y(i-nl), ..., y(i+nr) in the moving window, and
## then set g(i-nl+1) to the value of that polynomial at
## position x(i).
## Zdenek Remes, Mai 22, 1999

function [xnew,ynew]=savgol(x,y,nl, nr, M)
if max(diff(x,2))>100*eps
error("The x's must be equally spaced.")
endif
for i=-nl:nr
for j=0:M
A(i+nl+1,j+1)=i^j;
endfor
endfor
AA=inv(A'*A);
for i=-nl:nr
cc=0;
for m=0:M
cc=cc+AA(1,m+1)*i^m;
endfor
c(i+nl+1)=cc;
endfor

nx=length(x);
for i=nl:nx-nr-1
yy=0;
for j=-nl:nr
yy=yy+c(j+nl+1)*y(i+j+1);
endfor
xnew(i-nl+1)=x(i+1);
ynew(i-nl+1)=yy;
endfor
endfunction