hermite插值的上机实现及应用课程设计正文-大学论文.doc

学校代码： 10128
学 号： 20122090

课程设计说明书

题 目：Hermite插值的上机实现及应用
学生姓名：
学 院：理学院
班 级：
指导教师：任文秀 曹艳

2015年 1月 16日

目录
摘要 1
第一章Hermite插值的上机实现 2
§1.1 插值概述 2
§1.1.1插值问题的提出 2
§1.1.2插值的种类 2
§1.2 Hermite插值的问题 5
§1.2.1 Hermite插值的几种形式 5
§1.2.2 Hermite插值的几个重要定理 11
§1.2.3 Hermite插值的优点 12
§1.3 Hermite插值的源程序 12
§1.3.1 三次Hermite插值的C程序 12
§1.3.2 二重Hermite插值的matlab程序 13
第二章 Hermite插值的应用 14
§2.1 Hermite插值函数的工程应用 14
§2.2 应用Hermite插值作心电图基线漂移校正 16
参考文献 24
附录A 三次Hermite插值的C程序 25
附录B 二重Hermite插值的MATLAB程序 28

1

摘要
随着计算机技术的普及和应用的日益广泛,细分方法在近年来已经成为了计算机辅助设计(CAD)和计算机图形学(CG)领域内的一个国际性研究热点。通过近三十年的发展,细分方法日趋完善,多数经典的细分方法已经建立起了较为系统的理论知识体系。1992年Merrien首次提出了Hermite型的插值细分格式,随后Hermite插值型细分方法得到了迅速的发展,从一维区间上生成C1、C2细分曲线的格式到维矩形网格上生成光滑曲面的格式得以在短时间内展现,但是对于二维矩形上生成的光滑曲面在直观上与采样函数有不小的差距.
在构造插值时，对所构造的插值，不仅要求差值多项式节点的函数值与被插函数的函数值相同，还要求在节点处的插值函数与被插函数的一阶导数的值也相等对所构造的插值，不仅要求差值多项式节点的函数值与被插函数的函数值相同，还要求在节点处的插值函数与被插函数的一阶导数的值也相等.
关键词 Hermite插值；拉格朗日插值；Newton插值；余项；Hermite插值应用

第一章Hermite插值的上机实现
§1.1 插值概述
§1.1.1插值问题的提出
在许多实际问题及科学研究中，因素之间往往存在着函数关系，但这些关系的表达式不一定都知道，通常只是由观察或测试得到一些离散数值，所以只能从这些数据构造函数的近似表达式.有时，虽然给出了解析表达式，不过由于解析表达式过于复杂，使用或计算起来十分麻烦.这就需要建立某种近似表达，因此引入插值.
§1.1.2插值的种类
类型1 拉格朗日插值.
定义1.1 若函数y=f(x)在若干点的函数值=（i=0,1,,n），则另一个函数(x)：p()=,i=0,1,,n,则称p(x)为f(x)的插值函数，而f(x)为被插值函数.对于，且，用（x)的值作为f(x)的近似值或估计值，常称内插法.对于，用(x)的值去估计f(x)的值，又称外插法.
注解1.1 拉格朗日插值分为线性插值和n次插值.
注解1.2 拉格朗日插值的余项为

类型2 Newton插值
定义1.2任何一个不高于次多项式，都可以表示成函数的线性组合.既可以把满足插值条件的次插值多项式写成如下形式：

其中，为待定系数，这种形式的插值多项式称为牛顿插值多项式，记为.
注解1.3 设互不相同，则关于的阶差商为：
.
则一阶差商为
.
且二阶差商为
.
总结以上可得如下表1-1.
表1-1 差商表

一阶差商
二阶差商
三阶差商
n阶差商

类型3 分段插值
定义1.3 对给定区间做划分

