Endymion (bep)

One part of the Endymion program deals with solving the Bounded Extremal Problem. Given a set of data points, this finds an analytic function that fits the data. Several steps are needed, and a data structure, named bep, is used.

If you put foo in the file field, then foo_C will hold the result of conversion to the circle, foo_s will hold the splines, foo_F will hold the Fourier coefficients of the function, foo_G will hold the result of BEP, and foo_A will hold the input to the arl2 procedure. You can put foo[1/4] or foo[2][3] in the shared_input field, the same file names will be used. However the input file foo or foo.s3p is assumed to contain 4 or 3x3 functions and only the first or second one will be used (see parse_s2p for details). The function bep_execute executes all steps one after the others. For instance, it computes the value of fct_spline from fct_circle (found in foo_C) and stores it in foo_s; in some cases there is no need to recompute, and the content of the file foo_s is read. This computation can also be done explictly, case where nothing is written to or read from a file.

new_bep() (Symbolic function)

check_bep() (Symbolic function)

This function checks some fields (in particulat the file field). After that, the BEP objects prints differently: you will see all fields. An example is given below.

(1) X:=new_bep();

(1)                                

(2) set_bep("file","foo");

(2)                                  foo

(3) check_bep(X);

(3)  <bep2: file = foo, verbose = 0, data_on_circle = FALSE, re_inf = 0, 

im_inf = 0, inf_order = 1, symmetric = FALSE, fourier = 201, sr = 0, M = 1, 

precision = 0.0001, lambda = -0.5, gamma = 0, alpha = 0, beta = 0, nw0 = 0, 

ndw = 1, scale = 0, real_transform = FALSE, shrink_circle = FALSE, 

circle_len = 0, min = -1, max = 1, min_xn = 0, max_xn = 0, min_yn = 0, 

max_yn = 0, nb_points = 1000, nb_nyquist = 1000, nb_bode = 1000, 

data_min = 0, data_max = 0, otherflags = 0, continuity = FALSE, 

real = FALSE, cplot = FALSE, h = FALSE, nbdata = 0, have_data_min = 0, 

sobolev = FALSE, w0 = 0, dw = 1>

get_bep(B,F) (Symbolic function)

set_bep(B,F,V) (Symbolic function)

This function sets the value of the field F of the BEP B to V. The function returns false. Note that at BEP can be considered as an array so that B[a]:=B[b] is the same as set_bep(B,a, get_bep(B,b)). The value at infinity has real part re_inf and imaginary part im_inf. You can get (not set) it via inf_val. You should not set quantities

(4) X[inf_order]++;

(4)                                   2

(5) set_bep(X, "inf_order", 1+ get_bep(X,"inf_order"));

(5)                                 FALSE

(6)  X[inf_order];

(6)                                   3

handle_data_line(B) (Symbolic function)

This function takes one argument, a BEP, and realises step one. It modifiles fct_data, w0 dw and real_coef. The Endymion code is the function is the following.

Handle_data_line(B):=
  let([F=B["fct_data"],L=0,first=0,last=0,c=B["circle_len"],
      dm=B["data_min"],dM=B["data_max"], hdm=B["have_data_min"]],
    #:pw:sub_inf(F,B["inf_val"]),
    first := F["theta",0],
    last := F["theta",size(F)-1],
    if c<=0 then c:=1,
    c := c/(if B["symmetric"] then last else (last-first)/2),
    if hdm>0 then 
       { #:pw:skip_unwanted(F,dm,dM,hdm),
          first :=  F["theta",0],
          last := F["theta",size(F)-1]},
    if B["symmetric"] then { B["dw"]:= last,  B["w0"]:= 0 }
    else { B["dw"]:= (last-first)/2,  B["w0"]:= (last+first)/2 },
    if B["real_transform"] then c := c* B["w0"]/2,
    B["real_coef"] := c,
   done);

The following points are of interest. Endymion handles stricly proprer rational functions. If the function F is not strictly proper, and its value at infinity is known to be C, you can subtract C, otherwise you must guess the value. In some cases, you want to restrict the data interval to [dm,dM]. The data interval is [w0-dw,w0+dw]. Later on, it will be mapped to [-1,1] via a linear transformation, then to the half circle (see below). The rational approximation program finds a discrete-time approximation, and converts it back into a normalised contuinuous-time approximation. Here normalised means that the band width is [-1,1]. The quantity real_coef is the scaling factor to use. It is the inverse of half the length of the data interval (before truncation). A correction term is added if real_transform is true.

Assume that the data points are in the interval I. The interval is called symmetric in the continuous time case if it is of the form [-a,a]. In the discrete-time case, the condition becomes [pi-a,pi+a]. The completion algorithm works only for symmetric intervals. If F(s) is a real continuous-time transfert function, then F(-i*ω) is the complex conjugate of F(i*ω). If the symmetric flag is given, the function is assumed real, and the data file contains F(i*ω) for positive ω this explains the settings w0 and dw.

to_circle(B) (Symbolic function)

This function takes as argument a BEP structure, and performs the conversion from the line to the circle. We show here the Endymion equivalent of the C++ code. The objective is to compute the value R of fct_circle, from the value F of fct_data; these two objects are of type pwfunction, that is a sequence (xi,yi), representing a function that associates yi to the point xi. In the discrete-time case xi is a point on the circle (complex number z=cos(theta) +i* sin(theta)), and the data structure holds theta, a number between 0 and 2pi) and in the continuous-time case, it is a number between minus infinity and plus infinity (generally called ω).

