1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
|
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA - Serge Steer
// Copyright (C) ENPC -
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
function [best_h,best_w]=frep2tf(frq,repf,dg,dom,tols,weight)
// iterative use of frep2tf_b jpc fd 1997
// Copyright INRIA
[lhs,rhs]=argn(0);
if rhs <= 3 then dom="c" ; end
if rhs <= 4 then
rtol=1.e-2; atol=1.e-4, N=10;
else
rtol=tols(1);atol=tols(2);N=tols(3);
end
if dom==[] then dom="c";end
if dom=="d" then dom=1;end
if rhs <=5 then
[h,err]=frep2tf_b(frq,repf,dg,dom);
best_w = [];
else
[h,err]=frep2tf_b(frq,repf,dg,dom,weight);
best_w = weight;
end
best_h = h ;
for k=1:N
if dom=="c" then
// weight=(1)./abs(freq(h('den'),1,%i*frq*2*%pi));
weight=(1)./horner(h("den"),%i*frq*2*%pi);
else
weight=(1)./horner(h("den"),exp(dom*%i*frq*2*%pi));
end
[h,err1]=frep2tf_b(frq,repf,dg,dom,weight);
if ( (abs(err-err1) < rtol *err & err > err1 )| err1 < atol) then break;end;
if err1 < err then best_err = err1 ; best_h = h; end
err=err1;
mprintf(gettext("%s: Iteration %s, error=%s.\n"), "frep2tf", part(string(k+1),1:5), string(err1));
end
endfunction
function [h,err]=frep2tf_b(frq,repf,dg,dom,weight)
// steer, jpc, fd 1997 (Nov)
//============================
[lhs,rhs]=argn(0);
// test the system type 'c' 'd' or dt
if rhs <= 3 then dom="c" ; end
if rhs <= 4 then weight=[];end
if dom==[] then dom="c";end
if dom=="d" then dom=1;end
n=size(frq,"*");
if dom=="c" then
w=2*%i*%pi*matrix(frq,n,1);
else
w=exp(2*%i*%pi*dom*matrix(frq,n,1));
end
//initialization
m=2*dg
//We compute the linear system to be solved:
//w(k)=%i* frq(k)*2pi
//for k=1,n sum(a_i*(w(k))^i,i=1,dg)
// -repf(k)*sum(b_i*(w(k))^i,i=1,dg) = 0
//with sum x_i = 1
//building Van der monde matrix ( w_i^j ) i=1,n j=0:dg-1
a1=w.*.[ones(1,dg)];
//0.^0 is not accepted in Scilab....
a1=[ones(n,1),a1.^(ones(n,1).*.[1:(dg)])];
a2=a1; for k=1:n; a2(k,:)= -repf(k)*a2(k,:);end
a=[a1,a2];
// Computing constraints
// We impose N(i wk) - repfk D(i wk) =0 for k=imax
// as follows:
// N(i wk) = repfk*(1+%i*b)
// D(i wk) = 1+%i*b
// L*[x;b]=[repfk;1]
// Least squ. pb is min norm of [A,0] [x;b]
// under constraint L*[x;b]=[repfk;1]
[rmax,imax]=max(abs(repf))
L2=a(imax,1:dg+1);
L=[zeros(L2),L2,%i;
L2,zeros(L2),repf(imax)*%i];
BigL=[real(L);imag(L)]
c=[1;repf(imax)];
Bigc=[real(c);imag(c)];
[ww,dim]=rowcomp(BigL);
BigL=ww*BigL;Bigc=ww*Bigc;
BigL=BigL(1:dim,:);Bigc=Bigc(1:dim,:);
a=[a,zeros(size(a,1),1)];
// auto renormalization : if weight is not given
if dom == "c" then
if weight==[] then
nn= sqrt(sum(abs(a).^2,"c"))+ones(n,1);
a=a./(nn*ones(1,size(a,2)));
end
end
// user given renormalization
if weight<>[] then
if size(frq,"*")<>size(weight,"*") then
error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Same numbers of elements expected.\n"),"frep2tf",1,5))
end
w1=weight(:)*ones(1,size(a,2));
a= w1.*a;
end
BigA=[real(a);imag(a)];
// Constraints BigL x =Bigc
//
x=LSC(BigA,BigL,Bigc);
x=x(1:$-1);
h=syslin(dom,poly(x(1:dg+1),"s","c"),poly([x((dg+2):$)],"s","c"))
if lhs==2 then
repf1=repfreq(h,frq);
err = sum(abs(repf1(:)-repf(:)))/n;
end
endfunction
function [x]=LSC(A,L,c)
// Ax=0 Least sq. + Lx = c
[W,rk]=colcomp(L);
LW=L*W;
Anew=A*W
A1=Anew(:,1:($-rk))
A2=Anew(:,($-rk+1:$));
x2=inv(LW(:,$-rk+1:$))*c
b= -A2*x2
x1=A1\b
x=W*[x1;x2]
endfunction
|