在每个小区间上作以为节点的线性插值，记这个插值
，
， ()
把每一个区间的线性插值函数连接起来，得到的以为剖分节点的分段性函数.
注解1.4 分段插值的基本思想
将被插值函数的插值节点由小到大排序，然后在每对相邻的两个节点为端点的区间上用次多项式去近似.
类型4 Hermite插值
定义1.4 Hermite插值是利用未知函数在插值节点上的函数值及导数值来构造插值多项式；分为带导数的插值与不带导数的插值二类.
类型5 三次样条插值
样条插值是一种改进的分段插值.
定义1.5 函数
，
且在每个小区间上是三次多项式，其中

是给定节点，则称是节点上的三次样条函数.
若在节点上给定函数值
，
并且
，
则称为三次样条插值函数.
注解1.5 本文着重介绍Hermite插值
§1.2 Hermite插值的问题
§1.2.1 Hermite插值的几种形式
类型一 Hermite插值的一般形式
求一个次数不大于n+r+1的代数多项式,满足
, i=0,1,2,...,n . （1.1）

称以上的插值问题为Hermite插值问题
注解1.6 Hermite插值多项式的推导（即建立Hermite插值多项式的方法）
令
（1.2）
其中和都是次待定多项式，并且它们满足以下条件：
（1.3）
（1.4）
显然满足条件式（1.3）,（1.4）的多项式（1.2）的次数不大于次，且满足插值条件式（1.1）.
形式一 求解（不带导数的Hermite插值）
由条件式（1.3）知

是的二重零点.
且由条件式（1.3）知

是的零点.
当时具有如下形式：

（1.5）
其中，是待定系数.
由条件式（1.3）知

即

由上述两式解得
.
.

将A,B代入式（1.5），得
（1.6）
其中，
.
.
.
.

当时，具有如下形式
. （1.7）
由条件式（1.3）知
.
将C代入式（1.7），得
（1.8）
其中，
.
.
.

