大家好,又见面了,我是你们的朋友全栈君。
头文件
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账号...