常用插值算法程序

#include<stdio.h>#include<math.h>float f(float x) //计算ex的值{return (exp(x));}float g(float x) //计算根号x的值{return (pow(x,0.5));}void linerity () //线性插值{float px,x;float x0,x1;printf("请输入x0，x1的值\n");scanf("%f,%f",&x0,&x1);printf("请输入x的值: ");scanf("%f",&x);px=(x-x1)/(x0-x1)*f(x0)+(x-x0)/(x1-x0)*f(x1);printf("f(%f)=%f \n",x,px);}void second () //二次插值{float x0,x1,x2,x,px;x0=0;x1=0.5;x2=2;printf("请输入x的值:");scanf("%f",&x);px=((x-x1)*(x-x2))/((x0-x1)*(x0-x2))*f(x0)+((x-x0)*(x-x2))/((x1-x0)*(x1-x2))*f(x1)+((x-x0)* (x-x1))/((x2-x0)*(x2-x1))*f(x2);printf("f(%f)=%f\n",x,px);}void Hermite () //Hermite插值{int i,k,n=2;int flag1=0;printf("Hermite插值多项式H5(x)=");for(i=0;i<=n;i++){int flag=0;flag1++;if(flag1==1){printf("y%d[1-2(x-x%d)*(",i,i);}else{printf("+y%d[1-2(x-x%d)*(",i,i);}for(k=0;k<=n;k++){if(k!=i){flag++;if(flag==1){printf("(1/x%d-x%d)",i,k);}else{printf("+(1/x%d-x%d)",i,k);}}}printf(")]*（");for(k=0;k<=n;k++){if(i!=k){printf("[(x-x%d)/(x%d-x%d)]2",i,k,i);}}printf("）");}printf("\n");}void sectionl () //分段线性插值{float x[5]={2.0,2.1,2.2,2.3,2.4};float y;printf("请输入y：");scanf("%f",&y);if(y>=2.0&&y<2.1){float px;px=((y-x[1])/(x[0]-x[1]))*g (x[0])+((y-x[0])/(x[1]-x[0]))*g (x[1]);printf("f(%f)=%f\n",y,px);}else if(y>=2.1&&y<2.2){float px;px=((y-x[2])/(x[1]-x[2]))*g (x[1])+((y-x[1])/(x[2]-x[1]))*g (x[2]);printf("f(%f)=%f\n",y,px);}else if(y>=2.2&&y<2.3){float px;px=((y-x[3])/(x[2]-x[3]))*g (x[2])+((y-x[2])/(x[3]-x[2]))*g (x[3]);printf("f(%f)=%f\n",y,px);}else if(y>=2.3&&y<2.4){float px;px=((y-x[4])/(x[3]-x[4]))*g (x[3])+((y-x[3])/(x[4]-x[3]))*g (x[4]);printf("f(%f)=%f\n",y,px);}else if(y>2.4) printf("**********ERROR!******************\n"); }void sectionp (){int i;float a[5]={2.0,2.1,2.2,2.3,2.4};float x,y;printf("input the data: x?\n");scanf("%f",&x);if(x<a[1]){i=1;goto loop;}if(x>a[4]){i=4;goto loop;}i=1;loop1:i++;if(x>a[i])goto loop1;if(fabs(x-a[i-1])<=fabs(x-a[i]))i=i-1;loop:y=g(a[i-1])*(x-a[i])*(x-a[i+1])/((a[i-1]-a[i])*(a[i-1]-a[i+1]));y=y+g(a[i])*(x-a[i-1])*(x-a[i+1])/((a[i]-a[i-1])*(a[i]-a[i+1]));y=y+g(a[i+1])*(x-a[i-1])*(x-a[i])/((a[i+1]-a[i-1])*(a[i+1]-a[i]));printf("f(%f)=%f\n",x,y);}int main(){char flag1='y';while(flag1=='y'){int flag=0;printf("*******[1]:线性插值***************\n");printf("*******[2]:二次插值***************\n");printf("*******[3]:Hermite插值************\n");printf("*******[4]:分段线性插值***********\n");printf("*******[5]:分段抛物线插值*********\n");printf("请输入：");scanf("%d",&flag);switch(flag){case 1:linerity ();break;case 2:second ();break;case 3:Hermite ();break;case 4:sectionl ();break;case 5:sectionp ();break;default:printf("error!!\n");}printf("是否继续？y/n \n");getchar();scanf("%c",&flag1);}return 0;}。

常见的插值方法及其原理这一节无可避免要接触一些数学知识，为了让本文通俗易懂，我们尽量绕开讨厌的公式等。

为了进一步的简化难度，我们把讨论从二维图像降到一维上。

首先来看看最简单的‘最临近像素插值’。

A,B是原图上已经有的点，现在我们要知道其中间X位置处的像素值。

我们找出X位置和A,B位置之间的距离d1,d2，如图，d2要小于d1,所以我们就认为X处像素值的大小就等于B处像素值的大小。

显然，这种方法是非常苯的，同时会带来明显的失真。

在A,B中点处的像素值会突然出现一个跳跃，这就是为什么会出现马赛克和锯齿等明显走样的原因。

最临近插值法唯一的优点就是速度快。

图10，最临近法插值原理接下来是稍微复杂点的‘线性插值’（Linear）线性插值也很好理解，AB两点的像素值之间，我们认为是直线变化的，要求X点处的值，只需要找到对应位置直线上的一点即可。

换句话说，A,B间任意一点的值只跟A,B有关。

由于插值的结果是连续的，所以视觉上会比最小临近法要好一些。

线性插值速度稍微要慢一点，但是效果要好不少。

如果讲究速度，这是个不错的折衷。

图11，线性插值原理其他插值方法立方插值，样条插值等等，他们的目的是试图让插值的曲线显得更平滑，为了达到这个目的，他们不得不利用到周围若干范围内的点，这里的数学原理就不再详述了。

图12，高级的插值原理如图，要求B,C之间X的值，需要利用B,C周围A,B,C,D四个点的像素值，通过某种计算，得到光滑的曲线，从而算出X的值来。

计算量显然要比前两种大许多。

好了，以上就是基本知识。

所谓两次线性和两次立方实际上就是把刚才的分析拓展到二维空间上，在宽和高方向上作两次插值的意思。

在以上的基础上，有的软件还发展了更复杂的改进的插值方式譬如S-SPline, Turbo Photo等。

他们的目的是使边缘的表现更完美。

插值(Interpolation)，有时也称为“重置样本”，是在不生成像素的情况下增加图像像素大小的一种方法，在周围像素色彩的基础上用数学公式计算丢失像素的色彩。

题 7：一维函数插值算法课题内容：课题 7：一维函数插值算法课题内容：对函数||e-y x=，取[-5,5]之间步长为 1 的值*10作为粗值，以步长0.1 作为细值，编写程序实现插值算法：最邻近插值算法，线性插值算法和三次多项式函数插值算法，并对比插值效果。

课题要求：1、设计良好的人机交互 GUI 界面。

2、自己编写实现插值算法。

3、在同一个图形窗口显示对比最后的插值效果。

附录一、界面设计二、图像结果三、程序设计1、线性插值function pushbutton1_Callback(hObject, eventdata, handles) x=-5:5;y=10*exp(-abs(x));f1=[];for x1=-5:0.1:5a=(x1-floor(x1));%请读者认真逐一带入推导if x1==floor(x1)f1=[f1,y(floor(x1)+6)];elsef1=[f1,y(floor(x1)+6)+a*(y(floor(x1)+7)-y(floor(x1)+6))]; endendm=-5:0.1:5plot(m,f1,'-r',x,y,'+')axis([-5 5 0 10])legend('liner插值','原函数');xlabel('X');ylabel('Y');title('liner插值与原函数的对比');grid2、多项式插值x0=-5:1:-3;y0=10*exp(-abs(x0));x=-5:0.1:-3;n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;endaxis([-5 5 0 10])plot(x,y,'m',x0,y0,'+')legend('三次多项式插值','原函数');xlabel('X');ylabel('Y');title('三次多项式插值与原函数的对比');gridhold onx0=-3:1:-1;y0=10*exp(-abs(x0));x=-3:0.1:-1;n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;endaxis([-5 5 0 10])plot(x,y,'m',x0,y0,'+')legend('三次多项式插值','原函数');xlabel('X');ylabel('Y');title('三次多项式插值与原函数的对比');gridhold onx0=-1:1:1;y0=10*exp(-abs(x0));x=-1:0.1:1;n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;endaxis([-5 5 0 10])plot(x,y,'m',x0,y0,'+')legend('三次多项式插值','原函数');xlabel('X');ylabel('Y');title('三次多项式插值与原函数的对比');gridhold onx0=1:1:3;y0=10*exp(-abs(x0));x=1:0.1:3;n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;endaxis([-5 5 0 10])plot(x,y,'m',x0,y0,'+')legend('三次多项式插值','原函数');xlabel('X');ylabel('Y');title('三次多项式插值与原函数的对比');gridhold onx0=3:1:5;y0=10*exp(-abs(x0));x=3:0.1:5;n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;endaxis([-5 5 0 10])plot(x,y,'m',x0,y0,'+')legend('三次多项式插值','原函数');xlabel('X');ylabel('Y');title('三次多项式插值与原函数的对比');grid3、最邻近插值function pushbutton3_Callback(hObject, eventdata, handles) x=-5:5;y=10*exp(-abs(x));f2=[];for x1=-5:0.1:5if abs(x1-floor(x1))<0.5f2=[f2,y(floor(x1)+6)];elsef2=[f2,y(floor(x1)+7)];endendm=[-5:0.1:5];f4=10*exp(-abs(m));plot(m,f2,'-r',x,y,'+')axis([-5 5 0 10])legend('nearest插值','原函数');xlabel('X');ylabel('Y');title('nearest插值与原函数的对比');grid。

空间插值算法：1、距离倒数乘方法（Inverse Distance to a Power）距离倒数乘方格网化方法是一个加权平均插值法，可以进行确切的或者圆滑的方式插值。

方次参数控制着权系数如何随着离开一个格网结点距离的增加而下降。

对于一个较大的方次，较近的数据点被给定一个较高的权重份额，对于一个较小的方次，权重比较均匀地分配给各数据点。

计算一个格网结点时给予一个特定数据点的权值与指定方次的从结点到观测点的该结点被赋予距离倒数成比例。

当计算一个格网结点时，配给的权重是一个分数，所有权重的总和等于1.0。

当一个观测点与一个格网结点重合时，该观测点被给予一个实际为 1.0 的权重，所有其它观测点被给予一个几乎为0.0 的权重。

换言之，该结点被赋给与观测点一致的值。

这就是一个准确插值。

距离倒数法的特征之一是要在格网区域内产生围绕观测点位置的"牛眼"。

用距离倒数格网化时可以指定一个圆滑参数。

大于零的圆滑参数保证，对于一个特定的结点，没有哪个观测点被赋予全部的权值，即使观测点与该结点重合也是如此。

圆滑参数通过修匀已被插值的格网来降低"牛眼"影响。

2、克里金法（Kriging）克里金法是一种在许多领域都很有用的地质统计格网化方法。

克里金法试图那样表示隐含在你的数据中的趋势，例如，高点会是沿一个脊连接，而不是被牛眼形等值线所孤立。

克里金法中包含了几个因子：变化图模型，漂移类型和矿块效应。

3、最小曲率法（Minimum Curvature）最小曲率法广泛用于地球科学。

用最小曲率法生成的插值面类似于一个通过各个数据值的，具有最小弯曲量的长条形薄弹性片。

最小曲率法，试图在尽可能严格地尊重数据的同时，生成尽可能圆滑的曲面。

使用最小曲率法时要涉及到两个参数：最大残差参数和最大循环次数参数来控制最小曲率的收敛标准。

4、多元回归法（Polynomial Regression）多元回归被用来确定你的数据的大规模的趋势和图案。

第16章插值法Lagrange插值求n次Lagrange插值多项式算法1．输入n+1个插值点：(x i, y i),i=0,1,…,n2．计算插值基函数l0n(x), l1n(x) ,…, l nn(x)3．给出n次Lagrange插值多项式：L n(x)= y0 l0n(x)+ y1 l1n(x) +…+y n l nn(x)求Lagrange插值多项式程序Clear[lag,xi,x,yi];xi=Input["xi="]yi=Input["yi="]n=Length[xi]-1;p=Sum[yi[[i]]*(Product[(x-xi[[j]])/(xi[[i]]-xi[[j]]),{j,1,i-1}]*Product[(x-xi[[j]])/(xi[[i]]-xi[[j]]),{j,i+1,n+1}]),{i,1,n+1}];lag[x_]=Simplify[p]说明：本程序用于求n次Lagrange插值多项式。

程序执行后，按要求通过键盘输入插值基点xi:{x0 , x1, ... , x n }和对应函数值yi:{ y0 , y1, … , y n }后，程序即可给出对应的n次Lagrange插值多项式lag[x]。

程序中变量说明xi：存放插值基点{x0 , x1, ... , x n }yi: 存放对应函数值{y0 , y1, … , y n}lag[x]: 存放求出的n次Lagrange插值多项式Ln(x)注：语句lag[x_]=Simplify[p]用简化形式给出对应的n次Lagrange插值多项式。

例题与实验例1．给定数据表x 0 1 2 3y=f(x) 1 3 5 12用Lagrange插值法求三次插值多项式,并给出函数f(x)在x =1.4的近似值。

解：执行Lagrange插值程序后，在输入的两个窗口中按提示分别输入{0, 1, 2, 3}、{1, 3, 5, 12}，每次输入后用鼠标点击窗口的“OK”按扭，得如下插值函数。

