C++实现matlab中的interp1和interp2插值「建议收藏」

C++实现matlab中的interp1和interp2插值「建议收藏」头文件interpfun.h#ifndefINTERPFUN_H#defineINTERPFUN_H#include"math.h"#include<stdio.h>#include<string.h> #ifndefSafeDeleteVec#defineSafeDeleteVec(X){if((X))delete(X);(X)=NULL;}#endi…

大家好,又见面了,我是你们的朋友全栈君。

头文件

interpfun.h

#ifndef INTERPFUN_H
#define INTERPFUN_H
#include”math.h”
#include<stdio.h>
#include <string.h> 
#ifndef SafeDeleteVec
#define SafeDeleteVec(X) { if((X)) delete (X); (X)=NULL;}
#endif
class Interp
{

public:

char *strMethod1;
char *strMethod2;
char *strMethod3;


public:
Interp();
~Interp();
/*****以下两个函数为插值函数**/
void interp_onePoint(float *x,float *y,int n,float t,float *fval,char *method);//返回单个点的值
void interp_multiPoint(float *x,float *y,int n,float *t,float *fval,int m,char *method);//返回多个点的值,点的个数为m
float interp2_onePoint(float *x,float *y,float *z,int m,int n,float a,float b,char *method);//matlab interp2 实现 2D面数据插值
void interp2_multiPoint(float *x,float *y,float *z,int m,int n,float *a,float *b,int asize,int bsize,float *fval,char *method);//matlab interp2 实现 2D面向量插值


private:
//拉格朗日线性插值
float lgr(float *x,float *y,int n,float t);
//拉格朗日一元三点线性插值
float lg3(float *x,float *y,int n,float t);
//平滑插值 三次多项式插值
void spl(float *x,float *y,int n,int k,float t,float *s);

int findIndex(float *x,int t,int n);

};

实现:

interpfun.cpp

#include”interpfun.h”
Interp::Interp()
{

strMethod1 =”lgr”;
strMethod2 = “lg3”;
strMethod3 = “spl”;
}

Interp::~Interp()
{


}

float Interp::lgr(float *x,float *y,int n,float t)
{

int i,j,k,m;
float z,s;
z= 0.0f;
if(n<1)
return(z);
if(n==1)
{

z=y[0];
return(z);
}

if(n==2)
{

z=(y[0]*(t-x[1]) – y[1]*(t-x[0]))/(x[0]-x[1]);
return(z);

}
if(n>1 && t>x[n-1]){

z= y[n-1];
return(z);
}
if(n>1 && t<x[0]){

z= y[0];
return(z);
}

i =0;
while((x[i] <t) && (i<n)) {

i=i+1;
}

k = i-4;
if(k<0) 
k=0;
m = i+3;
if(m>n-1) 
m = n-1;
for( i = k;i<=m;i++)
{

s=1.0;
for( j =k;j<=m;j++){

if(j!=i) 
s= s*(t-x[j])/(x[i]-x[j]);
}
z = z+s*y[i];
}
return(z);

}

//拉格朗日一元三点插值
float Interp::lg3(float *x,float *y,int n,float t)
{

int i,j,k,m;
float z,s;
z= 0.0;
if(n<1) return (z);
if(n ==1) {z= y[0];return(z);}
if(n==2)
{

z = (y[0]*(t-x[1]) – y[1]*(t – x[0]))/(x[0] – x[1]);
return(z);
}
if(n>0&&t<x[0]){

z = y[0];
return(z);
}
if(n>0 && t>x[n-1]){

z = y[n-1];
return(z);
}

if(t<=x[1]) 
{

k =0;
m=2;
}
else if(t>=x[n-2])
{

k = n-3;
m = n-1;
}
else 
{

k =1;
m =n;
while ((m-k) !=1)
{

i = (k+m)/2;
if(t<x[i-1])
m =i;
else
k =i;
}
k = k-1;
m = m-1;
if(fabs(t-x[k]) <fabs(t- x[m]))
k = k-1;
else
m = m+1;
}
z = 0.0;
for(i =k;i<=m;i++)
{

s= 1.0;
for(j = k;j<=m;j++)
{

if(j!=i)
s = s*(t-x[j])/(x[i] – x[j]);
}
z = z+ s*y[i];

}
return (z);

}