The command to_circle converts (if necessary) continuous-time case into discrete-time case. If continuity is required, we subtract from the data the function az+b (where z=exp(i*θ)), where a and b are chosen so that the resulting function vanishes at the endpoints.

tocircle(B):=
 let([on_circle = B["data_on_circle"],
      ow0 = B["ow0"],
      F= B["fct_data"], res=0],
    res:=#copy(F),
    if(on_circle) then cvcircle0(res, F,ow0)
    else cvcircle1(res, F, ow0, real_trans, B),
    if B["continuity"] then #:pw:sub_line(res),
    B["fct_circle"]:=res)

This piece of code is called in the discrete-time case. The data interval is [w0-dw,w0+dw]. We shift the values so that w0 is replaced by pi. This makes tha data interval symmetric.

cvcircle0 (res, F, ow0):=
 let([n=size(F), mw0=%pi-ow0],
   for i in 0...n-1 do {
     res["theta",i] := mw0+F["theta",i],
     res[i] := F[i] });

The Moebius function associates to s=i*ω the value z=(s-1)/(s+1). Its inverse function φ(z) associates to z the value s=(1+z)/(1-z). We have Re s = (1-|z|^2)/(|1-z|).. As a consequence, φ(infinity)=1, and the real part of s is positive if and only if z is in the unit disk. If f is a function, we consider &sigma1(f)=g1, &sigma2(f)=g2 and &sigma3(f)=g3, defined as follows. g1(z)=f(phi(z)). If f(s) is rational and analytic in the half plane, then g1 is rational an analytic in the disk. We have g2(z)=g1(z)/(1-z) Assume f=p/q where p and q are polynomials, of degree at most n p=sum ak*s^k, q=sum bk*s^k. Then g1=p1/q1, where the coefficients of p1 are sum ak*(1+z)^k(1-z)^(n-k). Thus g1 is of degree at most n. The same holds, for g2, provided that f is strictly proper (the degree of p is less than the degree of q). Define tilde(g) to be cc,g(1/z)/z. Here cc stands for complex conjugate. In some cases, Endymion replaces all coefficients by their complex conjugate, in some other cases, this is not done. For simplicity, let us forget this for the moment being. We define g3 to be tilde(g2) . Hence g3(z)=f((z+1)/(z-1))/(z-1). If f is strictly proper, of degree at n, analytic in the right half plane, then g3 is of rational of degree at most n, analytic outside the unit disk. The denominator of g2 is sum bk*(z+1)^k*(z-1)*(n-k) and the numerator is sum ak*(z+1)^k*(z-1)*(n-k-1). The coefficient of z^n is the denominator is q(1), assumed to be non-zero. Thus g3 is strictly proper, of degree n. Note that, if the input is a matrix F, and if we compute a matrix G, via this technique, then F and G have the same McMillan degree (this generalises the degree of a rational function). If g2 =sum g_k *z^k, then g3 =sum g_k *z^(-k-1). The aim of the arl2 program is to find an approximation of g3 in the form P/Q (where the roots of Q are in the disk), and the BEP program computes g2.

This point is on the unit circle, hence is exp(i*θ). The following piece of code computes theta from omega. Note that this is a monotonically decreasing function, with range [0,2pi]. It maps the interval [-a,a] to [pi-b,pi+b], for some b; in particular b=pi/2 for the interval [-1,1]. Said otherwise, it maps a symmetric continuous-time interval to a symmetric discrete-time interval.

moebius(w):=
   let(t=#atan(w*2/(w^2-1)),
      if w<= -1 then t + %pi +%pi else if w<1 then t+%pi else t);

In the continous time case, a phase correction may be used. Here w is the frequency, w1 is the normalised frequency, and the phase correction is α*w1+β. These two quantities are zero by default, hence no correction is applied.

use_alpha_beta(val, w, alpha,beta,w0,dw) :=
 let(t=alpha*((w-w0)/dw)+beta,
   val*#exp(%%i*t));

The program starts by renormalising w, replacing w by a*(w+b). This produces a symmetric interval. The transformation might be w/w0 -w0/w, yielding a non-symmetric interval, case where make_symmetric will be called at the end. After that R(phi(s))=f*F(s) where the factor is 1/(1-z). This is also f=(1+iw)/2. In some cases, it may be useful to multiply by a power of f.

cvcircle1 (res, F,ow0, B):=
  let([n=size(F),t=0,l=0,f=0,w=0,k=B["inf_order"],r= B["real_transform"],
     alpha=B["alpha"],beta=B["beta"],nw0=B["nw0"],ndw=B["ndw"] ],
    l := n-1,
    for i=0...l do 
      let([ theta =  F["theta",i], val = F[i]],
          w:=theta,
          w:= B["real_coef"]*(if r then w/ow0 -ow0/w else w-ow0),
          if k>0
          then { f := (1+ %I*w)/2, 
                 for j=1...k do val := val*f },
          val:=use_alpha_beta(val,theta,B),
          res[l-i] := val,
          res["theta",l-i]:=moebius(w)),
    if r then #:pw:make_symmetric(res, B["shrink_circle"]));