I should be more busy to continue the series about problems solved with Euler Math Toolbox. This might be a good advertisement of EMT, and it might serve as a guide for users.

Suppose I am given some data yd. For a test, we generate and plot random data as follows.

>n=20; i=1:n; xd=i/n; yd=xd^3-xd+normal(n)*0.1;
>plot2d(xd,yd,>addpoints,thickness=2):

See below for the data and smoothing functions. Our first try (in red) is a quadratic polynomial.

>p2=polyfit(xd,yd,2); plot2d("polyval(p2,x)",color=red,>add):

For a more sophisticated solution, we can use weights depending on t for each t in the interval where we need our smooth function. The weights depend on the distance to t, of course. We can use the density of the normal distribution with mean value t to compute the weight at each xd.

>function map smooth (t;x,y,sigma=1) ...
$  w=qnormal(x,t,sigma); w=w/sum(w);
$  return sum(w*y);
$  endfunction
>plot2d("smooth(x,xd,yd,0.1)",color=green,>add):
>plot2d("smooth(x,xd,yd,0.3)",color=lightgreen,>add):

The function smooth() works only for scalar t. So we vectorize it with „map“. The semicolon in the parameter list keeps it from vectorizing to the vectors x and y. Also, note that we need to normalize our weights to the sum 1. We add two functions in green and light green with different spreads in the normal distribution. The light function incorporates more values.

There is another trick. We compute new data y from yd by minimizing

\(h(y) = \|y-y_d\|^2 + c \|\Delta^2 y\|^2.\)

The Delta operator is the operator taking twice the differences of the elements of the vector y. The norm of the vector of these double differences resembles the curvature of y. Since taking the differences is a linear process, the process has a matrix.

>function vdiag (v:vector,n:index) ...
$  k=cols(v); M=zeros(n,k+n-1);
$  for i=1 to n; M[i,i:i+k-1]=v; end;
$  return M;
$  endfunction
>shortest vdiag([1,-2,1],4)
           1        -2         1         0         0         0 
           0         1        -2         1         0         0 
           0         0         1        -2         1         0 
           0         0         0         1        -2         1

By the way, the coefficients of the Delta operator are connected to the binomial distribution by

\(\Delta^k = (S-id)^k = \sum\limits_{i=}^k (-1)^i \binom{k}{i} S^i\)

where S is the shift operator. This applies only to data on equidistant points. Otherwise, the matrix has to be computed line by line.

The function h(y) consists of a the distance of the points plus a punishment for the curvature. It can be minimized by

\(0 = \text{grad} \, h(y) = 2 (y-y_d)^T + 2 c y^T K K^T\)

Now we are ready to compute smoother y values. We connect the values with a spline.

>function smoothopt (y,c=1) ...
$  n=cols(y);
$  K=vdiag([1,-2,1],n);
$  return ((id(n)+c*K.K')\y')';
$  endfunction
>y=smoothopt(yd,10); s=spline(xd,y);
>plot2d("splineval(x,xd,y,s)",color=blue,>add);
>plot2d(xd,y,>points,style="+",color=blue,>add):