[matlab] data interpolation

Data interpolation

The difference between interpolation and fitting ：

1. Implementation method ： Interpolation requires the curve to pass through the sample points , The fitting does not need to pass through the sample points , Only the minimum overall error is required .
2. The form of the result ： Interpolation is a piecewise approximation of sample points , There is no same approximation function ; Function fitting uses a function to approximate , There is a complete expression .
3. Focus on ： Interpolation can be used to estimate the function values corresponding to some points in the interval ; Fitting can not only estimate the points in the interval , You can also predict points outside the interval .
4. applications ： Interpolation is often used for precise data sets ; Fitting is often used in statistical data sets .

The role of data interpolation ： In a given number of data , Simulate a curve through these known data , We can get the data information of unknown points according to the curve .

interp1 function

Y1=interp1(X,Y,X1,method)： One dimensional interpolation function . according to X、Y Value , The calculation function is in X1 Place the value of the . among ,X、Y It's two known vectors of equal length , Represents the sample point and the sample value respectively .X1 It's a vector or a scalar , Represents the point to be interpolated .

method For the string , Optional values are ：

1. linear, Default , linear interpolation , Broken line .
2. nearest, Nearest interpolation , Stepped .
3. pchip, Piecewise cubic polynomial interpolation , The first derivative of the interpolation function of adjacent segments at several nodes is equal , It's smooth .
4. spline, Cubic spline interpolation , Piecewise construct cubic polynomials , Satisfy that there are continuous first-order and second-order derivatives at each node , Very smooth .

Why? ’pchip’ and ’spline’ Both interpolation methods use 3 Degree polynomial instead of higher degree ？

Polynomial degree is not the higher the better . The more times , The more likely it is to oscillate and deviate from the original function , This phenomenon is called Runge （Runge ） The phenomenon .

%%  All interpolation methods 
xdata=0:pi/6:2*pi;
ydata=sin(xdata);
x=0:pi/20:2*pi;
subplot(4,2,1)
y=interp1(xdata,ydata,x,'nearest');
plot(xdata,ydata,'p',x,y,'k-')
title('nearest')
subplot(4,2,2)
y=interp1(xdata,ydata,x,'next');
plot(xdata,ydata,'p',x,y,'k-')
title('next')
subplot(4,2,3)
y=interp1(xdata,ydata,x,'prevoius');
plot(xdata,ydata,'p',x,y,'k-')
title('prevoius')
subplot(4,2,4)
y=interp1(xdata,ydata,x,'linear');
plot(xdata,ydata,'p',x,y,'k-')
title('linear')
subplot(4,2,5)
y=interp1(xdata,ydata,x,'spline');
plot(xdata,ydata,'p',x,y,'k-')
title('spline')
subplot(4,2,6)
y=interp1(xdata,ydata,x,'pchip');
plot(xdata,ydata,'p',x,y,'k-')
title('pchip')
subplot(4,2,7)
y=interp1(xdata,ydata,x,'cubic');
plot(xdata,ydata,'p',x,y,'k-')
title('cubic')

Add ： Newton interpolation code ：（ I don't know the application scenario ）

function s=newtoninterp(xdata,ydata,x) % xdata and ydata For known point data ,x To request point data 
%Newton interpolation at notes (xdata,ydata)
n=length(xdata);
s=ydata(1);xx=1;
for k=2:n
    xx=xx.*(x-xdata(k-1));
    for l=k:n
        ydata(l)=(ydata(l)-ydata(k-1))/(xdata(l)-xdata(k-1));
    end
    s=s+xx*ydata(k);
end
end

Other usages and similar function references

interp2 function

interp2(x, y, z, x1, y1, method)： Two dimensional interpolation .x、y It can be either grid data or vector ,z Corresponding to the grid data z coordinate ,x1、y1 To request data point information .

application

Braking distance of motor vehicles

