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