数值分析报告班级：专业：流水号：学号：姓名：常用的插值方法序言在离散数据的基础上补插连续函数，使得这条连续曲线通过全部给定的离散数据点。

插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。

早在6世纪，中国的刘焯已将等距二次插值用于天文计算。

17世纪之后，牛顿、拉格朗日分别讨论了等距和非等距的一般插值公式。

在近代，插值法仍然是数据处理和编制函数表的常用工具，又是数值积分、数值微分、非线性方程求根和微分方程数值解法的重要基础，许多求解计算公式都是以插值为基础导出的。

插值问题的提法是：假定区间[a，b〕上的实值函数f（x）在该区间上n＋1个互不相同点x0，x1……x n处的值是f（x0），……f（x n），要求估算f（x）在[a，b〕中某点的值。

其做法是：在事先选定的一个由简单函数构成的有n＋1个参数C0，C1，……C n的函数类Φ(C0，C1，……C n)中求出满足条件P（x i）＝f（x i）（i＝0，1，…… n）的函数P(x)，并以P(x)作为f(x)的估值。

此处f（x）称为被插值函数，x0，x1，……xn 称为插值结(节)点，Φ(C0，C1，……C n)称为插值函数类，上面等式称为插值条件，Φ（C0，……C n）中满足上式的函数称为插值函数，R（x）＝f（x）－P（x）称为插值余项。

求解这类问题，它有很多种插值法，其中以拉格朗日(Lagrange)插值和牛顿(Newton)插值为代表的多项式插值最有特点，常用的插值还有Hermit 插值，分段插值和样条插值。

一．拉格朗日插值1.问题提出：已知函数()y f x =在n+1个点01,,,n x x x 上的函数值01,,,n y y y ，求任意一点x '的函数值()f x '。

说明：函数()y f x =可能是未知的；也可能是已知的，但它比较复杂，很难计算其函数值()f x '。

5.2.3 n次拉格朗日插值
➢问题描述
插值基点：x0，x1，…，xn（n+1个点互异） 插值函数：不超过n次的多项式
插值条件：Ln(xi)=yi, i=0,1,2,…,n
➢基函数
li (x)
(x x0 ) (x xi1 )(x xi1 ) (x xn ) (xi x0 ) (xi xi1 )(xi xi1 ) (xi xn )
定义5-3
设H
是
n
不超过n次的多项式的全体的集合，
0
(
x)
,1
(
x),
,n (x)
是H n中n
1个线性无关的多项式，则0 (x),1 (x),
,
n
(
x)是H
的
n
一组基函数。
注意：基函数是不唯一的；
n
H n中的任一多项式pn (x)均可由基函数唯一表示，即pn (x) kii (x) i0
➢定理5-1 (插值函数的存在唯一性定理)
由于多项式有其优良的特性，所以通常都是用多项式作为 插值函数。还有其它类型的插值函数，如有理函数插值、 三角函数插值等
➢函数插值涉及的基本问题
插值函数的存在唯一性问题
插值函数的构造问题
截断误差估计与收敛性问题
➢ 代数多项式插值函数的构造方法
拉格朗日插值法 埃尔米特插值法
牛顿插值法
样条函数插值法
拉格朗日插值函数均可表示为一组基函数与函数值的线性组 合，这些基函数与被插函数无关，只需用插值基点有来构造。
5.2.1 拉格朗日线性插值L1(x) ➢线性插值及几何意义
n=1时的n次多项式L1(x) 称为线性插值。此时，有两个互异的 插值基点:x0，x1，插值条件为: L1(x0)=y0, L1(x1)=y1 。

opencv中resize函数五种插值算法;java -回复OpenCV中提供了多种图像处理函数和算法，其中之一就是resize函数，用于调整图片大小。

resize函数中有五种不同的插值算法，它们分别是：最近邻插值(Nearest-neighbor interpolation)、双线性插值(Bilinear interpolation)、双三次插值(Bicubic interpolation)、立方插值(Cubic interpolation)和区域插值(Area interpolation)。

这五种插值算法在不同的场景下有不同的效果和应用范围。

在本文中，我们将逐步解释这五种插值算法的原理和应用，然后给出一些使用Java实现的示例代码。

首先，让我们从最近邻插值算法开始。

最近邻插值算法是一种简单而直接的插值方法，它通过选择离目标像素最近的现有像素来进行图像的缩放。

这种插值算法的计算速度很快，但它会造成图像的锯齿状边缘，尤其是在图像的缩小过程中。

因此，最近邻插值算法主要适用于图像放大的情况，不过其精度相对较低。

接下来是双线性插值算法，它是一种基于线性插值的方法。

相比于最近邻插值算法，双线性插值算法通过对目标像素四个相邻像素的加权平均值来计算新的像素值。

这样可以得到更平滑的缩放效果，减少了锯齿状边缘的出现。

双线性插值算法的计算复杂度较低，是一种较常用且适用于大多数图像处理任务的插值方法。

双三次插值算法是一种更高级的插值方法，它通过在四个最近像素的基础上使用三次多项式来进行插值计算。

这种插值方法可以提供更高质量和更准确的图像缩放效果，同时也增加了计算复杂度。

双三次插值算法通常适用于需要高质量图像处理的场景，比如数字图像处理和医学图像处理等领域。

立方插值算法是一种更一般化的插值方法，它通过在二维空间中使用分数幂函数来对像素进行插值计算。

这种方法可以提供更平滑的缩放效果，并且能够处理更高阶的插值计算。

立方插值算法适用于一些对图像质量要求较高和需要更高阶插值计算的图像处理任务。

如何进行线性插值假设我们已知坐标 (x0, y0) 与(x1, y1)，要得到 [x0, x1] 区间内某一位置x在直线上的值。

根据图中所示，我们得到假设方程两边的值为α，那么这个值就是插值系数—从x0到x的距离与从x0到x1距离的比值。

由于x值已知，所以可以从公式得到α的值同样，这样，在代数上就可以表示为：或者，这样通过α就可以直接得到y。

实际上，即使x不在x0到x1之间并且α也不是介于 0 到 1 之间，这个公式也是成立的。

在这种情况下，这种方法叫作线性外插—参见外插值。

已知y求x的过程与以上过程相同，只是x与y要进行交换。

[编辑] 线性插值近似法线性插值经常用于已知函数f在两点的值要近似获得其它点数值的方法，这种近似方法的误差定义为其中p表示上面定义的线性插值多项式根据罗尔定理，我们可以证明：如果f有两个连续导数，那么误差范围是正如所看到的，函数上两点之间的近似随着所近似的函数的二阶导数的增大而逐渐变差。

从直观上来看也是这样：函数的曲率越大，简单线性插值近似的误差也越大。

[编辑] 应用线性插值经常用于补充表格中的间隔部分。

假设一个表格列出了一个国家1970年、1980年、1990年以及2000年的人口，那么如果需要估计1994年的人口的话，线性插值就是一种简便的方法。

两值之间的线性插值基本运算在计算机图形学中的应用非常普遍，以至于在计算机图形学领域的行话中人们将它称为lerp。

所有当今计算机图形处理器的硬件中都集成了线性插值运算，并且经常用来组成更为复杂的运算：例如，可以通过三步线性插值完成一次双线性插值运算。

由于这种运算成本较低，所以对于没有足够数量条目的光滑函数来说，它是实现精确快速查找表的一种非常好的方法。

[编辑] 历史线性插值从远古以来就一直用于补充表格中的间隔，它经常用于天文学数据。

人们相信公元前最後三个世纪的塞琉西王朝、公元前2世纪的希腊天文学家与数学家喜帕恰斯就曾经使用了这种方法。

最临近插值法原理：这种算法就是根据原图像和目标图像的尺寸，计算缩放的比例，然后根据缩放比例计算目标像素所依据的原像素，过程中自然会产生小数，这时就采用四舍五入，取与这个点最相近的点。

图解如下：如果(i+u, j+v)落在A区，即u<0.5, v<0.5，则将左上角象素的灰度值赋给待求象素，同理，落在B区则赋予右上角的象素灰度值，落在C区则赋予左下角象素的灰度值，落在D区则赋予右下角象素的灰度值。

具体Matlab源代码实现：clear all;A=imread('12.png'); %读取图像信息imshow(A); %显示原图title('原图128*128');Row = size(A,1); Col = size(A,2); %图像行数和列数nn=8; %放大倍数m = round(nn*Row); %求出变换后的坐标的最大值n = round(nn*Col);B = zeros(m,n,3); %定义变换后的图像for i = 1 : mfor j = 1 : nx = round(i/nn); y = round(j/nn); %最小临近法对图像进行插值if x==0 x = 1; endif y==0 y = 1; endif x>Row x = Row; endif y>Col y = Col;endB(i,j,:) = A(x,y,:);endendB = uint8(B); %将矩阵转换成8位无符号整数figure;imshow(B); %显示输出图片title('最邻近插值法放大8倍1024*1024');运行程序后，原图如图1所示：图1用最邻近插值法放大8倍后的图如图2所示：图2双线性内插值法：计算过程简单了解，如图，已知Q12，Q22，Q11，Q21，但是要插值的点为P 点，这就要用双线性插值了，首先在x轴方向上，对R1和R2两个点进行插值，这个很简单，然后根据R1和R2对P点进行插值，这就是所谓的双线性插值。

数学建模插值算法插值算法是数学建模中一种常用的技术，用于在已知数据点处的估计和未知数据点之间的预测。

插值算法可以帮助我们充分利用已知数据点的信息，获得更完整和连续的数据。

在数学建模中，插值算法有多种方法可选，常见的包括拉格朗日插值、牛顿插值、样条插值等。

拉格朗日插值是最常见和简单的插值方法之一、它的基本思想是通过构造一个n次多项式来近似通过已知数据点的曲线。

具体地说，我们可以根据已知数据点的横纵坐标，构造出n个满足这些坐标的插值基函数。

然后，将这些插值基函数分别与相应基函数在未知数据点处取值的乘积相加，得到插值多项式。

最后，利用这个多项式来估计未知数据点的纵坐标。

牛顿插值是另一种常用的插值方法。

它的基本思想是使用差商的概念来创建一个n次多项式。

差商是一个递归定义的概念，其基本思想是通过逐步添加一个已知数据点来计算多项式的高次项系数。

具体地说，我们可以根据已知数据点的横纵坐标，构造出n个差商。

然后，将这些差商与相应基函数在未知数据点处取值的乘积相加，得到插值多项式。

最后，利用这个多项式来估计未知数据点的纵坐标。

样条插值是一种更加复杂但更精确的插值方法。

它的基本思想是通过构造一组n次多项式的集合，使得每个多项式在相应数据点处完全符合已知数据。

具体地说，我们可以根据已知数据点的横纵坐标，构造出n个多项式，并设置它们在数据点处的约束条件。

然后，通过求解一个线性方程组来计算每个多项式的系数。

最后，利用这组多项式来估计未知数据点的纵坐标。

以上是数学建模中常用的几种插值算法，它们各有优缺点，在不同情景下有着不同的适用性。

插值算法在实际应用中具有广泛的用途，例如地图绘制、图像处理、信号处理等领域。

在进行插值计算时，要根据实际情况选择适当的算法，并合理处理计算误差，以提高插值结果的准确性和稳定性。

插值法一．实验目的•熟悉Matlab 编程； • 学习插值方法及程序设计算法二．实验内容1. 已知函数在下列各点的值为试用4次Newton 插值多项式4()P x 及分段线性插值函数()h I x 对数据进行插值.用图给出{}(,),0.20.08,0,1,,10i i i x y x i i =+=，4()P x 及()h I x .2. 在区间[1,1]-上分别取10n =、20用两组等距节点对龙格函数21()125f x x =+作多项式插值及分段线性插值，对每个n 值，分别画出插值函数及()f x 的图像.3. 下列数据点的插值可以得到平方根函数的近似，在区间[0,64]上作图.（1）用这9个点作8次Lagrange 插值8()L x .（2）用分段线性插值程序求()h I x .从得到结果看在[0,64]上，哪个插值更精确；在区间[0,1]上，两种插值哪个更精确？三．实验内容1.（1）牛顿插值程序function y0=new(x,y,x0)n=length(x);m=length(y);if n~=merror('the lengths of x and y must be equal.');end%计算均差表YY=zeros(n,n);Y(:,1)=y';for k=1:n-1for i=1:n-kif abs(x(i+k)-x(i))<epserror('the DATA is error .')endY(i+k,k+1)=(Y(i+k,k)-Y(i+k-1,k))./(x(i+k)-x(i));endend%计算牛顿插值公式y0=Y(1,1);for i=2:nz=1;for k=1:i-1z=z.*(x0-x(k));endy0=y0+Y(i,i).*z;end命令窗口输入命令>> x=[0.2:0.2:1.0];>> y=[0.98 0.92 0.81 0.64 0.38];>> x0=[0.2:0.08:1.0];>> P4=new(x,y,x0);（2）分段线性插值程序function y=divline(x0,y0,x)n=length(x0);m=length(x);for j=1:mfor i=1:n-1if (x(j)>=x0(i))&(x(j)<=x0(i+1))y(j)=(x(j)-x0(i+1))./(x0(i)-x0(i+1)).*y0(i)+(x(j)-x0(i))./(x0(i+1)-x0(i)).*y0(i+1);elsecontinue;endendend命令窗口输入命令>>x=[0.2:0.2:1.0];>> y=[0.98 0.92 0.81 0.64 0.38];>> x0=[0.2:0.08:1.0];>> I=divline(x,y,x0);实验结果（含实验图像）（1）牛顿插值画图命令及图像>> plot(x0,P4)（2）分段线性插值画图命令及图像>> plot(x0,I)2.function f=zuoye2_2(n)x0=linspace(-1,1,n+1);y0=1./(1+25*x0.^2);x=linspace(-1,1,100);y=1./(1+25*x.^2);y3=interp1(x0,y0,x,'spline');for k=1:length(x)im=ones(1,n+1);for i=1:n+1for j=1:n+1if x0(i)~=x0(j)im(i)=im(i)*(x(k)-x0(j))/(x0(i)-x0(j));endendendyp(k)=im*y0';endplot(x0,y0,'*',x,y,'b',x,y3,'-.',x,yp,'--');%画图legend('原插值点','原函数','三次样条插值','多项式插值')%加图例以zuoye2_2.m 保存文件。

常见图像插值算法只有3种么？电脑摄像头最高只有130万像素的，800万是通过软件修改的。

何为数码插值(软件插值)插值(Interpolation)，有时也称为“重置样本”，是在不生成像素的情况下增加图像像素大小的一种方法，在周围像素色彩的基础上用数学公式计算丢失像素的色彩。

简单地说，插值是根据中心像素点的颜色参数模拟出周边像素值的方法，是数码相机特有的放大数码照片的软件手段。

一、认识插值的算法“插值”最初是电脑术语，后来引用到数码图像上来。

图像放大时，像素也相应地增加，但这些增加的像素从何而来？这时插值就派上用场了。

插值就是在不生成像素的情况下增加图像像素大小的一种方法，在周围像素色彩的基础上用数学公式计算丢失像素的色彩(也有些相机使用插值，人为地增加图像的分辨率)。

所以在放大图像时，图像看上去会比较平滑、干净。

但必须注意的是插值并不能增加图像信息。

以图1为原图(见图1)，以下是经过不同插值算法处理的图片。

1.最近像素插值算法最近像素插值算法(Nearest Neighbour Interpolation)是最简单的一种插值算法，当图片放大时，缺少的像素通过直接使用与之最接近的原有像素的颜色生成，也就是说照搬旁边的像素，这样做的结果是产生了明显可见的锯齿(见图2)。

2.双线性插值算法双线性插值算法(Bilinear Interpolation)输出的图像的每个像素都是原图中四个像素(2×2)运算的结果，这种算法极大程度上消除了锯齿现象(见图3)。

3.双三次插值算法双三次插值算法(Bicubic Interpolation)是上一种算法的改进算法，它输出图像的每个像素都是原图16个像素(4×4)运算的结果(见图4)。

这种算法是一种很常见的算法，普遍用在图像编辑软件、打印机驱动和数码相机上。

4.分形算法分形算法(Fractal Interpolation)是Altamira Group提出的一种算法，这种算法得到的图像跟其他算法相比更清晰、更锐利(见图5)。

常见插值算法--拉格朗⽇插值、三次卷积插值、三次样条插值、兰克索斯插值写在前⾯本⽂简单介绍了⼏种常见的插值算法并附带了相应的python代码，本⽂公式使⽤latex编写，如有错误欢迎评论指出，如果谁知道如何修改latex字号也欢迎留⾔关于⼀维、⼆维和多维插值三次卷积插值、拉格朗⽇两点插值（线性插值）、兰克索斯插值在⼆维插值时改变x和y⽅向的计算顺序不影响最终结果，这三个也是图像缩放插值时常⽤的插值算法，⽽其他插值在改变计算顺序时会产⽣明显差异，多维的情况笔者没有尝试，读者可以⾃⾏尝试或推导最近邻插值法（Nearest Neighbour Interpolation）在待求像素的四邻像素中，将距离待求像素最近的像素值赋给待求像素p_{11}p_{12}pp_{21}p_{22}python代码1def NN_interpolation(srcImg, dstH, dstW):2 scrH, scrW, _ = srcImg.shape3 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)4for i in range(dstH - 1):5for j in range(dstW - 1):6 scrX = round(i * (scrH / dstH))7 scrY = round(j * (scrW / dstW))8 dstImg[i, j] = srcImg[scrX, scrY]9return dstImg拉格朗⽇插值（Lagrange Interpolation）拉格朗⽇插值法需要找到k个p_i(x)函数，使得每个函数分别在在x_i处取值为1，其余点取值为0，则y_ip_i(x)可以保证在x_i处取值为y_i，在其余点取值为0，因此L_k(x)能恰好经过所有点，这样的多项式被称为拉格朗⽇插值多项式，记为L_k(x)=\sum_{i=1}^ky_ip_i(x)p_i(x)=\prod_{j \neq i}^{1 \leq j \leq k}\frac{x-x_j}{x_i-x_j}以四点即三次图像插值为例，因为横坐标间隔为1，则设四个点横坐标为-1、0、1和2，可得p_1(x)、p_2(x)、p_3(x)和p_4(x)假设y_1、y_2、y_3和y_4分别为1、2、-1、4，则可得拉格朗⽇函数如下图所⽰，待插值点横坐标范围为[0,1]在K=2时在k=2时，也被称为线性插值通⽤公式p_1=\frac{x-x_2}{x_1-x_2}p_2=\frac{x-x_1}{x_2-x_1}\begin{align} L_2x &= p_1y_1+p_2y_2 \nonumber \\ &= \frac{x-x_2}{x_1-x_2}y_1 + \frac{x-x_1}{x_2-x_1}y_2 \nonumber \end{align}图像插值像素分布如图所⽰p_{11}p_{12}pp_{21}p_{22}即当x_{i+1}=x_i+1时，设p与p_{11}的横纵坐标差分别为dx和dy\begin{align} L_2x &= \frac{x-x_2}{x_1-x_2}y_1 + \frac{x-x_1}{x_2-x_1}y_2 \nonumber \\ &= (x_2-x)y_1+(x-x_1)y_2 \nonumber \\ &= (1-dx)y_1+dxy_2 \nonumber \\ &= (y_2-y_1)dx+y_1 \nonumber \end{align}L_2'x=y_2-y_1在K=3时通⽤公式p_1=\frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}p_2=\frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}p_3=\frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\begin{align} L_3x &= p_1y_1+p_2y_2+p_3y_3 \nonumber \\ &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}y_1+\frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}y_2+\frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}y_3 \nonumber \end{align}图像插值像素分布如图所⽰p_{11}p_{12}p_{13}p_{21}p_{22}p_{23}pp_{31}p_{32}p_{33}即当x_{i+1}=x_i+1时，设p与p_{11}的横纵坐标差分别为dx和dy\begin{align} L_3x &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}y_1 + \frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}y_2 + \frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}y_3 \nonumber \\ &= \frac{-dx(1-dx)}{(-1)\cdot(-2)}y_1 + \frac{-(1+dx)(1-dx)}{1\cdot(-1)}y_2 + \frac{(1+dx)dx}{2\cdot 1}y_3 \nonumber \\ &= (\frac{1}{2}d^2x-\frac{1}{2}dx)y_1 - (d^2x-1)y_2 + (\frac{1}{2}d^2x+\frac{1}{2}dx)y_3 \nonumber \\ &= d^2x(\frac{1}{2}y_1-y_2+\frac{1}{2}y_3)+dx(-\frac{1}{2}y_1+\frac{1}{2}y_3)+y_2 \nonumber \end{align}L_3'x=dx(y_1-2y_2+y_3)+(\frac{1}{2}y_3-\frac{1}{2}y_1)在K=4时通⽤公式p_1=\frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}\frac{x-x_4}{x_1-x_4}p_2=\frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}\frac{x-x_4}{x_2-x_4}p_3=\frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\frac{x-x_4}{x_3-x_4}p_4=\frac{x-x_1}{x_4-x_1}\frac{x-x_2}{x_4-x_2}\frac{x-x_3}{x_4-x_3}\begin{align} L_4x &= p_1y_1+p_2y_2+p_3y_3+p_4y_4 \nonumber \\ &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}\frac{x-x_4}{x_1-x_4}y_1 + \frac{x-x_1}{x_2-x_1}\frac{x-x_3} {x_2-x_3}\frac{x-x_4}{x_2-x_4}y_2 + \frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\frac{x-x_4}{x_3-x_4}y_3 + \frac{x-x_1}{x_4-x_1}\frac{x-x_2}{x_4-x_2}\frac{x-x_3}{x_4-x_3}y_4\nonumber \end{align}图像插值p_{11}p_{12}p_{13}p_{14}p_{21}p_{22}p_{23}p_{24}pp_{31}p_{32}p_{33}p_{34}p_{41}p_{42}p_{43}p_{44}即当x_{i+1}=x_i+1时，设p与p_{11}的横纵坐标差分别为dx和dy\begin{align} L_4x &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}\frac{x-x_4}{x_1-x_4}y_1 + \frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}\frac{x-x_4}{x_2-x_4}y_2 + \frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\frac{x-x_4}{x_3-x_4}y_3 + \frac{x-x_1}{x_4-x_1}\frac{x-x_2}{x_4-x_2}\frac{x-x_3}{x_4-x_3}y_4 \nonumber \\ &= \frac{dx[-(1-dx)][-(2-dx)]}{(-1)\cdot(-2)\cdot(-3)}y_1 + \frac{(1+dx)[-(1-dx)][-(2-dx)]}{1\cdot(-1)\cdot(-2)}y_2 + \frac{(1+dx)dx[-(2-dx)]}{2\cdot 1\cdot(-1)}y_3 + \frac{(1+dx)dx[-(1-dx)]}{3\cdot 2\cdot 1}y_4 \nonumber \\ &= \frac{d^3x-3d^2x+2dx}{-6}y1 + \frac{d^3x-2d^2x-dx+2}{2}y_2 + \frac{d^3x-d^2x-2dx}{-2}y_3 + \frac{d^3x-dx}{6}y_4 \nonumber \\ &= d^3x(-\frac{1}{6}y_1+\frac{1}{2}y_2-\frac{1} {2}y_3+\frac{1}{6}y_4)+d^2x(\frac{1}{2}y_1-y_2+\frac{1}{2}y_3)+dx(-\frac{1}{3}y_1-\frac{1}{2}y_2+y_3-\frac{1}{6}y_4)+y_2 \nonumber \end{align}\begin{align} L_4'x &= d^2x(-\frac{1}{2}y_1+\frac{3}{2}y_2-\frac{3}{2}y_3+\frac{1}{2}y_4)+dx(y_1-2y_2+y_3)+(-\frac{1}{3}y_1-\frac{1}{2}y_2+y_3-\frac{1}{6}y_4) \nonumber \\ &= -[\frac{1}{2}d^2x(y_1-3y_2+3y_3-y_4)-dx(y_1-2y_2+y_3)+\frac{1}{6}(2y_1+3y_2-6y_3+y_4)] \nonumber \end{align}python代码插值核计算的时候乘法和加减法计算的顺序不同可能会导致结果存在细微的差异，读者可以⾃⾏研究⼀下1class BiLagrangeInterpolation:2 @staticmethod3def LagrangeInterpolation2(x, y1, y2):4 f1 = 1 - x5 f2 = x6 result = y1 * f1 + y2 * f27return result89 @staticmethod10def LagrangeInterpolation3(x, y1, y2, y3):11 f1 = (x ** 2 - x) / 2.012 f2 = 1 - x ** 213 f3 = (x ** 2 + x) / 2.014 result = y1 * f1 + y2 * f2 + y3 * f315return result1617 @staticmethod18def LagrangeInterpolation4(x, y1, y2, y3, y4):19 f1 = - (x ** 3 - 3 * x ** 2 + 2 * x) / 6.020 f2 = (x ** 3 - 2 * x ** 2 - x + 2) / 2.021 f3 = - (x ** 3 - x ** 2 - 2 * x) / 2.022 f4 = (x ** 3 - x) / 6.023 result = y1 * f1 + y2 * f2 + y3 * f3 + y4 * f424return result2526def biLag2_2(self, srcImg, dstH, dstW):27 dstH, dstW = int(dstH), int(dstW)28 srcH, srcW, _ = srcImg.shape29 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')30 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)31for dstY in range(dstH):32for dstX in range(dstW):33for channel in [0, 1, 2]:34# p11 p1235# p36# p21 p2237# 储存为 p(y, x)38 p = [dstY * srcH / dstH, dstX * srcW / dstW]39 p11 = [math.floor(p[0]), math.floor(p[1])]40 p12 = [p11[0], p11[1] + 1]4142 p21 = [p11[0] + 1, p11[1]]43 p22 = [p21[0], p12[1]]4445 diff_y, diff_x = p[0] - p11[0], p[1] - p11[1]46 r1 = grangeInterpolation2(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel])47 r2 = grangeInterpolation2(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel])4849 c = grangeInterpolation2(diff_y, r1, r2)5051 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)52return dstImg5354def biLag3_3(self, srcImg, dstH, dstW):55 dstH, dstW = int(dstH), int(dstW)56 srcH, srcW, _ = srcImg.shape57 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')58 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)59for dstY in range(dstH):60for dstX in range(dstW):61for channel in [0, 1, 2]:62# p11 p12 p1363#64# p21 p22 p2365# p66# p31 p32 p3367# 储存为 p(y, x)68 p = [dstY * srcH / dstH, dstX * srcW / dstW]69 p22 = [math.floor(p[0]), math.floor(p[1])]70 p21 = [p22[0], p22[1] - 1]71 p23 = [p22[0], p22[1] + 1]7273 p11 = [p21[0] - 1, p21[1]]74 p12 = [p11[0], p22[1]]75 p13 = [p11[0], p23[1]]7677 p31 = [p21[0] + 1, p21[1]]78 p32 = [p31[0], p22[1]]79 p33 = [p31[0], p23[1]]8081 diff_y, diff_x = p[0] - p22[0], p[1] - p22[1]82 r1 = grangeInterpolation3(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel], srcImg[p13[0], p13[1], channel])83 r2 = grangeInterpolation3(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel], srcImg[p23[0], p23[1], channel])84 r3 = grangeInterpolation3(diff_x, srcImg[p31[0], p31[1], channel], srcImg[p32[0], p32[1], channel], srcImg[p33[0], p33[1], channel]) 8586 c = grangeInterpolation3(diff_y, r1, r2, r3)8788 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)89return dstImg9091def biLag4_4(self, srcImg, dstH, dstW):92 dstH, dstW = int(dstH), int(dstW)93 srcH, srcW, _ = srcImg.shape94 srcImg = np.pad(srcImg, ((1, 2), (1, 2), (0, 0)), 'edge')95 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)96for dstY in range(dstH):97for dstX in range(dstW):98for channel in [0, 1, 2]:99# p11 p12 p13 p14100#101# p21 p22 p23 p24102# p103# p31 p32 p33 p34104#105# p41 p42 p43 p44106# 储存为 p(y, x)107 p = [dstY * srcH / dstH, dstX * srcW / dstW]108 p22 = [math.floor(p[0]), math.floor(p[1])]109 p21 = [p22[0], p22[1] - 1]110 p23 = [p22[0], p22[1] + 1]111 p24 = [p22[0], p22[1] + 2]112113 p11 = [p21[0] - 1, p21[1]]114 p12 = [p11[0], p22[1]]115 p13 = [p11[0], p23[1]]116 p14 = [p11[0], p24[1]]117118 p31 = [p21[0] + 1, p21[1]]119 p32 = [p31[0], p22[1]]120 p33 = [p31[0], p23[1]]121 p34 = [p31[0], p24[1]]122123 p41 = [p21[0] + 2, p21[1]]124 p42 = [p41[0], p22[1]]125 p43 = [p41[0], p23[1]]126 p44 = [p41[0], p24[1]]127128 diff_y, diff_x = p[0] - p22[0], p[1] - p22[1]129 r1 = grangeInterpolation4(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel], srcImg[p13[0], p13[1], channel], srcImg[p14[0], p14[1], channel]) 130 r2 = grangeInterpolation4(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel], srcImg[p23[0], p23[1], channel], srcImg[p24[0], p24[1], channel]) 131 r3 = grangeInterpolation4(diff_x, srcImg[p31[0], p31[1], channel], srcImg[p32[0], p32[1], channel], srcImg[p33[0], p33[1], channel], srcImg[p34[0], p34[1], channel]) 132 r4 = grangeInterpolation4(diff_x, srcImg[p41[0], p41[1], channel], srcImg[p42[0], p42[1], channel], srcImg[p43[0], p43[1], channel], srcImg[p44[0], p44[1], channel]) 133134 c = grangeInterpolation4(diff_y, r1, r2, r3, r4)135136 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)137return dstImg三次卷积插值法（Cubic Convolution Interpolation）使⽤上图中的卷积核进⾏加权平均计算，卷积核为u(s)，四个等距（距离为1）的采样点记为x_0、x_1、x_2和x_3，采样数值记为y_0、y_1、y_2和y_3，且保证四个点均在[-2,2]区间上，计算得到g(x)，假设y_1、y_2、y_3和y_4分别为1、2、-1、4，则可得三次卷积插值函数如下图所⽰，待插值点横坐标范围为[0,1]公式推导设u(s)=\begin{cases} A_1|s|^3+B_1|s|^2+C_1|s|+D_1, &0<|s|<1 \\ A_2|s|^3+B_2|s|^2+C_2|s|+D_2, &1<|s|<2 \\ 1, &s=0 \\ 0, &otherwise \end{cases}\because函数在s=0,1,2处连续\therefore\begin{cases} 1=u(0^+)=D_1 \\ 0=u(1^-)=A_1+B_1+C_1+D_1 \\ 0=u(1^+)=A_2+B_2+C_2+D_2 \\ 0=u(2^-)=8A_2+4B_2+2C_2+D_2 \end{cases} (1)\because函数在s=0,1,2处导函数连续\therefore\begin{cases} u'(0^-)=u'(0+) \\ u'(1^-)=u'(1+) \\ u'(2^-)=u'(2+)\end{cases} \Rightarrow \begin{cases} -C_1=C_1 \\ 3A_1+2B_1+C_1=3A_2+2B_2+C_2\\ 12A_2+4B_2+C+2=0 \end{cases} ~~~~ (2)联⽴⽅程组(1)(2)，设A_2=a，解得\begin{cases} A_1=a+2 \\ B_1=-(a+3) \\ C_1=0 \\ D_1=1 \\ A_2=a \\ B_2=-5a \\ C_2=8a \\ D_2=-4a \end{cases}\Rightarrow u(s)=\begin{cases} (a+2)|s|^3-(a+3)|s|^2+1, &0<|s|<1 \\ A_2|s|^3+B_2|s|^2+C_2|s|+D_2, &1<|s|<2\\ 1, &s=0 \\ 0, &otherwise \end{cases}\because g(x)=\sum_kC_ku(s+j-k), ~~~~k=j-1,j, j+1,j+2且0<s<1⼜\because \begin{cases}\begin{align} u(s+1)&=as^3-2as^2+as \nonumber \\ u(s)&=(a+2)s^3-(a+3)s^2+1 \nonumber \\ u(s-1)&=-(a+2)s^3+(2a+3)s^2-as \nonumber \\ u(s-2)&=-as^3+as^2 \nonumber \end{align}\end{cases}\begin{align} \therefore g(x) &= C_{j-1}u(s+1)+C_{j}u(s)+C_{j+1}u(s-1)+C_{j+2}u(s-2) \nonumber \\ &= C_{j-1}(as^3-2as^2+as)+C_j[(a+2)s^3-(a+3)s^2+1]+C_{j+1}[-(a+2)s^3+ (2a+3)s^2-as]+C_{j+2}[-a^3+as^2] \nonumber \\ &= s^3[aC_{j-1}+(a+2)C_j-(a+2)C_{j+1}-aC_{j+2}]+s^2[-2aC_{j-1}-(a+3)C_j+(2a+3)C_{j+1}+aC_{j+2}]+s[aC_{j-1}-aC_{j+1}]+C_j \nonumber \end{align} ~~(3)f在x_j处泰勒展开得到f(x)=f(x_j)+f'(x_j)(x-x_j)+\frac{1}{2}f''(x_j)(x-x_j)^2+\cdots\therefore \begin{cases} f(x_{j+1})=f(x_j)+f'(x_j)(x_{j+1}-x_j)+\frac{1}{2}f''(x_j)(x_{j+1}-x_j)^2+\cdots \\ f(x_{j+2})=f(x_j)+f'(x_j)(x_{j+2}-x_j)+\frac{1}{2}f''(x_j)(x_{j+2}-x_j)^2+\cdots \\ f(x_{j-1})=f(x_j)+f'(x_j)(x_{j-1}-x_j)+\frac{1}{2}f''(x_j)(x_{j-1}-x_j)^2+\cdots \end{cases}令x_{j+1}-x_j=h\because x_{i+1}=x_i+1\therefore x_{j+2}-x_j=2h，x_{j-1}-x_j=-h\therefore \begin{cases} f(x_{j+2})=f(x_j)+2f'(x_j)h+2f''(x_j)h^2+\cdots \\ f(x_{j+1})=f(x_j)+f'(x_j)h+\frac{1}{2}f''(x_j)h^2+\cdots \\ f(x_{j-1})=f(x_j)-f'(x_j)h+\frac{1}{2}f''(x_j)h^2+\cdots \end{cases}\therefore \begin{cases} c_{j-1}=f(x_j)-f'(x_j)h+\frac{1}{2}f''(x_j)h^2+o(h^3) \\ c_j=f(x_j) \\ c_{j+1}=f(x_j)+f'(x_j)h+\frac{1}{2}f''(x_j)h^2+o(h^3)\\ c_{j+2}=f(x_j)+2f'(x_j)h+2f''(x_j)h^2+o(h^3) \end{cases} ~~ (4)将(4)代⼊(3)，得g(x)=-(2a+1)[2hf'(x_j)+h^2f''(x_j)]s^3+[(6a+3)hf'(x_j)+\frac{4a+3}{2}h^2f''(x_j)]s^2-2ahf'(x_j)s+f(x_j)+o(h^3)\because h=1,s=x-x_J\therefore sh=x-x_j\begin{align}\therefore f(x)&= f(x_j)+f'(x_j)(x-x_j)+\frac{1}{2}f''(x_j)(x-x_j)^2+o(h^3) \nonumber \\ &= f(x_j)+f'(x_j)sh+\frac{1}{2}f''(x_j)s^2h^2+o(h^3) \nonumber \end{align}\therefore f(x)-g(x)=(2a+1)[2hf'(x_j)+h^2f''(x_j)]s^3-(2a+1)[3hf'(x_j)+h^2f''(x_j)]s^2+[(2a+1)hf'(x_j)]s+o(h^3)\because 期望f(x)-g(x)趋于0\therefore 2a+1=0 \Rightarrow a=-\frac{1}{2}\therefore u(s)=\begin{cases} \frac{3}{2}|s|^3-\frac{5}{2}|s|^2+1, &0<|s|<1 \\ -\frac{1}{2}|s|^3+\frac{5}{2}|s|^2-4|s|+2, &1<|s|<2 \\ 1, &s=0 \\ 0, &otherwise \end{cases}\therefore g(s)=s^3[-\frac{1}{2}c_{j-1}+\frac{3}{2}c_j-\frac{3}{2}c_{j+1}+\frac{1}{2}c_{j+2}]+s^2[c_{j-1}-\frac{5}{2}c_j+2c_{j+1}-\frac{1}{2}c_{j+2}]+s[-\frac{1}{2}c_{j-1}+\frac{1} {2}c_{j+1}]+c_j图像插值p_{11}p_{12}p_{13}p_{14}p_{21}p_{22}p_{23}p_{24}pp_{31}p_{32}p_{33}p_{34}p_{41}p_{42}p_{43}p_{44}python代码1class BiCubicConvInterpolation:2 @staticmethod3def CubicConvInterpolation1(p0, p1, p2, p3, s):4# ⽤g（s）公式计算，已经将四个u（s）计算完毕并整理5# as^3 + bs^2 + cs + d6 a = 0.5 * (-p0 + 3.0 * p1 - 3.0 * p2 + p3)7 b = 0.5 * (2.0 * p0 - 5.0 * p1 + 4.0 * p2 - p3)8 c = 0.5 * (-p0 + p2)9 d = p110return d + s * (c + s * (b + s * a))1112 @staticmethod13def CubicConvInterpolation2(s):14# ⽤u（s）公式计算15 s = abs(s)16if s <= 1:17return 1.5 * s ** 3 - 2.5 * s ** 2 + 118elif s <= 2:19return -0.5 * s ** 3 + 2.5 * s ** 2 - 4 * s + 220else:21return 02223def biCubic1(self, srcImg, dstH, dstW):24# p11 p12 p13 p1425#26# p21 p22 p23 p2427# p28# p31 p32 p33 p3429#30# p41 p42 p43 p4431 dstH, dstW = int(dstH), int(dstW)32 scrH, scrW, _ = srcImg.shape33 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')34 dstImg = np.zeros((dstH, dstW, 1), dtype=np.uint8)35for dstY in range(dstH):36for dstX in range(dstW):37for channel in [0]:38 y = dstY * scrH / dstH39 x = dstX * scrW / dstW40 y1 = math.floor(y)41 x1 = math.floor(x)4243 array = []44for i in [-1, 0, 1, 2]:45 temp = self.CubicConvInterpolation1(srcImg[y1 + i, x1 - 1, channel],46 srcImg[y1 + i, x1, channel],47 srcImg[y1 + i, x1 + 1, channel],48 srcImg[y1 + i, x1 + 2, channel],49 x - x1)50 array.append(temp)5152 temp = self.CubicConvInterpolation1(array[0], array[1], array[2], array[3], y - y1)53 dstImg[dstY, dstX, channel] = np.clip(temp, 0, 255)5455return dstImg5657def biCubic2(self, srcImg, dstH, dstW):58# p11 p12 p13 p1459#60# p21 p22 p23 p2461# p62# p31 p32 p33 p3463#64# p41 p42 p43 p4465 dstH, dstW = int(dstH), int(dstW)66 scrH, scrW, _ = srcImg.shape67 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')68 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)69for dstY in range(dstH):70for dstX in range(dstW):71for channel in [0, 1, 2]:72 y = dstY * scrH / dstH73 x = dstX * scrW / dstW74 y1 = math.floor(y)75 x1 = math.floor(x)7677 array = []78for i in [-1, 0, 1, 2]:79 temp = 080for j in [-1, 0, 1, 2]:81 temp += srcImg[y1 + i, x1 + j, channel] * self.CubicConvInterpolation2(x - (x1 + j))82 array.append(temp)8384 temp = 085for i in [-1, 0, 1, 2]:86 temp += array[i + 1] * self.CubicConvInterpolation2(y - (y1 + i))87 dstImg[dstY, dstX, channel] = np.clip(temp, 0, 255)8889return dstImg三次样条插值在n-1个区间上寻找n-1个三次曲线，使其满⾜相邻曲线在端点处值相等、⼀阶导数相等，⼆阶导数相等，在加以边界条件后可得每个曲线的⽅程，然后沿x轴依次偏移对应的距离即可得到插值结果，如仅需要特定范围内的结果，则可以⼤幅减少计算量公式推导设S_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^3, ~~~~i=0,1,...,n-1则 \begin{cases} S_i'(x)=b_i+2c_i(x-x_i)+3d_i(x-x_i)^2\\ S_i''(x)=2c_i+6d_i(x-x_i)\\ S_i'''(x)=6d_i\\ \end{cases} ~~~~i=0,1,...,n-1设h_i(x)=x_{i+1}-x_i，可得\begin{cases} S_i(x)=a_i+b_ih_i+c_ih_i^2+d_ih_i^3\\ S_i'(x)=b_i+2c_ih_i+3d_ih_i^2\\ S_i''(x)=2c_i+6d_ih_i\\ S_i'''(x)=6d_i\\ \end{cases} ~~~~i=0,1,...,n-1\because S_i(x)过点(x_i,y_i)\therefore S_i(x)=a_i=y+i, ~~~~i=0,1,...,n-1 ~~~~~~(1)\because S_i(x)与S_{i+1}(x)在X_{i+1}处相等\therefore S_i(x_{i+1})=S_{i+1}(x_{i+1})\Rightarrow a_i+b_ih_i+c_ih_i^2+d_ih_i^3=y_{i+1}, ~~~~i=0,1,...,n-2~~~~~~(2)\because S_i'(x)与S_{i+1}'(x)在X_{i+1}处相等\therefore S_i'(x)-S_{i+1}'(x)=0\Rightarrow b_i+2c_ih_i+3d_ih_i^2-b_{i+1}=0~~~~~~(3)\because S_i''(x)与S_{i+1}''(x)在X_{i+1}处相等\therefore S_i''(x)-S_{i+1}''(x)=0\Rightarrow 2c_i+6d_ih_i-2c_{i+1}=0, ~~~~i=0,1,...,n-2~~~~~~(4)设m_i=S_i(x_i)=2c_i，即c_i=\frac{1}{2}m_i, ~~~~i=0,1,...,n-1~~~~~~(5)将(5)代⼊(4)，得2c_i+6d_ih_i-2c_{i+1}=0\Rightarrow m_i+6h_id_i-m_{i+1}=0\Rightarrow d_i=\frac{m_{i+1}-m_i}{6h_i}, ~~~~i=0,1,...,n-2~~~~~~(6)将(1)(5)(6)代⼊(2)，得\begin{align} &a_i+b_ih_i+c_ih_i^2+d_ih_i^3=y_{i+1} \nonumber \\ \Rightarrow&y_i+b_ih_i+\frac{1}{2}m_ih_i^2+\frac{m_{i+1}-m_i}{6h_i}h_i^3=y_{i+1} \nonumber \\\Rightarrow&b_i=\frac{y_{i+1}-y_i}{h_i}-\frac{1}{2}m_ih_i-\frac{1}{6}(m_{i+1}-m_i)h_i \nonumber \\ \Rightarrow&b_i=\frac{y_{i+1}-y_i}{h_i}-\frac{1}{3}m_ih_i-\frac{1}{6}m_{i+1}h_i, ~~~~i=0,1,...,n-2~~~~~~(7) \nonumber \end{align}将(5)(6)(7)代⼊(3)，得\begin{align} &\frac{y_{i+1}-y{i}}{h_i}-\frac{1}{3}m_ih_i-\frac{1}{6}m_{i+1}h_i+2\cdot\frac{1}{2}m_ih_i+3\frac{m_{i+1}-m_i}{6h_i}h_i^2-(\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{1}{3}m_{i+1}h_{i+1}-\frac{1}{6}m_{i+2}h_{i+1})=0 \nonumber \\ \Rightarrow&\frac{y_{i+1}-y{i}}{h_i}-\frac{1}{3}m_ih_i-\frac{1}{6}m_{i+1}h_i+m_ih_i+\frac{1}{2}(m_{i+1}-m_i)h_i-\frac{y_{i+2}-y_{i+1}}{h_{i+1}}+\frac{1}{3}m_{i+1}h_{i+1}+\frac{1}{6}m_{i+2}h_{i+1}=0 \nonumber \\ \Rightarrow&m_ih_i(-\frac{1}{3}+1-\frac{1}{2})+m_{i+1}h_i(-\frac{1}{6}+\frac{1} {2})+\frac{1}{3}m_{i+1}h_{i+1}+\frac{1}{6}m_{i+2}h_{i+1}=\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_{i}}{h_{i}} \nonumber \\ \Rightarrow&\frac{1}{6}(m_ih_i+2m_{i+1}h_i+2m_{i+1}h_{i+1}+m_{i+2}h_{i+1})=\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_{i}}{h_{i}} \nonumber \\ \Rightarrow&m_ih_i+2m_{i+1}(h_i+h_{i+1})+m_{i+2}h_{i+1}=6(\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_{i}}{h_{i}}), ~~~~i=0,1,...,n-2~~~~~~(8) \nonumber \end{align}由(8)可知i=0,1,...,n-2，则有m_0,m_1,...,m_n，需要两个额外条件⽅程组才有解⾃然边界(Natural)m_0=0，m_n=0\begin{bmatrix} \tiny 1 & 0 & 0 & 0 & 0 & \cdots & 0\\ h_0 & 2(h_0+h_1) & h_1 & 0 & 0 & \cdots & 0\\ 0 & h_1 & 2(h_1+h_2) & h_2 & 0 & \cdots & 0\\ 0 & 0 & h_2 & 2(h_2+h_3) & h_3 & \cdots & 0\\ \vdots& & & \ddots & \ddots & \ddots & \vdots\\ 0 & \cdots & & & h_{n-2} & 2(h_{n-2}+h_{n-1}) & h_{n-1}\\ 0 & \cdots & & & 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3\\\vdots\\m_{n-1}\\m_n \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ \frac{y_4-y_3}{h_3}-\frac{y_3-y_2}{h_2}\\ \vdots\\ \frac{y_n-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\\ 0 \end{bmatrix}固定边界(Clamped)\begin{align} &\begin{cases} S_0'(x_0)=A\\ S_{n-1}'(x_n)=B \end{cases} \nonumber \\ \Rightarrow&\begin{cases} b_0=A\\ b_{n-1}+2c_{n-1}h_{n-1}+3d_{n-1}h_{n-1}^2=B\end{cases} \nonumber \\ \Rightarrow&\begin{cases} A=\frac{y_1-y_0}{h_0}-\frac{h_0}{2}m_0-\frac{h_0}{6}(m_1-m_0)\\ B=\frac{y_n-y_{n-1}}{h_{n-1}}-\frac{1}{3}m_{n-1}h_{n-1}+m_{n-1}h_{n-1}+\frac{1}{2}m_nh_{n-1}-\frac{1}{2}m_{n-1}h_{n-1} \end{cases} \nonumber \\ \Rightarrow&\begin{cases} 2h_0m_0+h_0m_1=6(\frac{y_1-y_0}{h_0}-A)\\ h_{n-1}m_{n-1}+2h_{n-1}m_{n}=6(B-\frac{y_n-y_{n-1}}{h_{n-1}}) \end{cases} \nonumber \\ \end{align}\begin{bmatrix} 2 & 1 & 0 & 0 & 0 & \cdots & 0\\ h_0 & 2(h_0+h_1) & h_1 & 0 & 0 & \cdots & 0\\ 0 & h_1 & 2(h_1+h_2) & h_2 & 0 & \cdots & 0\\ 0 & 0 & h_2 & 2(h_2+h_3) & h_3 & \cdots & 0\\ \vdots& & & \ddots & \ddots & \ddots & \vdots\\ 0 & \cdots & & & h_{n-2} & 2(h_{n-2}+h_{n-1}) & h_{n-1}\\ 0 & \cdots & & & 0 & 1 & 2 \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3\\\vdots\\m_{n-1}\\m_n \end{bmatrix}=6\begin{bmatrix} \frac{y_1-y_0}{h_0}-A\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ \frac{y_4-y_3}{h_3}-\frac{y_3-y_2}{h_2}\\ \vdots\\\frac{y_n-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\\ B-\frac{y_n-y_{n-1}}{h_{n-1}} \end{bmatrix}⾮节点边界(Not-A-Knot)\begin{align} &\begin{cases} S_0'''(x_1)=S_1'''(x_1)\\ S_{n-2}'''(x_{n-1})=S_{n-1}'''(x_{n-1}) \end{cases} \nonumber \\ \Rightarrow&\begin{cases} 6\cdot\frac{m_1-m_0}{6h_0}=6\cdot\frac{m_2-m_1}{6h_1}\\ 6\cdot\frac{m_{n-1}-m_{n-2}}{6h_{n-2}}=6\cdot\frac{m_n-m_{n-1}}{6h_{n-1}} \end{cases} \nonumber \\ \Rightarrow&\begin{cases} h_1(m_1-m_0)=h_0(m_2-m_1)\\ h_{n-1}(m_{n-1}-m_{n-2})=h_{n-2}(m_n-m_{n-1}) \end{cases} \nonumber \\ \Rightarrow&\begin{cases} -h_1m_0+(h_1+h_0)m_1-h_0m_2=0\\ -h_{n-1}m_{n-2}+(h_{n-1}+h_{n-2})m_{n-1}-h_{n-2}m_n=0 \end{cases} \nonumber \\ \end{align}\begin{bmatrix} -h_1 & h_0+h_1 & -h_0 & 0 & 0 & \cdots & 0\\ h_0 & 2(h_0+h_1) & h_1 & 0 & 0 & \cdots & 0\\ 0 & h_1 & 2(h_1+h_2) & h_2 & 0 & \cdots & 0\\ 0 & 0 & h_2 &2(h_2+h_3) & h_3 & \cdots & 0\\ \vdots& & & \ddots & \ddots & \ddots & \vdots\\ 0 & \cdots & & & h_{n-2} & 2(h_{n-2}+h_{n-1}) & h_{n-1}\\ 0 & \cdots & & & -h_{n-1} & h_{n-1}+h_{n-2} & -h_{n-2} \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3\\\vdots\\m_{n-1}\\m_n \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ \frac{y_4-y_3}{h_3}-\frac{y_3-y_2}{h_2}\\ \vdots\\ \frac{y_n-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\\ 0 \end{bmatrix}在n=4时通⽤公式⾃然边界\begin{bmatrix} 1 & 0 & 0 & 0 \\ h_0 & 2(h_0+h_1) & h_1 & 0 \\ 0 & h_1 & 2(h_1+h_2) & h_2 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ 0 \end{bmatrix}固定边界\begin{bmatrix} 2 & 1 & 0 & 0 \\ h_0 & 2(h_0+h_1) & h_1 & 0 \\ 0 & h_1 & 2(h_1+h_2) & h_2 \\ 0 & 0 & 1 & 2 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} \frac{y_1-y_0}{h_0}-A\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ B-\frac{y_3-y_2}{h_2} \end{bmatrix}⾮节点边界\begin{bmatrix} -h_1 & h_0+h_1 & -h_0 & 0 \\ h_0 & 2(h_0+h_1) & h_1 & 0 \\ 0 & h_1 & 2(h_1+h_2) & h_2 \\ 0 & -h_2 & h_1+h_2 & -h_1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ 0 \end{bmatrix}图像插值x_{i+1}-x_i=1 \Rightarrow h_i(x)=1在n=4时，即四个点时如下所⽰p_{11}p_{12}p_{13}p_{14}p_{21}p_{22}p_{23}p_{24}pp_{31}p_{32}p_{33}p_{34}p_{41}p_{42}p_{43}p_{44}⾃然边界（可⽤TDMA或化简计算）\begin{bmatrix} 1 & 0 & 0 & 0 \\ 1 & 4 & 1 & 0 \\ 0 & 1 & 4 & 1 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ y_0+y_2-2y_1\\ y_1+y_3-2y_2\\ 0 \end{bmatrix}固定边界（只能⽤TDMA计算）\begin{bmatrix} 2 & 1 & 0 & 0 \\ 1 & 4 & 1 & 0 \\ 0 & 1 & 4 & 1 \\ 0 & 0 & 1 & 2 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} y_1-y_0-A\\ y_0+y_2-2y_1\\ y_1+y_3-2y_2\\ y_2-y_3+B \end{bmatrix}⾮节点边界（只能化简计算）\begin{bmatrix} -1 & 2 & -1 & 0 \\ 1 & 4 & 1 & 0 \\ 0 & 1 & 4 & 1 \\ 0 & -1 & 2 & -1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ y_0+y_2-2y_1\\ y_1+y_3-2y_2\\ 0 \end{bmatrix}python代码1class BiSplineInterpolation:2 @staticmethod3 def TDMA(a, b, c, d):4 n = len(d)56 c[0] = c[0] / b[0]7 d[0] = d[0] / b[0]89for i in range(1, n):10 coef = 1.0 / (b[i] - a[i] * c[i - 1])11 c[i] = coef * c[i]12 d[i] = coef * (d[i] - a[i] * d[i - 1])1314for i in range(n - 2, -1, -1):15 d[i] = d[i] - c[i] * d[i + 1]1617return d1819 @staticmethod20 def Simplified_Natural4(y1, y2, y3, y4):21 # 四点⾃然边界化简公式22 d1 = y1 + y3 - 2 * y223 d2 = y2 + y4 - 2 * y32425 k0 = 026 k1 = (4 * d1 - d2) * 0.427 k2 = (4 * d2 - d1) * 0.428 k3 = 02930return [k0, k1, k2, k3]3132 @staticmethod33 def Simplified_Not_A_Knot4(y1, y2, y3, y4):34 # 四点⾮节点边界化简公式35 d1 = y1 + y3 - 2 * y236 d2 = y2 + y4 - 2 * y33738 k0 = 2 * d1 - d239 k1 = d140 k2 = d241 k3 = 2 * d2 - d14243return [k0, k1, k2, k3]4445 # TDMA矩阵说明46 # a0 和 c3 没有实际意义，占位⽤47 # a0 [b0 c0 00 ] [x0] [d0]48 # [a1 b1 c1 0 ] [x1] = [d1]49 # [0 a2 b2 c2] [x2] [d2]50 # [00 a3 b3] c3 [x3] [d3]5152 def SplineInterpolationNatural4(self, x, y1, y2, y3, y4):53 # ⽤TDMA计算54 # matrix_a = [0, 1, 1, 0]55 # matrix_b = [1, 4, 4, 1]56 # matrix_c = [0, 1, 1, 0]57 # matrix_d = [0, 6 * (y1 + y3 - 2 * y2), 6 * (y2 + y4 - 2 * y3), 0]58 # matrix_x = self.TDMA(matrix_a, matrix_b, matrix_c, matrix_d)5960 # 化简计算61 matrix_x = self.Simplified_Natural4(y1, y2, y3, y4)6263 a = y264 b = y3 - y2 - matrix_x[1] / 3.0 - matrix_x[2] / 6.065 c = matrix_x[1] / 2.066 d = (matrix_x[2] - matrix_x[1]) / 6.06768 s = a + b * x + c * x * x + d * x * x * x69return s7071 def SplineInterpolationClamped4(self, x, y1, y2, y3, y4):72 # 仅有TDMA计算，⽆法化简73 A, B = 1, 17475 matrix_a = [0, 1, 1, 1]76 matrix_b = [2, 4, 4, 2]77 matrix_c = [1, 1, 1, 0]78 matrix_d = [6 * (y2 - y1 - A), 6 * (y1 + y3 - 2 * y2), 6 * (y2 + y4 - 2 * y3), 6 * (B - y4 + y3)]79 matrix_x = self.TDMA(matrix_a, matrix_b, matrix_c, matrix_d)8081 a = y282 b = y3 - y2 - matrix_x[1] / 3.0 - matrix_x[2] / 6.083 c = matrix_x[1] / 2.084 d = (matrix_x[2] - matrix_x[1]) / 6.08586 s = a + b * x + c * x * x + d * x * x * x87return s8889 def SplineInterpolationNotAKnot4(self, x, y1, y2, y3, y4):90 # ⽆法使⽤TDMA计算91 matrix_x = self.Simplified_Not_A_Knot4(y1, y2, y3, y4)9293 a = y294 b = y3 - y2 - matrix_x[1] / 3.0 - matrix_x[2] / 6.095 c = matrix_x[1] / 2.096 d = (matrix_x[2] - matrix_x[1]) / 6.09798 s = a + b * x + c * x * x + d * x * x * x99return s100101 def biSpline4(self, srcImg, dstH, dstW):102 dstH, dstW = int(dstH), int(dstW)103 srcH, srcW, _ = srcImg.shape104 srcImg = np.pad(srcImg, ((1, 2), (1, 2), (0, 0)), 'edge')105 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)106for dstY in range(dstH):107for dstX in range(dstW):108for channel in [0, 1, 2]:109 # p11 p12 p13 p14110 #111 # p21 p22 p23 p24112 # p113 # p31 p32 p33 p34114 #115 # p41 p42 p43 p44116 # 储存为 p(y, x)117 p = [dstY * srcH / dstH, dstX * srcW / dstW]118 p22 = [math.floor(p[0]), math.floor(p[1])]119 p21 = [p22[0], p22[1] - 1]120 p23 = [p22[0], p22[1] + 1]121 p24 = [p22[0], p22[1] + 2]122123 p11 = [p21[0] - 1, p21[1]]124 p12 = [p11[0], p22[1]]125 p13 = [p11[0], p23[1]]126 p14 = [p11[0], p24[1]]127128 p31 = [p21[0] + 1, p21[1]]129 p32 = [p31[0], p22[1]]130 p33 = [p31[0], p23[1]]131 p34 = [p31[0], p24[1]]132133 p41 = [p21[0] + 2, p21[1]]134 p42 = [p41[0], p22[1]]135 p43 = [p41[0], p23[1]]136 p44 = [p41[0], p24[1]]137138 diff_y, diff_x = p[0] - p22[0], p[1] - p22[1]139 r1 = self.SplineInterpolationNatural4(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel], srcImg[p13[0], p13[1], channel], srcImg[p14[0], p14[1], channel]) 140 r2 = self.SplineInterpolationNatural4(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel], srcImg[p23[0], p23[1], channel], srcImg[p24[0], p24[1], channel]) 141 r3 = self.SplineInterpolationNatural4(diff_x, srcImg[p31[0], p31[1], channel], srcImg[p32[0], p32[1], channel], srcImg[p33[0], p33[1], channel], srcImg[p34[0], p34[1], channel]) 142 r4 = self.SplineInterpolationNatural4(diff_x, srcImg[p41[0], p41[1], channel], srcImg[p42[0], p42[1], channel], srcImg[p43[0], p43[1], channel], srcImg[p44[0], p44[1], channel]) 143144 c = self.SplineInterpolationNatural4(diff_y, r1, r2, r3, r4)145146 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)。