Driving , Start when the driver sees the obstacle , The shortest distance required to make a judgment and take braking measures to stop is called stopping sight distance . Stopping sight distance consists of three parts ： One is the driving distance within the driver's reaction time （ Reaction distance ); The second is the distance from the start of braking to the complete stop of the vehicle （ Braking distance ）; The third is the safe distance between the vehicle and the obstacles when the vehicle stops . among , The braking distance is mainly related to the driving speed and the type of road surface . According to the test , When a certain type of vehicle runs on Asphalt Pavement in wet weather , Its driving speed （ Company : km/h） Distance from braking （ Company ：m） The relationship is shown in the following table .

Speed	20	30	40	50	60	70	80	90	100	110	120	130	140	150
Braking distance	3.15	7.08	12.59	19.68	28.34	38.57	50.4	63.75	78.71	95.22	113.29	132.93	154.12	176.87

Suppose the reaction time of the driver is 10s, The safe distance is 10m. Excuse me, ：

① According to the actual vision and visual habits of a driver , The effective sight distance when driving is 120m, When it is driving on the road , The maximum speed per hour should not exceed （ Result rounding ）?

② If the data in the table is used as a reference , Design a maximum speed of 125km/h The highway , Then the designer should ensure that the visible distance of the driver at any point on the highway is how many meters ？

---

Set the speed to $v$, The stopping sight distance is $d$, The reaction distance is $d_1$, The braking distance is $d_2$, The safe distance is $d_3$, The reaction time is $a_s$, be
$$d=d_1+d_2+d_3$$
among ,$d_1=a_{s}v$,$d_2$ by $v$ Function of ,$d_3$ It is known that .

First question ： According to the actual vision and visual habits of a driver , When driving Effective sight distance （ That is, the maximum distance ）120m, When it is driving on the road , The maximum speed per hour should not exceed （ Result rounding ）?

The known reaction time is 10s, The safe distance is 10m, May adopt Method of solving equation :
$$10v+d_2+10=120$$
The problem is ,$d_2$ yes $v$ Function of , But the functional relationship is unknown , The equation is not solvable .

Now consider Data interpolation method , Take the data in the table as a sample , Perform data interpolation , To calculate the 120m The speed index corresponding to the stopping sight distance .

Programming ideas ：

First step ： Establish speed and stopping sight distance vectors .

The second step ： With 1 In units of , Interpolate all velocities in the sampling interval , Calculate the corresponding stopping sight distance .

The third step ： Calculate the stopping sight distance 120 Corresponding speed .

Step four ： Graphic display .

How to stop according to sight distance 120 Find the corresponding speed ？

First step ： Let the vector representing the stopping sight distance $d_i$ subtract 120, Then take the absolute value , Get a new vector x.
The second step ： take x In ascending order , And record the serial number of the smallest element , This serial number is the stopping sight distance 120 The corresponding velocity data is in the vector $v_i$ The serial number in .
The third step ： Obtain speed data according to serial number .

v=20:10:150;
vs=v*(1000/3600); %  Unit conversion  km/s m/s
d1=10*vs;
d2=[3.15,7.08,12.59,19.68,28.34,38.57,50.4,63.75,78.71,95.22,113.29,132.93,154.12,176.87];
d3=10;
d=d1+d2+d3;
vi=20:150;
di=interp1(v,d,vi,'spline'); %  interpolation 
%%  The following is a selection equal to 120 Corresponding abscissa method , There's not only one way ！
%%  Select and 120 The abscissa of the number with the smallest difference 
x=abs(di-120);
[y i]=sort(x);
vi(i(1)) %  The first is the smallest  % 36
plot(vi,di,vi(i(1)),di(i(1)),'rp') %  Points that meet the requirements 

Second questions ： Design a maximum speed of 125km/h The highway , Then the designer should ensure that the visible distance of the driver at any point on the highway is how many meters ？

j=find(vi==125); %  Find a top speed of 125 The index of the corresponding distance 
di(j)
plot(vi,di,125,di(j),'rp')