//平滑插值 三次多项式插值
void Interp::spl(float *x,float *y,int n,int k,float t,float *s)
{

int kk, m,lc;

float u[5], p,q;

s[4] = 0.0f;
s[0] = 0.0f;
s[1] = 0.0f;
s[2] = 0.0f;
s[3] = 0.0f;

if(n<1)
return;
if(n == 1)
{

s[0] = y[0];
s[4] = y[0];
return;
}
if(n == 2)
{

s[0] = y[0];
s[1] = (y[1] – y[0])/(x[1] – x[0]);
if(k<0)
s[4] = (y[0]*(t – x[1]) – y[1]*(t – x[0]))/(x[0] – x[1]);
return;
}

if(k<0 && n>0 && t<x[0] )
{

s[4] = y[0];
return;
}

if(k<0 && n>0 && t>x[n-1] )
{

s[4] = y[n-1];
return;
}

if(k<0)
{

if(t<=x[1]) 
kk =0;
else if(t>=x[n-1])
kk = n-2;
else
{

kk=1;
m = n;
while (((kk -m) !=1) && ((kk – m) != -1))
{

lc=(kk +m)/2;
if(t<x[lc-1]) 
m = lc;
else
kk =lc;

}
kk = kk -1;
}
}
else
kk = k;

if( kk>n-1)
kk = n-2;

u[2] = (y[kk +1] – y[kk])/(x[kk +1] – x[kk]);
if (n == 3)
{

if(kk == 0)
{

u[3] = (y[2] – y[1])/(x[2]-x[1]);
u[4] = 2.0*u[3] – u[2];
u[1] = 2.0*u[2] – u[3];
u[0] = 2.0*u[1] – u[2];
}
else
{

u[1] = (y[1] – y[0])/(x[1] – x[0]);
u[0] = 2.0*u[1] – u[2];
u[3] = 2.0*u[2] – u[1];
u[4] = 2.0*u[3] – u[2];
}

}
else
{

if(kk <=1)
{

u[3] = (y[kk+2] – y[kk +1])/(x[kk+2]-x[kk+1]);
if(kk ==1)
{

u[1] = (y[1] -y[0])/(x[1] – x[0]);
u[0] = 2.0*u[1] – u[2];
if(n == 4)
u[4] = 2.0*u[3] – u[2];
else
u[4] = (y[4]-y[3])/(x[4] – x[3]);
}
else
{

u[1] = 2.0*u[2] – u[3];
u[0] = 2.0*u[1] – u[2];
u[4] = (y[3]-y[2])/(x[3] – x[2]);
}

}
else if(kk >=(n-3))
{

u[1] = (y[kk]-y[kk-1])/(x[kk] – x[kk-1]);
if(kk == (n-3))
{

u[3] = (y[n-1] – y[n-2])/(x[n-1] – x[n-2]);
u[4] = 2.0*u[3] – u[2];
if( n== 4)
u[0] = 2.0*u[1] – u[2];
else 
u[0] = (y[kk-1] – y[kk-2])/(x[kk-1] – x[kk-2]);

}
else
{

u[3] = 2.0*u[2] – u[1];
u[4] = 2.0*u[3] – u[2];
u[0] = (y[kk -1] – y[kk-2])/(x[kk-1] – x[kk -2]);
}

}
else
{

u[1] = (y[kk] – y[kk -1])/(x[kk] – x[kk-1]);
u[0] = (y[kk-1] – y[kk -2])/(x[kk-1] – x[kk-2]);
u[3] = (y[kk+2] – y[kk +1])/(x[kk+2] – x[kk +1]);
u[4] = (y[kk+3] – y[kk +2])/(x[kk+3] – x[kk +2]);
}

}

s[0] = fabs(u[3] – u[2]);
s[1] = fabs(u[0] – u[1]);
if( (s[0] + 1.0 == 1.0) && (s[1] + 1.0 == 1.0))
p = (u[1] + u[2])/2.0;
else
p = (s[0]*u[1] + s[1]*u[2])/(s[0] + s[1]);
s[0] = fabs(u[3] – u[4]);
s[1] = fabs(u[2] – u[1]);
if((s[0] + 1.0 == 1.0) && (s[1] + 1.0 == 1.0))
q = (u[2]+u[3])/2.0;
else 
q= (s[0]*u[2] + s[1]*u[3])/(s[0] + s[1]);

s[0] = y[kk];
s[1] = p;
s[3] = x[kk +1] – x[kk];
s[2] = (3.0*u[2] – 2.0*p – q)/s[3];
s[3] = (p + q – 2.0*u[2])/(s[3]*s[3]);
if (k<0)
{

p = t-x[kk];
s[4] = s[0] + s[1]*p + s[2]*p*p + s[3]*p*p*p;

}
return;

}