综合式（1.1）、（1.2）可以得到，即式（1.6）、（1.8）
形式二 求解(即带导数的Hermite插值）
由条件式（1.4）知是的二重零点，且由条件式（1.4）知
是的零点.
当时，具有如下形式：
（1.9）
由条件式（1.4）知

将D代入式（1.9），得
. (1.10)
其中，
.
.
由式（1.2）,（1.6）,（1.8）,（1.10）所表示的多项式称为Hermite插值多项式，其中由式（1.6）,（1.8），（1.10）所表示的多项式称为Hermite插值基函数.
Hermite插值多项式的余项为
(x)=(x).
类型二 二重Hermite插值多项式
一般的Hermite插值为m=2 的情况，即给定的插值节点

均为二重节点，更具体些
，
及插值节点，若有

满足
.
.
就称为关于节点的二重Hermite插值多项式.
类型三 三次Hermite插值
设是区间[a,b]上的实函数, 是[a,b]上相异两点, 且 ,
在上的函数值和一阶导数值分别为和, 求三次多项式, 使其满足：
.
称为三次埃尔米特插值多项式.
注解1.7 误差估计
定理1.1 设f(x)在包含、的区间[a,b]内存在四阶导数，则当x∈[a,b]时有余项
（且与x有关）
设

则当

时,余项有如下估计式
.
类型四 两点三次Hermite插值
设在节点处的函数值为、，在节点、处的一阶导数值为、，两个节点最高可以用3次Hermite多项式，作为插值函数应满足插值条件
.
.
应用四个插值基函数表示，设的插值基函数为，

如果希望插值系数与Lagrange插值一样简单，那么重新假设

其中

可知是的二重零点，即可假设
.
由

可得
.

... ...
Lagrange插值基函数如下式所示

类似可得

将以上结果代入

得两个节点的三次Hermite插值公式

注解1.8 二点三次Hermite插值的余项
(x)=
§1.2.2 Hermite插值的几个重要定理
定理1.2 误差定理
若
，
则关于上节点的二重Hermite插值多项式误差为

定理1.3 唯一性定理
Hermite插值问题的表达式
.
.
的解H(x) 存在而且唯一.
定理1.4 Hermite插值余项定理
Hermite插值公式的余项为
.
其中,是插值区间内的某一点.
§1.2.3 Hermite插值的优点
分段线性插值的算法简单，计算量小，然而从整体上看，逼近函数不够光滑，在节点处，逼近函数的左右导数不相等.
Hermite插值的逼近函数与被逼近函数不仅在插值节点上取相同的函数值，而且逼近函数与被逼近函数在插值节点上去相同的若干阶导数值.
Hermite插值法结合了函数的导数值，使得插值的精度更为提高.
Hermite插值具有少节点得到高次插值多项式的特点.
Hermite插值插值多项式灵活多样.
Hermite插值在节点一定的条件下，可以多种构造插值条件.
§1.3 Hermite插值的源程序
§1.3.1 三次Hermite插值的C程序
例1.1 已知函数
y=1/(1+x2)
在区间[0,3]上取等距插值节点，求区间[0,3]上的分段三次埃尔米特插值函数，并利用它求出f(1.5)的近似值（0.3075）
表1-2 例题数据表

0
1
2

1
0.5
0.2

0
-0.5
-0.16

注解1.9 本例题程序流程图及C程序详见附录A.
§1.3.2 二重Hermite插值的matlab程序
注解1.10 程序及程序演示详见附录B

第二章 Hermite插值的应用
§2.1 Hermite插值函数的工程应用
对于同一个问题运用不同的方法或许都能得到相同的结果，但是每一个方法都有其得天独厚的优势以及劣势.特别是在现在这个现代化的信息时代，计算已经变得越来越重要，对计算结果的要求也十分苛刻.插值方法在实际问题中有着广泛的应用它能使一个有着大量数据的问题变得简单明了、易于观察，因此，地位自然不喻.Hermite插值为使插值函数能更好地和原来的函数重合，不但要求二者在节点上函数值相等，而且还要求相切，对应的导数值也相等，甚至要求高阶导数也相等,凭借其精度高，计算严谨被大多数人应用了起来.
算例分析
在土方工程中，土的最大干密度与最优含水量是确保路基压实质量的两个关键指标,利用埃尔米特插值函数求得的干密度、含水量能更好地逼近试验得到的 一 曲线，求解精度较高.
通过绘制干密度与含水量的相关曲线，即pd 一曲线，求得最大干密度与最优含水量的方法为图解法。图解法因简便直观而在实际工作中被广泛采用，但图解法随意性大，易产生人为误差。目前，数解法主要有两类：一是利用曲线拟合法求解，二是利用代数插值求解.用上述方法分别对试验的工程实例进行了求解，发现所得结果 的差值较大，其中最大干密度差值达0.01～0.06，最佳含水量差值达 0.5％ ～1.4％.在本研究中利用埃尔米特插值问题，试图更加精确地求解最大干密度与最优含水量.
例2.1 某公路工程路基填七的一组室内标准击实试验结果见表2-1，由表2-1可知，其最火干密度应在含水量11.6％附近.
表2-1 室内标准击实试验结果
试验序号
1
2
3
4
5
含水量%
5.8
7.4
11.6
15.5
17.6
干密度（g)
1.77
1.80
1.85
1.82
1.78
根据图解法,将最大干密度定为1.85 g/,对应的最优含水最为 l1.6％.而根据曲线图，最优含水量在 12％附近更为恰当.下面利用埃尔米特插值函数求解最大干密度与最优含水量.取分别为7.4、l1.6、 15.5， 对应的分 别 为 1.80、1.85、1.82 ，得 到
= 0．01l 905,
=-0.0024194.
步骤一 建立干密度、含水量的埃尔米特插值函数.利用式

建立干密度、含水量的埃尔米特插值函数为
H ( ) = 1.8+0.0l1905 (-7.4 )-0.002 419 4(-7.4) (-11.6 ) +A(-7.4 ) (-11.6 )(-15.5)，
利用式子
.
可得
A = -0.06105 + 0.00010644.
步骤二 求解最大干密度与最优含水量.取= 5.8％， 对应的
=1.77g/,