Sand table making problem

A ground force is divided into red and blue sides in a designated unfamiliar area （ Flat areas [0,2000]×[0,2000] Inside , Company ：m） Conduct operational exercises . During the exercise , The red investigation unit has measured the elevations of some places as shown in the following table .

① According to the data in the table , Making a military sand table .

② Within the scope of the exercise , The party occupying the largest highland will gain the advantage of condescending . Which area should the red side seize first .

---

Their thinking ：

First question ： Using two-dimensional interpolation to estimate data , To facilitate the production of military sand table .

Second questions ： On the basis of interpolation , Draw contour map , Find the biggest highland .

x=0:200:1800;
y=x;
z = [2000,2000,2001,1992,1954,1938,1972,1995,1999,1999;
     2000,2002,2006,1908,1533,1381,1728,1959,1998,2000;
     2000,2005,2043,1921, 977, 897,1310,1930,2003,2000; 
     1997,1978,2009,2463,2374,1445,1931,2209,2050,2003;
     1992,1892,1566,1971,2768,2111,2653,2610,2121,2007;
     1991,1875,1511,1556,2221,1986,2660,2601,2119,2007;
     1996,1950,1797,2057,2849,2798,2608,2303,2052,2003;
     1999,1999,2079,2685,3390,3384,2781,2165,2016,2000;
     2000,2002,2043,2271,2668,2668,2277,2049,2003,2000;
     2000,2000,2004,2027,2067,2067,2027,2004,2000,2000];
surf(x,y,z);
colormap(mycolor)

%%  After interpolation 
x1=0:100:1800; %  With 100 interval 
y1=x1';
z1=interp2(x,y,z,x1,y1,'spline');
figure(2)
surf(x1,y1,z1);
colormap(mycolor)

x1=0:50:1800; %  With 50 interval 
y1=x1';
z1=interp2(x,y,z,x1,y1,'spline');
figure(3)
surf(x1,y1,z1);
colormap(mycolor)

%%  Draw contours 
figure(4)
[c h] = contourf(x1, y1, z1, 10)
clabel(c, h) %  Displays the height represented by contour lines 
colorbar %  Display height color table 
colormap(mycolor)

%%  The following is just my favorite color map 
mycolorpoint=[[0 0 16];
    [8 69 99];
    [57 174 156];
    [198 243 99];
    [222 251 123];
    [239 255 190]];
mycolorposition=[1 11 33 50 57 64];
mycolormap_r=interp1(mycolorposition,mycolorpoint(:,1),1:64,'linear','extrap');
mycolormap_g=interp1(mycolorposition,mycolorpoint(:,2),1:64,'linear','extrap');
mycolormap_b=interp1(mycolorposition,mycolorpoint(:,3),1:64,'linear','extrap');
mycolor=[mycolormap_r',mycolormap_g',mycolormap_b']/255;
mycolor=round(mycolor*10^4)/10^4;% Retain 4 Decimal place 

contour Functions and contourf Function explanation

Here is only the basic usage ：

[c,h] = contourf(X,Y,Z,t); X、Y For grid data or vectors （ And mesh The parameter usage in is consistent ）,Z Is the function value ,t For the spacing between contour lines and the range of contour lines （ It can be understood that the larger the contour line, the denser the contour line ）; The return value does not know the meaning .

clabel(c, h) ： Dimension the height of the contour line （ Although I don't know the meaning of the return value , But it can be used directly ）.

colorbar： Color label .

contour Function and contourf Function usage is consistent , It's just contour Lines are displayed , There is no color fill between lines ; and contourf There is a color fill between the lines .

figure('NumberTitle', 'off', 'Name', 'contourf'); %  Modify the drawing window name 
contourf(x1, y1, z1, 10)
title('contourf')
colormap(mycolor)
figure('NumberTitle', 'off', 'Name', 'contour');
contour(x1, y1, z1, 10)
title('contour')
colormap(mycolor)