matlab interp1 linear interpolation/extrapolation to C code

最新推荐文章于 2024-01-01 11:37:52 发布

chongchong.ding

最新推荐文章于 2024-01-01 11:37:52 发布

阅读量442

点赞数

分类专栏： Note 文章标签： matlab c语言开发语言

本文链接：https://blog.csdn.net/weixin_41277132/article/details/127869569

版权

Note 专栏收录该内容

2 篇文章 0 订阅

订阅专栏

该程序使用C语言实现了一维数据的插值与外推功能。通过对给定的数据集进行线性插值或外推，计算出新的插值点的值。程序读取二进制文件中的数据，调用`Interp1`函数处理数据，并将结果与已计算的值进行比较。适用于数据处理和科学计算场景。

摘要由CSDN通过智能技术生成

#include <stdio.h>
#include <windows.h>
#include "float.h"
void Interp1(double * x, double * y,int n, double *xn,double *yn,int len)
{
    double max,min;
    //to find the max&min value of x-series
    if(x[0]>x[n-1]){
        max=x[0];
        min=x[n-1];
    }else{
        min=x[0];
        max=x[n-1];
    }


    for (int i = 0; i < len; i++) {
        int j;
        //in case the interoplation point alraedy exists in the array
        for (j = 0; j < n; j++)
            if(xn[i]==x[j])
                break;
        if(j!=n)
            yn[i] = x[j];
        else{
            // we need linear interpolation/extrapolation.
            if((xn[i]>=min)&&(xn[i]<=max)){

                for (int idx = 0; idx < n; idx++)
                {
                    if (xn[i] >= x[idx] && xn[i] <= x[idx+1])
                    {
                        //y = y0 + (y1-y0)*(x-x0)/(x1-x0);
                        double x0 = x[idx];
                        double x1 = x[idx + 1];
                        double y0 = y[idx];
                        double y1 = y[idx + 1];
                        yn[i] = y0 + (y1 - y0)*(xn[i] - x0) / (x1 - x0);
                    }
                }

            }else{//extrapolate
                int idx;
                //y = y0 + ((x - x0) / (x1 - x0)) * (y1 - y0)
                if (xn[i] >= x[n-1])
                {
                    idx = n-1;
                }
                else
                {
                    idx = 0;
                }
                double y0 = y[idx];
                double y1 = y[idx + 1];
                double x0 = x[idx];
                double x1 = x[idx + 1];
                yn[i] = y0 + ((xn[i] - x0) / (x1 - x0)) * (y1 - y0);
            }
        }

    }

}
double array[188<<2];
double yn_c[188];
int main()
{
    FILE *fid;
    fid = fopen("F:\\osphotonics\\matlab\\sync\\interp_data.bin","rb");
    fread(array,8,188<<2,fid);
    fclose(fid);
    double *x;
    double *y;
    double *xn;
    double *yn;
    x = &array[188*0];
    y = &array[188*1];
    xn = &array[188*2];
    yn = &array[188*3];
    Interp1(x,y,188,xn,yn_c,188);
    for (int i = 0; i < 188; ++i) {
        printf("[%d]x=%.8llf y=%.8llf xn=%.8llf yn_mat=%.8llf yn_calc=%.8llf (diff)=%.8llf\n",i,x[i],y[i],xn[i],yn[i],yn_c[i],yn[i]-yn_c[i]);
    }

    return 0;
}