1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
|
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA - Farid BELAHCENE
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
function yi=interp1(varargin)
// yi=interp1(x,y,xi[,method[,interpolation)
// This function performs the yi values corresponding to xi by interpolation defined by x and y.
// Inputs :
// x , y : given data, x is a reals vector, y is a vector, matrix, or hypermatrix of reals
// if y is a vector, the length of x must be equal to the length of y,
// else the size of the first dimension of y must be equal to length of x.
// xi : a vector, matrix, or hypermatrix of reals
// Output
// yi : reals vector, matrix or hypermatrix, the values corresponding to xi by interpolation defined by x and y
// if size(y)=[C,N1,N2,N3,....] and size(xi)=[M1,M2,M3,M4] then size(xi)=[M1,M2,M3,M4,N1,N2,N3,N4,..], and length of x must be equal to C.
// Several kind of intepolations may be computed by selecting the appropriate method parameter:
// The methods are:
// linear : this is the default method (using the interp Scilab function)
// spline : this is the cubic spline interpolation (using interpln and splin Scilab functions)
// nearest : yi take the values corresponding to the nearest neighbor of xi
//
// Several kind of extrapolations may be computed :
// 'extrap' : the extrapolation points is performed by the defined method
// real value : you can choose a real value for extrapolation, in this way yp(i) takes this value for xp(i) not in [x1,xn] interval, for example 0 (but also nan or inf).
// by default the extrapolation is performed by the defined method (for spline method), and by nan for linear and nearest method.
// F.B
rhs=size(varargin)
// 2 <= Number of input arguments <= 5
if rhs<2 | rhs>5 then
error(msprintf(gettext("%s: Wrong number of input arguments: Must be between %d and %d.\n"),"interp1",2,5));
end
//if yi=interp1(x,y,xi,..) not change
//if yi=interp1(y,xi,...) replace input argument by yi=interp1(x,y,xi,..),whith x=1:size(y,1) by default
if rhs==2 | (rhs>2 & type(varargin(3))==10) then
if isvector(varargin(1)) then
X=1:size(varargin(1),"*")
elseif size(size(varargin(1)),"*")==2 then
if (size(varargin(1),1)>1 & size(varargin(1),2)>1) then
X=1:size(varargin(1),1)
end
else
error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector or matrix expected.\n"),"interp1",2));
end
for i=rhs:-1:1
varargin(i+1)=varargin(i)
end
varargin(1)=X
end
//............................
//ininialisation of xi
//............................
//xi components must be reals
xi=varargin(3)
if type(xi)<>1 then
if typeof(xi)=="hypermat" then
if type(xi.entries)<>1 then
error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of reals expected.\n"),"interp1",3));
end
else
error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of reals expected.\n"),"interp1",3));
end
end
//delete the dimension of xi equal to one after the second dimension
//or the first dimension
xisize=size(xi);
while size(xisize,"*")>=2 & xisize($)==1
xisize=xisize(1:$-1);
end
xisizetemp=xisize
if size(xisize,"*")>=2 then
if xisize(1)==1 then
xisize=xisize(2:$);
end
end
//-------------------------
//Initialisation of x, y
//-------------------------
x=varargin(1);
y=varargin(2);
//x must be real vector
if type(x)<>1 then
error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of reals expected.\n"),"interp1",1));
end
//y components must be reals
if type(y)<>1 then
if typeof(y)=="hypermat" then
if type(y.entries)<>1 then
error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of reals expected.\n"),"interp1",2));
end
else
error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of reals expected.\n"),"interp1",2));
end
end
//verification of x,y line/column
if isvector(x) then
if find(isnan(x))<>[] then
error(msprintf(gettext("%s: Wrong value for input argument #%d: Reals expected but some NaN found.\n"),"interp1",1));
end
if isvector(y) then
if size(x,"*")<>size(y,"*") then
error(msprintf(gettext("%s: Wrong size for input arguments #%d and #%d: Same size expected.\n"),"interp1",1,2));
end
elseif size(size(y),"*")>=2 then
if size(x,"*")<>size(y,1) then
error(msprintf(gettext("%s: Wrong size for input arguments #%d and #%d: Same size expected.\n"),"interp1",1,2));
end
else
error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector or matrix expected.\n"),"interp1",2));
end
else
error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector expected.\n"),"interp1",1));
end
// xi : increase order sorting (for xi)
[xtemp,p]=gsort(matrix(x,1,-1),"c","i")
x=matrix(xtemp,size(x))
x=matrix(x,1,-1)
if isvector(y) then
y=y(p)
elseif size(size(y),"*") then
for l=1:size(y,"*")/size(y,1)
y(:,l)=y(p,l)
end
else
error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector or matrix expected.\n"),"interp1",2));
end
//-------------------------------------------------
// CASE : 3 inputs arguments : yi=interp1(x,y,xi)
//-------------------------------------------------
//default method : linear method is used
if size(varargin)==3 then
yi=interp1(x,y,xi,"linear",%nan)
end
//--------------------------------------------------
// CASE : 4 inputs arguments : yi=interp1(x,y,xi,method)
//--------------------------------------------------
if size(varargin)==4 then
select part(varargin(4),1)
//-------------------------------------------
// Linear method : yi=linear(x,y,xi,'linear')
//-------------------------------------------
// the values of extrapolation points are nan for linear method
case "l"
yi=interp1(x,y,xi,"linear",%nan)
//-------------------------------------------
// Spline method yi=interp1(x,y,xi,'spline')
//-------------------------------------------
// the extrapolation used the spline method
case "s"
if xi==[] then
yi=[]
return
end
yi=interp1(x,y,xi,"spline","extrap")
//----------------------------------------------
// Nearest method yi=interp1(x,y,xi,'nearest')
//----------------------------------------------
// the values of extrapolation points are nan for nearest method
case "n"
if xi==[] then
yi=[]
return
end
yi=interp1(x,y,xi,"nearest",%nan)
else
error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or ''%s'' or ''%s'' expected.\n"),"interp1",4,"linear","nearest"));
end
end
//-------------------------------------------------------------------------------------------------
// CASE : 5 inputs arguments : yi=interp1(x,y,xi,method,'extrap') or yi=interp1(x,y,xi,method,extrapval)
//-------------------------------------------------------------------------------------------------
if size(varargin)==5 then
select part(varargin(4),1)
//-----------------------------------------------------------------------------------
// Linear method : yi=linear(x,y,xi,'linear','extrap') or yi=interp1(x,y,xi,method,extrapval)
//------------------------------------------------------------------------------------
case "l"
xitemp=matrix(xi,-1,1)
// y is a vector
if isvector(y) then
yi=hypermat(size(xitemp))
[x,ind]=gsort(matrix(x,1,-1),"c","i")
if varargin(5)==%nan then
yi=linear_interpn(xitemp,x,y(ind),"by_nan");
end
if type(varargin(5))==10 then
if varargin(5)<>"extrap" then
error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
else
yi=linear_interpn(xitemp,x,y(ind),"natural");
end
elseif type(varargin(5))==1 then
yi=linear_interpn(xitemp,x,y(ind),"by_nan");
if ~isnan(varargin(5)) then
k=find(xitemp>max(x)|xitemp<min(x))
yi(k)=varargin(5)
end
end
if size(xisize,"*")>=2
yi=matrix(yi,xisize)
else
yi=matrix(yi,xisizetemp)
end
// y is matrix or hypermatrix
elseif size(size(y),"*")>=2 then
ysize=size(y)
ky=ysize(2:$)
yi=hypermat([size(xitemp),ky])
[x,ind]=gsort(matrix(x,1,-1),"c","i")
//extrapolation
if type(varargin(5))==10 then
if varargin(5)<>"extrap" then
error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
else
if xitemp==[] then
yi=[]
return
end
for l=1:size(y,"*")/size(y,1)
ytemp=y(:,l)
yi(:,l)=matrix(linear_interpn(xitemp,x,ytemp(ind),"natural"),size(xitemp))
end
end
elseif type(varargin(5))==1 then
if xitemp==[] then
yi=[]
return
end
for l=1:size(y,"*")/size(y,1)
ytemp=y(:,l)
yi(:,l)=matrix(linear_interpn(xitemp,x,ytemp(ind),"by_nan"),size(xitemp))
end
if ~isnan(varargin(5)) then
k=find(xitemp>max(x)|xitemp<min(x))
yi(k,:)=varargin(5)
end
end
yi=matrix(yi,[xisize,ky])
else
error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector or matrix expected.\n"),"interp1",2));
end
//-------------------------------------------------------------------------------------
// Spline method yi=interp1(x,y,xi,'spline','extrap') or yi=interp1(x,y,xi,'spline',extrapval)
//-------------------------------------------------------------------------------------
case "s"
if xi==[] then
if varargin(5)=="extrap"|type(varargin(5))==1 then
yi=[]
return
else
error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
end
end
xitemp=matrix(xi,-1,1)
//y is a vector
if isvector(y) then
yi=hypermat(size(xitemp))
yi=interp(xitemp,matrix(x,1,-1),matrix(y,1,-1),splin(matrix(x,1,-1),matrix(y,1,-1)),"natural");
if type(varargin(5))==10 then
if varargin(5)<>"extrap" then
error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
end
elseif type(varargin(5))==1 then
k=find(xitemp>max(x)|xitemp<min(x))
yi(k)=varargin(5)
end
if size(xisize,"*")>=2
yi=matrix(yi,xisize)
else
yi=matrix(yi,xisizetemp)
end
//y is a matrix or a hypermatrix
elseif size(size(y),"*")>=2 then
ky=size(y)
ky=ky(2:$)
yi=hypermat([size(xitemp),ky])
for l=1:size(y,"*")/size(y,1)
yi(:,l)=matrix(interp(matrix(xi,-1,1),matrix(x,-1,1),y(:,l),splin(matrix(x,-1,1),y(:,l)),"natural"),size(xitemp))//les composante de yi
end
//extrapolation
if type(varargin(5))==10 then
if varargin(5)<>"extrap" then
error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
end
elseif type(varargin(5))==1 then
k=find(xitemp>max(x)|xitemp<min(x))
yi(k,:)=varargin(5)
end
yi=matrix(yi,[xisize,ky])
else
error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector or matrix expected.\n"),"interp1",2));
end
//---------------------------------------------------------------------------------------
// Nearest method yi=interp1(x,y,xi,'nearest','extrap') or yi=interp1(x,y,xi,'nearest',extrapval)
//---------------------------------------------------------------------------------------
case "n"
//if all xi values are nan, retuns nan values for yi
if size(find(isnan(xi)),"*")==size(xi,"*") then
if varargin(5)=="extrap"|type(varargin(5))==1 then
yi=xi
return
else
error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
end
end
if xi==[] then
if varargin(5)=="extrap"|type(varargin(5))==1 then
yi=[]
return
else
error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
end
end
//y is vector
if isvector(y) then
xitemp=matrix(xi,1,-1)
knan=find(isnan(xitemp))
knotnan=find(~isnan(xitemp))
[xitemp,p]=gsort(matrix(xitemp(knotnan),1,-1),"c","i")
yi=matrix(xi,1,-1)
k=zeros(xitemp)
x_size=size(x,"*")
j=size(xitemp,"*")
i=x_size
while j>=1 & i>=1
if xitemp(j)>=x(i) then
if i<>x_size then
k(j)=i
end
j=j-1
else
i=i-1
end
end
k(xitemp<x(1)) = 1;
k(xitemp>=x(x_size)) = x_size-1;
i = find(xitemp >= matrix((x(k)+x(k+1))/2,size(k)));
k(i) = k(i)+1;
yi=y(k)
yi=matrix(yi,1,-1)
//extrapolation
if type(varargin(5))==10 then
if varargin(5)<>"extrap" then
error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
end
elseif type(varargin(5))==1 then
n=find(xitemp>max(x)|xitemp<min(x)|isnan(xitemp)|isnan(xitemp))
yi(n)=varargin(5)
else
error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
end
yitemp=yi
yi(p)=yitemp
ytemp=yi
yi=matrix(xi,1,-1)
yi(knan)=%nan
yi(knotnan)=ytemp
if size(xisize,"*")>=2
yi=matrix(yi,xisize)
else
yi=matrix(yi,xisizetemp)
end
//y is a matrix or a hypermatrix
elseif size(size(y),"*")>=2 then
xitemp=matrix(xi,1,-1)
knan=find(isnan(xitemp))
knotnan=find(~isnan(xitemp))
[xitemp,p]=gsort(xitemp(knotnan),"c","i")
ind=size(y)
ind=ind(2:$)
yi=hypermat([size(xitemp,"*"),ind])
k=zeros(xitemp)
x_size=size(x,"*")
j=size(xitemp,"*")
i=x_size
while j>=1 & i>=1
if xitemp(j)>=x(i) then
if i<>x_size then
k(j)=i
end
j=j-1
else
i=i-1
end
end
k(xitemp<x(1)) = 1;
k(xitemp>=x(x_size)) = x_size-1;
i = find(xitemp >= matrix((x(k)+x(k+1))/2,size(k)));
k(i) = k(i)+1;
for l=1:size(y,"*")/size(y,1)
ytemp=matrix(y(:,l),1,-1)
yi(:,l) =ytemp(k)
end
//extrapolation
if type(varargin(5))==10 then
if varargin(5)<>"extrap" then
error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
end
elseif type(varargin(5))==1 then
n=find(xitemp>max(x)|xitemp<min(x))
for l=1:size(y,"*")/size(y,1)
yi(n,l)=varargin(5)
end
else
error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
end
yitemp=yi
for l=1:size(y,"*")/size(y,1)
yi(p,l)=yitemp(:,l)
end
yitemp=yi
yi=hypermat([size(xi,"*"),ind])
for l=1:size(y,"*")/size(y,1)
yi(knan,l)=%nan
yi(knotnan,l)=yitemp(:,l)
end
yi=matrix(yi,[xisize,ind])
else
error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector or matrix expected.\n"),"interp1",2));
end
else
error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
end
end
endfunction
|