根据插值条件
,
代入式
A = -0.06105 + 0.00010644,
令
,
得
,
解此方程得最优含水量为12.3％.得最大干密度为 1．85 g／.
步骤三 误差分析：由表2-1中试验数据可得
，
由

和

可得

=
=4.751,
根据误差分析可知，此法求解最大干密度与最优含水量的精度较高，能更好地逼近试验中得到的曲线.
模糊矩阵综合评价得：

=
=
以上计算结果表明，I级水的隶属度为0.3787，II级水的隶属度0.1336，III级水的隶属度为0.136 8，I V 级水的隶属度为0.2996，V级水的隶属度为0.4313.由于V 级水的隶属度最大，因此鉴湖水体综合评价等级应为V级.
总结
应用模糊数学原理综合评判鉴湖水质等级，比采用单因子极值评价更为合理.评判结果表明，鉴湖所在地区由于经济社会的快速发展，已经造成了严重的水体污染，因此水质等级很快由Ⅲ类变成V类.水体的污染引起的一系列问题应该引起足够的重视，如果这样发展下去，鉴湖将失去它原来的价值，因而政府应该采取措施，防止和减轻水污染，努力提高鉴湖水质等级，从而使之能发挥更好的作用.
§2.2 应用Hermite插值作心电图基线漂移校正
消除心电图的基线漂移是个重要向题.采用分段三次Hermite插值来作基线漂移 校正,提出了当心率变化引起插值区间信号长度变化时,插值墓函数的线性变换规 则.由此可以保持拟合的高精度,又减少计算量.有可能用于实时心电监护.如果监护仪 中的CPU能力有限,本文还提出了一种计算 Hermite插值函数硬件电路,使每一点的计算时间缩短为12微秒 .
心电图(ECG)信号的计算机处理历来国内外十分重视.国内外其临床应用主要分为二大类：一是ECG计算机辅助诊断,主要用于医院的心电分析中心,常为离线分析,使用的计算机也多为中小型机甚至大型机；二是作ECG实时监护,主要用于临床危重 病人、手术病人的监护,强调实时性要求,计算机多是由ECG等集成片构成,计算能力与存贮容量均受到限制.尽管ECG计算机分析已有二十多年的历史,国内外已做了大量的工作,但是仍然存在不少困难问题未予妥善解决.例如：消除ECG基线漂移是实 时监护中的一个重要而又困难的问题.
导致ECG基线漂移的主要因素有：电极的极化电位的变化,心电放大器的直流偏 置漂移,人体由于呼吸或其它肌肉、体位的缓慢移动等.尽管可以努力消除产生基线漂 移的原因例如努力使病人静卧不动，改善电极材料与导电膏的性能，改善心电放大器的特性等，但基线漂移仍然是不可避免的，因而会造成诊断疾病的困难.
消除基线漂移的困难在于基飘的频率很低，其范围为0.05Hz至1Hz，主要分量在0.1Hz左右，如图2.1所示，而ST段的频率成分也很低，其最大分量在0.6Hz-0.7Hz左右，它们的频谱非常接近.所以若使用高通频率滤波的方法以消除基飘，即使采用线性相位滤波器，仍会引起ST段的严重失真，而ST段在临床上有重要的价值.

图2.1 基线漂移与ST段的频谱
目前解决基线漂移的方法，除了高通滤波外，常采用某种数学函数校正法，如分段直线校正，三次样条函数校正，二次函数校正及三次函数校正法.在每个心电周期中选取1-2个零电位点作为插值结点，俩结点之间的心电漂移，以消除基飘.若采用直线进行逼近，是为直线校正法，这种方法计算量小，可实时实现，对慢变化的基线漂移效果尚好，对变化较快的基飘误差就严重.应用三次样条函数插值，可以获得较高的精度，本次报告就三次样条函数插值进行谈论.
今设二个相邻结点为和，并已知这二个结点的函数值和一阶导数值为：
，则三次Hermite插值函数为：
.
满足下述条件：

