PCHIP armadillo函数

Question

我正在尝试将一些MATLAB代码转换成C++，使用的MATLAB函数之一是interp1。“足够简单”，我想，Armadillo有interp1。“线性应该足够好了”，好吧，我错了。所以我搜索interp1.m源代码，找到了Octave源代码，它使用的是pchip.m，我找到了octave源代码，但pchip.m使用的是pchip_deriv.cc，它似乎使用的是fortran fcn。

所以，在我真正深入到pchip的转换之前，有没有其他的库或源代码可以让我使用呢？

Answer 1

为了防止其他人需要这个，我在here上查找了方程式。对于大数据集，它似乎不能正常工作，因此使用矩阵公式，我创建了一个滑动窗口PCHIP interp。仍然与MATLAB不匹配，所以我用线性interp平均，它相当接近。这是正确的方法吗？可能不会..。它起作用了吗？是啊！

void pchip_slide(const arma::vec& x, const arma::vec& y, const arma::vec& x_new, arma::vec& y_custom,const int M=4){
vec y_interp=zeros<vec>(size(x_new));
interp1(x,y,x_new,y_interp);
int I=0;
int start_interp=0;
int end_interp=0;
for(int ii=0;ii<x_new.n_elem;ii++){
    
    I=index_min(abs(x-x_new(ii)));
    start_interp=std::max((I-2),0);
    end_interp=std::min((start_interp+M-1),int(numel(x)-1));
    
    vec x_mini=x(span(start_interp,end_interp));
    vec y_mini=y(span(start_interp,end_interp));
    mat x_mat=ones<mat>(x_mini.n_elem,x_mini.n_elem);
    for(int ll=0;ll<x_mini.n_elem-1;ll++){
        x_mat.row(ll)=pow((x_mini.t()),(x_mini.n_elem-ll));
    }
    vec c_mini=solve( x_mat, y_mini,solve_opts::fast + solve_opts::no_approx);
    rowvec x_pchip_mini=ones<rowvec>(x_mini.n_elem);
    for (int ll=0;ll<x_mini.n_elem-1;ll++){
        x_pchip_mini(ll)=pow((x_new(ii)),(x_mini.n_elem-ll));
    }
    y_custom(ii)=conv_to<double>::from(x_pchip_mini*c_mini);
    if ((x_new(ii)>=x(0))&& (x_new(ii)<=x(x.n_elem-1))){
        y_custom(ii)=(y_custom(ii)*1/M+y_interp(ii)*3/M);
    }
}
return;
}