轉(zhuǎn)自：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