.
上述式子中

(2.1)
是为插值基函数.由于实际心电信号的心率是不断随机变化的，所以不能按照等间隔计算，即值将随心率的变化而变化的.由于这个变化，将使得上述4个插值基函数随之改变.因而必须重新计算新的插值基函数，因此用一种简化插值基函数的计算方法，令，T为采样周期，k=0,1...m,t=kT,则可将式2.1插值基函数写成离散形式：
(2.2)
如若将

代入式（2.2），可得插值基函数为：

(2.3)

比较式（2.2）与（2.3），可得

由此可见，当变化时，插值基函数的幅值不变，只是时间轴发生线性变化.而的幅值也将发生线性变化.因而可得变换公式如下.

这样的变换可以使计算简化，图所示为各插值基函数随着，的变化而变化的图形.

图2.2 当变化时的插值基函数
确定插值结点，即选择心电信号的零电位点.一般常可选TP段，为此可先估计T波的终了点T1.根据临床心电图学Bazett（巴泽特）公式：

式中QT表示的是Q波起点Q1至T1的间期.RR是二R波的间隔.若确定了Q1则由QT值可得T1点.所以可先检测R波峰值 ，再往后定H点.该H点约在T1之后30至70ms处，便可作插值结点图2.3.这样可以吧连续二个心拍的H点作为二个插值节点，进行三次Hermite计算，然后作基线漂移校正 . 具体步骤如下.

图2.3 确定插值结点
Step1 确定信号长度m.如下确定了几个典型的m值，如表2-2所示.
表2-2 几个典型的m值

心率（次/分)

压缩方式d
修正系数
375
40

1.000
335
45
17
0.893
295
50
7
0.787
255
60
4
0.680
215
70
3
0.573
175
85
2
0.467
135
115
1
0.360
Step2 压缩方式是指由于m压缩后，基函数的点数也需作相应减少.
Step3 计算插值的边界条件.
实验是用8组不同心率的心电信号，迭加上不同频率的基线漂移（0.1Hz,0.2Hz10.3Hz)来进行基漂校正.图2.4所示为其中一例，心率为105次/分.迭加三种不同频率的正弦波作为基线漂移.图中 为原始信号,为由插值函数计算所得的基漂，为经校正后的心电信号.

图2.4 实时基漂校正结果.HR=105次/分.
原始ECG，由插值拟合的基漂，校正后的ECG.
而实际临床情况，心率一般均在60次/分以上，基漂频率为0.17Hz至0.33Hz之间，所以基漂校正的仿真结果误差都在1.0%以下，可以满足临床要求.ST段的计算也是令人满意的.
硬件电路实现
虽然上述插值方法经过改进与简化，计算量已有很大减少，但对小型实时心电监护仪来说，CPU还可能不能承担.因此又用专用硬件电路实现了上述的插值计算，并且还构成了一个ST段检测仪.其框图如图2.5所示.其中插值基函数电路是将Hermite插值基函数，其中（k=12m,m=375).计算并量化后写入EPROM片.再在乘法控制线的控制下可依次读出.插值条件寄存电路则由由CPU送入的每段插值结点的边界条件,它们也可在乘法控制线的控制下依次读出.这样每当由插值基函数电路端口读出函数值时，乘法控制线变回产生含有4个负脉冲的脉冲序列，乘法电路就依次产生4个对应的乘积,这四个乘积经累加电路累加后送至输出端口，完成一次基漂拟合值的计算.由此连续运行k=0至m，即可完成一个周期的基漂校正.这个电路具有高速性能，插值基函数的确定、乘法运算、累加、翻转、技数、清零、重复等操作均是在乘法控制线的控制下同步进行的，有一部分操作室并行进行的，最大限度的减少了运算时间，提高了运算速度，可在12微秒内完成一个点的插值计算，时钟脉冲频率为100MHz.且整个电路的成本也很低.此监测仪对于心率在40次/分以上，基线漂移频率在0.4Hz以下的ECG基线漂移能相当好地进行校正.对ST段的检测，在一般情况下也能满足临床要求.