float Interp::interp2_onePoint(float *x,float *y,float *z,int m,int n,float a,float b,char *method)
{

//find a,b的位置
int tempi,tempj;
float w1,w2,w;
float tempx[2] = {0,0};
float tempz[2] = {0,0};
/*if(a>x[m-1] || b>y[n-1] || a<x[0] || b<y[0])
{

w = 0.0f;

}
else
{*/

tempi = findIndex(x,a,m);
tempj = findIndex(y,b,n);
if(tempi<0)
tempi =0;
if(tempi>m-2)
tempi = m-2;
if(tempj<0)
tempj = 0;
if(tempj>n-2)
tempj = n-2;

tempx[0] = x[tempi];
tempx[1] = x[tempi+1];
tempz[0] = z[tempj*m + tempi];
tempz[1] = z[tempj*m + tempi + 1];

interp_onePoint(tempx,tempz,2,a,&w1,method);
//update tempz
tempz[0] = z[(tempj+1)*m + tempi];
tempz[1] = z[(tempj+1)*m + tempi + 1];
interp_onePoint(tempx,tempz,2,a,&w2,method);

//update tempx and tempz
tempx[0] = y[tempj];
tempx[1] = y[tempj+1];

tempz[0] = w1;
tempz[1] = w2;

interp_onePoint(tempx,tempz,2,b,&w,method);
//}
return(w);
}

void Interp::interp2_multiPoint(float *x,float *y,float *z,int m,int n,float *a,float *b,int asize,int bsize,float *fval,char *method)
{

//grid nets
for(int i=0;i<asize;i++)
{

for(int j = 0;j<bsize;j++)
*(fval+j*asize+i) = interp2_onePoint(x,y,z,m,n,a[i],b[j],method);
}

}

/*
void Interp::interp2_multiPoint(float *x,float *y,float *z,int m,int n,float *gridx,float *gridy,int asize,int bsize,float *fval,char *method)
{

//grid nets
for(int i=0;i<asize;i++)
{

for(int j = 0;j<bsize;j++)
*(fval+j*asize+i) = interp2_onePoint(x,y,z,m,n,gridx[j*asize + i],gridx[j*asize + i],method);
}

}*/

int Interp::findIndex(float *x,int t,int n)
{

int kk,m,lc;
/*if(t<x[0])
return( 0);
if(t>=x[n-1])
return(n-1);*/

if(t<=x[1]) 
kk =0;
else if(t>=x[n-1])
kk = n-2;
else
{

kk=1;
m = n;
while (((kk -m) !=1) && ((kk – m) != -1))
{

lc=(kk +m)/2;
if(t<x[lc-1]) 
m = lc;
else
kk =lc;

}
kk = kk -1;
}

return(kk);

}

void Interp::interp_onePoint(float *x,float *y,int n,float t,float *fval,char *method)
{


if(nullptr == method)
*fval = lgr(x,y,n,t);
else if(_strcmpi(method,strMethod1) == 0)
*fval = lgr(x,y,n,t);
else if(_strcmpi(method,strMethod2) == 0)
*fval = lg3(x,y,n,t);
else if(_strcmpi(method,strMethod3) == 0)
{

static float s[5];
spl(x,y,n,-1,t,s);
*fval = s[4];
}
else
*fval = lgr(x,y,n,t);
}

