当前位置：首页 > 编程日记 > 正文

Spline interpolation and Savitzki-Golay smoothing

编程日记 2024-12-30 06:10:00

转自：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

https://www.dkcj.cn/info/34053.html

Spline interpolation and Savitzki-Golay smoothing

相关文章：

SpringBoot实现热部署(修改class不需要重启)

Oracle中查看表空间的使用率的脚本

vue实现多个元素或多个组件之间动画效果

干货 | 20个教程，掌握时间序列的特征分析（附代码）

利用OpenCV求取图像的重心

C++中Ansi、Unicode、UTF8字符串之间的转换和写入

MyBatis的扩展点(plugins)

linux中使用CST时间

到「黄埔学院」去：打造AI首席架构师，第二期限量招募！

linux守护进程的创建

tesseract3.01的训练和使用

旷视推出鼻纹识别，用AI寻找丢失宠物

Qt浅谈之一：内存泄露（总结）

LINUX新手入门-1.装系统

Log4cplus1.04的使用

FRVT赛程全纪录：格灵深瞳全球排名前五

反转比特位(文章最后有干货)【转】

过关斩将打进Kaggle竞赛Top 0.3%，我是这样做的

JavaBean规范

vigra1.8.0的使用

17个Python小窍门

用产品思路建设中台，这走得通吗？| 白话中台

利用cvMinAreaRect2求取轮廓最小外接矩形

电脑开机显示Invalidsystemdisk

Windows7 64位下vs2008配置OpenCV2.3.1

腾讯拥抱开源：首次公布开源路线图，技术研发向共享、复用和开源迈进

Babylon.js 3.3发布：更强大的粒子系统和WebVR支持

基于Erlang语言的视频相似推荐系统 | 深度

MinGW简介

Confluence 6 创建小组的公众空间