第一章 插值方法求得近似函数的方法通常有两种：一种是插值法，另一种是曲线拟合法.本章主要讨论插值法用多项式来逼近列表函数问题，即具有唯一插值函数的多项式插值,对其中的多项式插值主要讨论n 次多项式插值的方法，即给定n+1各点处的函数值后，怎样构造一个n 次插值多项式的方法.§1.1插值问题的由来在许多实际问题及科学研究中，因素之间往往存在函数关系.然而这种关系经常很难有明显的解析表达式，通常只是由观察与测试得到一些离散数值.有时由于表达式过于复杂，不仅使用不便，而且不易于进行计算与理论分析.解决此类问题的方法有两种：一种是插值法，另一种是拟合法.插值法是一种古老的数学方法，他来自生产实践，早在一千多年前，我国科学家在研究历法上就应用了线性插值与二次插值，但它的基本理论却是在微积分产生之后才逐渐完善的，其应用也逐步增多.特别是在计算机软件中，许三多库函数.等的计算实际上归结于它的逼近函数的计算.在工程的实际问题中，由于航空、造船、精密机械加工等实际问题的需要，插值方法在实践上和理论上显得更为重要，并得到了进一步的发展和广泛的应用.因此我们希望能够得到一个“简单函数”逼近被计算函数的函数值.于是就得到了一种方法.叫做插值逼近或者插值法.§1.2 几个基本概念§1.2.1 插值函数的概念及相关性质定义 1.1 设连续函数)(x f y -在区间],[b a 上有定义，已知在1+n 个互异的点n x x x ,,,10 上取值分别为n y y y ,,10（设b x x x a n ≤<<<≤ 10）.若在函数类中存在以简单函数)(x p ，使得i i y x p =)(（n i ,1,0=）,则称)(x p 为)(x f 的插值函数. 称n x x x ,,,10 为插值节点，称],[b a 为插值区间. 定义1.2 设如下表i x0x 1xn x iyy1yny求一个)(x p ，使其满足： ○1)()(x IP x p n ∈○2i i i y x f x P ==)()(， n i ,,1,0 =则称)(x p 为)(x f 的n 次代数插值多项式，该问题称为n 次代数插值问题.定理1.1 n 次代数插值问题的解存在且唯一 .§1.2.2 Lagrange 插值函数的概念及相关性质定义1.3 若n 次多项式)(x l k （n k ,1,0=）在1+n 个互异节点n x x x ,,,10 上满足),,1,0,(,0,1)(n j k kj kj x l kj j k =⎩⎨⎧≠===δ则称)(x l k （n k ,,1,0 =）为节点n x x x ,,,10 上的n 次插值基函数. n 次差值多项式)(x L n 为)()(0x l y x L k nk k n ∑==称之为拉格朗日多项式.定理1.2 设],[1b a C f n +∈，则余项为：)()!1()()(11x w n f x R n n n +++=ξ其中)())(()(101n n x x x x x x x w ---=+ .§1.2.3 差商的概念及相关性质定义1.4 称函数)(x f 于点j i x x ,（j i x x ≠）上的平均变化率为)(x f 在j i x x ,处的一阶差商，记作],[j i x x f ，即ij i j j i x x x f x f x x f --=)()(],[，j i x x ≠.一阶差商的差商为)(x f 在点k j i x x x ,,处的二阶差商.记作],,[k j i x x x f ，即ik j i k j k j i x x x x f x x f x x x f --=],[],[],,[一般地，在求出)(x f 的1-m 阶差商后就可以构造)(x f 的m 阶差商.称10110],,[],[],,,[i im im i im i im i i x x x x f x x f x x x f --=-为)(x f 在im i i x x x ,,10处的m 阶差商.特别地，规定零阶差商)(][i i x f x f =.性质1.1 任意调换插值节点的次序其值不变 即],,[],,[],,[120201210x x x f x x x f x x x f ==性质1.2 差商是一个数，可表示为某个函数的线性组合.性质1.3 当)(x f 的m 阶导数)()(x f m 在包含节点im i i x x x ,,10的某个区间上存在时，在 im i i x x x ,,10之间至少存在一点ξ使!)(],,,[)(10m f x x x f m im i i ξ=注1.1 由性质2可以推出，n 次多项式的一阶差商是1-n 多项式.表1-1差商表k x )(kx f 一阶差商 二阶差商三阶差商四阶差商x)(kx f1x )(kx f ],[10x x f2x)(kxf],[21x x f ],,[21x x x f3x )(k x f ],[32x x f ],,[321x x x f ],,,[321x x x x f4x)(kx f],[43x x f],,[432x x x f],,[321x x x f],,,,[4321x x x x x f§1.2.3 Newton 插值函数与差分的概念及相关性质 定义1.5 我们定义牛顿插值公式为)())(](,,,[))(](,,[)](,[)()(1111211----++--+-+=n nnx x x x x x x x x f x x x x x x x f x x x x f x f x N （1.1）定义1.6 设函数)(x f y =在等距节点kh x x k +=0),,1,0(n k =上的函数值)(k k x f y = ),,1,0(n k =，其中h 为常数，称为步长.称函数在每个小区间上的增量k k y y -+1为函数)(x f 在k x 处以h 为步长的一阶向前差分，记作k y ∆，即k y ∆k k y y -=+1称一阶差分的差分1+∆k y -k y ∆为二阶差分，记作k y 2∆，即k k k k k k y y y y y y +-=∆-∆=∆+++12122由此可得1-m 阶差分的差分k m k m k m y y y 111-+-∆-∆=∆为)(x f 在k x 点处的m 阶差分.称1--k k y y 为函数)(x f 在k x 处以h 为步长的一阶向后差分，记作k y ∇.即1--=∇k k k y y y称111---∇-∇k m k m y y 为函数)(x f 在k x 处以h 为步长的m 阶向后差分，记作k m y ∇. 即111---∇-∇=∇k m k m k m y y y表1-2 向前差分表ix i y i y ∆y 2∆y 3∆y 3∆x1y1x 1y 0y ∆2x2y1y ∆02y ∆3x 3y 2y ∆ 12y ∆03y ∆4x4y2y ∆22y ∆13y ∆3y ∆表1-3向后差分表ix i y i y ∇y 2∇y 3∇y 3∇x1y1x 1y 0y ∇2x2y1y ∇02y ∇3x 3y 2y ∇ 12y ∇ 03y ∇4x4y2y ∇22y ∇13y ∇3y ∇性质1.4 各阶差分可以表示成函数值的线性组合.即(1)k n j knnk k j ny C y -+=∑-=∆0)1( (2)k j knn k k j ny C y -=∑-=∇0)1( 性质1.5 差商与差分有如下关系n n nn n h n f h n f x x x f !!],,,[0010∇=∆=定义1.7 设1+n 个等距插值节点的顺序为n n x x x x ,,,110- ，则由差商表示的牛顿插值公式为)()](,,,[))(](,,[)](,[][)(111211---++--+-+=n nnx x x x x x x f x x x x x x x f x x x x f x f x N又设插值节点th x x +=0 )0(n t <<由性质1.5差商与差分关系kk k h k f x x x f !],,,[010∆=可得∑=∆+--=∆+--++∆-+∆+=nk knn f k k t t t f n n t t t f t t f t f x N 00200!)1()1(!)1()1(!2)1()( （1.2） 称为牛顿前插公式.其余项公式为,)()1()!1()()()(110++--+=+=n n n n h n t t t n f th x R x R ξ],[0n x x ∈ξ定义1.8 设1+n 个等距插值节点的顺序为011,,,x x x x n n -，则由差商表示的牛顿插值公式为)()](,,,[))()(](,,[)](,[][)(10121211x x x x x x x f x x x x x x x x x f x x x x f x f x N n n n n n n n n n n n n n n --++---+-+=------又设插值节点th x x n -= )0(n t <<由性质1.5差商与差分关系knk k n n n hk f x x x f !],,,[1∇=-- 可得∑=∆+---=∆+---+-∇-+∇-=nk nkk nnn n n n n f k k t t t f n n t t t f t t f t f x N 02!)1()1()1(!)1()1()1(!2)1()( （1.3） 称为牛顿后插公式.其余项公式为,)1)(()1()!1()()()(111+++---+=-=n n n n n n h n t t t n f th x R x R ξ],[0n x x ∈ξ第二章 Newton 插值法的程序设计与应用利用插值基函数很容易得到Lagrange 插值多项式，公式结构紧凑，在理论分析中也很方便，但当插值节点增减时全部插值基函数)(x l i n i ,,1,0 =都要随之变化，整个公式也将发生变化，这在实际计算中是很不方便的，于是我们导出插值多项式的又一种表示形式——Newton 插值公式.远远减少了工作量.但由于通过相同结点的差值多项式是惟一的，因此，拉格朗日差值公式与牛顿公式本质上是一个多项式，只不过形式不同而已.为了更方便的应用Newton 插法,本章主要讨论Newton 插值法的程序设计,并对其应用进行讨论.§2.1 Newton 插值法的算法与程序设计本节主要讨论n 次多项式插值的方法，即给定n+1各点处的函数值后，怎样构造一个n 次插值多项式的方法.为了更方便的计算并与计算机结合，我们便作出了相应的程序.并利用Matlab 进行绘图.§2.1.1 Newton 插值法的应用步骤步骤 首先我们按照表1-1，求得各点的差商.然后利用Newton 前差或后差公式即（公式1.3)把数值带入.即可以求得n 次多项式. 它在计算机上的应用步骤如下：步骤1 输入所要求的牛顿多项式的次数n ，并依次输入1+n 个节点),(i i y x .for : i=0 to n+1 {scanf("%f",&x[i]) scanf("%f",&y[0][i]); }步骤2 计算各阶差商 for : i=1 to n+1 {for : j=i to n+1 if(i>1)y[i][j]=(y[i-1][j]-y[i-1][j-1])/(x[j]-x[j-i]); elsey[i][j]=(y[i-1][j]-y[i-1][j-1])/(x[j]-x[j-1]);}步骤3 代入牛顿插值公式,可计算得出结果. for ：i=1 to n+1 {temp=temp*(xx-x[i-1]);newton=newton+y[i][i]*temp;}由于篇幅所限,流程图、C 语言程序和Matlab 程序详见附录（附录A 、B 、C ）.§2.1.2 利用Newton 插值法程序的实例为了更方便的应用牛顿插值法，我们进行了与计算机的结合，下面我们将展示几个例子.例2.1 已知x y sin =的一组数据为x6π4π3π2πsinx21 22 23 1（1）构造牛顿插值函数并作图分析. （2）并分别利用程序估计12sinπ，1514sinπ的估计值.分析 首先我们可以通过程序求出差商表：带入定义1.5中可求得如下牛顿插值多项式如下)6(2086076.095492926.0)(2π--=x x x x N (2.1)kx)(kx f一阶差商 二阶差商三阶差商四阶差商6π 21 954929268.04π 22 791089631.0086076.2-3π 23607024424.0 35153865.0- 13648909.0-2π125587263.044710035.0-09125470.0-028797106.0)4)(6(13648909.0)6(2086076.095492926.0)(3πππ-----=x x x x x x x N (2.2))3)(4)(6(028797106.0)4)(6(13648909.0)6(2086076.095492926.0)(4ππππππ---+-----=x x x x x x x x x x x N (2.3)第二步 我们分别利用C++程序和Matlab 分别计算12sinπ和1514sinπ的值.步骤如下：○1利用C++程序：首先我们输入所要求的牛顿多项式的次数n ，然后输入n+1个节点的值.即可以得出12sinπ和1514sin π的值为0.2586和0.0804； ○2利用Matlab 程序：首先我们在Matlab 中新建一个文件，文件内容即为附录C 中Matlab 程序代码，然后输入以下步骤：x=[0,pi/6,pi/4,pi/3,pi/2];y=[0,0.5,sqrt(2)./2,sqrt(3)./2,1]; x0=pi/12;Newton(x,y,x0] ans =0.2586x=[0,pi/6,pi/4,pi/3,pi/2];y=[0,0.5,sqrt(2)./2,sqrt(3)./2,1]; x0=14.*pi/15; Newton(x,y,x0) ans =0.0804即结果为0.2586和0.0804利用Matlab 分别对求得的牛顿差值多项式作图(代码详见附录D)图形如下.图2-1例2.2 设)1ln()(x x f +=的函数表如下：x0.25 0.30 0.36 0.39 0.45 )(x f0.2231440.2623640.3074850.3293040.371564试计算)275.1ln(，3.2ln )(=x f分析：同上题步骤我们先求差商表，进而代入公式可得)39.0)(36.0)(30.0)(25.0(057720032.0)36.0)(30.0)(25.0(1411736357.0)30.0)(25.0(294393939.0)25.0(7844.022314.0)(4----+---+----+=x x x x x x x x x x x N （ 2.4）)36.0)(30.0)(25.0(1411736357.0)30.0)(25.0(294393939.0)25.0(7844.022314.0)(3---+----+=x x x x x x x N （2.5）)30.0)(25.0(294393939.0)25.0(7844.022314.0)(2----+=x x x x N （2.6）利用C++程序我们可以得到计算结果 ：242946.0)275.0(4=N ，242945.0)275.0(3=N ，242938.0)275.0(2=N 928823.0)3.1(4=N ，877002.0)3.1(3=N ，737653.0)275.0(2=N 我们利用Matlab 进行绘图.图2-2从上述两个例子我们可以看出,多项式（2.1）在区间]4,0[π周围与原函数逼近的较好.离这个区间越远与原函数的误差越大在5.1=x 处时，该图像就已经开始背离图像了.所以该多项式只能在一个小的区间里可以逼近原函数,不适合作为原函数的逼近函数.也可以看出多项式（2.2）在]3,0[π区间的周围逼近的较好，但是2=x 处时，该图像就离原图像误差较大.多项式（2.3）在区间[0,2.5]都逼近的挺好，从图中我们看出在远离这个区间的图像误差相对较大，但是在这三个多项式中是逼近最好的,同样在例2.2中也是如此.于是我们我们可以得出节点越多，函数逼近的相对较好.在节点附近逼近的越好，越远离节点误差越大，所以公式适用于计算节点附近的值于是为了减小误差，在下一节的等距节点下的插值公式根据所求的点的函数值的不同分别采取了前插和后插公式.§2.2 等距节点下的Newton 插值算法与程序设计前面我们讲述了一般节点下的牛顿插值公式，为了计算方便于是有了对等距节点下的牛顿多项式的研究，本节将对等距节点下的插值多项式进行总结讨论.§2.2.1 等距节点下的Newton 插值法的程序算法步骤步骤 我们按照表1-2或1-3，求差分.然后利用Newton 前插公式（公式1.2)或Newton 后插公式（公式1.3)并把数值带入.即可以求得n 次多项式. 它在计算机上的应用步骤如下：步骤1 输入所要求的牛顿多项式的次数n ，步长h ，并依次输入1+n 个节点),(i i y x .for : i=0 to n+1 {scanf("%f",&x[i]) scanf("%f",&y[0][i]); }步骤2 求得各界差分 for : i=1 to n+1{ for : j=i to n+1y[i][j]=(y[i-1][j+1]-y[i-1][j-1])/(x[j]-x[j]); }//求向前差分for : i=1 to n+1{ for : j=i to n+1y[i][j]=(y[i-1][j]-y[i-1][j-1])/(x[j]-x[j-i]); }//求向后差分步骤3 代入牛顿插值公式,可计算得出结果 for(i=1;i<n+1;i++) {temp=temp*((t-i+1)/i);newton=newton+y[i][i]*temp; }printf("求得的结果为:N(%.4f)=%9f\n",xx,newton); }//求得运用前插公式的值for(i=n;i<0;i--) {temp=temp*((t-i+1)/i);newton=newton+y[n-i][n-i]*temp; }printf("求得的结果为:N(%.4f)=%9f\n",xx,newton); }//求得运用后插公式的值由于篇幅所限 C 语言程序(详见附录D 、E).§2.2.2 等距节点下Newton 插值的实例例2.4 已知x tan 的值列表如下： x 1. 3 1.31 1.32 1.33 x tan 3.6021 3.7471 3.90334.0723近似计算325.1tan ,305.1tan .采用Newton 向后插公式.为此，做差分表i x i f i f ∇ i f 2∇ i f 3∇ 1.3 3.6021 0.1450 0.0112 0.0016 1.31 3.7471 0.1562 0.0128 1.32 3.9033 0.1690 1.334.0723从而，有).2)(1(!30016.0)1(!20128.01690.00723.4)(3++++++=t t t t t t x p 将5.001.033.1325.1-=-=t 代入上式，得9869.3)325.1(325.1tan 3=≈p .将5.201.033.1305.1-=-=t 带入后插多项式中可以得到0797.3)305.1(305.1tan 3=≈p现在我们利用C++程序 首先输入所求插值的次数3和步长0.01.然后输入各个节点，并输入所要求的点1.325既可以求出该点的函数值.即9869.3)325.1(=f . 0797.3)26.1()26.1(3=≈p f下面我们可得图图2-4 例2.4 设函数在各节点的取值如下x0 0.2 0.4 0.6 0.8 1.0 )(x f 1.0 0.818731 0.670320 0.548812 0.449329 0.367879 试利用插值公式求)3.0(f的值.采用Newton向前插公式，同上题我们先做差分表，然后相应带入到差分公式.中求得后插公式.利用C语言程序步骤如下:首先输入所求插值的次数5和步长0.2.然后输入各个节点，并输入所要求的点0.3既可以求出该点的函数值.即7880f.)3.0(.0由以上例子我们看到例1用了牛顿后插公式，例2用了牛顿前插公式，我们该怎样选取.这个经过验证得出，如果所要求的点较靠近节点x，则采用前插公式；如果1x，则采用牛顿后插公式.靠近n附录A ：牛顿插值法的流程图N YNYY NNY开始 输入ni=0i=i+1分别输入n+1个输入各节点 j=j+1 i=1;t=1;N=00y j=ii<n j<n? i>1ij j ij x x y --=1-j 1,-i 1j -i y -y 11-j 1,-i 1j -i y -y --=j j ij x x y 输入ij y i<n+1 t=t*(x-1-i x )N=N+i i ty i=i+1输出n 的值附录B：牛顿插值法的程序#include<stdio.h>void main(){float x[11],y[11][11],xx,temp,newton;int i,j,n;printf("Newton插值:\n请输入要运算的值:x=");scanf("%f",&xx);printf("请输入插值的次数(n<11):n=");scanf("%d",&n);printf("请输入%d组值:\n",n+1);for(i=0;i<n+1;i++){ printf("x%d=",i);scanf("%f",&x[i]);printf("y%d=",i);scanf("%f",&y[0][i]);}for(i=1;i<n+1;i++)for(j=i;j<n+1;j++){ if(i>1)y[i][j]=(y[i-1][j]-y[i-1][j-1])/(x[j]-x[j-i]);elsey[i][j]=(y[i-1][j]-y[i-1][j-1])/(x[j]-x[j-1]); printf("%f\n",y[i][i]);}temp=1;newton=y[0][0];for(i=1;i<n+1;i++){ temp=temp*(xx-x[i-1]);newton=newton+y[i][i]*temp;}printf("求得的结果为:N(%.4f)=%9f\n",xx,newton);}附录C：牛顿插值法的Matlab程序function f = Newton(x,y,x0)syms t;if(length(x) == length(y))n = length(x);c(1:n) = 0.0;elsedisp('x和y的维数不相等！');return;endf = y(1);y1 = 0;l = 1;for(i=1:n-1)for(j=i+1:n)y1(j) = (y(j)-y(i))/(x(j)-x(i));endc(i) = y1(i+1);l = l*(t-x(i));f = f + c(i)*l;simplify(f);y = y1;if(i==n-1)if(nargin == 3)f = subs(f,'t',x0);elsef = collect(f); %将插值多项式展开 f = vpa(f, 6);endendend附录D：等距节点下的前插公式的C语言程序#include<stdio.h>void main(){float x[11],y[11][11],xx,temp,newton,t,h;int i,j,n;printf("Newton插值:\n请输入要运算的值:x=");scanf("%f",&xx);printf("请输入插值的次数(n<11):n=");scanf("%d",&n);printf("步长:\n请输入要运算的值:h=");scanf("%f",&h);printf("请输入%d组值:\n",n+1);for(i=0;i<n+1;i++){printf("x%d=",i);scanf("%f",&x[i]);printf("y%d=",i);scanf("%f",&y[0][i]);}t=(xx-x[0])/h;for(i=1;i<n+1;i++)for(j=i;j<n+1;j++){y[i][j]=(y[i-1][j+1]-y[i-1][j]);printf("%f\n",y[i][i]);}temp=1;newton=y[0][0];for(i=1;i<n+1;i++){temp=temp*((t-i+1)/i);newton=newton+y[i][i]*temp;}printf("求得的结果为:N(%.4f)=%9f\n",xx,newton);}附录E：等距节点下的后插公式的C语言程序#include<stdio.h>void main(){float x[11],y[11][11],xx,temp,newton,t,h;int i,j,n;printf("Newton插值:\n请输入要运算的值:x=");scanf("%f",&xx);printf("请输入插值的次数(n<11):n=");scanf("%d",&n);printf("步长:\n请输入要运算的值:h=");scanf("%f",&h);printf("请输入%d组值:\n",n+1);for(i=0;i<n+1;i++){printf("x%d=",i);scanf("%f",&x[i]);printf("y%d=",i);scanf("%f",&y[0][i]);}t=(xx-x[n-1])/h;for(i=1;i<n+1;i++)for(j=i;j<n+1;j++){y[i][j]=(y[i-1][j]-y[i-1][j-1]);printf("%f\n",y[i][i]);}temp=1;newton=y[n-1][n-1];for(i=n;i<0;i--){temp=temp*((t-i+1)/i);newton=newton+y[n-i][n-i]*temp;}printf("求得的结果为:N(%.4f)=%9f\n",xx,newton);}。

01,,n1,,n1,,)n x及数值分析各算法流程图一、插值1、 拉格朗日插值流程图：( 相应程序：lagrintp(x,y,xx))2,,n ,,j n 1,2,,n 1,,)n 2、 牛顿插值流程图(1)产生差商表的算法流程图(相应程序：divdiff(x,y))注：1、另一程序divdiff1(x,y)，输出的矩阵包含了节点向量。