void Interp::interp_multiPoint(float *x,float *y,int n,float *t,float *fval,int m,char *method)
{

if(t == NULL)
return;
float tempVal = 0.0;;
for(int k = 0;k<m;k++){


interp_onePoint(x,y,n,t[k],&tempVal,method);
fval[k] = tempVal;
}

}




测试:

#include””

int main()

{

 float t ,z;
static float x1[10] = {0.10,0.15,0.25,0.4,0.5,0.57,0.7,0.85,0.93,1.00};
static float y1[10] = {0.904837,0.860708,0.778801,0.670320,0.606531,
0.565525,0.496585,0.427415,0.394554,0.367879};

char *method = “spl”;
float t[3] = {0.63,0.71,0.84};
float multval[3] = {0};
float val;
interp.interp_onePoint(x1,y1,10,1.5,&val,method);
cout<<“t=10,”<<“z = “<<val<<endl;

getchar();


interp.interp_multiPoint(x1,y1,10,t,multval,3,method);

cout<<“t[0] = “<<t[0]<<“,”<<“val[0]=”<<multval[0]<<endl;
cout<<“t[1] = “<<t[1]<<“,”<<“val[1]=”<<multval[1]<<endl;
cout<<“t[2] = “<<t[2]<<“,”<<“val[2]=”<<multval[2]<<endl;
getchar();

 return 0;

}


版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/144931.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(1)


相关推荐

  • 查询数据库空间使用情况的函数_查看当前数据库

    查询数据库空间使用情况的函数_查看当前数据库sp_spaceused[[@objname=]'objname'][,[@updateusage=]'updateusage'][@objname

  • 文件无法保存(文件夹已损坏无法读取怎么办)

    要使用EmguCV,首先需要去官网下载安装包(除了下载安装包,还有其他几种方式)。安装完成后,得配置环境变量。这里不再赘述。网上的教程很多。我这里分享一个官网的下载链接:https://sourceforge.net/projects/emgucv/files/emgucv/这里使用窗体程序进行说明,首先新建个窗体1.引用EmguCV的DLL:usingEmgu.CV;usingEmg…

  • 小程序怎么开发自己的小程序_微信小程序建议使用

    小程序怎么开发自己的小程序_微信小程序建议使用微信小程序入门前言随着科技的不断进步,很多功能将会开放,那么很多需求也将会因为现实而得到满足,这是一种不需要下载和安装就可以使用的应用软件。用户只需扫描和搜索就可以打开应用程序。它很容易使用,而且很容易实现。小程序带来了巨大的流量,吸引了很多大的人和企业家前来追踪,也显示出它强大的生命力。小程序是下一个被确定为互联网新品种的程序,信已经成为不可缺少的交流工具,小程序依附于微信,用户搜索起来也会更方便,其实这就是小程序发展的前景和优势。提示:以下是本篇文章正文内容,下面案例可供参考一、小程序的概

  • 程序员在哪里可以接私活(程序员私活一个月能拿多少)

    最近很郁闷,一个女粉丝找我要开发一个系统。说是很着急。但是费用不高,说自己没钱。我平时事也很多,又不是很想接。说了一大堆苦情的话,然后说自己的要求不高,没有ui要求,我就接了。但是后来,越扯越严重……上升到600块要开发一个知乎的高度。这让我想起了预算茅草屋的价格,想要别墅的效果。扯皮扯的很累,项目我已经远程调试好了,也讲解了。最后全部退款了。关键是我还整理了很多讲解的说明:既然退款,项目就不是你的了,我开发的就是我的,项目的源码就开源吧。…

  • CentOS7 network.service loaded failed 处理技巧

    CentOS7 network.service loaded failed 处理技巧

    2021年10月18日
  • SplitContainer控件设置固定比例

    SplitContainer控件设置固定比例SplitContainer控件2个panel如何设置均等大小.(竖直拆分)先随意设置空间高的大小,然后如果想要均等显示,则直接设置SplitterDistance为高的一半,即可!随后任意改变控件大小,都不会改变均等显示比例。同理3:7,4:6…

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号