精华内容
下载资源
问答
  • MATLAB插值与拟合§1曲线拟合实例:温度曲线问题气象部门观测到一天某些时刻的温度变化数据为:t012345678910T1315171416192624262729试描绘出温度变化曲线。曲线拟合就是计算出两组数据之间的一种函数关系,由此可...

    MATLAB插值与拟合

    §1曲线拟合

    实例:温度曲线问题

    气象部门观测到一天某些时刻的温度变化数据为:

    t

    0

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    T

    13

    15

    17

    14

    16

    19

    26

    24

    26

    27

    29

    试描绘出温度变化曲线。

    曲线拟合就是计算出两组数据之间的一种函数关系,由此可描绘其变化曲线及估计非采集数据对应的变量信息。

    曲线拟合有多种方式,下面是一元函数采用最小二乘法对给定数据进行多项式曲线拟合,最后给出拟合的多项式系数。

    1.线性拟合函数:regress()

    调用格式: b=regress(y,X)

    [b,bint,r,rint,stats]= regress(y,X)

    [b,bint,r,rint,stats]= regress(y,X,alpha)

    说明:b=regress(y,X)返回X与y的最小二乘拟合值,及线性模型的参数值β、ε。该函数求解线性模型:

    y=Xβ+ε

    β是p´1的参数向量;ε是服从标准正态分布的随机干扰的n´1的向量;y为n´1的向量;X为n´p矩阵。

    bint返回β的95%的置信区间。r中为形状残差,rint中返回每一个残差的95%置信区间。Stats向量包含R2统计量、回归的F值和p值。

    例1:设y的值为给定的x的线性函数加服从标准正态分布的随机干扰值得到。即y=10+x+ε

    ;求线性拟合方程系数。

    程序: x=[ones(10,1) (1:10)'];

    y=x*[10;1]+normrnd(0,0.1,10,1);

    [b,bint]=regress(y,x,0.05)

    结果: x =

    1 1

    1 2

    1 3

    1 4

    1 5

    1 6

    1 7

    1 8

    1 9

    1 10

    y =

    10.9567

    11.8334

    13.0125

    14.0288

    14.8854

    16.1191

    17.1189

    17.9962

    19.0327

    20.0175

    b =

    9.9213

    1.0143

    bint =

    9.7889 10.0537

    0.9930 1.0357

    即回归方程为:y=9.9213+1.0143x

    2.多项式曲线拟合函数:polyfit( )

    调用格式: p=polyfit(x,y,n)

    [p,s]= polyfit(x,y,n)

    说明:x,y为数据点,n为多项式阶数,返回p为幂次从高到低的多项式系数向量p。矩阵s用于生成预测值的误差估计。(见下一函数polyval)

    例2:由离散数据

    x

    0

    .1

    .2

    .3

    .4

    .5

    .6

    .7

    .8

    .9

    1

    y

    .3

    .5

    1

    1.4

    1.6

    1.9

    .6

    .4

    .8

    1.5

    2

    拟合出多项式。

    程序:

    x=0:.1:1;

    y=[.3

    .5 1 1.4 1.6 1.9 .6 .4 .8 1.5 2];

    n=3;

    p=polyfit(x,y,n)

    xi=linspace(0,1,100);

    z=polyval(p,xi);

    %多项式求值

    plot(x,y,'o',xi,z,'k:',x,y,'b')

    legend('原始数据','3阶曲线')

    结果:

    p =

    16.7832 -25.7459 10.9802 -0.0035

    多项式为:16.7832x3-25.7459x2+10.9802x-0.0035

    曲线拟合图形:

    a4c26d1e5885305701be709a3d33442f.png

    如果是n=6,则如下图:

    a4c26d1e5885305701be709a3d33442f.png

    也可由函数给出数据。

    例3:x=1:20,y=x+3*sin(x)

    程序:

    x=1:20;

    y=x+3*sin(x);

    p=polyfit(x,y,6)

    xi=linspace(1,20,100);

    z=polyval(p,xi);

    %多项式求值函数

    plot(x,y,'o',xi,z,'k:',x,y,'b')

    legend('原始数据','6阶曲线')

    结果:

    p =

    0.0000 -0.0021 0.0505 -0.5971 3.6472

    -9.7295 11.3304

    a4c26d1e5885305701be709a3d33442f.png

    再用10阶多项式拟合

    程序:x=1:20;

    y=x+3*sin(x);

    p=polyfit(x,y,10)

    xi=linspace(1,20,100);

    z=polyval(p,xi);

    plot(x,y,'o',xi,z,'k:',x,y,'b')

    legend('原始数据','10阶多项式')

    结果:p =

    Columns 1 through 7

    0.0000 -0.0000 0.0004 -0.0114 0.1814 -1.8065 11.2360

    Columns 8 through 11

    -42.0861 88.5907 -92.8155 40.2671

    a4c26d1e5885305701be709a3d33442f.png

    可用不同阶的多项式来拟合数据,但也不是阶数越高拟合的越好。

    3.

    多项式曲线求值函数:polyval( )

    调用格式: y=polyval(p,x)

    [y,DELTA]=polyval(p,x,s)

    说明:y=polyval(p,x)为返回对应自变量x在给定系数P的多项式的值。

    [y,DELTA]=polyval(p,x,s) 使用polyfit函数的选项输出s得出误差估计Y

    DELTA。它假设polyfit函数数据输入的误差是独立正态的,并且方差为常数。则Y DELTA将至少包含50%的预测值。

    (未完)

    展开全文
  • matlab 插值算法

    2013-01-23 10:05:22
    根据沈祖仪的《水轮机调节》》一书及其他相关文献建立的水轮机传递系数建模与插值模型
  • 在信号处理中,滤波器的系数我们往往都是通过MATLAB来设计,只要我们知道滤波器的通带截止频率和阻带起始频率,就可以通过MATLAB中的fdatool(在MATLAB2020中使用filterDesigner)来设计滤波器了。 我们使用归一化的...

    90f24d355ca1ff4ea876c06cdaa3e703.png

    在信号处理中,滤波器的系数我们往往都是通过MATLAB来设计,只要我们知道滤波器的通带截止频率和阻带起始频率,就可以通过MATLAB中的fdatool(在MATLAB2020中使用filterDesigner)来设计滤波器了。

    我们使用归一化的参数来设计,通带截止频率是025,阻带起始频率是0.3,通带内纹波是0.2,阻带衰减是60dB,参数设置如下:

    d054cd0283343f086501810eb8640706.png

    那么问题来了,对于插值滤波器,如何确定通带和阻带的频率呢?这就涉及到我们刚开始学习数字信号处理时的插值和抽取理论。当信号抽取时,在数字频率上,信号的频谱是展宽的,当信号插值时,在数字频率上,信号的频谱是压缩的。这里我们强调数字频率,不是模拟频率,因为100MHz的采样率去采中频10MHz、带宽1MHz的信号,那么抽取2倍后,这个信号的频率还是10MHz,带宽还是1MHz,那为什么说数字频率上频谱展宽了呢?因为数字频率的2pi对应采样率。

    我们以信号处理书上这个经典的例子为例,原始信号的带宽是2pi/3,采样率是2pi,经过3倍抽取后,采样率由fs变为fs/3;而抽取后信号的采样率依旧对于数字域的2pi,因此原先的fs就对应6pi,信号带宽也就变成了2pi。

    ea4b3c152390118aa5a3b7724cf1a6af.png

    而插值滤波器则刚好相反,对于3倍的插值滤波器,信号带宽在数字频率上,缩小了1/3。也就是原来0~pi的区间缩小到0~pi/3,因此信号的截止频率就是pi/3,我们在设计滤波器时,直接指定截止频率是pi/3即可,至于阻带起始频率,我们可以设计的比通带截止频率稍大一些即可,同时还要考虑滤波器阶数,如果过渡带太窄了,滤波器阶数会太高。像我们上面设计的那个滤波器,正好可以适用于4倍插值滤波器。

    这里我们再提供一种解决方案,这种方法也是我强烈推荐的,就是当我们对一种设计没有头绪时,可以参考mathworks给出的设计。从哪参考呢?当然是MATLAB程序。我们知道Matlab的一个强大之处在于给我们提供了很多API可以调用,为我们节省了不少时间,而且大多数的函数我们都是可以看到源码的。比如我们今天所说的插值滤波器,可以直接使用resample函数,比如要对向量sig插值4倍,就可以直接使用sig2 = resample(sig, 4, 1)。这次我们再打开resample这个函数,可以看到:

    4aa3c88a7052100d2f94e5e129a19083.png

    这里的N是10,也就是说,如果是p倍插值,Matlab给出的插值滤波器阶数是2x10xp,也就是4倍插值滤波器对应阶数是80阶。再用firls来设计滤波器,最后再给滤波器加个kaiser窗。其实我们也可以直接使用fir2函数来设计,就是把图中高亮的行换成:

    h = fir2(L - 1, [0 2*fc 2*fc 1], [1 1 0 0]);
    

    fir2函数默认就是加了hamming窗的。使用这种方法设计的滤波器频响如下:

    66aa2319b46960539c385a24ed0824ee.png

    微信公众号:Quant_Times

    展开全文
  • Matlab插值与拟合

    2009-05-22 12:43:31
    弹簧在力F的作用下伸长x,一定范围内服从胡克定律:F与x成正比,即F=kx,k为弹性系数。现在得到下面一组x,F数据,并在(x, F)坐标下作图。可以看出,当F大到一定数值(如x=9以后)后,就不服从这个定律了。试由...
  • MATLAB实验报告 题目 第二次实验报告 学生姓名 学院 专业班级 学号 年月 MATLAB第二次实验报告 插 与 合 插 即在离散数据的基 上 插 函数使得 条 曲 通 全部 定的离散数据点插 是离散函数逼近的重要方法利用它可通 ...
  • robustfit参考资料:§1曲线拟合实例:温度曲线问题气象部门观测到一天某些时刻的温度变化数据为: ...曲线拟合有多种方式,下面是一元函数采用最小二乘法对给定数据进行多项式曲线拟合,最后给出拟合的多项式系数。...

    robustfit参考资料:

    §1曲线拟合

    实例:温度曲线问题

    气象部门观测到一天某些时刻的温度变化数据为: t

    0

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    T

    13

    15

    17

    14

    16

    19

    26

    24

    26

    27

    29

    试描绘出温度变化曲线。

    曲线拟合就是计算出两组数据之间的一种函数关系,由此可描绘其变化曲线及估计非采集数据对应的变量信息。

    曲线拟合有多种方式,下面是一元函数采用最小二乘法对给定数据进行多项式曲线拟合,最后给出拟合的多项式系数。

    1. 线性拟合函数:regress()

    调用格式: b=regress(y,X)

    [b,bint,r,rint,stats]= regress(y,X)

    [b,bint,r,rint,stats]= regress(y,X,alpha)

    说明:b=regress(y,X)返回X处y的最小二乘拟合值。该函数求解线性模型:

    y=Xβ+ε

    β是p?1的参数向量;ε是服从标准正态分布的随机干扰的n?1的向量;y为n?1的向量;X为n?p矩阵。

    bint返回β的95%的置信区间。r中为形状残差,rint中返回每一个残差的95%置信区间。Stats向量包含R2统计量、回归的F值和p值。

    例1:设y的值为给定的x的线性函数加服从标准正态分布的随机干扰值得到。即y=10+x+ε ;求线性拟合方程系数。

    程序: x=[ones(10,1) (1:10)’]

    y=x*[10;1]+normrnd(0,0.1,10,1)

    [b,bint]=regress(y,x,0.05)

    结果: x =

    1 1

    1 2

    1 3

    1 4

    1 5

    1 6

    1 7

    1 8

    1 9

    1 10

    y =

    10.9567

    11.8334

    13.0125

    14.0288

    14.8854

    16.1191

    17.1189

    17.9962

    19.0327

    20.0175

    b =

    9.9213

    1.0143

    bint =

    9.7889 10.0537

    0.9930 1.0357

    即回归方程为:y=9.9213+1.0143x

    2. 多项式曲线拟合函数:polyfit( )

    调用格式: p=polyfit(x,y,n)

    [p,s]= polyfit(x,y,n)

    说明:x,y为数据点,n为多项式阶数,返回p为幂次从高到低的多项式系数向量p。矩阵s用于生成预测值的误差估计。(见下一函数polyval)

    例2:由离散数据 x

    0

    .1

    .2

    .3

    .4

    .5

    .6

    .7

    .8

    .9

    1

    y

    .3

    .5

    1

    1.4

    1.6

    1.9

    .6

    .4

    .8

    1.5

    2

    拟合出多项式。

    程序:

    x=0:.1:1;

    y=[.3 .5 1 1.4 1.6 1.9 .6 .4 .8 1.5 2]

    n=3;

    p=polyfit(x,y,n)

    xi=linspace(0,1,100);

    z=polyval(p,xi); %多项式求值

    plot(x,y,’o’,xi,z,’k:’,x,y,’b’)

    legend(‘原始数据’,’3阶曲线’)

    结果:

    p =

    16.7832 -25.7459 10.9802 -0.0035

    多项式为:16.7832x3-25.7459x2+10.9802x-0.0035

    曲线拟合图形:

    也可由函数给出数据。

    例3:x=1:20,y=x+3*sin(x)

    程序:

    x=1:20;

    y=x+3*sin(x);

    p=polyfit(x,y,6)

    xi=1inspace(1,20,100);

    z=poyval(p,xi); %多项式求值函数

    plot(x,y,’o’,xi,z,’k:’,x,y,’b’)

    legend(‘原始数据’,’6阶曲线’)

    结果:

    p =

    0.0000 -0.0021 0.0505 -0.5971 3.6472 -9.7295 11.3304

    再用10阶多项式拟合

    程序:x=1:20;

    y=x+3*sin(x);

    p=polyfit(x,y,10)

    xi=linspace(1,20,100);

    z=polyval(p,xi);

    plot(x,y,'o',xi,z,'k:',x,y,'b')

    legend('原始数据','10阶多项式')

    结果:p =

    Columns 1 through 7

    0.0000 -0.0000 0.0004 -0.0114 0.1814 -1.8065 11.2360

    Columns 8 through 11

    -42.0861 88.5907 -92.8155 40.2671

    可用不同阶的多项式来拟合数据,但也不是阶数越高拟合的越好。

    3. 多项式曲线求值函数:polyval( )

    调用格式: y=polyval(p,x)

    [y,DELTA]=polyval(p,x,s)

    说明:y=polyval(p,x)为返回对应自变量x在给定系数P的多项式的值。

    [y,DELTA]=polyval(p,x,s) 使用polyfit函数的选项输出s得出误差估计Y

    DELTA。它假设polyfit函数数据输入的误差是独立正态的,并且方差为常数。则Y DELTA将至少包含50%的预测值。

    4. 多项式曲线拟合的评价和置信区间函数:polyconf( )

    调用格式: [Y,DELTA]=polyconf(p,x,s)

    [Y,DELTA]=polyconf(p,x,s,alpha)

    说明:[Y,DELTA]=polyconf(p,x,s)使用polyfit函数的选项输出s给出Y的95%置信区间Y

    DELTA。它假设polyfit函数数据输入的误差是独立正态的,并且方差为常数。1-alpha为置信度。

    例4:给出上面例1的预测值及置信度为90%的置信区间。

    程序: x=0:.1:1;

    y=[.3 .5 1 1.4 1.6 1.9 .6 .4 .8 1.5 2]

    n=3;

    [p,s]=polyfit(x,y,n)

    alpha=0.05;

    [Y,DELTA]=polyconf(p,x,s,alpha)

    结果: p =

    16.7832 -25.7459 10.9802 -0.0035

    s =

    R: [4x4

    double]

    df: 7

    normr: 1.1406

    Y =

    Columns 1 through 7

    -0.0035 0.8538 1.2970 1.4266 1.3434 1.1480 0.9413

    Columns 8 through 11

    0.8238 0.8963 1.2594 2.0140

    DELTA =

    Columns 1 through 7

    1.3639 1.1563 1.1563 1.1589 1.1352 1.1202 1.1352

    Columns 8 through 11

    1.1589 1.1563 1.1563 1.3639

    5. 稳健回归函数:robust( )

    稳健回归是指此回归方法相对于其他回归方法而言,受异常值的影响较小。

    调用格式: b=robustfit(x,y)

    [b,stats]=robustfit(x,y)

    [b,stats]=robustfit(x,y,’wfun’,tune,’const’)

    说明:b返回系数估计向量;stats返回各种参数估计;’wfun’指定一个加权函数;tune为调协常数;’const’的值为’on’(默认值)时添加一个常数项;为’off

    ’时忽略常数项。

    例5:演示一个异常数据点如何影响最小二乘拟合值与稳健拟合。首先利用函数y=10-2x加上一些随机干扰的项生成数据集,然后改变一个y的值形成异常值。调用不同的拟合函数,通过图形观查影响程度。

    程序:x=(1:10)’;

    y=10-2*x+randn(10,1);

    y(10)=0;

    bls=regress(y,[ones(10,1) x]) %线性拟合

    brob=robustfit(x,y) %稳健拟合

    scatter(x,y)

    hold on

    plot(x,bls(1)+bls(2)*x,’:’)

    plot(x,brob(1)+brob(2)*x,’r‘)

    结果 : bls =

    8.4452

    -1.4784

    brob =

    10.2934

    -2.0006

    分析:稳健拟合(实线)对数据的拟合程度好些,忽略了异常值。最小二乘拟合(点线)则受到异常值的影响,向异常值偏移。

    6. 向自定义函数拟合

    对于给定的数据,根据经验拟合为带有待定常数的自定义函数。

    所用函数:nlinfit( )

    调用格式: [beta,r,J]=nlinfit(X,y,’fun’,betao)

    说明:beta返回函数’fun’中的待定常数;r表示残差;J表示雅可比矩阵。X,y为数据;‘fun’自定义函数;beta0待定常数初值。

    例6:在化工生产中获得的氯气的级分y随生产时间x下降,假定在x≥8时,y与x之间有如下形式的非线性模型:

    现收集了44组数据,利用该数据通过拟合确定非线性模型中的待定常数。

    x y x y x y

    8 0.49 16 0.43 28 0.41

    8 0.49 18 0.46 28 0.40

    10 0.48 18 0.45 30 0.40

    10 0.47 20 0.42 30 0.40

    10 0.48 20 0.42 30 0.38

    10 0.47 20 0.43 32 0.41

    12 0.46 20 0.41 32 0.40

    12 0.46 22 0.41 34 0.40

    12 0.45 22 0.40 36 0.41

    12 0.43 24 0.42 36 0.36

    14 0.45 24 0.40 38 0.40

    14 0.43 24 0.40 38 0.40

    14 0.43 26 0.41 40 0.36

    16 0.44 26 0.40 42 0.39

    16 0.43 26 0.41

    首先定义非线性函数的m文件:fff6.m

    function yy=model(beta0,x)

    a=beta0(1);

    b=beta0(2);

    yy=a+(0.49-a)*exp(-b*(x-8));

    程序:

    x=[8.00 8.00 10.00 10.00 10.00 10.00 12.00 12.00 12.00 14.00

    14.00 14.00...

    16.00 16.00 16.00 18.00 18.00 20.00 20.00 20.00 20.00 22.00 22.00

    24.00...

    24.00 24.00 26.00 26.00 26.00 28.00 28.00 30.00 30.00 30.00 32.00

    32.00...

    34.00 36.00 36.00 38.00 38.00 40.00 42.00]';

    y=[0.49 0.49 0.48 0.47 0.48

    0.47 0.46 0.46 0.45 0.43 0.45 0.43 0.43 0.44 0.43...

    0.43 0.46 0.42 0.42 0.43 0.41 0.41 0.40 0.42 0.40 0.40 0.41 0.40

    0.41 0.41...

    0.40 0.40 0.40 0.38 0.41 0.40 0.40 0.41 0.38 0.40 0.40 0.39

    0.39]';

    beta0=[0.30 0.02];

    betafit = nlinfit(x,y,'sta67_1m',beta0)

    结果:betafit =

    0.3896

    0.1011

    即:a=0.3896 ,b=0.1011 拟合函数为:

    =================================================================================

    http://zhan.renren.com/matlab2012?gid=3602888498024090959&checked=true

    regress函数和regstats函数利用普通最小二乘法估计模型中的参数,参数的估计值受异常值的影响比较大。robustfit函数采用加权最小二乘法估计模型中的参数,受异常值的影响就比较小。robustfit函数用来作稳健的多重线性或广义线性回归分析,下面介绍robustfit函数的用法。

    1.4.1.robustfit函数的用法

    robustfit函数有以下几种调用方式:

    b = robustfit(X,y)

    b = robustfit(X,y,wfun,tune)

    b = robustfit(X,y,wfun,tune,const)

    [b,stats] = robustfit(…)

    (1)b = robustfit(X,y)

    返回多重线性回归方程中系数向量β的估计值b,这里的b为一个1p×的向量。输入参数X为自变量观测值矩阵(或设计矩阵),它是的矩阵。与regress函数不同的是,默认情况下,robustfit函数自动在X第1列元素的左边加入一列1,不需要用户自己添加。输入参数y为因变量的观测值向量,是的列向量。robustfit函数把y或X中不确定数据NaN作为缺失数据而忽略它们。np×1n×

    (2)b = robustfit(X,y,wfun,tune)

    用参数wfun指定加权函数,用参数tune 指定调节常数。wfun为字符串,其可能的取值如表1-3所示。

    表1-3 robustfit函数支持的加权函数

    加权函数(wfun)

    函数表达式

    默认调节常数值

    'andrews' sin(||)rwIrrπ=⋅<

    1.339

    'bisquare'(默认值)

    22(1)(||1)wrIr=−⋅<

    4.685

    'cauchy' 21(1)wr=+

    2.385

    'fair' 1(1||)wr=+

    1.400

    'huber' 1max(1, ||)wr=

    1.345

    'logistic' tanh()wr=

    1.205

    'ols'

    普通最小二乘,无加权函数

    'talwar'

    (||1)wIr=<

    2.795

    'welsch'

    2rwe−=

    2.985

    若调用时没有指定调节常数tune,则用表1-3中列出的默认调节常数值进行计算。表1-3中加权函数中的r通过下式计算residr

    =tunes1-h××

    其中resid为上一步迭代的残差向量,tune为调节常数,h是由最小二乘拟合得到的中心化杠杆值向量,s为误差项的标准差的估计。s的计算公式为:s

    =

    MAD/0.6745,其中MAD为残差绝对值的中位数,在正态分布下,这个估计是无偏的。若X中有p列,计算MAD时,将残差绝对值向量的前p个最小值舍去。

    用户可以定义自己的权重函数,函数的输入必须是残差向量,输出是权重向量。在调用robustfit函数时,把自定义权重函数的句柄(形如@myfun)作为wfun参数传递给robustfit函数,此时必须指定tune参数。

    (3)b = robustfit(X,y,wfun,tune,const)

    用参数const来控制模型中是否包含常数项。若const取值为 'on'

    或1,则模型中包含常数项,此时自动在X第1列的左边加入一列1,若const取值为 'off'

    或0,则模型中不包含常数项,此时不改变X的值。

    (4)[b,stats] = robustfit(…)

    返回一个结构体变量stats,它的字段包含了用于模型诊断的统计量。stats有以下字段:

    • stats.ols_s — 普通最小二乘法得出的σ的估计(RMSE);

    • stats.robust_s — σ的稳健估计;

    • stats.mad_s — 用残差绝对值的中位数计算σ的估计;

    • stats.s — σ的最终估计,是ols_s 和robust_s的加权平均与robust_s中的最大值;

    • stats.se — 系数估计的标准误差;

    • stats.t — b与stats.se的比值;

    • stats.p — t检验的p值;

    • stats.covb — 系数向量的协方差矩阵的估计;

    • stats.coeffcorr — 系数向量的相关系数矩阵的估计;

    • stats.w — 稳健拟合的权重向量;

    • stats.h — 最小二乘拟合的中心化杠杆值向量;

    • stats.R — 矩阵X的QR分解中的R因子

    =================================================================================

    §2 插值问题

    在应用领域中,由有限个已知数据点,构造一个解析表达式,由此计算数据点之间的函数值,称之为插值。

    实例:海底探测问题

    某公司用声纳对海底进行测试,在5×5海里的坐标点上测得海底深度的值,希望通过这些有限的数据了解更多处的海底情况。并绘出较细致的海底曲面图。

    一、一元插值

    一元插值是对一元数据点(xi,yi)进行插值。

    1. 线性插值:由已知数据点连成一条折线,认为相临两个数据点之间的函数值就在这两点之间的连线上。一般来说,数据点数越多,线性插值就越精确。

    调用格式:yi=interp1(x,y,xi,’linear’) %线性插值

    zi=interp1(x,y,xi,’spline’) %三次样条插值

    wi=interp1(x,y,xi,’cubic’) %三次多项式插值

    说明:yi、zi、wi为对应xi的不同类型的插值。x、y为已知数据点。

    例1:已知数据: x

    0

    .1

    .2

    .3

    .4

    .5

    .6

    .7

    .8

    .9

    1

    y

    .3

    .5

    1

    1.4

    1.6

    1.9

    .6

    .4

    .8

    1.5

    2

    求当xi=0.25时的yi的值。

    程序:

    x=0:.1:1;

    y=[.3 .5 1 1.4 1.6 1 .6 .4 .8 1.5 2];

    yi0=interp1(x,y,0.025,'linear')

    xi=0:.02:1;

    yi=interp1(x,y,xi,'linear');

    zi=interp1(x,y,xi,'spline');

    wi=interp1(x,y,xi,'cubic');

    plot(x,y,'o',xi,yi,'r+',xi,zi,'g*',xi,wi,'k.-')

    legend('原始点','线性点','三次样条','三次多项式')

    结果:yi0 = 0.3500

    要得到给定的几个点的对应函数值,可用:

    xi =[ 0.2500 0.3500 0.4500]

    yi=interp1(x,y,xi,'spline')

    结果:

    yi =1.2088 1.5802 1.3454

    (二) 二元插值

    二元插值与一元插值的基本思想一致,对原始数据点(x,y,z)构造见世面函数求出插值点数据(xi,yi,zi)。

    一、单调节点插值函数,即x,y向量是单调的。

    调用格式1:zi=interp2(x,y,z,xi,yi,’linear’)

    ‘liner’ 是双线性插值 (缺省)

    调用格式2:zi=interp2(x,y,z,xi,yi,’nearest’)

    ’nearest’ 是最近邻域插值

    调用格式3:zi=interp2(x,y,z,xi,yi,’spline’)

    ‘spline’是三次样条插值

    说明:这里x和y是两个独立的向量,它们必须是单调的。z是矩阵,是由x和y确定的点上的值。z和x,y之间的关系是z(i,:)=f(x,y(i))

    z(:,j)=f(x(j),y)

    即:当x变化时,z的第i行与y的第i个元素相关,当y变化时z的第j列与x的第j个元素相关。如果没有对x,y赋值,则默认x=1:n,

    y=1:m。n和m分别是矩阵z的行数和列数。

    例2:已知某处山区地形选点测量坐标数据为:

    x=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

    y=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6

    海拔高度数据为:

    z=89 90 87 85 92 91 96 93 90 87 82

    92 96 98 99 95 91 89 86 84

    82 84

    96 98 95 92 90 88 85 84 83

    81 85

    80 81 82 89 95 96 93 92 89

    86 86

    82 85 87 98 99 96 97 88 85

    82 83

    82 85 89 94 95 93 92 91 86

    84 88

    88 92 93 94 95 89 87 86 83

    81 92

    92 96 97 98 96 93 95 84 82

    81 84

    85 85 81 82 80 80 81 85 90

    93 95

    84 86 81 98 99 98 97 96 95

    84 87

    80 81 85 82 83 84 87 90 95

    86 88

    80 82 81 84 85 86 83 82 81

    80 82

    87 88 89 98 99 97 96 98 94

    92 87

    其地貌图为:

    对数据插值加密形成地貌图。

    程序:

    x=0:.5:5;

    y=0:.5:6;

    z=[89 90 87 85 92 91 96 93 90 87 82

    92 96 98 99 95 91 89 86 84

    82 84

    96 98 95 92 90 88 85 84 83

    81 85

    80 81 82 89 95 96 93 92 89

    86 86

    82 85 87 98 99 96 97 88 85

    82 83

    82 85 89 94 95 93 92 91 86

    84 88

    88 92 93 94 95 89 87 86 83

    81 92

    92 96 97 98 96 93 95 84 82

    81 84

    85 85 81 82 80 80 81 85 90

    93 95

    84 86 81 98 99 98 97 96 95

    84 87

    80 81 85 82 83 84 87 90 95

    86 88

    80 82 81 84 85 86 83 82 81

    80 82

    87 88 89 98 99 97 96 98 94

    92 87];

    mesh(x,y,z) %绘原始数据图

    xi=linspace(0,5,50); %加密横坐标数据到50个

    yi=linspace(0,6,80); %加密纵坐标数据到60个

    [xii,yii]=meshgrid(xi,yi); %生成网格数据

    zii=interp2(x,y,z,xii,yii,'cubic'); %插值

    mesh(xii,yii,zii) %加密后的地貌图

    hold

    on % 保持图形

    [xx,yy]=meshgrid(x,y); %生成网格数据

    plot3(xx,yy,z+0.1,’ob’) %原始数据用‘O’绘出

    2、二元非等距插值

    调用格式:zi=griddata(x,y,z,xi,yi,’指定插值方法’)

    插值方法有:

    linear % 线性插值 (默认)

    bilinear % 双线性插值

    cubic % 三次插值

    bicubic % 双三次插值

    nearest % 最近邻域插值

    例:用随机数据生成地貌图再进行插值

    程序:

    x=rand(100,1)*4-2;

    y=rand(100,1)*4-2;

    z=x.*exp(-x.^2-y.^2);

    ti=-2:.25:2;

    [xi,yi]=meshgrid(ti,ti); % 加密数据

    zi=griddata(x,y,z,xi,yi);% 线性插值

    mesh(xi,yi,zi)

    hold on

    plot3(x,y,z,'o')

    该例中使用的数据是随机形成的,故函数griddata可以处理无规则的数据。

    =================================================================================

    展开全文
  • polynomialSize总是等长(outputConv),它总是等于n(因为你有n个数据点,第n-1个多项式有n个系数),所以最后一个for循环和下一个行也可以用简单的L(k, ???? =乘数* outputConv;所以我在http://en.wikipedia.org/wik...

    是的,一些建议(在下面的版本1中实现):if循环可以与上面的组合(只需通过下面的jr(jr~ = j)使索引跳过k); polynomialSize总是等长(outputConv),它总是等于n(因为你有n个数据点,第n-1个多项式有n个系数),所以最后一个for循环和下一个行也可以用简单的L(k, 🙂 =乘数* outputConv;

    所以我在http://en.wikipedia.org/wiki/Lagrange_polynomial上复制了这个例子(并采用了他们的j-m表示法,但对我来说j是1:n而m是1:n和m~ = j),因此我的初始化看起来像

    clear; clc;

    X=[-9 -4 -1 7]; %example taken from http://en.wikipedia.org/wiki/Lagrange_polynomial

    Y=[ 5 2 -2 9];

    n=length(X); %Lagrange basis polinomials are (n-1)th order, have n coefficients

    lj = zeros(1,n); %storage for numerator of Lagrange basis polyns - each w/ n coeff

    Lj = zeros(n); %matrix of Lagrange basis polyns coeffs (lj(x))

    L = zeros(1,n); %the Lagrange polynomial coefficients (L(x))

    然后v 1.0看起来像

    jr=1:n; %j-range: 1<=j<=n

    for j=jr %my j is your k

    multiplier = 1;

    outputConv = 1; %numerator of lj(x)

    mr=jr(jr~=j); %m-range: 1<=m<=n, m~=j

    for m = mr %my m is your index

    outputConv = conv(outputConv,[1 -X(m)]);

    multiplier = multiplier * ((X(j) - X(m))^-1);

    end

    Lj(j,:) = multiplier * outputConv; %jth Lagrange basis polinomial lj(x)

    end

    L = Y*Lj; %coefficients of Lagrange polinomial L(x)

    如果你意识到l_j(x)的分子只是一个具有特定根的多项式,那么可以进一步简化 – 因为在matlab中有一个很好的命令 – poly.类似地,分母只是在X(j)处评估的多边形 – 因为存在多边形.因此,v 1.9:

    jr=1:n; %j-range: 1<=j<=n

    for j=jr

    mr=jr(jr~=j); %m-range: 1<=m<=n, m~=j

    lj=poly(X(mr)); %numerator of lj(x)

    mult=1/polyval(lj,X(j)); %denominator of lj(x)

    Lj(j,:) = mult * lj; %jth Lagrange basis polinomial lj(x)

    end

    L = Y*Lj; %coefficients of Lagrange polinomial L(x)

    为什么版本1.9而不是2.0?好吧,有可能有办法摆脱这最后的循环,并将其全部写入1行,但我现在想不到它 – 它是v 2.0的todo 🙂

    而且,对于甜点,如果你想获得与维基百科相同的图片:

    figure(1);clf

    x=-10:.1:10;

    hold on

    plot(x,polyval(Y(1)*Lj(1,:),x),'r','linewidth',2)

    plot(x,polyval(Y(2)*Lj(2,:),x),'b','linewidth',2)

    plot(x,polyval(Y(3)*Lj(3,:),x),'g','linewidth',2)

    plot(x,polyval(Y(4)*Lj(4,:),x),'y','linewidth',2)

    plot(x,polyval(L,x),'k','linewidth',2)

    plot(X,Y,'ro','linewidth',2,'markersize',10)

    hold off

    xlim([-10 10])

    ylim([-10 10])

    set(gca,'XTick',-10:10)

    set(gca,'YTick',-10:10)

    grid on

    产生

    享受并随意重用/改进

    展开全文
  • 克里金插值matlab程序

    热门讨论 2012-09-10 22:08:17
    克里金插值matlab程序 克里金(Kriging)插值法又称空间自协方差最佳插值法,它是以南非矿业工程师D.G.Krige的名字命名的一种最优内插法。克里金法广泛地应用于地下水模拟、土壤制图等领域,是一种很有用的地质统计...
  • Matlab三次样条插值

    2019-08-16 09:22:41
    先看这个,内容很详细: 添加链接描述 插值实操的一些补充 如何看插值多项式的系数 以csape为例: p_x = csape(t_t,x); p_x_xishu = p_x.coefs; 插值多项式长什么样子 ...
  • 拉格朗日插值算法Matlab实现

    千次阅读 2017-12-04 16:39:43
    拉格朗日插值算法是什么网上有更多详细的解释,这里附上自己所写的Matlab源代码:代码1(计算拉格朗日多项式系数):function [ result ] = lagr( X,Y ) %lagr 返回所求拉格朗日多项式的系数 % X为每个点的横坐标,Y...
  • 称 Nn (x) 为牛顿Newton均差插值多项式系数 ak 就是书本表 5-1 中第一条斜线上对应的数值 式1-2为插值余项由插值多项式唯一性可知它与书本式5.1.19是等价的事实上利用均差与导数关系式可由 式1-2推出书本式5.1.19但...
  • Matlab学习手记】拉格朗日插值

    千次阅读 2018-08-11 10:00:36
    源码 function [y0, coef] = Interpolation_Lagrange(x, y, x0) %{ 函数功能:拉格朗日插值法; 输入: x:已知点横坐标;... coef:插值系数; 示例: clear; clc; x = 0 : 0.2 : 2; y = sin(x); x0 = ...
  • % 放大、缩小 function J = myResize(I,Kx,Ky) %放大尺寸 % Kx = 1.8,Ky=1.3; % 获取原图像的尺寸 % 这里只做rgb图像缩放 if Kx<... error('请输入正确的缩放系数'); end if size(I,3)~=3 er...
  • MATLAB三次样条插值之三弯矩法首先说这个程序并不完善,为了实现通用(1,2,…,n)格式解题,以及为调用追赶法程序,没有针对节点数在三个以下的情况进行分类讨论。希望能有朋友给出更好的方法。首先,通过函数...
  • MATLAB回归、插值、逼近、拟合总结

    万次阅读 多人点赞 2017-06-29 11:11:58
    在求回归前,已经假设所有型值点同时满足某一曲线方程,计算只要求出该方程的系数 2、多项式插值:用一个多项式来近似代替数据列表函数,并要求多项式通过列表函数中给定的数据点。(插值曲线要经过型值点。) 3、...
  • 个人收集整理 仅供参考学习 % 高斯赛德尔迭代解线 性方程组 a=input' 请输入系数矩阵 a) %输入系数矩阵 a disp(a) [m,n]=size(a) if m~=n %若 a 不是方阵 则显示错误 disp'false) end de=det(a; if de==0 % 判断矩阵...
  • %高斯赛德尔迭代解线 性方程组 a=i nput'请输入系数矩阵 a) %输入系数矩阵a disp(a) [m,n]=size(a) if m~=n %若a不是方阵 则显示错误 disp'false) end de=det(a; if de==0 %判断矩阵a是 否为奇异矩阵 disp'矩阵奇异)...
  • 敲了一下午终于把代码敲完了,留念。 文章目录题目:三次样条插值程序1:Gauss列主元消去法程序2:三次样条插值程序3:主函数及输入数据 题目:三次样条插值 程序1:Gauss列主元消去法 ... % 比例系数矩阵 X =
  • % 图像缩放变换 双线性插值 % 输入: % img 灰白图像(彩色要多一个color维度或转化为灰度图) % s_x x方向上的比例系数 % s_y y方向上的比例系数 % 输出: % 缩放后的新图像 %% 双线性插值注意点: % 双线性插值...
  • matlab编写的三次样条插值源代码,最终画出拟合图像,并输出拟合函数的系数
  • %图像缩放变换 最近邻插值 % 输入: % img 灰白图像(若是彩色需多一个color维度或转化为灰度图) % s_x x方向上的比例系数 % s_y y方向上的比例系数 % 输出: % 缩放后的新图像 %% 思考? % 1. 为什么放大,新图能...
  • 在求回归前,已经假设所有型值点同时满足某一曲线方程,计算只要求出该方程的系数 2、多项式插值:用一个多项式来近似代替数据列表函数,并要求多项式通过列表函数中给定的数据点。(插值曲线要经过型值点。) 3、...
  • matlab开发-瑞利杭古拉伏尔散射系数插值大气瑞利角体积散射系数,A.Bucholtz,1995
  • 关于多项式Matlab命令一个多项式的幂级数形式可表示为:也可表示为嵌套形式:或因子形式幂系数:在matlab里,多项式用行向量表示,其元素未多项式的系数,并从左至右按降幂排列。Roots:求多项式的零点。%%例如:y=2...
  • MATLAB对一组数据进行插值的方法

    万次阅读 2017-09-25 20:07:09
    用多项式函数(10.2)作为插值函数时,希望通过解方程组(10.3)而得到待定系数 function y=lagrange(x0,y0,x); n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j...
  • Description: 计算4-3-4分段多项式的系数 Input: 向量t(递增序列), 向量p, 起始速度vs, 结束速度ve, 起始加速度as, 结束加速度ae, 向量维数n Output: 4-3-4分段多项式的系数a Author: Marc Pony(marc_pony@163.com) ...
  • 将数据点按多项式的形式进行拟合,使用最小二乘法,可以确定多项式的系数。多项式拟合有指令语句和图形窗口两种方法:1、多项式拟合指令polyfit(x,y,n) :多项式拟合,返回降幂排列的多项式系数。polyval(p,xi) :...
  • 一、算法原理 设曲线F(x),寻找到其极值区间[x1,x2],使其满足f(x1)>...f(x)=ax^2+bx+c,a b c 为系数。 ax1^2+bx1+c=f(x1) ax2^2+bx2+c=f(x2) ax3^2+bx3+c=f(x3) 写成矩阵形式...

空空如也

空空如也

1 2 3 4 5 6
收藏数 103
精华内容 41
关键字:

matlab插值系数

matlab 订阅