而divdiff(x,y)不含节点向量。

2、另一程序tableofdd(x,y,m)，输出的是表格形式，添加了表头。

1,,),,n m 及1,,m (2)非等距节点的牛顿插值流程图(相应程序：newtint11(x,y,xx,m)) 、注：1、虽然程序newtint11(x,y,xx,m)考虑了多种情形，看上去很复杂，但基本流程结构还是如上图所示。

2、程序中调用的子程序是divdiff 。

若调用的子程序是divdiff1的话，流程图中的第三，第四，第五步要相应的改一下数字。

2,3,,1m +1,,j1,2,,n=1,2,,)n m 及(3)求差分表的流程图(相应程序：difference(y,m))注：1、difference 输出的是矩阵D 。

而另一程序tableofd(y,m)，输出的是带有表头的差分表。

n x m1,,),,1,,m注：1、程序newtforward1(x,y,xx,m))的结构与上述流程图一致，xx可以是数组。

2、另一程序newtforward(x,y,xx,m))先求出插值多项式，再求插值多项式在插值点的函数值。

基本结构还是和上面的流程图一样。

n x m1,,),,-x x1,,m注：1、程序newtbackward1(x,y,xx,m))的结构与上述流程图一致，xx可以是数组。

2、另一程序newtbackward(x,y,xx,m))先求出插值多项式，再求插值多项式在插值点的函数值。

