* `unique` is enabled for any 2D sparse arrays, in simple, 'c' and 'r' modes.
* `%chars` constant added, to easily access to some selected sets of unicode symbols.
* Lists are displayed in a more compact and comprehensive way.
+* `interp1` is upgraded:
+ - complex numbers `y` now supported: the real and imaginary parts are interpolated separately.
+ - extrapolation option extended: `edgevalue` mode added for all interpolations; `periodic` mode added for linear and spline interpolations. `linear` mode added for the spline interpolations.
Help pages:
-----------
* [#16553](https://bugzilla.scilab.org/16553): `unique(["" ""])` returned `["" ""]`.
* [#16557](https://bugzilla.scilab.org/16557): `macr2tree` + `tree2code` translated `e={2}` into `"e=1"` and `e={2,"ab"}` into `"e=[2,"ab"]"`.
* [#16559](https://bugzilla.scilab.org/16553): `isempty(A)` was true for sparse matrix of dimension 2^16 or larger.
+<<<<<<< HEAD
* [#16622](https://bugzilla.scilab.org/16622): `inv` could no longer be overloaded for hypermatrices of decimal or complex numbers.
* [#16623](https://bugzilla.scilab.org/16623): `rand(2,2,2)^2` yielded a wrong result instead of trying to call the `%s_p_s` overload for input hypermatrices.
* [#16644](https://bugzilla.scilab.org/16644): `input("message:")` yielded a wrong error message about `mprintf` in case of non-interpretable input.
* [#16654](https://bugzilla.scilab.org/16654): `interp` was leaking memory.
+=======
+* [#16629](https://bugzilla.scilab.org/16629): `interp1`'s documentation did not tell the spline edges conditions ; extrapolation modes were poorly explained ; the description of the result's size was completely wrong ; x as an option was not documented. A wrong extrapolation value could silently return a wrong result. There was some dead code like `if varargin(5)==%nan`. A bugged error message yielded its own error. When x is implicit, the argument index in error messages could be wrong. `periodic` and `edgevalue` extrapolation modes were not available. `linear` extrapolation was not available for splines. When `xp` is an hypermatrix with `size(xp,1)==1`, the size of the result was irregular/wrong.
+>>>>>>> ae66d0ee461... * Bug 16629 fixed: interp1() fixed + complex + extended extrapolations
### Bugs fixed in 6.1.0:
hsvcolormap_1.png=11918d88bcc793200af0b9f1b58b0554
iir_1.png=e675c2755f68ddc4611c849895b63012
intdec_1.png=da3896e7d2e8468dd33edf57ccbe4480
-interp1_1.png=0e9a3f4319b2818ce4921b9bc3008d80
+interp1_1.png=e76a8e0a6070fbe078060844707c49fe
interp2d_1.png=f4af61bc3faf493d778169ec7decc7ae
interp2d_2.png=a62363849a3a9a3571dc1942fcebfbd3
interp2d_3.png=a9f2a0b0e59ac8cedae3ad22e0b2cc46
<!--
* Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
* Copyright (C) 2008 - INRIA - Farid BELAHCENE
- *
* Copyright (C) 2012 - 2016 - Scilab Enterprises
+ * Copyright (C) 2021 - Samuel GOUGEON - Le Mans Université
*
* This file is hereby licensed under the terms of the GNU GPL v2.0,
* pursuant to article 5.3.4 of the CeCILL v.2.1.
xmlns:scilab="http://www.scilab.org" xml:id="interp1" xml:lang="en">
<refnamediv>
<refname>interp1</refname>
- <refpurpose>one_dimension interpolation function</refpurpose>
+ <refpurpose>1D interpolation in nearest, linear or spline mode</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Syntax</title>
<synopsis>
+ yp = interp1(y, xp)
yp = interp1(x, y, xp)
- yp = interp1(x, y, xp, method)
- yp = interp1(x, y, xp, method, extrapolation)
+ yp = interp1(.., xp, method)
+ yp = interp1(.., xp, method, extrapolation)
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Arguments</title>
<variablelist>
<varlistentry>
- <term>xp</term>
+ <term>x</term>
<listitem>
- <para>reals scalar, vector or matrix (or hypermatrix)</para>
+ vector of at least 2 real numbers: Abscissas of known interpolation nodes,
+ without duplicates. By default,
+ <itemizedlist>
+ <listitem>
+ if <varname>y</varname> is a vector: <literal>x=1:length(y)</literal>.
+ </listitem>
+ <listitem>
+ if <varname>y</varname> is a matrix or an hypermatrix:
+ <literal>x=1:size(y,1)</literal>.
+ </listitem>
+ </itemizedlist>
+ <para/>
</listitem>
</varlistentry>
<varlistentry>
- <term>x</term>
+ <term>y</term>
<listitem>
- <para>reals vector</para>
+ vector, matrix or hypermatrix of real or complex numbers: values at known
+ interpolation nodes, at the corresponding <varname>x</varname> abscissas.
+ <itemizedlist>
+ <listitem>
+ if <varname>y</varname> is a vector, <varname>x</varname> and
+ <varname>y</varname> must have the same length.
+ </listitem>
+ <listitem>
+ if <varname>y</varname> is a matrix or an hypermatrix, we must have
+ <literal>length(x)==size(y,1)</literal>. Each column of
+ <varname>y</varname> is then interpolated versus the same
+ <varname>x</varname> abscissas, for the given <varname>xp</varname>.
+ </listitem>
+ </itemizedlist>
+ <para/>
</listitem>
</varlistentry>
<varlistentry>
- <term>method</term>
+ <term>xp</term>
<listitem>
- <para>(optional) string defining the interpolation method</para>
+ scalar, vector, matrix or hypermatrix or decimal numbers: abscissas of
+ points whose values <varname>yp</varname> must be computed according
+ to data of interpolation nodes.
+ <para/>
</listitem>
</varlistentry>
<varlistentry>
- <term>extrapolation</term>
+ <term>yp</term>
<listitem>
- <para>(optional) string, or real value defining the yp(j) components
- for xp(j) values outside [x1,xn] interval.
- </para>
+ vector, matrix, or hypermatrix of numbers: interpolated <varname>y</varname>
+ values at the given <varname>xp</varname>.
+ <itemizedlist>
+ <listitem>
+ if <varname>y</varname> is a vector: <varname>yp</varname> has the
+ size of <varname>xp</varname>.
+ <para/>
+ </listitem>
+ <listitem>
+ if <varname>y</varname> is a matrix or an hypermatrix:
+ <itemizedlist>
+ <listitem>
+ if <varname>xp</varname> is a scalar or a vector:
+ <literal>size(yp)</literal> is <literal>[length(xp) size(y)(2:$)]</literal>
+ </listitem>
+ <listitem>
+ if <varname>xp</varname> is a matrix or an hypermatrix:
+ <literal>size(yp)</literal> is <literal>[size(xp) size(y)(2:$)]</literal>
+ </listitem>
+ </itemizedlist>
+ </listitem>
+ </itemizedlist>
+ <para/>
</listitem>
</varlistentry>
<varlistentry>
- <term>yp</term>
+ <term>method</term>
<listitem>
- <para>vector, or matrix (or hypermatrix)</para>
+ string defining the interpolation method. Possible values and processing are:
+ <para/>
+ <table>
+ <tr>
+ <td><emphasis role="bold">"linear"</emphasis>:</td>
+ <td>
+ linear interpolation between consecutive nodes, used by default.
+ </td>
+ </tr>
+ <tr>
+ <td><emphasis role="bold">"spline"</emphasis>:</td>
+ <td>
+ interpolation by cubic splines
+ </td>
+ </tr>
+ <tr>
+ <td><emphasis role="bold">"nearest"</emphasis>:</td>
+ <td>
+ <para>for each value xp(j), yp(j) takes the value or y(i)
+ corresponding to x(i) the nearest neighbor of xp(j)
+ </para>
+ </td>
+ </tr>
+ </table>
+ <para/>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>extrapolation</term>
+ <listitem>
+ string or number defining the yp(j) components for xp(j) values outside the
+ [x(1)=min(x),x($)=max(x)] interval. We suppose here-below that <varname>x</varname>
+ and <varname>y</varname> have already been sorted accordingly.
+ <para/>
+ <table>
+ <tr>
+ <td><emphasis role="bold">"extrap"</emphasis>:</td>
+ <td>
+ <literal>interp1(x,y,xp, method, "extrap")</literal> is equivalent
+ to <literal>interp1(x,y,xp, method, method)</literal>.
+ </td>
+ </tr>
+ <tr>
+ <td><emphasis role="bold">"linear"</emphasis>:</td>
+ <td>
+ Can be used with the <literal>"spline"</literal> (and obviously
+ <literal>"linear"</literal>) interpolation methods.
+ </td>
+ </tr>
+ <tr>
+ <td><emphasis role="bold">"periodic"</emphasis>:</td>
+ <td>
+ This extrapolation type can be used with the
+ <literal>"linear"</literal> or <literal>"spline"</literal>
+ interpolation methods. Then: if <varname>y</varname> is a
+ vector, <literal>y(1)==y($)</literal> is required ; otherwise
+ <literal>y(1,:)==y($,:)</literal> is required.
+ </td>
+ </tr>
+ <tr>
+ <td><emphasis role="bold">"edgevalue"</emphasis>:</td>
+ <td>
+ Then <literal>yp(i)=y(1)</literal> for every
+ <literal>xp(i)<x(1)</literal>, and
+ <literal>yp(i)=y($)</literal> for every <literal>xp(i)>x($)</literal>.
+ </td>
+ </tr>
+ <tr>
+ <td><emphasis role="bold">padding</emphasis>:</td>
+ <td>
+ <varname>padding</varname> is a decimal or complex number
+ used to set <literal>yp(i)=padding</literal> for every
+ <literal>xp(i) ∉ [min(x),max(x)]</literal>. Example:
+ <literal>yi=interp1(x,y,xp,method, 0)</literal>.
+ </td>
+ </tr>
+ <tr>
+ <td><emphasis role="bold">(none)</emphasis>:</td>
+ <td>
+ By default, the extrapolation is performed by splines when
+ splines are used for the interpolation, and by padding
+ with %nan when the interpolation is linear or by "nearest" node.
+ </td>
+ </tr>
+ </table>
</listitem>
</varlistentry>
</variablelist>
<refsection>
<title>Description</title>
<para>
- Given <literal>(x,y,xp)</literal>, this function performs the yp
- components corresponding to xp by the interpolation(linear by default)
- defined by x and y.
+ Given <literal>(x,y,xp)</literal>, this function computes the yp
+ components corresponding to xp by the interpolation between known
+ data provided by (x,y) nodes.
</para>
<para>
- If <literal>yp</literal> is a vector then the length of xp is equal
- to the length of <literal>yp,</literal> if <literal>yp</literal> is a
- matrix then <literal>xp</literal> have the same length than the length of
- each columns of yp, if <literal>yp</literal> is a hypermatrix then the
- length of <literal>xp</literal> have the same length than the first
- dimension of <literal>yp. </literal>
- </para>
- <para>If size(y)=[C,M1,M2,M3,...,Mj] and size(xp)=[N1,N2,N3,...,Nk] then
- size(yp)=[N1,N2,..,Nk,M1,M2,...Mj] and length of x must be equal to
- size(y,1)
+ <varname>x</varname> is priorly sorted in ascending order, and
+ <varname>y</varname> values or per column are then sorted accordingly.
</para>
<para>
- The <literal>method</literal> parameter sets the evaluation rule for
- interpolation
+ <emphasis role="bold">Interpolation of complex values</emphasis>: When
+ <varname>y</varname> is complex, its real and imaginary parts are interpolated
+ separately, and then added to build the complex <varname>yp</varname>.
</para>
- <variablelist>
- <varlistentry>
- <term>"linear"</term>
- <listitem>
- <para>
- the interpolation is defined by linear method (see <link linkend="interpln">interpln</link>)
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>"spline"</term>
- <listitem>
- <para>the interpolation is defined by cubic spline interpolation (
- see <link linkend="splin">splin</link> , <link linkend="interp">interp</link>)
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>"nearest"</term>
- <listitem>
- <para>for each value xp(j), yp(j) takes the value or y(i)
- corresponding to x(i) the nearest neighbor of xp(j)
- </para>
- </listitem>
- </varlistentry>
- </variablelist>
<para>
- The <literal>extrapolation</literal> parameter sets the evaluation
- rule for extrapolation, i.e for <literal>xp(i) </literal>not in [x1,xn]
- interval
+ <emphasis role="bold">interp1(x,y,xp,"nearest")</emphasis>: For any xp
+ at the middle of an <literal>[x(i),x(i+1)]</literal> interval, the upper
+ bound <literal>x(i+1)</literal> is considered as the nearest
+ <varname>x</varname> value, and <literal>yp=y(i+1)</literal> is assigned.
</para>
- <variablelist>
- <varlistentry>
- <term>"extrap"</term>
- <listitem>
- <para>the extrapolation is performed by the defined method.
- yp=interp1(x,y,xp,method,"extrap")
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>real value</term>
- <listitem>
- <para>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). yi=interp1(x,y,xp,method, 0)
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>by default</term>
- <listitem>
- <para>the extrapolation is performed by the defined method (for
- spline method), and by nan for linear and nearest methods.
- yp=interp1(x,y,xp,method)
- </para>
- </listitem>
- </varlistentry>
- </variablelist>
+ <refsect3>
+ <title>linear interpolations</title>
+ They are performed through the <literal>linear_interpn(..)</literal> function,
+ with the corresponding <literal>"edgevalue"→"C0"</literal>,
+ <literal>"linear"→"natural"</literal>, <literal>"periodic"→"periodic"</literal>
+ extrapolation option.
+ </refsect3>
+ <refsect3>
+ <title>spline interpolations</title>
+ <para>
+ <emphasis role="bold">interp1(..,xp,"spline")</emphasis> or
+ <emphasis role="bold">interp1(..,xp,"spline","spline")</emphasis> or
+ <emphasis role="bold">interp1(..,xp,"spline","extrap")</emphasis>
+ use <literal>not_a_knot</literal> edges conditions. Extrapolation is performed
+ by using both spline polynomials computed at the <literal>(x,y)</literal> edges.
+ </para>
+ <para>
+ <emphasis role="bold">interp1(..,xp,"spline","edgevalue")</emphasis>
+ uses <literal>not_a_knot</literal> edges conditions and then calls
+ <literal>interp(..,"C0")</literal> to perform the actual interpolation
+ and extrapolation.
+ </para>
+ <para>
+ <emphasis role="bold">interp1(..,xp,"spline","periodic")</emphasis>
+ calls both <literal>splin(..)</literal> and then <literal>interp(..)</literal>
+ with their <literal>"periodic"</literal> option.
+ </para>
+ <para>
+ <emphasis role="bold">interp1(..,xp,"spline","linear")</emphasis>
+ calls <literal>splin(..,"natural")</literal> for linear edges conditions,
+ and then feeds <literal>interp(..,"linear")</literal>.
+ </para>
+ </refsect3>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
-x = linspace(0, 3, 20);
-y = x.^2;
-xx = linspace(0, 3, 100);
-yy1 = interp1(x, y, xx, 'linear');
-yy2 = interp1(x, y, xx, 'spline');
-yy3 = interp1(x, y, xx, 'nearest');
-plot(xx', [yy1; yy2; yy3]', x, y, '*')
-xtitle('interpolation of square function')
-legend(['linear','spline','nearest'], "in_upper_left");
+x = linspace(0, 10, 11)';
+y = sin(x);
+xx = linspace(0,10,1000)';
+yy2 = interp1(x, y, xx, 'linear');
+yy1 = interp1(x, y, xx, 'nearest');
+yy3 = interp1(x, y, xx, 'spline');
+clf
+h = plot(xx, [yy1 yy2 yy3], x, y, '.')
+h(1).mark_size = 8;
+title "Interpolation of a poorly sampled sin() function" fontsize 3
+legend(['nearest','linear','spline','nodes'], "in_lower_left");
]]></programlisting>
<scilab:image>
- x=linspace(0,3,20);
- y=x.^2;
- xx=linspace(0,3,100);
- yy1=interp1(x,y,xx,'linear');
- yy2=interp1(x,y,xx,'spline');
- yy3=interp1(x,y,xx,'nearest');
- plot(xx',[yy1;yy2;yy3]',x,y,'*')
- xtitle('interpolation of square function')
- legend(['linear','spline','nearest'],a=2)
+ x = linspace(0, 10, 11)';
+ y = sin(x);
+ xx = linspace(0,10,1000)';
+ yy2 = interp1(x, y, xx, 'linear');
+ yy1 = interp1(x, y, xx, 'nearest');
+ yy3 = interp1(x, y, xx, 'spline');
+ clf
+ h = plot(xx, [yy1 yy2 yy3], x, y, '.')
+ h(1).mark_size = 8;
+ title "Interpolation of a poorly sampled sin() function" fontsize 3
+ legend(['nearest','linear','spline','nodes'], "in_lower_left");
+ gcf().axes_size = [600 400];
</scilab:image>
</refsection>
<refsection role="see also">
<link linkend="interp">interp</link>
</member>
<member>
- <link linkend="interpln">interpln</link>
+ <link linkend="splin">splin</link>
</member>
<member>
- <link linkend="splin">splin</link>
+ <link linkend="linear_interpn">linear_interpn</link>
</member>
</simplelist>
</refsection>
<title>History</title>
<revhistory>
<revision>
- <revnumber>5.4.0</revnumber>
- <revremark>previously, imaginary part of input arguments were implicitly ignored.</revremark>
+ <revnumber>6.1.1</revnumber>
+ <revremark>
+ <itemizedlist>
+ <listitem>
+ For complex <varname>y</varname> values, <literal>imag(y)</literal> is
+ no longer ignored: <literal>real(y)</literal> and <literal>imag(y)</literal>
+ parts are separately interpolated.
+ </listitem>
+ <listitem>
+ <literal>"periodic"</literal> extrapolation added for the linear and
+ spline interpolations.
+ </listitem>
+ <listitem>
+ <literal>"edgevalue"</literal> extrapolation added for all nearest,
+ linear and spline interpolations.
+ </listitem>
+ <listitem>
+ <literal>"linear"</literal> extrapolation added for the spline
+ interpolation.
+ </listitem>
+ <listitem>
+ When <varname>xp</varname> is an hypermatrix and
+ <literal>size(xp,1)==1</literal>, <literal>size(yp)</literal>
+ is now always <literal>[size(xp) size(y)(2,$)</literal> instead of
+ <literal>[size(xp,2:$), size(y)(2,$)</literal>.
+ </listitem>
+ </itemizedlist>
+ </revremark>
</revision>
</revhistory>
</refsection>
+++ /dev/null
-<?xml version="1.0" encoding="UTF-8"?>
-<!--
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2008 - INRIA - Farid BELAHCENE
- *
- * Copyright (C) 2012 - 2016 - Scilab Enterprises
- *
- * This file is hereby licensed under the terms of the GNU GPL v2.0,
- * pursuant to article 5.3.4 of the CeCILL v.2.1.
- * This file was originally licensed under the terms of the CeCILL v2.1,
- * and continues to be available under such terms.
- * For more information, see the COPYING file which you should have received
- * along with this program.
- *
- -->
-<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink"
- xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns5="http://www.w3.org/1999/xhtml"
- xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook"
- xmlns:scilab="http://www.scilab.org" xml:id="interp1" xml:lang="ja">
- <refnamediv>
- <refname>interp1</refname>
- <refpurpose>一次元補間関数</refpurpose>
- </refnamediv>
- <refsynopsisdiv>
- <title>呼び出し手順</title>
- <synopsis>
- yp = interp1(x, y, xp)
- yp = interp1(x, y, xp, method)
- yp = interp1(x, y, xp, method, extrapolation)
- </synopsis>
- </refsynopsisdiv>
- <refsection>
- <title>パラメータ</title>
- <variablelist>
- <varlistentry>
- <term>xp</term>
- <listitem>
- <para>実数スカラー, ベクトルまたは行列 (またはハイパー行列)</para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>x</term>
- <listitem>
- <para>実数ベクトル</para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>method</term>
- <listitem>
- <para>(オプションの) 補間方法を定義する文字列</para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>extrapolation</term>
- <listitem>
- <para>(オプションの) 文字列, または領域 [x1,xn] の外側の xp(j)の値について
- yp(j) の要素を定義する実数値.
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>yp</term>
- <listitem>
- <para>ベクトルまたは行列 (またはハイパー行列)</para>
- </listitem>
- </varlistentry>
- </variablelist>
- </refsection>
- <refsection>
- <title>説明</title>
- <para>
- <literal>(x,y,xp)</literal>を指定すると, この関数は
- xp に対応する yp の要素を x および y により定義された
- (デフォルトでは線形)補間により計算します.
- </para>
- <para>
- <literal>yp</literal> がベクトルの場合,
- xp の長さは <literal>yp</literal> の長さに等しくなります.
- <literal>yp</literal> が行列の場合,
- <literal>xp</literal> は yp の各列の長さと同じ長さとなります.
- <literal>yp</literal> がハイパー行列の場合,
- <literal>xp</literal> の長さは<literal>yp</literal> の最初の次元と同じ
- になります.
- </para>
- <para>size(y)=[C,M1,M2,M3,...,Mj] かつ size(xp)=[N1,N2,N3,...,Nk] の場合,
- size(yp)=[N1,N2,..,Nk,M1,M2,...Mj] となり, x の長さは size(y,1) に等しくなります.
- </para>
- <para>
- <literal>method</literal> パラメータは補間の計算手法を設定します.
- </para>
- <variablelist>
- <varlistentry>
- <term>"linear"</term>
- <listitem>
- <para>
- 補間は線形手法により定義されます ( <link linkend="interpln">interpln</link>参照)
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>"spline"</term>
- <listitem>
- <para>補間は3次スプライン補間により定義されます (
- <link linkend="splin">splin</link> , <link linkend="interp">interp</link>参照)
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>"nearest"</term>
- <listitem>
- <para>各 xp(j), yp(j) の値は,
- xp(j)の最近傍にある x(i) に対応する y(i) の値を選択します.
- </para>
- </listitem>
- </varlistentry>
- </variablelist>
- <para>
- <literal>extrapolation</literal> パラメータは捕外用,すなわち
- <literal>xp(i) </literal> が領域[x1,xn]内にない時
- の計算手法を設定します.
- </para>
- <variablelist>
- <varlistentry>
- <term>"extrap"</term>
- <listitem>
- <para>捕外が定義された手法により実行されます.
- yp=interp1(x,y,xp,method,"extrap")
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>実数値</term>
- <listitem>
- <para>捕外用の実数値を選択できます.
- この場合, 区間 [x1,xn] の中にない xp(i) について yp(i) は
- この値とします.
- yi=interp1(x,y,xp,method, 0)
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>デフォルトの動作</term>
- <listitem>
- <para>捕外は(スプライン手法の場合)定義された手法により行われます
- 線形または最近傍手法の場合は nan となります.
- yp=interp1(x,y,xp,method)
- </para>
- </listitem>
- </varlistentry>
- </variablelist>
- </refsection>
- <refsection>
- <title>例</title>
- <programlisting role="example"><![CDATA[
-x = linspace(0, 3, 20);
-y = x^2;
-xx = linspace(0, 3, 100);
-yy1 = interp1(x, y, xx, 'linear');
-yy2 = interp1(x, y, xx, 'spline');
-yy3 = interp1(x, y, xx, 'nearest');
-plot(xx, [yy1; yy2; yy3], x, y, '*')
-xtitle('interpolation of square function')
-legend(['linear','spline','nearest'], "in_upper_left");
- ]]></programlisting>
- <scilab:image>
- x=linspace(0,3,20);
- y=x.^2;
- xx=linspace(0,3,100);
- yy1=interp1(x,y,xx,'linear');
- yy2=interp1(x,y,xx,'spline');
- yy3=interp1(x,y,xx,'nearest');
- plot(xx',[yy1;yy2;yy3]',x,y,'*')
- xtitle('interpolation of square function')
- legend(['linear','spline','nearest'],a=2)
- </scilab:image>
- </refsection>
- <refsection role="see also">
- <title>参照</title>
- <simplelist type="inline">
- <member>
- <link linkend="interp">interp</link>
- </member>
- <member>
- <link linkend="interpln">interpln</link>
- </member>
- <member>
- <link linkend="splin">splin</link>
- </member>
- </simplelist>
- </refsection>
- <refsection>
- <title>履歴</title>
- <revhistory>
- <revision>
- <revnumber>5.4.0</revnumber>
- <revremark>以前では, 入力引数の虚部は暗黙のうちに無視されていました.</revremark>
- </revision>
- </revhistory>
- </refsection>
-</refentry>
+++ /dev/null
-<?xml version="1.0" encoding="UTF-8"?>
-<!--
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2008 - INRIA - Farid BELAHCENE
- *
- * Copyright (C) 2012 - 2016 - Scilab Enterprises
- *
- * This file is hereby licensed under the terms of the GNU GPL v2.0,
- * pursuant to article 5.3.4 of the CeCILL v.2.1.
- * This file was originally licensed under the terms of the CeCILL v2.1,
- * and continues to be available under such terms.
- * For more information, see the COPYING file which you should have received
- * along with this program.
- *
- -->
-<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink"
- xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns5="http://www.w3.org/1999/xhtml"
- xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook"
- xmlns:scilab="http://www.scilab.org" xml:id="interp1" xml:lang="pt">
- <refnamediv>
- <refname>interp1</refname>
- <refpurpose>função de interpolação 1d</refpurpose>
- </refnamediv>
- <refsynopsisdiv>
- <title>Seqüência de Chamamento</title>
- <synopsis>
- yp = interp1(x, y, xp)
- yp = interp1(x, y, xp, method)
- yp = interp1(x, y, xp, method, extrapolation)
- </synopsis>
- </refsynopsisdiv>
- <refsection>
- <title>Parâmetros</title>
- <variablelist>
- <varlistentry>
- <term>xp</term>
- <listitem>
- <para>escalar real, vetor ou matriz (ou hipermatriz) de reais
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>x</term>
- <listitem>
- <para>vetor de reais</para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>method</term>
- <listitem>
- <para>(opcional) string definindo o método de interpolação </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>extrapolation</term>
- <listitem>
- <para>(opcional) string ou valor real defindo os componentes yp(j)
- para os valores xp(j) fora do intervalo [x1,xn].
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>yp</term>
- <listitem>
- <para>vetor ou matriz (ou hipermatriz) </para>
- </listitem>
- </varlistentry>
- </variablelist>
- </refsection>
- <refsection>
- <title>Descrição</title>
- <para>
- Dados <literal>(x,y,xp)</literal>, esta função faz corresponder os
- yp componentes aos xp por interpolação (linear por padrão) definida por x
- e y.
- </para>
- <para>
- Se <literal>yp</literal> é um vetor, então o comprimento de xp é
- igual ao comprimento de <literal>yp</literal>. Se <literal>yp</literal> é
- uma matriz, então <literal>xp</literal> tem o mesmo comprimento que cada
- uma das colunas de yp. Se <literal>yp</literal> é uma hipermatriz, então o
- comprimento de <literal>xp</literal> é o mesmo da primeira dimensão de
- <literal>yp. </literal>
- </para>
- <para>Se size(y)=[C,M1,M2,M3,...,Mj] e size(xp)= [N1,N2,N3,...,Nk] então
- size(yp)= [N1,N2,..,Nk,M1,M2,...Mj] e o comprimento de x deve ser igual a
- size(y,1)
- </para>
- <para>
- O parâmetro <literal>method</literal> ajusta a regra avaliação para
- interpolação.
- </para>
- <variablelist>
- <varlistentry>
- <term>"linear"</term>
- <listitem>
- <para>
- a interpolação é definida pelo método linear (ver <link linkend="interpln">interpln</link>)
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>"spline"</term>
- <listitem>
- <para>
- definição de interpolação por spline cúbico (ver <link linkend="splin">splin</link> , <link linkend="interp">interp</link>)
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>"nearest"</term>
- <listitem>
- <para>para cada valor xp(j), yp(j) toma o valor ou y(i)
- correspondente ao x(i), o vízinho mais próximo de xp(j)
- </para>
- <para/>
- </listitem>
- </varlistentry>
- </variablelist>
- <para>
- O parâmetro <literal>extrapolation</literal> ajusta a regra de
- avaliação para extrapolação, i.e para <literal>xp(i) </literal>fora do
- intervalo [x1,xn]
- </para>
- <variablelist>
- <varlistentry>
- <term>"extrap"</term>
- <listitem>
- <para>a extrapolação é realizada pelo método definido.
- yp=interp1(x,y,xp,method,"extrap")
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>valor real</term>
- <listitem>
- <para>você pode escolher um valor real para extrapolação. Deste
- modo, yp(i) toma este valor para xp(i) fora do intervalo [x1,xn],
- por exemplo 0 (mas também nan ou inf). yi=interp1(x,y,xp,method, 0)
- </para>
- </listitem>
- </varlistentry>
- <varlistentry>
- <term>por padrão</term>
- <listitem>
- <para>a extrapolação é realizada pelo método definido (para método
- spline), e por nan para os métodos "linear" e "nearest".
- yp=interp1(x,y,xp,method)
- </para>
- </listitem>
- </varlistentry>
- </variablelist>
- </refsection>
- <refsection>
- <title>Exemplos</title>
- <programlisting role="example"><![CDATA[
-x = linspace(0, 3, 20);
-y = x^2;
-xx = linspace(0, 3, 100);
-yy1 = interp1(x, y, xx, 'linear');
-yy2 = interp1(x, y, xx, 'spline');
-yy3 = interp1(x, y, xx, 'nearest');
-plot(xx, [yy1; yy2; yy3], x, y, '*')
-xtitle('interpolation of square function')
-legend(['linear','spline','nearest'], "in_upper_left");
- ]]></programlisting>
- </refsection>
- <refsection role="see also">
- <title>Ver Também</title>
- <simplelist type="inline">
- <member>
- <link linkend="interp">interp</link>
- </member>
- <member>
- <link linkend="interpln">interpln</link>
- </member>
- <member>
- <link linkend="splin">splin</link>
- </member>
- </simplelist>
- </refsection>
-</refentry>
--- /dev/null
+<?xml version="1.0" encoding="UTF-8"?>
+<!--
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) 2008 - INRIA - Farid BELAHCENE
+ * Copyright (C) 2012 - 2016 - Scilab Enterprises
+ * Copyright (C) 2021 - Samuel GOUGEON - Le Mans Université
+ *
+ * This file is hereby licensed under the terms of the GNU GPL v2.0,
+ * pursuant to article 5.3.4 of the CeCILL v.2.1.
+ * This file was originally licensed under the terms of the CeCILL v2.1,
+ * and continues to be available under such terms.
+ * For more information, see the COPYING file which you should have received
+ * along with this program.
+ *
+ -->
+<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink"
+ xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns5="http://www.w3.org/1999/xhtml"
+ xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook"
+ xmlns:scilab="http://www.scilab.org" xml:id="interp1" xml:lang="ru">
+ <refnamediv>
+ <refname>interp1</refname>
+ <refpurpose>одномерная интерполяция в режиме nearest, linear или spline</refpurpose>
+ </refnamediv>
+ <refsynopsisdiv>
+ <title>Синтаксис</title>
+ <synopsis>
+ yp = interp1(y, xp)
+ yp = interp1(x, y, xp)
+ yp = interp1(.., xp, method)
+ yp = interp1(.., xp, method, extrapolation)
+ </synopsis>
+ </refsynopsisdiv>
+ <refsection>
+ <title>Аргументы</title>
+ <variablelist>
+ <varlistentry>
+ <term>x</term>
+ <listitem>
+ вектор по меньшей мере из двух вещественных чисел: ось абсцисс известных узлов
+ интерполяции, без повторов. По умолчанию
+ <itemizedlist>
+ <listitem>
+ если <varname>y</varname> является вектором: <literal>x=1:length(y)</literal>.
+ </listitem>
+ <listitem>
+ если <varname>y</varname> является матрицей или гиперматрицей:
+ <literal>x=1:size(y,1)</literal>.
+ </listitem>
+ </itemizedlist>
+ <para/>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>y</term>
+ <listitem>
+ вектор, матрица или гиперматрица вещественных или комплексных
+ чисел: значения в узлах интерполяции в соответствующих значениях
+ оси абсцисс <varname>x</varname>.
+ <itemizedlist>
+ <listitem>
+ если <varname>y</varname> является вектором, то <varname>x</varname>
+ и <varname>y</varname> должны быть одной длины.
+ </listitem>
+ <listitem>
+ если <varname>y</varname> является матрицей или гиперматрицей, то мы
+ должны иметь <literal>length(x)==size(y,1)</literal>. Каждый столбец
+ <varname>y</varname> тогда интерполирован в зависимости от тех же
+ значений оси абсцисс <varname>x</varname>, для указанного
+ <varname>xp</varname>.
+ </listitem>
+ </itemizedlist>
+ <para/>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>xp</term>
+ <listitem>
+ скаляр, вектор, матрица или гиперматрица десятичных чисел:
+ ось абсцисс точек, чьи значения <varname>yp</varname> должны быть
+ вычислены в соответствии с данными узлов интерполяции.
+ <para/>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>yp</term>
+ <listitem>
+ вектор, матрица или гиперматрица чисел: интерполированные значения
+ <varname>y</varname> в указанном <varname>xp</varname>.
+ <itemizedlist>
+ <listitem>
+ если <varname>y</varname> является вектором:
+ <varname>yp</varname> имеет размер <varname>xp</varname>.
+ <para/>
+ </listitem>
+ <listitem>
+ если <varname>y</varname> является матрицей или гиперматрицей:
+ <itemizedlist>
+ <listitem>
+ если <varname>xp</varname> является скаляром или вектором:
+ <literal>size(yp)</literal> равен
+ <literal>[length(xp) size(y)(2:$)]</literal>
+ </listitem>
+ <listitem>
+ если <varname>xp</varname> является матрицей или
+ гиперматрицей: <literal>size(yp)</literal> равен
+ <literal>[size(xp) size(y)(2:$)]</literal>
+ </listitem>
+ </itemizedlist>
+ </listitem>
+ </itemizedlist>
+ <para/>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>method</term>
+ <listitem>
+ строковое значение, определяющее метод интерполяции. Возможными
+ значениями и обработкой являются:
+ <para/>
+ <table>
+ <tr>
+ <td><emphasis role="bold">"linear"</emphasis>:</td>
+ <td>
+ линейная интерполяция между последовательными узлами, используется
+ по умолчанию.
+ </td>
+ </tr>
+ <tr>
+ <td><emphasis role="bold">"spline"</emphasis>:</td>
+ <td>
+ интерполяция кубическими сплайнами
+ </td>
+ </tr>
+ <tr>
+ <td><emphasis role="bold">"nearest"</emphasis>:</td>
+ <td>
+ <para>
+ для каждого значения <literal>xp(j)</literal> <literal>yp(j)</literal>
+ принимает значение <literal>y(i)</literal>, соответствующее
+ <literal>x(i)</literal>, ближайшего соседа
+ <literal>xp(j)</literal>
+ </para>
+ </td>
+ </tr>
+ </table>
+ <para/>
+ </listitem>
+ </varlistentry>
+ <varlistentry>
+ <term>extrapolation</term>
+ <listitem>
+ строковое значение или число, определяющее элементы <literal>yp(j)</literal>
+ для значений <literal>xp(j)</literal> за пределами интервала
+ <literal>[x(1)=min(x),x($)=max(x)]</literal>. Мы полагаем здесь и далее,
+ что <varname>x</varname> и <varname>y</varname> уже соответственно
+ отсортированы.
+ <para/>
+ <table>
+ <tr>
+ <td><emphasis role="bold">"extrap"</emphasis>:</td>
+ <td>
+ <literal>interp1(x,y,xp, method, "extrap")</literal> эквивалентен
+ <literal>interp1(x,y,xp, method, method)</literal>.
+ </td>
+ </tr>
+ <tr>
+ <td><emphasis role="bold">"linear"</emphasis>:</td>
+ <td>
+ Может использоваться с методами интерполяции
+ <literal>"spline"</literal> (и, очевидно,
+ <literal>"linear"</literal>).
+ </td>
+ </tr>
+ <tr>
+ <td><emphasis role="bold">"periodic"</emphasis>:</td>
+ <td>
+ Этот тип экстраполяции может использоваться с методами
+ интерполяции <literal>"linear"</literal> или
+ <literal>"spline"</literal>. Тогда: если <varname>y</varname>
+ является вектором, то требуется, чтобы
+ <literal>y(1)==y($)</literal>; в противном случае требуется,
+ чтобы <literal>y(1,:)==y($,:)</literal>.
+ </td>
+ </tr>
+ <tr>
+ <td><emphasis role="bold">"edgevalue"</emphasis>:</td>
+ <td>
+ Тогда <literal>yp(i)=y(1)</literal> для каждого
+ <literal>xp(i)<x(1)</literal>, и
+ <literal>yp(i)=y($)</literal> для каждого
+ <literal>xp(i)>x($)</literal>.
+ </td>
+ </tr>
+ <tr>
+ <td><emphasis role="bold">padding</emphasis>:</td>
+ <td>
+ <varname>padding</varname> является десятичным или
+ комплексным числом, используемым для установки
+ <literal>yp(i)=padding</literal> для каждого
+ <literal>xp(i) ∉ [min(x),max(x)]</literal>. Например:
+ <literal>yi=interp1(x,y,xp,method, 0)</literal>.
+ </td>
+ </tr>
+ <tr>
+ <td><emphasis role="bold">(none)</emphasis>:</td>
+ <td>
+ По умолчанию экстраполяция выполняется сплайнами, когда
+ сплайны используются для интерполяции и дополнением
+ значениями <literal>%nan</literal>, если интерполяция
+ линейна или по "ближайшему" узлу.
+ </td>
+ </tr>
+ </table>
+ </listitem>
+ </varlistentry>
+ </variablelist>
+ </refsection>
+ <refsection>
+ <title>Описание</title>
+ <para>
+ Указывая <literal>(x,y,xp)</literal>, данная функция вычисляет элементы
+ <varname>yp</varname>, соответствующие <varname>xp</varname> с помощью интерполяции
+ между известными данными, указанными в узлах <literal>(x,y)</literal>.
+ </para>
+ <para>
+ <varname>x</varname> предварительно отсортирован в порядке возрастания, а значения
+ <varname>y</varname> либо по столбцам тогда сортируются соответственно.
+ </para>
+ <para>
+ <emphasis role="bold">Интерполяция комплексных значений</emphasis>:
+ Если <varname>y</varname> является комплексным, то его вещественная и мнимая
+ части интерполируются отдельно, и потом суммируются для формирования
+ комплексного <varname>yp</varname>.
+ </para>
+ <para>
+ <emphasis role="bold">interp1(x,y,xp,"nearest")</emphasis>:
+ Для любого <varname>xp</varname> в середине интервала
+ <literal>[x(i),x(i+1)]</literal> верхняя граница <literal>x(i+1)</literal>
+ рассматривается как ближайшее значение <varname>x</varname>, и присваивается
+ <literal>yp=y(i+1)</literal>.
+ </para>
+ <refsect3>
+ <title>линейные интерполяции</title>
+ Они выполняются через функцию <literal>linear_interpn(..)</literal>,
+ с соответствующей опцией интерполяции <literal>"edgevalue"→"C0"</literal>,
+ <literal>"linear"→"natural"</literal>, <literal>"periodic"→"periodic"</literal>.
+ </refsect3>
+ <refsect3>
+ <title>интерполяции сплайнами</title>
+ <para>
+ <emphasis role="bold">interp1(..,xp,"spline")</emphasis> или
+ <emphasis role="bold">interp1(..,xp,"spline","spline")</emphasis> или
+ <emphasis role="bold">interp1(..,xp,"spline","extrap")</emphasis>
+ используют условия граней <literal>не_узел</literal>.
+ Экстраполяция выполняется с помощью обоих сплайновых полиномов,
+ вычисленных на гранях <literal>(x,y)</literal>.
+ </para>
+ <para>
+ <emphasis role="bold">interp1(..,xp,"spline","edgevalue")</emphasis>
+ использует условия граней <literal>не_узел</literal>, а затем вызывает
+ <literal>interp(..,"C0")</literal>, чтобы выполнить фактическую
+ интерполяцию и экстраполяцию.
+ </para>
+ <para>
+ <emphasis role="bold">interp1(..,xp,"spline","periodic")</emphasis>
+ вызывает оба <literal>splin(..)</literal>, а затем <literal>interp(..)</literal>
+ с их опцией <literal>"periodic"</literal>.
+ </para>
+ <para>
+ <emphasis role="bold">interp1(..,xp,"spline","linear")</emphasis>
+ вызывает <literal>splin(..,"natural")</literal> для условий линейных граней, а затем
+ передаётся в <literal>interp(..,"linear")</literal>.
+ </para>
+ </refsect3>
+ </refsection>
+ <refsection>
+ <title>Примеры</title>
+ <programlisting role="example"><![CDATA[
+x = linspace(0, 10, 11)';
+y = sin(x);
+xx = linspace(0,10,1000)';
+yy2 = interp1(x, y, xx, 'linear');
+yy1 = interp1(x, y, xx, 'nearest');
+yy3 = interp1(x, y, xx, 'spline');
+clf
+h = plot(xx, [yy1 yy2 yy3], x, y, '.')
+h(1).mark_size = 8;
+title "Interpolation of a poorly sampled sin() function" fontsize 3
+legend(['nearest','linear','spline','nodes'], "in_lower_left");
+ ]]></programlisting>
+ <scilab:image>
+ x = linspace(0, 10, 11)';
+ y = sin(x);
+ xx = linspace(0,10,1000)';
+ yy2 = interp1(x, y, xx, 'linear');
+ yy1 = interp1(x, y, xx, 'nearest');
+ yy3 = interp1(x, y, xx, 'spline');
+ clf
+ h=plot(xx, [yy1 yy2 yy3],x, y, '.')
+ h(1).mark_size = 8;
+ title "Interpolation of a poorly sampled sin() function" fontsize 3
+ legend(['nearest','linear','spline','nodes'], "in_lower_left");
+ gcf().axes_size = [600 400];
+ </scilab:image>
+ </refsection>
+ <refsection role="see also">
+ <title>Смотрите также</title>
+ <simplelist type="inline">
+ <member>
+ <link linkend="interp">interp</link>
+ </member>
+ <member>
+ <link linkend="splin">splin</link>
+ </member>
+ <member>
+ <link linkend="linear_interpn">linear_interpn</link>
+ </member>
+ </simplelist>
+ </refsection>
+ <refsection>
+ <title>История</title>
+ <revhistory>
+ <revision>
+ <revnumber>6.1.1</revnumber>
+ <revremark>
+ <itemizedlist>
+ <listitem>
+ Для комплексных значений <varname>y</varname>,
+ <literal>imag(y)</literal> более не игнорируется:
+ <literal>real(y)</literal> и <literal>imag(y)</literal>
+ части интерполируются раздельно.
+ </listitem>
+ <listitem>
+ Экстраполяция <literal>"periodic"</literal> добавлена для
+ линейной и сплайновой интерполяций.
+ </listitem>
+ <listitem>
+ Экстраполяция <literal>"edgevalue"</literal> добавлена для
+ всех интерполяций ближайшего соседа, линейной и сплайновой.
+ </listitem>
+ <listitem>
+ Экстраполяция <literal>"linear"</literal> добавлена для
+ сплайновой интерполяции.
+ </listitem>
+ <listitem>
+ Когда <varname>xp</varname> является гиперматрицей и
+ <literal>size(xp,1)==1</literal>, <literal>size(yp)</literal>
+ теперь всегда <literal>[size(xp) size(y)(2,$)</literal> вместо
+ <literal>[size(xp,2:$), size(y)(2,$)</literal>.
+ </listitem>
+ </itemizedlist>
+ </revremark>
+ </revision>
+ </revhistory>
+ </refsection>
+</refentry>
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA - Farid BELAHCENE
-//
// Copyright (C) 2012 - 2016 - Scilab Enterprises
+// Copyright (C) 2021 - Samuel GOUGEON - Le Mans Université
//
// This file is hereby licensed under the terms of the GNU GPL v2.0,
// pursuant to article 5.3.4 of the CeCILL v.2.1.
// For more information, see the COPYING file which you should have received
// along with this program.
-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
+function yi = interp1(varargin)
+ // yi = interp1(x, y, xi)
+ // yi = interp1(x, y, xi, method)
+ // yi = interp1(x, y, xi, method, extrapolation)
+ // yi = interp1(y, xi)
+ // yi = interp1(y, xi, method)
+ // yi = interp1(y, xi, method, extrapolation)
//
- // 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
+ // x: real vector
+ // y: real or complex vector of length(x) or array of size(y,1)==length(x)
+ // xi: real scalar, vector, matrix or hypermatrix of any size.
+ // method = "linear" (default) | "nearest" | "spline"
+ // extrapolation = "extrap" (= method) | real or complex padding scalar |
+ // * if method = "linear": "linear" | "periodic" | "edgevalue"
+ // * if method = "spline": "linear" | "periodic" | "edgevalue" | "spline"
+ // * if method = "nearest": "nearest" = "edgevalue"
- rhs=size(varargin)
+ // ========================
+ // CHECKING INPUT ARGUMENTS
+ // ========================
+ 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));
+ msg = _("%s: Wrong number of input arguments: Must be between %d and %d.\n")
+ error(msprintf(msg, "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
+ // Checking x
+ // ----------
+ iy = 2 // index of y argin, for error messages
+ if rhs >= 3 & type(varargin(1))<>0 & type(varargin(3))<>10 then
+ // x is explicit: interp1(x, y, xi,..)
+ x = varargin(1)
+ if type(x)<>1 | ~isreal(x) then
+ msg = _("%s: Argument #%d: Real numbers expected.\n")
+ error(msprintf(msg,"interp1",1))
+ elseif ~isvector(x)
+ msg = _("%s: Argument #%d: Vector expected.\n")
+ error(msprintf(msg, "interp1", 1))
+ elseif or(isnan(x)) then
+ msg = _("%s: Argument #%d: Nan value forbidden.\n")
+ error(msprintf(msg, "interp1", 1))
+ end
+ x = matrix(x,1,-1)
+ y = varargin(2)
+ else // Default x:
+ if type(varargin(1))<>0 // interp1(y, xi,..)
+ iy = 1
+ y = varargin(1)
+ else // interp1(, y, xi,..)
+ y = varargin(2)
+ end
+ if isvector(y) then
+ x = 1:size(y,"*")
+ y = y(:)
else
- error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector or matrix expected.\n"),"interp1",2));
+ x = 1:size(y,1)
end
- for i=rhs:-1:1
- varargin(i+1)=varargin(i)
+ end
+ if length(x) <= 1
+ if iy==2
+ msg = _("%s: Argument #%d: At least 2 elements expected.\n")
+ error(msprintf(msg, "interp1", 1))
+ else
+ msg = _("%s: Argument #%d: vector or at least 2 rows expected.\n")
+ error(msprintf(msg, "interp1", 2))
end
- varargin(1)=X
end
+ // Sorting x in increasing order:
+ [?, k] = gsort(x,"g","i")
+ x = x(k) // should anyway be a row
- //............................
- //ininialisation of xi
- //............................
- //xi components must be reals
- xi=varargin(3)
- if type(xi)<>1 then
- error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of reals expected.\n"),"interp1",3));
- 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);
+ // Checking y
+ // ----------
+ if type(y) <> 1 then
+ msg = _("%s: Argument #%d: Decimal or complex numbers expected.\n")
+ error(msprintf(msg, "interp1", iy));
end
- xisizetemp=xisize
- if size(xisize,"*")>=2 then
- if xisize(1)==1 then
- xisize=xisize(2:$);
+ if isvector(y) then
+ y = y(:)
+ if length(y) <> length(x)
+ msg = _("%s: Arguments #%d and #%d: Same numbers of elements expected.\n")
+ error(msprintf(msg, "interp1", 1, 2))
+ end
+ else
+ if size(y,1) <> length(x) then
+ msg = _("%s: Arguments #%d and #%d: Same numbers of rows expected.\n")
+ error(msprintf(msg, "interp1", 1, 2))
end
end
+ // Sorting y according to x:
+ cmd = "y = y(k," + strcat(repmat(":",[1 ndims(y)-1]),",")+")"
+ execstr(cmd)
- //-------------------------
- //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));
+ // Checking xi
+ // -----------
+ if type(varargin(iy+1))<>1 | ~isreal(varargin(iy+1)) then
+ msg = _("%s: Argument #%d: Decimal numbers expected.\n")
+ error(msprintf(msg, "interp1", iy+1))
end
- //y components must be reals
- if type(y)<>1 then
- error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of reals expected.\n"),"interp1",2));
+ xi = varargin(iy+1)
+ if xi == [] then
+ yi = []
+ return
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));
+ xisize = size(xi)
+ if length(size(y))==2 & or(size(y)==1) then // y = 2D vector
+ if iscolumn(xisize) then
+ xisize = length(xi)
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)
+ if length(size(xi))==2 & or(size(xi)==1) then // xi = scalar or 2D vector
+ xisize = length(xi)
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"));
+ // Checking method
+ // ---------------
+ if rhs > iy+1 & type(varargin(iy+2))<>0 then
+ method = varargin(iy+2)
+ if size(method,"*") <> 1
+ msg = _("%s: Argument #%d: Scalar (1 element) expected.\n")
+ error(msprintf(msg, "interp1", iy+2))
+ end
+ if type(method) <> 10 | ~or(convstr(method)==["linear","nearest","spline"])
+ msg = _("%s: Argument #%d: Must be in the set {%s}.\n")
+ error(msprintf(msg, "interp1", iy+2, """linear"",""spline"",""nearest"""))
end
+ method = convstr(method)
+ else
+ method = "linear"
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 = resize_matrix(0, 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");
+ // Checking extrapolation method
+ // -----------------------------
+ extraVal = []
+ if rhs > iy+2 then
+ extra = varargin(iy+3)
+ err = %f
+ if type(extra)==10 & size(extra,"*")==1
+ extra = convstr(extra)
+ if extra=="extrap"
+ extra = "natural"
+ elseif or(method==["linear" "spline"])
+ if extra=="linear" then
+ if method=="linear"
+ extra = "natural"
+ end
+ elseif extra=="edgevalue"
+ extra = "C0"
+ elseif extra=="periodic"
+ msg = []
+ if isvector(y) & y($)<>y(1)
+ msg = _("%s: Argument #%d: periodicity y($)==y(1) expected.\n")
+ elseif ~isvector(y) & or(y($,:)<>y(1,:))
+ msg = _("%s: Argument #%d: periodicity y($,:)==y(1,:) expected.\n")
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)
+ if msg <> []
+ error(msprintf(msg, "interp1",iy))
end
+ elseif extra<>"natural"
+ err = %t
end
- if size(xisize,"*")>=2
- yi=matrix(yi,xisize)
+ else // nearest
+ if or(extra==["edgevalue" "nearest" "natural"])
+ extra = "natural"
else
- yi=matrix(yi,xisizetemp)
- end
-
- // y is matrix or hypermatrix
- elseif size(size(y),"*")>=2 then
- ysize=size(y)
- ky=ysize(2:$)
- yi = resize_matrix(0, [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
+ err = %t
end
- yi=matrix(yi,[xisize,ky])
+ end
+ else
+ if type(extra)<>1 | length(extra)<>1
+ err = %t
else
- error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector or matrix expected.\n"),"interp1",2));
+ extraVal = extra
+ extra = "natural"
end
+ end
+ if err
+ if method=="nearest"
+ tmp = """extrap"",""edgevalue"",""nearest"""
+ else
+ tmp = """extrap"",""linear"",""periodic"",""edgevalue"""
+ end
+ msg = _("%s: Argument #%d: scalar number or one of {%s} expected.\n")
+ error(msprintf(msg, "interp1", iy+3, tmp))
+ end
+ else
+ extra = "natural"
+ if method <> "spline"
+ extraVal = %nan
+ end
+ 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"));
+ // ==========
+ // PROCESSING
+ // ==========
+ if ~isreal(y) then
+ if extraVal==[]
+ [extraR, extraI] = (extra, extra);
+ else
+ [extraR, extraI] = (real(extraVal), imag(extraVal));
+ end
+ yir = interp1(x, real(y), xi, method, extraR)
+ yii = interp1(x, imag(y), xi, method, extraI)
+ yi = complex(yir, yii)
+ return
+ end
+ // REAL CASE
+ // =========
+ select part(method, 1)
+ case "l"
+ // LINEAR
+ // ------
+ xitemp = matrix(xi,-1,1)
+ // y is a vector
+ if isvector(y) then
+ yi = resize_matrix(0, size(xitemp))
+ if extraVal <> [] then
+ yi = linear_interpn(xitemp, x, y, "by_nan");
+ if ~isnan(extraVal) then
+ k = find(xitemp<min(x) | xitemp>max(x))
+ yi(k) = extraVal
end
+ else
+ yi = linear_interpn(xitemp, x, y, extra);
end
- xitemp=matrix(xi,-1,1)
- //y is a vector
- if isvector(y) then
- yi = resize_matrix(0, 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)
+ yi = matrix(yi, xisize)
+
+ else // y is matrix or hypermatrix
+ ky = size(y)(2:$)
+ ncolumns = prod(ky)
+ yi = resize_matrix(0, [size(xitemp) ky])
+ if extraVal <> [] then
+ for ic = 1:ncolumns
+ tmp = linear_interpn(xitemp, x, y(:,ic), "by_nan")
+ yi(:,ic) = matrix(tmp, size(xitemp))
end
- if size(xisize,"*")>=2
- yi=matrix(yi,xisize)
- else
- yi=matrix(yi,xisizetemp)
+ if ~isnan(extraVal) then
+ k = find(xitemp < min(x) | xitemp > max(x))
+ yi(k,:) = extraVal
end
- //y is a matrix or a hypermatrix
- elseif size(size(y),"*")>=2 then
- ky=size(y)
- ky=ky(2:$)
- yi = resize_matrix(0, [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
+ else // extra=="linear"=="natural" | "periodic" | "edgevalue"=="C0"
+ for ic = 1:ncolumns
+ tmp = linear_interpn(xitemp, x, y(:,ic), extra)
+ yi(:,ic) = matrix(tmp, size(xitemp))
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])
+ end
+
+ case "s"
+ // SPLINE
+ // ------
+ xitemp = matrix(xi,-1,1)
+ // y is a vector
+ if isvector(y) then
+ yi = resize_matrix(0, size(xitemp))
+ yrow = matrix(y,1,-1)
+ if extra=="linear"
+ yi = interp(xitemp, x, yrow, splin(x, yrow, "natural"), extra)
+ elseif extra=="periodic"
+ yi = interp(xitemp, x, yrow, splin(x, yrow, extra), extra)
+ else // default not_a_knot
+ yi = interp(xitemp, x, yrow, splin(x, yrow), extra)
+ end
+ if extraVal <> [] then
+ k = find(xitemp<min(x) | xitemp>max(x))
+ yi(k) = extraVal
+ end
+ yi = matrix(yi, xisize)
+
+ else // y is a matrix or a hypermatrix
+ ky = size(y)(2:$)
+ ncolumns = prod(ky)
+ yi = resize_matrix(0, [size(xitemp) ky])
+ x = matrix(x,-1,1)
+ xi = matrix(xi,-1,1)
+ for ic = 1:ncolumns
+ tmp = y(:,ic)
+ tmp = interp(xi, x, tmp, splin(x, tmp), extra)
+ yi(:,ic) = matrix(tmp, size(xitemp))
+ end
+ // extrapolation
+ if extraVal <> [] then
+ k = find(xitemp<min(x) | xitemp>max(x))
+ yi(k,:) = extraVal
+ end
+ yi = matrix(yi,[xisize,ky])
+ end
+
+ case "n"
+ // NEAREST
+ // -------
+ // if all xi values are nan, retuns nan values for yi
+ if and(isnan(xi)) then
+ yi = xi
+ return
+ end
+
+ xitemp = matrix(xi,1,-1)
+ knan = isnan(xitemp)
+ [xitemp,p] = gsort(xitemp(~knan),"g","i")
+ // For each xi, we search for the k index of the x nearest to xi:
+ k = zeros(xitemp);
+ nx = length(x);
+ j = length(xitemp);
+ i = nx;
+ while j>=1 & i>=1
+ if xitemp(j) >= x(i) then
+ if i <> nx then
+ k(j) = i;
end
- yi=matrix(yi,[xisize,ky])
+ j = j - 1;
else
- error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector or matrix expected.\n"),"interp1",2));
+ i = i - 1;
end
+ end
+ k(xitemp < x(1)) = 1;
+ k(xitemp >= x(nx)) = nx - 1;
+ i = find(xitemp >= matrix((x(k)+x(k+1))/2,size(k)));
+ if i~=[]
+ k(i) = k(i)+1;
+ 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
+ // y is vector
+ if isvector(y) then
+ yi = y(k)
+ yi = matrix(yi,1,-1)
+ // extrapolation
+ if extraVal <> [] then
+ n = find(xitemp<min(x) | xitemp>max(x) | isnan(xitemp))
+ yi(n) = extraVal
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
+ yitemp = yi
+ yi(p) = yitemp
+ ytemp = yi
+ yi = matrix(xi,1,-1)
+ yi(knan) = %nan
+ yi(~knan) = ytemp
+ yi = matrix(yi, xisize)
+
+ else // y is a matrix or a hypermatrix
+ ky = size(y)(2:$)
+ yi = resize_matrix(0, [size(xitemp,"*") ky])
+ ncolumns = prod(ky)
+ for ic = 1:ncolumns
+ yi(:,ic) = y(k,ic)
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)));
- if i~=[]
- k(i) = k(i)+1;
- end
- 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 = resize_matrix(0, [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
+ // extrapolation
+ if extraVal <> [] then
+ n = find(xitemp < min(x) | xitemp > max(x))
+ for ic = 1:ncolumns
+ yi(n,ic) = extraVal
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)));
- if i~=[]
- k(i) = k(i)+1;
- end
- 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 = resize_matrix(0, [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"));
+ yitemp = yi
+ for ic = 1:ncolumns
+ yi(p,ic) = yitemp(:,ic)
+ end
+ yitemp = yi
+ yi = resize_matrix(0, [size(xi,"*") ky])
+ for ic = 1:ncolumns
+ yi(knan,ic) = %nan
+ yi(~knan,ic) = yitemp(:,ic)
+ end
+ yi = matrix(yi,[xisize,ky])
end
end
endfunction
--- /dev/null
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2021 - Samuel GOUGEON - Le Mans Université
+//
+// This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- CLI SHELL MODE -->
+// <-- NO CHECK REF -->
+
+// Checking the size of the result
+// ===============================
+for m = ["linear", "spline", "nearest"]
+ for x = list(1:0.5:4, (1:0.5:4)')
+ // Y vector => size(result) == size(xi)
+ Xi = list(%pi, 1.2:0.2:2, (1.2:0.2:2)', 1+rand(2,5)*4, ..
+ 1+rand(1,2,5)*4, 1+rand(2,1,5)*4, 1+rand(2,5,6)*4);
+ Y = list(5:11, (5:11)');
+ for y = Y
+ for xi = Xi
+ r = interp1(x, y, xi, m);
+ assert_checkequal(size(r), size(xi));
+ end
+ end
+
+ // y = matrix or hypermatrix, xi = scalar or vector
+ Y = list(rand(7,3), rand(7,1,3), rand(7,3,4));
+ Xi = list(%pi, 1.2:0.2:2, (1.2:0.2:2)');
+ for y = Y
+ for xi = Xi
+ r = interp1(x, y, xi, m);
+ assert_checkequal(size(r), [length(xi) size(y)(2:$)]);
+ end
+ end
+
+ // Y and xi are matrices or hypermatrices
+ Y = list(rand(7,3), rand(7,1,3), rand(7,3,4));
+ Xi = list(1+rand(2,5)*4, 1+rand(1,2,5)*4, 1+rand(2,1,5)*4, 1+rand(2,5,6)*4);
+ for y = Y
+ for xi = Xi
+ r = interp1(x, y, xi, m);
+ assert_checkequal(size(r), [size(xi) size(y)(2:$)]);
+ end
+ end
+ end
+end
+
+
+// =============
+// y is a vector
+// =============
+x = 1:4;
+// method = "linear" (default)
+// ---------------------------
+// extrapolation "by_nan" (default)
+for y = list(2:5, (2:5)')
+ assert_checkequal(interp1(x,y, 0), %nan);
+ assert_checkequal(interp1(x,y,1.5), 2.5);
+ assert_checkequal(interp1(x,y,[0.5:5]), [%nan 2.5 3.5 4.5 %nan]);
+ assert_checkequal(interp1(x,y,[0.5:5]'), [%nan 2.5 3.5 4.5 %nan]');
+ assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5]), [%nan 2.5 3.5;4.5 %nan %nan]);
+ Ref = cat(3, [%nan 2.5 3.5], [4.5 %nan %nan]);
+ assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5])), Ref);
+end
+// extrapolation "extrap" = "linear"
+for y = list(2:5, (2:5)')
+ assert_checkequal(interp1(x,y, 0, ,"extrap"), 1);
+ assert_checkequal(interp1(x,y,1.5, ,"extrap"), 2.5);
+ assert_checkequal(interp1(x,y,[0.5:5], ,"extrap"), [1.5 2.5 3.5 4.5 5.5]);
+ assert_checkequal(interp1(x,y,[0.5:5]', ,"extrap"), [1.5 2.5 3.5 4.5 5.5]');
+ assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], ,"extrap"), ..
+ [1.5 2.5 3.5;4.5 5.5 6.5]);
+ Ref = cat(3, [1.5 2.5 3.5], [4.5 5.5 6.5]);
+ assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), ,"extrap"), Ref);
+end
+// extrapolation by padding
+for y = list(2:5, (2:5)')
+ assert_checkequal(interp1(x,y, 0, ,%pi), %pi);
+ assert_checkequal(interp1(x,y,1.5, ,%pi), 2.5);
+ assert_checkequal(interp1(x,y,[0.5:5], ,%pi), [%pi 2.5 3.5 4.5 %pi]);
+ assert_checkequal(interp1(x,y,[0.5:5]', ,%pi), [%pi 2.5 3.5 4.5 %pi]');
+ assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], ,%pi), ..
+ [%pi 2.5 3.5;4.5 %pi %pi]);
+ Ref = cat(3, [%pi 2.5 3.5], [4.5 %pi %pi]);
+ assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), ,%pi), Ref);
+end
+// periodic extrapolation
+for y = list([2 3 4 2], [2 3 4 2]')
+ assert_checkequal(interp1(x,y, 0, ,"periodic"), 4);
+ assert_checkequal(interp1(x,y,1.5, ,"periodic"), 2.5);
+ assert_checkequal(interp1(x,y,[0.5:5], ,"periodic"), [3 2.5 3.5 3 2.5]);
+ assert_checkequal(interp1(x,y,[0.5:5]', ,"periodic"), [3 2.5 3.5 3 2.5]');
+ assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], ,"periodic"), ..
+ [3 2.5 3.5 ; 3 2.5 3.5]);
+ Ref = cat(3, [3 2.5 3.5], [3 2.5 3.5]);
+ assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), ,"periodic"), Ref);
+end
+// extrapolation by edgevalue
+e = "edgevalue";
+for y = list(2:5, (2:5)')
+ assert_checkequal(interp1(x,y, 0, , e), 2);
+ assert_checkequal(interp1(x,y,1.5, , e), 2.5);
+ assert_checkequal(interp1(x,y,[0.5:5], , e), [2 2.5 3.5 4.5 5]);
+ assert_checkequal(interp1(x,y,[0.5:5]', , e), [2 2.5 3.5 4.5 5]');
+ assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], , e), [2 2.5 3.5 ; 4.5 5 5]);
+ Ref = cat(3, [2 2.5 3.5], [4.5 5 5]);
+ assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), , e), Ref);
+end
+
+// method = "nearest"
+// ------------------
+m = "nearest";
+// extrapolation "by_nan" (default)
+for y = list(2:5, (2:5)')
+ assert_checkequal(interp1(x,y, [0 5.4], m), [%nan %nan]);
+ assert_checkequal(interp1(x,y,1.5, m), 3);
+ assert_checkequal(interp1(x,y,[0.5:5], m), [%nan 3 4 5 %nan]);
+ assert_checkequal(interp1(x,y,[0.5:5]', m), [%nan 3 4 5 %nan]');
+ assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m), [%nan 3 4;5 %nan %nan]);
+ Ref = cat(3, [%nan 3 4], [5 %nan %nan]);
+ assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m), Ref);
+end
+// extrapolation "extrap" = "nearest" = "edgevalue"
+for e = ["extrap" "edgevalue" "nearest"]
+ for y = list(2:5, (2:5)')
+ assert_checkequal(interp1(x,y, [0 5.4], m, e), [2 5]);
+ assert_checkequal(interp1(x,y,1.5, m, e), 3);
+ assert_checkequal(interp1(x,y,[0.5:5], m, e), [2 3 4 5 5]);
+ assert_checkequal(interp1(x,y,[0.5:5]', m, e), [2 3 4 5 5]');
+ assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, e), ..
+ [2 3 4 ; 5 5 5]);
+ Ref = cat(3, [2 3 4], [5 5 5]);
+ assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, e), Ref);
+ end
+end
+// extrapolation by padding
+for y = list(2:5, (2:5)')
+ assert_checkequal(interp1(x,y, [0 5.4], m, %pi), [%pi %pi]);
+ assert_checkequal(interp1(x,y, 1.5, m, %pi), 3);
+ assert_checkequal(interp1(x,y, 0.5:5, m, %pi), [%pi 3 4 5 %pi]);
+ assert_checkequal(interp1(x,y,[0.5:5]', m, %pi), [%pi 3 4 5 %pi]');
+ assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, %pi), ..
+ [%pi 3 4 ; 5 %pi %pi]);
+ Ref = cat(3, [%pi 3 4], [5 %pi %pi]);
+ assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, %pi), Ref);
+end
+
+// method = "spline"
+// -----------------
+// extrapolation default == "spline" (=="natural")
+m = "spline";
+for y = list(2:5, (2:5)')
+ assert_checkalmostequal(interp1(x,y, [0 5.4], m), [1 6.4], 5*%eps);
+ assert_checkequal(interp1(x,y,1.5, m), 2.5);
+ assert_checkalmostequal(interp1(x,y,[0.5:5], m), 1.5:5.5, %eps);
+ assert_checkalmostequal(interp1(x,y,[0.5:5]', m), (1.5:5.5)', %eps);
+ assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m), [1.5 2.5 3.5;4.5 5.5 6.5], 5*%eps);
+ Ref = cat(3, [1.5 2.5 3.5], [4.5 5.5 6.5]);
+ assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m), Ref, 5*%eps);
+end
+// extrapolation "extrap" == "spline"
+for y = list(2:5, (2:5)')
+ assert_checkalmostequal(interp1(x,y, [0 5.4], m, "extrap"), [1 6.4], 5*%eps);
+ assert_checkequal(interp1(x,y,1.5, m, "extrap"), 2.5);
+ assert_checkalmostequal(interp1(x,y,[0.5:5], m, "extrap"), 1.5:5.5, %eps);
+ assert_checkalmostequal(interp1(x,y,[0.5:5]', m, "extrap"), (1.5:5.5)', %eps);
+ Ref = [1.5 2.5 3.5 ; 4.5 5.5 6.5];
+ assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, "extrap"), Ref, 5*%eps);
+ Ref = cat(3, [1.5 2.5 3.5], [4.5 5.5 6.5]);
+ assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, "extrap"), Ref, 5*%eps);
+end
+// extrapolation by padding
+for y = list(2:5, (2:5)')
+ assert_checkequal(interp1(x,y, [0 5.4], m, %pi), [%pi %pi]);
+ assert_checkequal(interp1(x,y, 1.5, m, %pi), 2.5);
+ assert_checkequal(interp1(x,y, 0.5:5, m, %pi), [%pi 2.5 3.5 4.5 %pi]);
+ assert_checkequal(interp1(x,y,[0.5:5]', m, %pi), [%pi 2.5 3.5 4.5 %pi]');
+ assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, %pi), ..
+ [%pi 2.5 3.5 ; 4.5 %pi %pi]);
+ Ref = cat(3, [%pi 2.5 3.5], [4.5 %pi %pi]);
+ assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, %pi), Ref);
+end
+// periodic extrapolation
+for y = list([2 3 4 2], [2 3 4 2]')
+ assert_checkalmostequal(interp1(x,y, [0 5.4], m, "periodic"), [4 3.736], %eps);
+ assert_checkequal(interp1(x,y, 1.5, m, "periodic"), 2.125);
+ assert_checkalmostequal(interp1(x,y, 0.5:5, m, "periodic"), [3 2.125 3.875 3 2.125], %eps);
+ assert_checkalmostequal(interp1(x,y,[0.5:5]', m, "periodic"), [3 2.125 3.875 3 2.125]', %eps);
+ assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, "periodic"), ..
+ [3 2.125 3.875 ; 3 2.125 3.875], %eps);
+ Ref = cat(3, [3 2.125 3.875], [3 2.125 3.875]);
+ assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, "periodic"), Ref, %eps);
+end
+// extrapolation by edgevalue
+e = "edgevalue";
+for y = list(2:5, (2:5)')
+ assert_checkequal(interp1(x,y, [0 5.4], m, e), [2 5]);
+ assert_checkequal(interp1(x,y, 1.5, m, e), 2.5);
+ assert_checkequal(interp1(x,y, 0.5:5, m, e), [2 2.5 3.5 4.5 5]);
+ assert_checkequal(interp1(x,y,[0.5:5]', m, e), [2 2.5 3.5 4.5 5]');
+ assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, e), ..
+ [2 2.5 3.5 ; 4.5 5 5]);
+ Ref = cat(3, [2 2.5 3.5], [4.5 5 5]);
+ assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, e), Ref);
+end
+
+// =============
+// y is a matrix
+// =============
+x = 1:4;
+y = [0 2 4 6 ; 0.5 1.5 2.5 3.5]';
+// method = "linear" (default)
+// ---------------------------
+// extrapolation "by_nan" (default)
+assert_checkequal(interp1(x,y, 0), [%nan %nan]);
+assert_checkequal(interp1(x,y,1.5), [1 1]);
+assert_checkequal(interp1(x,y, 0.5:5 ), [%nan 1 3 5 %nan; %nan 1 2 3 %nan]');
+assert_checkequal(interp1(x,y,(0.5:5)'), [%nan 1 3 5 %nan; %nan 1 2 3 %nan]');
+Ref = cat(3,[%nan 1 3 ; 5 %nan %nan], [%nan 1 2 ; 3 %nan %nan]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5]), Ref);
+Ref = matrix([%nan 1 3 5 %nan %nan %nan 1 2 3 %nan %nan],[1 3 2 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5])), Ref);
+
+// extrapolation "extrap" = "linear"
+assert_checkequal(interp1(x,y, 0, ,"extrap"), [-2 -0.5]);
+assert_checkequal(interp1(x,y,1.5, ,"extrap"), [1 1]);
+assert_checkequal(interp1(x,y, 0.5:4, ,"extrap"), [-1 1 3 5 ; 0 1 2 3]');
+assert_checkequal(interp1(x,y,(0.5:4)', ,"extrap"), [-1 1 3 5 ; 0 1 2 3]');
+Ref = cat(3, [-1 1 3 ; 5 7 9], [0 1 2 ; 3 4 5]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], ,"extrap"), Ref);
+Ref = matrix([-1 1 3 5 7 9 0 1 2 3 4 5],[1 3 2 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), ,"extrap"), Ref);
+
+// extrapolation by padding
+assert_checkequal(interp1(x,y, 0, ,%pi), [%pi %pi]);
+assert_checkequal(interp1(x,y,1.5, ,%pi), [1 1]);
+assert_checkequal(interp1(x,y, 0.5:5, ,%pi ), [%pi 1 3 5 %pi; %pi 1 2 3 %pi]');
+assert_checkequal(interp1(x,y,(0.5:5)', ,%pi), [%pi 1 3 5 %pi; %pi 1 2 3 %pi]');
+Ref = cat(3,[%pi 1 3 ; 5 %pi %pi], [%pi 1 2 ; 3 %pi %pi]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], ,%pi), Ref);
+Ref = matrix([%pi 1 3 5 %pi %pi %pi 1 2 3 %pi %pi],[1 3 2 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), ,%pi), Ref);
+
+// periodic extrapolation
+y = [0 2 4 0 ; 0.5 1.5 2.5 0.5]';
+assert_checkalmostequal(interp1(x,y, 5, ,"periodic"), [2 1.5], 2*%eps);
+assert_checkalmostequal(interp1(x,y,1.5, ,"periodic"), [1 1], 2*%eps);
+assert_checkalmostequal(interp1(x,y, 0.5:5, ,"periodic" ), [2 1 3 2 1; 1.5 1 2 1.5 1]', 2*%eps);
+assert_checkalmostequal(interp1(x,y,(0.5:5)', ,"periodic"), [2 1 3 2 1; 1.5 1 2 1.5 1]', 2*%eps);
+Ref = cat(3,[2 1 3 ; 2 1 3], [1.5 1 2 ; 1.5 1 2]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], ,"periodic"), Ref, 2*%eps);
+Ref = matrix([2 1 3 2 1 3 1.5 1 2 1.5 1 2],[1 3 2 2]);
+assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), ,"periodic"), Ref, 2*%eps);
+
+// extrapolation by edgevalue
+e = "edgevalue";
+y = [0 2 4 6 ; 0.5 1.5 2.5 3.5]';
+assert_checkequal(interp1(x,y, 0, , e), [0 0.5]);
+assert_checkequal(interp1(x,y,1.5, , e), [1 1]);
+assert_checkequal(interp1(x,y, 0.5:5, , e ), [0 1 3 5 6 ; 0.5 1 2 3 3.5]');
+assert_checkequal(interp1(x,y,(0.5:5)', , e), [0 1 3 5 6 ; 0.5 1 2 3 3.5]');
+Ref = cat(3,[0 1 3 ; 5 6 6], [0.5 1 2 ; 3 3.5 3.5]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], , e), Ref);
+Ref = matrix([0 1 3 5 6 6 0.5 1 2 3 3.5 3.5],[1 3 2 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), ,e), Ref);
+
+// method = "nearest"
+// ------------------
+m = "nearest";
+// extrapolation "by_nan" (default)
+assert_checkequal(interp1(x,y, [0 5.4], m), [%nan %nan ; %nan %nan]);
+assert_checkequal(interp1(x,y,1.5, m), [2 1.5]);
+assert_checkequal(interp1(x,y,[0.5:4], m), [%nan %nan ; 2 1.5 ; 4 2.5 ; 6 3.5]);
+assert_checkequal(interp1(x,y,[0.5:4]', m),[%nan %nan ; 2 1.5 ; 4 2.5 ; 6 3.5]);
+Ref = cat(3, [%nan 2 4 ; 6 %nan %nan], [%nan 1.5 2.5 ; 3.5 %nan %nan]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m), Ref); // OErr
+Ref = matrix([%nan 2 4 6 %nan %nan %nan 1.5 2.5 3.5 %nan %nan],[1 3 2 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m), Ref);
+
+// extrapolation "extrap" = "nearest" = "edgevalue"
+for e = ["extrap" "nearest" "edgevalue"]
+ assert_checkequal(interp1(x,y, [0 5.4], m, e), [0 0.5 ; 6 3.5]);
+ assert_checkequal(interp1(x,y,1.5, m, e), [2 1.5]);
+ assert_checkequal(interp1(x,y,[0.5:4], m, e), y);
+ assert_checkequal(interp1(x,y,[0.5:4]', m, e), y);
+ Ref = cat(3, [0 2 4 ; 6 6 6], [0.5 1.5 2.5 ; 3.5 3.5 3.5]);
+ assert_checkequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m, e), Ref);
+ Ref = matrix([0 2 4 6 6 6 0.5 1.5 2.5 3.5 3.5 3.5],[1 3 2 2]);
+ assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, e), Ref);
+end
+
+// extrapolation by padding
+assert_checkequal(interp1(x,y, [0 5.4], m, %pi), [%pi %pi ; %pi %pi]);
+assert_checkequal(interp1(x,y,1.5, m, %pi), [2 1.5]);
+assert_checkequal(interp1(x,y,[0.5:4], m, %pi), [%pi %pi ; 2 1.5 ; 4 2.5 ; 6 3.5]);
+assert_checkequal(interp1(x,y,[0.5:4]', m, %pi),[%pi %pi ; 2 1.5 ; 4 2.5 ; 6 3.5]);
+Ref = cat(3, [%pi 2 4 ; 6 %pi %pi], [%pi 1.5 2.5 ; 3.5 %pi %pi]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m, %pi), Ref); // OErr
+Ref = matrix([%pi 2 4 6 %pi %pi %pi 1.5 2.5 3.5 %pi %pi],[1 3 2 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, %pi), Ref); // OErr
+
+// method = "spline"
+// -----------------
+m = "spline";
+// extrapolation default == "spline" (=="natural")
+assert_checkalmostequal(interp1(x,y, [0 5.4], m), [-2 -0.5 ; 8.8 4.9], 5*%eps);
+assert_checkequal(interp1(x,y,1.5, m), [1 1]);
+assert_checkequal(interp1(x,y,[0.5:4], m), [-1 1 3 5 ; 0 1 2 3]');
+assert_checkequal(interp1(x,y,[0.5:4]', m), [-1 1 3 5 ; 0 1 2 3]');
+Ref = cat(3,[-1 1 3 ; 5 7 9], [0 1 2 ; 3 4 5]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m), Ref, 5*%eps);
+Ref = matrix([-1 1 3 5 7 9 0 1 2 3 4 5 ],[1 3 2 2]);
+assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m), Ref, 5*%eps);
+
+// extrapolation "extrap" == "spline"
+assert_checkalmostequal(interp1(x,y, [0 5.4], m, "extrap"), [-2 -0.5 ; 8.8 4.9], 5*%eps);
+assert_checkequal(interp1(x,y,1.5, m, "extrap"), [1 1]);
+assert_checkequal(interp1(x,y,[0.5:4], m, "extrap"), [-1 1 3 5 ; 0 1 2 3]');
+assert_checkequal(interp1(x,y,[0.5:4]', m, "extrap"), [-1 1 3 5 ; 0 1 2 3]');
+Ref = cat(3,[-1 1 3 ; 5 7 9], [0 1 2 ; 3 4 5]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, "extrap"), Ref, 5*%eps);
+Ref = matrix([-1 1 3 5 7 9 0 1 2 3 4 5 ],[1 3 2 2]);
+assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, "extrap"), Ref, 5*%eps);
+
+// extrapolation by padding
+assert_checkalmostequal(interp1(x,y, [0 5.4], m, %pi), [%pi %pi ; %pi %pi], 5*%eps);
+assert_checkequal(interp1(x,y,1.5, m, %pi), [1 1]);
+assert_checkequal(interp1(x,y,[0.5:4], m, %pi), [%pi 1 3 5 ; %pi 1 2 3]');
+assert_checkequal(interp1(x,y,[0.5:4]', m, %pi), [%pi 1 3 5 ; %pi 1 2 3]');
+Ref = cat(3,[%pi 1 3 ; 5 %pi %pi], [%pi 1 2 ; 3 %pi %pi]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, %pi), Ref, 5*%eps);
+Ref = matrix([%pi 1 3 5 %pi %pi %pi 1 2 3 %pi %pi],[1 3 2 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, %pi), Ref);
+
+// extrapolation by edgevalue
+e = "edgevalue";
+assert_checkequal(interp1(x,y, [0 5.4], m, e), [0 0.5 ; 6 3.5]);
+assert_checkequal(interp1(x,y,1.5, m, e), [1 1]);
+assert_checkequal(interp1(x,y,[0.5:4], m, e), [0 1 3 5 ; 0.5 1 2 3]');
+assert_checkequal(interp1(x,y,[0.5:4]', m, e), [0 1 3 5 ; 0.5 1 2 3]');
+Ref = cat(3, [0 1 3 ; 5 6 6], [0.5 1 2 ; 3 3.5 3.5]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, e), Ref);
+Ref = matrix([0 1 3 5 6 6 0.5 1 2 3 3.5 3.5], [1 3 2 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, e), Ref);
+
+// periodic extrapolation
+y = [0 2 4 0 ; 0.5 1.5 2.5 0.5]';
+assert_checkalmostequal(interp1(x,y, [0 5.4], m, "periodic"), [4 2.5 ; 3.136 2.068], 2*%eps);
+assert_checkalmostequal(interp1(x,y,1.5, m, "periodic"), [0.625 0.8125], 5*%eps);
+Ref = [3.125 2.0625 ; 0.625 0.8125 ; 3.375 2.1875 ; 3.125 2.0625];
+assert_checkequal(interp1(x,y,[0.5:4], m, "periodic"), Ref);
+assert_checkequal(interp1(x,y,[0.5:4]', m, "periodic"), Ref);
+Ref = cat(3, [3.125 0.625 3.375 ; 3.125 0.625 3.375], ..
+ [2.0625 0.8125 2.1875 ; 2.0625 0.8125 2.1875]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, "periodic"), Ref, 5*%eps);
+Ref = matrix([3.125 0.625 3.375 3.125 0.625 3.375 2.0625 0.8125 2.1875 2.0625 0.8125 2.1875],[1 3 2 2]);
+assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, "periodic"), Ref, 5*%eps);
+
+
+// ===================
+// y is an hypermatrix
+// ===================
+y = cat(4,[0 2 4 6 ; 0.5 1.5 2.5 3.5]', [1:4 ; 2:5]');
+// method = "linear" (default)
+// ---------------------------
+// extrapolation "by_nan" (default)
+assert_checkequal(interp1(x,y, 0), cat(4,[%nan %nan],[%nan %nan])); // OWrong size and values
+assert_checkequal(interp1(x,y,1.5), cat(4,[1 1],[1.5 2.5])); // OWrong size
+Ref = cat(4,[%nan 1 3 5 ; %nan 1 2 3]', [%nan 1.5:3.5 ; %nan 2.5:4.5]');
+assert_checkequal(interp1(x,y, 0.5:4), Ref);
+assert_checkequal(interp1(x,y,(0.5:4)'), Ref);
+Ref = matrix([%nan 5 1 %nan 3 %nan %nan 3 1 %nan 2 %nan %nan 3.5 1.5 %nan ..
+ 2.5 %nan %nan 4.5 2.5 %nan 3.5 %nan],[2 3 2 1 2]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5]), Ref);
+Ref = matrix([%nan 1 3 5 %nan %nan %nan 1 2 3 %nan %nan %nan 1.5 2.5 3.5 ..
+ %nan %nan %nan 2.5 3.5 4.5 %nan %nan],[1 3 2 2 1 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5])), Ref);
+
+// extrapolation "extrap" = "linear"
+assert_checkequal(interp1(x,y, 0, ,"extrap"), cat(4, [-2 -0.5],[0 1])); // OWrong size
+assert_checkequal(interp1(x,y,1.5, ,"extrap"), cat(4, [1 1],[1.5 2.5])); // OWrong size
+Ref = cat(4, [-1 1 3 5 ; 0 1 2 3]',[0.5:3.5 ; 1.5:4.5]');
+assert_checkequal(interp1(x,y, 0.5:4, ,"extrap"), Ref);
+assert_checkequal(interp1(x,y,(0.5:4)', ,"extrap"), Ref);
+
+Ref = matrix([-1 5 1 7 3 9 0 3 1 4 2 5 0.5 3.5 1.5 4.5 2.5 5.5 1.5 4.5 2.5 5.5 3.5 6.5],..
+ [2 3 2 1 2]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], ,"extrap"), Ref);
+Ref = matrix([-1 1 3 5 7 9 0 1 2 3 4 5 0.5 1.5 2.5 3.5 4.5 5.5 1.5 2.5 3.5 4.5 5.5 6.5],..
+ [1 3 2 2 1 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), ,"extrap"), Ref);
+
+// extrapolation by padding
+assert_checkequal(interp1(x,y, 0, , %pi), cat(4,[%pi %pi],[%pi %pi])); // OWrong size and values
+assert_checkequal(interp1(x,y,1.5, , %pi), cat(4,[1 1],[1.5 2.5]));
+Ref = cat(4,[%pi 1 3 5 ; %pi 1 2 3]', [%pi 1.5:3.5 ; %pi 2.5:4.5]');
+assert_checkequal(interp1(x,y, 0.5:4, , %pi), Ref);
+assert_checkequal(interp1(x,y,(0.5:4)', , %pi), Ref);
+Ref = matrix([%pi 5 1 %pi 3 %pi %pi 3 1 %pi 2 %pi %pi 3.5 1.5 %pi ..
+ 2.5 %pi %pi 4.5 2.5 %pi 3.5 %pi],[2 3 2 1 2]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], , %pi), Ref);
+Ref = matrix([%pi 1 3 5 %pi %pi %pi 1 2 3 %pi %pi %pi 1.5 2.5 3.5 ..
+ %pi %pi %pi 2.5 3.5 4.5 %pi %pi],[1 3 2 2 1 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), , %pi), Ref);
+
+// periodic extrapolation
+y = cat(4,[0 2 4 0 ; 0.5 1.5 2.5 0.5]', [1 2 3 1 ; 2 3 4 2]');
+assert_checkequal(interp1(x,y, 0, , "periodic"), cat(4,[4 2.5],[3 4]));
+assert_checkequal(interp1(x,y,1.5, , "periodic"), cat(4,[1 1],[1.5 2.5]));
+Ref = cat(4,[2 1 3 2 ; 1.5 1 2 1.5]', [2 1.5 2.5 2 ; 3 2.5 3.5 3]');
+assert_checkequal(interp1(x,y, 0.5:4, , "periodic"), Ref);
+assert_checkequal(interp1(x,y,(0.5:4)', , "periodic"), Ref);
+Ref = matrix([2 2 1 1 3 3 1.5 1.5 1 1 2 2 2 2 1.5 1.5 ..
+ 2.5 2.5 3 3 2.5 2.5 3.5 3.5],[2 3 2 1 2]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], , "periodic"), Ref, 2*%eps);
+Ref = matrix([2 1 3 2 1 3 1.5 1 2 1.5 1 2 2 1.5 2.5 2 ..
+ 1.5 2.5 3 2.5 3.5 3 2.5 3.5],[1 3 2 2 1 2]);
+assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), , "periodic"), Ref, 2*%eps);
+
+// extrapolation by edgevalue
+e = "edgevalue";
+y = cat(4,[0 2 4 6 ; 0.5 1.5 2.5 3.5]', [1:4 ; 2:5]');
+assert_checkequal(interp1(x,y, 0, , e) , cat(4,[0 0.5],[1 2]));
+assert_checkequal(interp1(x,y,1.5, , e), cat(4,[1 1],[1.5 2.5]));
+Ref = cat(4,[0 1 3 5 ; 0.5 1 2 3]', [1 1.5:3.5 ; 2 2.5:4.5]');
+assert_checkequal(interp1(x,y, 0.5:4, , e), Ref);
+assert_checkequal(interp1(x,y,(0.5:4)', , e), Ref);
+Ref = matrix([0 5 1 6 3 6 0.5 3 1 3.5 2 3.5 1 3.5 1.5 4 ..
+ 2.5 4 2 4.5 2.5 5 3.5 5],[2 3 2 1 2]);
+assert_checkequal(interp1(x, y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], , e), Ref);
+Ref = matrix([0 1 3 5 6 6 0.5 1 2 3 3.5 3.5 1 1.5 2.5 3.5 ..
+ 4 4 2 2.5 3.5 4.5 5 5],[1 3 2 2 1 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), , e), Ref);
+
+// method = "nearest"
+// ------------------
+m = "nearest";
+// extrapolation "by_nan" (default)
+assert_checkequal(interp1(x,y, [0 5.4], m), %nan*ones(2,2,1,2));
+assert_checkequal(interp1(x,y,1.5, m), cat(4,[2 1.5],[2 3]));
+Ref = cat(4, [%nan %nan ; 2 1.5 ; 4 2.5 ; 6 3.5], [%nan %nan ; 2 3 ; 3 4 ; 4 5]);
+assert_checkequal(interp1(x,y,[0.5:4], m), Ref);
+assert_checkequal(interp1(x,y,[0.5:4]', m), Ref);
+Ref = matrix([%nan 6 2 %nan 4 %nan %nan 3.5 1.5 %nan 2.5 %nan %nan 4 2 %nan 3 %nan %nan 5 3 %nan 4 %nan],[2 3 2 1 2]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m), Ref); // OErr
+Ref = matrix([%nan 2 4 6 %nan %nan %nan 1.5 2.5 3.5 %nan %nan %nan 2 3 4 %nan %nan %nan 3 4 5 %nan%nan],[1 3 2 2 1 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m), Ref); // OErr
+
+// extrapolation "extrap" = "nearest"
+for e = ["extrap" "nearest" "edgevalue"]
+ assert_checkequal(interp1(x,y, [0 5.4], m, e), cat(4,[0 0.5 ; 6 3.5],[1 2 ; 4 5]));
+ assert_checkequal(interp1(x,y,1.5, m, e), cat(4,[2 1.5],[2 3])); // OWrong values
+ assert_checkequal(interp1(x,y,[0.5:4], m, e), y);
+ assert_checkequal(interp1(x,y,[0.5:4]', m, e), y);
+ Ref = matrix([0 6 2 6 4 6 0.5 3.5 1.5 3.5 2.5 3.5 1 4 2 4 3 4 2 5 3 5 4 5], [2 3 2 1 2]);
+ assert_checkequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m, e), Ref);
+ Ref = matrix([0 2 4 6 6 6 0.5 1.5 2.5 3.5 3.5 3.5 1 2 3 4 4 4 2 3 4 5 5 5],[1 3 2 2 1 2]);
+ assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, e), Ref);
+end
+
+// extrapolation by padding
+assert_checkequal(interp1(x,y, [0 5.4], m, %pi), %pi*ones(2,2,1,2));
+assert_checkequal(interp1(x,y,1.5, m, %pi), cat(4,[2 1.5],[2 3]));
+Ref = cat(4, [%pi %pi ; 2 1.5 ; 4 2.5 ; 6 3.5], [%pi %pi ; 2 3 ; 3 4 ; 4 5]);
+assert_checkequal(interp1(x,y,[0.5:4], m, %pi), Ref);
+assert_checkequal(interp1(x,y,[0.5:4]', m, %pi), Ref);
+Ref = matrix([%pi 6 2 %pi 4 %pi %pi 3.5 1.5 %pi 2.5 %pi %pi 4 2 %pi 3 %pi %pi 5 3 %pi 4 %pi],[2 3 2 1 2]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m, %pi), Ref); // OErr
+Ref = matrix([%pi 2 4 6 %pi %pi %pi 1.5 2.5 3.5 %pi %pi %pi 2 3 4 %pi %pi %pi 3 4 5 %pi%pi],[1 3 2 2 1 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, %pi), Ref); // OErr
+
+// method = "spline"
+// -----------------
+m = "spline";
+// extrapolation default == "spline" (=="natural")
+Ref = cat(4, [-2 -0.5 ; 8.8 4.9], [0 1 ; 5.4 6.4]);
+assert_checkalmostequal(interp1(x,y, [0 5.4], m), Ref, 5*%eps);
+assert_checkequal(interp1(x,y,1.5, m), cat(4,[1 1],[1.5 2.5])); // OWrong values
+Ref = cat(4,[-1 1 3 5 ; 0 1 2 3]',[0.5:3.5 ; 1.5:4.5]');
+assert_checkequal(interp1(x,y,[0.5:4], m), Ref);
+assert_checkequal(interp1(x,y,[0.5:4]', m), Ref);
+Ref = matrix([-1 5 1 7 3 9 0 3 1 4 2 5 0.5 3.5 1.5 4.5 2.5 5.5 1.5 4.5 2.5 5.5 3.5 6.5],..
+ [2 3 2 1 2]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m), Ref, 5*%eps);
+Ref = matrix([-1 1 3 5 7 9 0 1 2 3 4 5 0.5 1.5 2.5 3.5 4.5 5.5 1.5 2.5 3.5 4.5 5.5 6.5],..
+ [1 3 2 2 1 2]);
+assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m), Ref, 5*%eps);
+
+// extrapolation "extrap" == "spline"
+Ref = cat(4, [-2 -0.5 ; 8.8 4.9], [0 1 ; 5.4 6.4]);
+assert_checkalmostequal(interp1(x,y, [0 5.4], m, "extrap"), Ref, 5*%eps); // OWrong size
+assert_checkequal(interp1(x,y,1.5, m, "extrap"), cat(4, [1 1], [1.5 2.5])); // OWrong size
+Ref = cat(4, [-1 1 3 5; 0:3]', [0.5:3.5 ; 1.5:4.5]');
+assert_checkequal(interp1(x,y,[0.5:4], m, "extrap"), Ref);
+assert_checkequal(interp1(x,y,[0.5:4]', m, "extrap"), Ref);
+Ref = matrix([-1 5 1 7 3 9 0 3 1 4 2 5 0.5 3.5 1.5 4.5 2.5 5.5 1.5 4.5 2.5 5.5 3.5 6.5],..
+ [2 3 2 1 2]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m, "extrap"), Ref, 5*%eps); // OWrong size
+Ref = matrix([-1 1 3 5 7 9 0 1 2 3 4 5 0.5 1.5 2.5 3.5 4.5 5.5 1.5 2.5 3.5 4.5 5.5 6.5],..
+ [1 3 2 2 1 2]);
+assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, "extrap"), Ref, 5*%eps);
+
+// extrapolation by padding
+assert_checkalmostequal(interp1(x,y, [0 5.4], m, %pi), %pi*ones(2,2,1,2), 5*%eps);
+assert_checkequal(interp1(x,y,1.5, m, %pi), cat(4,[1 1],[1.5 2.5])); // OWrong values
+Ref = cat(4,[%pi 1 3 5 ; %pi 1 2 3]',[%pi 1.5 2.5 3.5 ; %pi 2.5 3.5 4.5]');
+assert_checkequal(interp1(x,y,[0.5:4], m, %pi), Ref);
+assert_checkequal(interp1(x,y,[0.5:4]', m, %pi), Ref);
+Ref = matrix([%pi 5 1 %pi 3 %pi %pi 3 1 %pi 2 %pi %pi 3.5 1.5 %pi 2.5 %pi %pi 4.5 2.5 %pi 3.5 %pi], [2 3 2 1 2]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m, %pi), Ref, 5*%eps);
+Ref = matrix([%pi 1 3 5 %pi %pi %pi 1 2 3 %pi %pi %pi 1.5 2.5 3.5 %pi %pi %pi 2.5 3.5 4.5 %pi %pi], [1 3 2 2 1 2]);
+assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, %pi), Ref, 5*%eps);
+
+// extrapolation by edgevalue
+e = "edgevalue";
+assert_checkequal(interp1(x,y, [0 5.4], m, e), cat(4, [0 0.5 ; 6 3.5],[1 2 ; 4 5]));
+assert_checkequal(interp1(x,y,1.5, m, e), cat(4,[1 1],[1.5 2.5]));
+Ref = cat(4, [0 1 3 5 ; 0.5 1 2 3]', [1 1.5 2.5 3.5 ; 2 2.5 3.5 4.5]');
+assert_checkequal(interp1(x,y,[0.5:4], m, e), Ref);
+assert_checkequal(interp1(x,y,[0.5:4]', m, e), Ref);
+Ref = matrix([0 5 1 6 3 6 0.5 3 1 3.5 2 3.5 1 3.5 1.5 4 2.5 4 2 4.5 2.5 5 3.5 5], [2 3 2 1 2]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m, e), Ref);
+Ref = matrix([0 1 3 5 6 6 0.5 1 2 3 3.5 3.5 1 1.5 2.5 3.5 4 4 2 2.5 3.5 4.5 5 5], [1 3 2 2 1 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, e), Ref);
+
+// periodic extrapolation
+y = cat(4,[0 2 4 0 ; 0.5 1.5 2.5 0.5]', [1 2 3 1 ; 2 3 4 2]');
+Ref = cat(4, [4 2.5 ; 3.136 2.068], [3 4 ; 2.568 3.568]);
+assert_checkalmostequal(interp1(x,y, [0 5.4], m, "periodic"), Ref, 5*%eps);
+assert_checkequal(interp1(x,y,1.5, m, "periodic"), cat(4,[0.625 0.8125],[1.3125 2.3125]));
+Ref = cat(4, [3.125 2.0625 ; 0.625 0.8125 ; 3.375 2.1875 ; 3.125 2.0625], ..
+ [2.5625 3.5625; 1.3125 2.3125 ; 2.6875 3.6875 ; 2.5625 3.5625]);
+assert_checkequal(interp1(x,y,[0.5:4], m, "periodic"), Ref);
+assert_checkequal(interp1(x,y,[0.5:4]', m, "periodic"), Ref);
+Ref = matrix([3.125 3.125 0.625 0.625 3.375 3.375 2.0625 2.0625 0.8125 0.8125 ..
+ 2.1875 2.1875 2.5625 2.5625 1.3125 1.3125 2.6875 2.6875 3.5625 3.5625 ..
+ 2.3125 2.3125 3.6875 3.6875], [2 3 2 1 2]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m, "periodic"), Ref, 5*%eps);
+Ref = matrix([3.125 0.625 3.375 3.125 0.625 3.375 2.0625 0.8125 2.1875 2.0625 ..
+ 0.8125 2.1875 2.5625 1.3125 2.6875 2.5625 1.3125 2.6875 3.5625 2.3125 ..
+ 3.6875 3.5625 2.3125 3.6875], [1 3 2 2 1 2]);
+assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, "periodic"), Ref, 5*%eps);
+
+
+// y is complex
+// ============
+// method = "linear" (default)
+// ---------------------------
+assert_checkequal(interp1(0:3, (1:4)+%i, 0.5:2.5), (1.5:3.5)+%i);
+// extrapolation linear
+assert_checkequal(interp1(0:3,(1:4)+%i, 0.5:2:5, , "extrap"), [1.5 3.5 5.5]+%i);
+// Default padding with %nan
+assert_checkequal(interp1((1:4)+%i, 0.5:2.5), [%nan, 1.5+%i, 2.5+%i]);
+// Padding real
+assert_checkequal(interp1(1:3, (2:4)+%i, 0.5:2.5,,7), [7, 2.5+%i, 3.5+%i]);
+// Padding complex
+assert_checkequal(interp1(1:3, (2:4)+%i, [0.5 2.5 4.5],,7-%i), [7-%i, 3.5+%i, 7-%i]);
+
+// method = "nearest"
+// ------------------
+assert_checkequal(interp1(0:3, (1:4)+%i, 0.5:2.5, "nearest"), (2:4)+%i);
+// Default padding with %nan
+assert_checkequal(interp1((1:4)+%i, 0.5:2.5, "nearest"), [%nan, 2+%i, 3+%i]);
+// extrapolation with nearest
+assert_checkequal(interp1((1:4)+%i, 0.5:2:5, "nearest", "extrap"), [1 3 4]+%i);
+// Padding real
+assert_checkequal(interp1(1:3,(2:4)+%i, 0.5:2.5, "nearest", 7), [7, 3+%i, 4+%i]);
+// Padding complex
+assert_checkequal(interp1(1:3,(2:4)+%i, [0.5 2.5 4.5], "nearest",7-%i), [7-%i, 4+%i, 7-%i]);
+
+// method = "spline"
+// -----------------
+assert_checkequal(interp1(0:3, (1:4)+%i, 0.5:2.5, "spline"), (1.5:3.5)+%i);
+// Default padding with spline
+assert_checkequal(interp1((1:4)+%i, 0.5:2.5, "spline"), (0.5:2.5)+%i);
+assert_checkequal(interp1((1:4)+%i, 0.5:2.5, "spline", "extrap"), (0.5:2.5)+%i);
+// Padding real
+assert_checkequal(interp1(1:3,(2:4)+%i, 0.5:2.5, "spline", 7), [7, 2.5+%i, 3.5+%i]);
+// Padding complex
+assert_checkequal(interp1(1:3,(2:4)+%i, [0.5 2.5 4.5], "spline",7-%i), [7-%i, 3.5+%i, 7-%i]);
+
+
+// ==============
+// Error messages
+// ==============
+msg = _("%s: Wrong number of input arguments: Must be between %d and %d.\n")
+msg = msprintf(msg, "interp1", 2, 5);
+assert_checkerror("interp1()", msg);
+assert_checkerror("interp1(1)", msg);
+assert_checkerror("interp1(1,2,3,4,5,6)", msg);
+// x
+msg = msprintf(_("%s: Argument #%d: Real numbers expected.\n"),"interp1",1);
+assert_checkerror("interp1([1 %i],[1 2],1.5)", msg);
+msg = msprintf(_("%s: Argument #%d: Vector expected.\n"), "interp1", 1);
+assert_checkerror("interp1([1 2; 3 4],[1 2],1.5)", msg);
+assert_checkerror("interp1(1,[1 2],1.5)", msg);
+msg = msprintf(_("%s: Argument #%d: Nan value forbidden.\n"), "interp1", 1);
+assert_checkerror("interp1([1 2 %nan],[1 2],1.5)", msg);
+// y
+msg = msprintf(_("%s: Argument #%d: vector or at least 2 rows expected.\n"), "interp1", 2);
+assert_checkerror("interp1(1,1.5)", msg);
+msg = msprintf(_("%s: Argument #%d: Decimal or complex numbers expected.\n"), "interp1",1);
+assert_checkerror("interp1([1 %z],1.5)", msg);
+assert_checkerror("interp1([%t %f],1.5)", msg);
+//
+msg = msprintf(_("%s: Arguments #%d and #%d: Same numbers of elements expected.\n"),"interp1", 1, 2);
+assert_checkerror("interp1([1 2],1:3,3)", msg);
+assert_checkerror("interp1(1:3,1:2,3)", msg);
+msg = msprintf(_("%s: Arguments #%d and #%d: Same numbers of rows expected.\n"), "interp1", 1, 2);
+assert_checkerror("interp1(1:3,rand(2,3),3)", msg);
+ // periodicity
+msg = msprintf(_("%s: Argument #%d: periodicity y($)==y(1) expected.\n"), "interp1", 2);
+assert_checkerror("interp1(1:4,1:4,3,,""periodic"")", msg);
+msg = msprintf(_("%s: Argument #%d: periodicity y($)==y(1) expected.\n"), "interp1", 1);
+assert_checkerror("interp1(1:4,3,""linear"",""periodic"")", msg);
+msg = msprintf(_("%s: Argument #%d: periodicity y($,:)==y(1,:) expected.\n"), "interp1", 2);
+assert_checkerror("interp1(1:3,rand(3,2),3,,""periodic"")", msg);
+msg = msprintf(_("%s: Argument #%d: periodicity y($,:)==y(1,:) expected.\n"), "interp1", 1);
+assert_checkerror("interp1(rand(3,2),3,""linear"",""periodic"")", msg);
+
+// xp
+msg = msprintf(_("%s: Argument #%d: Decimal numbers expected.\n"), "interp1",2);
+assert_checkerror("interp1(1:3,[%i 1 1])", msg);
+msg = msprintf(_("%s: Argument #%d: Decimal numbers expected.\n"), "interp1",3);
+assert_checkerror("interp1(1:3,rand(3,2),%i)", msg);
+assert_checkerror("interp1(1:3,2:4,,""linear"")", msg);
+// method
+msg = msprintf(_("%s: Argument #%d: Scalar (1 element) expected.\n"), "interp1", 3);
+assert_checkerror("interp1(1:3,2,[""linear"" ""spline""])", msg);
+msg = msprintf(_("%s: Argument #%d: Scalar (1 element) expected.\n"), "interp1", 4);
+assert_checkerror("interp1(1:3,2:4,5,[""linear"" ""spline""])", msg);
+Set = """linear"",""spline"",""nearest""";
+msg = msprintf(_("%s: Argument #%d: Must be in the set {%s}.\n"), "interp1", 3, Set);
+assert_checkerror("interp1(1:3,5,""splines"")", msg);
+msg = msprintf(_("%s: Argument #%d: Must be in the set {%s}.\n"), "interp1", 4, Set);
+assert_checkerror("interp1(2:4,1:3,5,""splines"")", msg);
+assert_checkerror("interp1(2:4,1:3,5,%z)", msg);
+// extrapolation
+msg0 = _("%s: Argument #%d: scalar number or one of {%s} expected.\n");
+msg = msprintf(msg0, "interp1", 5, """extrap"",""linear"",""periodic"",""edgevalue""");
+assert_checkerror("interp1(1:3,2:4,5,,""abc"")", msg);
+assert_checkerror("interp1(1:3,2:4,5,,%z)", msg);
+assert_checkerror("interp1(1:3,2:4,5,,1:2)", msg);
+assert_checkerror("interp1(1:3,2:4,5,,int8(3))", msg);
+msg = msprintf(msg0, "interp1", 5, """extrap"",""edgevalue"",""nearest""");
+assert_checkerror("interp1(1:3,2:4,5,""nearest"",""linear"")", msg);
+assert_checkerror("interp1(1:3,2:4,5,""nearest"",""periodic"")", msg);
+msg = msprintf(msg0, "interp1", 4, """extrap"",""linear"",""periodic"",""edgevalue""");
+assert_checkerror("interp1(1:3,2:4,""linear"",""abc"")", msg);
+assert_checkerror("interp1(1:3,2:4,""linear"",%z)", msg);
+assert_checkerror("interp1(1:3,2:4,""linear"",1:2)", msg);
+assert_checkerror("interp1(1:3,2:4,""linear"",int8(3))", msg);
+msg = msprintf(msg0, "interp1", 4, """extrap"",""edgevalue"",""nearest""");
+assert_checkerror("interp1(1:3,2:4,""nearest"",""linear"")", msg);
+assert_checkerror("interp1(1:3,2:4,""nearest"",""periodic"")", msg);
projects / scilab.git / commitdiff