图2.5 插值计算硬件电路框图

参考文献
[1] 文畅平.人民黄河.湖南：邵阳学院，2006.
[2] 李庆阳，王能超，易大义.数值分析[M].北京：清华大学出版社，2008.
[3] 白峰杉.数值计算引论.北京：高等教育出版社，2004.
[4] 李庆阳.计算科学方法基础.北京：清华大学出版社，2006.
[5] 冯康等.数值计算方法.北京：国防工业大学，1978.
[6] 张雪敏.MATLAB基础及应用.北京：中国电力出版社，2009.

附录A 三次Hermite插值的C程序
1. 流程图

N
Y
开始
输入
y=0, j=0
t=1
i=0,…j-1,j+1,…n

j=n?
输出y
结束
j=j+1

2.C程序代码
#include<stdio.h>
#include<math.h>
float f0(float x)
{
return((x-1)*(x-1)*(2*x+1));
}
float f1(float x)
{
return(x*x*(-2*x+3));
}
float g0(float x)
{
return(x*(x-1)*(x-1));
}
float g1(float x)
{
return(x*x*(x-1));
}
void main()
{float x0,x1,x,y0,y1,yy0,yy1,h,p;
printf("输入x0,x1,x,y0,y1和yy0,yy1的取值");
scanf("%f%f%f%f%f%f%f",&x0,&x1,&x,&y0,&y1,&yy0,&yy1);
h=x1-x0;
p=y0*f0((x-x0)/h)+y1*f1((x-x0)/h)+h*yy0*g0((x-x0)/h)+h*yy1*g1((x-

x0)/h);
printf("%f\n",p);
}
3. 运行结果截图

图1 三次Hermite插值的C程序结果截图

附录B 二重Hermite插值的MATLAB程序

1 MATLAB程序代码
function [f,f0] = Hermite1(x,y,y_1)
syms t；
f = 0.0；
if(length(x) == length(y))
if(length(y) == length(y_1))
n = length(x)；
else
disp('y和y的导数的维数不相等')；
return；
end
else
disp('x和y的维数不相等！ ')；
Return；
end
for i=1 n
h = 1.0;
a = 0.0;
for j=1 n
if( j ~= i)
h = h*(t-x(j))^2/((x(i)-x(j))^2);
a = a + 1/(x(i)-x(j));
end
end
f = f + h*((x(i)-t)*(2*a*y(i)-y_1(i))+y(i));
end
f0 = subs(f,'t');

程序输入值：
x=1 0.2 1.8;
y_1=[0.5 0.4564 0.4226 0.3953 0.3727];
y=[1 1.0954 1.1832 1.2649 1.3416];
>> [f,f0] = Hermite1(x,y,y_1)；
可得如下运行结果：
运行结果如下：
f =
(390625*((337232823972231*t)/35184372088832-114418258618928321/10995116277760000)*(t - 1)^2*(t - 7/5)^2*(t - 8/5)^2*(t - 9/5)^2)/36 - (390625*(t - 1)^2*(t-6/5)^2*(t-7/5)^2*(t - 9/5)^2*((2855713758717179*t)/281474976710656 - 384779664999124623/21990232555520000))/36 + (390625*((64*t)/3 - 61/3)*(t -6/5)^2*(t-7/5)^2*(t-8/5)^2*(t-9/5)^2)/576+(390625*((7612884810106783*t)/18014398509481984 + 6660373488918492043/11258999068426240000)*(t - 1)^2*(t -6/5)^2*(t-8/5)^2*(t-9/5)^2)/16-(390625*((7762319875242775*t)/281474976710656 - 8968626627620006931/175921860444160000)*(t - 1)^2*(t - 6/5)^2*(t - 7/5)^2*(t - 8/5)^2)/576
2. 程序演示图

图2 二次Hermite插值的MATLAB图形演示

28