matlab的interp1函数用C++实现

好的，以下是一个简单的C++实现`interp1`函数的示例代码：

#include <iostream>
#include <vector>
#include <cmath>

using namespace std;

enum InterpMethod {
    InterpLinear,
    InterpCubic
};

// 线性插值
double InterpLinear(double x, const vector<double>& X, const vector<double>& Y)
{
    int n = X.size();
    int k = 0;
    for (int i = 1; i < n; ++i) {
        if (x > X[i]) {
            k = i;
        } else {
            break;
        }
    }

    double h = X[k+1] - X[k];
    double a = (X[k+1] - x) / h;
    double b = (x - X[k]) / h;

    double result = a * Y[k] + b * Y[k+1];

    return result;
}

// 三次样条插值
double InterpCubic(double x, const vector<double>& X, const vector<double>& Y, const vector<double>& M)
{
    int n = X.size();
    int k = 0;
    for (int i = 1; i < n; ++i) {
        if (x > X[i]) {
            k = i;
        } else {
            break;
        }
    }

    double h = X[k+1] - X[k];
    double a = (X[k+1] - x) / h;
    double b = (x - X[k]) / h;

    double result = a * Y[k] + b * Y[k+1] + ((a*a*a - a) * M[k] + (b*b*b - b) * M[k+1]) * h * h / 6.0;

    return result;
}

// interp1函数
vector<double> interp1(const vector<double>& X, const vector<double>& Y, const vector<double>& Xq, InterpMethod method = InterpLinear)
{
    int n = X.size();
    int nq = Xq.size();
    vector<double> Yq(nq);

    // 计算二阶导数
    vector<double> M(n);
    if (method == InterpCubic) {
        vector<double> h(n);
        for (int i = 0; i < n-1; ++i) {
            h[i] = X[i+1] - X[i];
        }

        vector<double> alpha(n);
        for (int i = 1; i < n-1; ++i) {
            alpha[i] = 3.0 / h[i] * (Y[i+1] - Y[i]) - 3.0 / h[i-1] * (Y[i] - Y[i-1]);
        }

        vector<double> l(n), mu(n), z(n);
        l[0] = 1.0;
        mu[0] = z[0] = 0.0;

        for (int i = 1; i < n-1; ++i) {
            l[i] = 2.0 * (X[i+1] - X[i-1]) - h[i-1] * mu[i-1];
            mu[i] = h[i] / l[i];
            z[i] = (alpha[i] - h[i-1] * z[i-1]) / l[i];
        }

        l[n-1] = 1.0;
        z[n-1] = 0.0;
        M[n-1] = 0.0;

        for (int i = n-2; i >= 0; --i) {
            M[i] = z[i] - mu[i] * M[i+1];
        }
    }

    // 插值
    for (int i = 0; i < nq; ++i) {
        double xq = Xq[i];
        if (xq < X[0]) {
            Yq[i] = Y[0];
        } else if (xq > X[n-1]) {
            Yq[i] = Y[n-1];
        } else {
            switch (method) {
                case InterpLinear:
                    Yq[i] = InterpLinear(xq, X, Y);
                    break;
                case InterpCubic:
                    Yq[i] = InterpCubic(xq, X, Y, M);
                    break;
                default:
                    break;
            }
        }
    }

    return Yq;
}

int main()
{
    // 定义数据点和插值点
    vector<double> X = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0};
    vector<double> Y = {0.0, 1.0, 8.0, 27.0, 64.0, 125.0};
    vector<double> Xq = {0.5, 1.5, 2.5, 3.5, 4.5};

    // 线性插值
    vector<double> Yq = interp1(X, Y, Xq, InterpLinear);
    cout << "线性插值：" << endl;
    for (int i = 0; i < Yq.size(); ++i) {
        cout << "Yq[" << i << "] = " << Yq[i] << endl;
    }

    // 三次样条插值
    Yq = interp1(X, Y, Xq, InterpCubic);
    cout << "三次样条插值：" << endl;
    for (int i = 0; i < Yq.size(); ++i) {
        cout << "Yq[" << i << "] = " << Yq[i] << endl;
    }

    return 0;
}

在这个示例代码中，我们通过 `interp1` 函数对一组数据点进行插值，可以选择线性插值或三次样条插值。对于线性插值，我们直接调用 `InterpLinear` 函数进行插值；对于三次样条插值，我们先计算二阶导数，然后调用 `InterpCubic` 函数进行插值。最后，我们可以将插值结果保存到一个向量 `Yq` 中，并输出结果。