基本结构还是和上面的流程图一样。

1,2,,n1,2,,n ,2,,)n x及3、Hermite 插值流程图(1) 已知条件中一阶导数的个数与插值节点的个数相等时的Hermite 插值流程图。

1.#include "stdafx.h"2.#include "Fft.h"3.#include "math.h"4.5.const double MOD_MAX = 65535.0;6.const double DISP_MAX = 1.0/255.0;7.extern FILE *fp;8.9.10.11.12.13.//FFT运算必须参数14.int fft_point,fft_order,fft_divide,fft_window,fft_scale;15.bool fft_cover;16.float filter[7];//FIR滤波参数17.18.extern FILE *fpIandQ;19.extern bool m_bIqWrite;20.double prFilter[256],piFilter[256];21.22.23./**************************************************************************************24. 0 RectangleWindow矩形窗25. FFT变换结果为对称型，矩形窗是使信号突然截断，旁瓣会很大，且衰减较慢，旁瓣的第一个负26.峰值为主瓣的21％，第一个正峰值为主瓣的12.6％，第二个负峰值为主瓣的9％，效果一般，泄27.漏较大。

28.**************************************************************************************/29.double WINAPI RectangleWindow(int t)30.{31.double wt;32. wt=1.0;33.return wt;34.}35.36./**************************************************************************************37. 1 TriangleWindow三角窗,也称费杰（Fejer）窗，Bartlett38.**************************************************************************************/39.double WINAPI TriangleWindow(int t)40.{41.42.double wt;43. wt=1-t/fft_point;44.return wt;45.}46.47./**************************************************************************************48. 2 HanningWindwo汉宁窗，即升余弦窗49.**************************************************************************************/50.double WINAPI HanningWindow(int t)51.{52.double wt;53. wt=(1-cos(2*PI*t/fft_point))/2;54.return wt;55.}56.57./**************************************************************************************58. 3 HammingWindow海明窗，即改进的升余弦窗59.**************************************************************************************/60.double WINAPI HammingWindow(int t)61.{62.double wt;63. wt=0.54-0.46*cos(2*PI*t/fft_point);64.return wt;65.}66.67./**************************************************************************************68. 4 BlackmanWindow布来克曼窗，即三阶升余弦窗69.**************************************************************************************/70.double WINAPI BlackmanWindow(int t)71.{72.double wt;73. wt=0.42-0.5*cos(2*PI*t/fft_point)+0.08*cos(4*PI*t/fft_point);74.return wt;75.}76.***********78. 5 CosgradeWindow余弦坡度窗79.**************************************************************************************/80.double WINAPI CosgradeWindow(int t)81.{82.double wt;83.if(t= int(4*fft_point/5) )84. wt=1.0;85.else86. wt=(1+cos(5*PI*t/fft_point))/2;87.return wt;88.}89.90.91./**************************************************************************************92. 6 ParzenWindow帕曾窗93.**************************************************************************************/94.double WINAPI ParzenWindow(int t)95.{96.double wt;97.if(t= int(fft_point/2) )98. wt=1-6*pow(t/fft_point,2)+6*pow(t/fft_point,3);99.else100. wt=2*pow(1-t/fft_point,3);101.return wt;102.}103.104./*************************************************************************** ***********105. 7 ExponentWindow指数窗106.**************************************************************************** **********/107.d ouble WINAPI ExponentWindow(int t)108.{109.double wt,alf;110. alf=0.0001;111. wt=exp((-1.0)*alf*t);112.return wt;113.}114.***********116. 8 GaussWindow高斯窗117.**************************************************************************** **********/118.d ouble WINAPI GaussWindow(int t)119.{120.double wt,alf;121. alf=0.0001;122. wt=exp(-1*alf*t*t);123.return wt;124.}125.126./*************************************************************************** ***********127. 9 KaiserWindow凯泽窗,128.凯泽在1966（1974）发现，利用第一类零阶修正（变形）贝赛尔函数可以构成一种近似最佳的窗129.最主要参数：beat－可同时调整主瓣宽度和旁瓣，beat越大，窗越窄，频谱旁瓣越小，但主瓣相应增加130. beta=0相当与矩形窗131.**************************************************************************** **********/132.d ouble WINAPI KaiserWindow(int t)133.{134.double wt,alfa,beta;135. alfa=(fft_point-1)/2;136. beta=0;137.//Bessel零阶第一类贝塞尔函数138. wt=Jim_Bessel0_R(0,beta*sqrt(1-pow((t-alfa)/alfa,2)))/Jim_Bessel0_R(0,be ta);139.return wt;140.}141.142./*************************************************************************** ***********143. 10 ChebWindow契比雪夫窗144.**************************************************************************** **********/145.d ouble WINAPI ChebWindow(int t)146.{147.double wt;148. wt=1.0;149.return wt;150.}151.152./*************************************************************************** ***********153. 11 BartlettWindow巴特利特窗154.**************************************************************************** **********/155.d ouble WINAPI BartlettWindow(int t)156.{157.double wt;158. wt=1.0;159.return wt;160.}161.162.b ool WINAPI fftInit(int point,int order,int divide,bool cover,int scale) 163.{164.165. fft_point=point;166. fft_order=order;167. fft_divide=divide;168. fft_cover=cover;169. fft_scale=scale;170.//采用8阶的FIR滤波171. Jim_FirFilter(7,3,50,7800/2.0,7800.0,1,filter);172.return TRUE;173.}174.175.176.B OOL WINAPI fftEnd()177.{178.return TRUE;179.}180.181.182.183.184.//傅立叶变换185.B OOL WINAPI fftTransform(VOID *inBuf,VOID *outBuf,int windows,int scale,int nProbe,int nWork)186.{187.BYTE *srcPtr = (BYTE*)inBuf;//数据源1024字节188.BYTE *dstPtr = (BYTE*)outBuf;//数据果256字节189.//double *testPtr =(double *)testBuf;190.WORD Ivalue,Qvalue;191.//WORD IvalueRev,QvalueRev;192. unsigned char flagQ,flagI;193.//double alfa,e;194.int i,j,Iorg=0,Qorg=0;195.double mod = 0;196.if(nProbe==LEFT_PROBE)197. {198. flagQ=0x30,flagI=0x20;199. }200.else201. {202. flagQ=0x10,flagI=0x0;203. }204.205.//每次取4个字节数据，分离I/Q分量，判断db(12)，0表示I分量，1表示Q分量206.if (((srcPtr[1]&0xf0) ==flagQ) && ((srcPtr[3]&0xf0) ==flagI)) 207. {208. Iorg = 2;209. Qorg = 0;210. }211.if (((srcPtr[1]&0xf0) ==flagI) && ((srcPtr[3]&0xf0) ==flagQ)) 212. {213. Iorg = 0;214. Qorg = 2;215. }216.double pr[256],pi[256],w;217.218.219.//计算自相关函数和互相关函数220.//double Rii=0.0,Rqq=0.0,Riq=0.0;221.222.for ( i=0;i<fft_point;i++ )223. {224.225.//更改存储的偏移地址分离I/Q，226. Ivalue = *((short*)(srcPtr + i*4 + Iorg));227. Qvalue = *((short*)(srcPtr + i*4 + Qorg));228.//~位非运算，低12位位有效数据229. pr[i]=double(Ivalue & 0xfff);230. pi[i]=double(Qvalue & 0xfff);231.232./*Rii=Rii+double(pr[i]*pr[i]);233. Rqq=Rqq+double(pi[i]*pi[i]);234. Riq=Riq+double(pr[i]*pi[i]);*/235. }236.237./*Rii=Rii/fft_point;238. Rqq=Rqq/fft_point;239. Riq=Riq/fft_point;240.241. alfa=asin(Riq/sqrt(Rii*Rqq));242. e=sqrt(Rqq/Rii)-1;*/243.244.for ( i=0;i<fft_point;i++ )245. {246.//防止谱泄漏，进行加窗处理247.switch(windows)248. {249.case WND_RECTANGLE: w=RectangleWindow(i); break; 250.case WND_TRIANGLE: w=TriangleWindow(i); break; 251.case WND_HANNING: w=HanningWindow(i); break; 252.case WND_HAMMING: w=HammingWindow(i); break; 253.case WND_BLACKMAN: w=BlackmanWindow(i); break; 254.case WND_COSGRADE: w=CosgradeWindow(i); break; 255.case WND_PARZEN: w=ParzenWindow(i); break; 256.case WND_EXPONENT: w=ExponentWindow(i); break; 257.case WND_GAUSS: w=GaussWindow(i); break; 258.case WND_KAISER: w=KaiserWindow(i); break; 259.case WND_CHEB: w=ChebWindow(i); break; 260.case WND_BARTLETT: w=BartlettWindow(i); break; 261.default: w=RectangleWindow(i); break; 262. }263.264./*pr[i]=((1+e)*cos(alfa)*pr[i])/((1+e)*cos(alfa))*w; 265. pi[i]=(pi[i]-(1+e)*sin(alfa)*pr[i])/((1+e)*cos(alfa))*w;*/ 266.267. pr[i]=pr[i]*w;268. pi[i]=pi[i]*w;269. prFilter[i]=piFilter[i]=0.0;270.for(j=0;j<7;j++)271. {272.if(i+j>255)273. {274. prFilter[i]=prFilter[i];275. piFilter[i]=piFilter[i];276. }277.else278. {279. prFilter[i]=prFilter[i]+pr[i+j]*w*filter[j];280. piFilter[i]=piFilter[i]+pi[i+j]*w*filter[j];281. }282. }283.//fprintf(fp,"%d,%f\n",i,filter[i]);284.//pr[i]=pr[i]*filter[i];285.//pi[i]=pi[i]*filter[i];286. }287.//fclose(fpIandQ);288.289.//fft处理290.//Jim_FFT(pr,pi,fft_point,fft_order,0,1);291. Jim_FFT(prFilter,piFilter,fft_point,fft_order,0,1);292./*for(int t=0;t<256;t++)293. {294. fprintf(fpIandQ,"%d\n",int(prFilter[t]));295. }*/296.//归一化297.//Jim_Unitary(pr,dstPtr,scale,nProbe);298. Jim_Unitary(prFilter,dstPtr,scale,nWork);299.300.301.return TRUE;302.}303.304./*************************************************************************** ***********305. I/Q相位校正算法 by Jim Fang at 2007306.**************************************************************************** **********/307.v oid WINAPI Jim_Pharev(void *pIvalue,void *pQvalue,void *pRev)308.{309.double *iPtr =(double *)pIvalue;310.double *qPtr =(double *)pQvalue;311.double *rPtr =(double *)pRev;312.//计算自相关函数和互相关函数313.double Rii=0.0,Rqq=0.0,Riq=0.0;314.int i;315.for ( i=0;i<fft_point;i++)316. {317. Rii=Rii+double(iPtr[i]*iPtr[i]);318. Rqq=Rqq+double(qPtr[i]*qPtr[i]);319. Riq=Riq+double(iPtr[i]*qPtr[i]);320. }321. rPtr[0]=sqrt(Rqq/Rii)-1;322. rPtr[1]=asin(Riq/sqrt(Rii*Rqq));323.}324.325./*************************************************************************** ***********326.能量归一化算法 by Jim Fang at 2007327. FFT的源数据为WORD型，范围在±32768之间，用于显示能量的色阶范围为0-255328. 1.消除奇异点329. 2.找出峰值330. 3.归一化331.**************************************************************************** **********/332.v oid WINAPI Jim_Unitary(void *pSrcData,BYTE *pDstData,int nScale,int nWork)333.{334.//平滑处理后取最大值作为255，进行归一化。