Mercurial > hg > ltpda

comparison m-toolbox/classes/@ssm/doBode.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % DOBODE makes a bode computation from the given inputs to outputs.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: makes a bode computation from all the inputs to all outputs.
5 %
6 % CALL: varargout = doBode(a, b, c, d, w, Ts)
7 %
8 % INPUTS:
9 % 'sys' - ssm object
10 % 'inputnames, statenames, outputnames' - 3 cellstr
11 %
12 % OUTPUTS:
13 %
14 % 'as' - array of output TFs containing the requested responses.
15 %
16 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
17
18 function varargout = doBode(a, b, c, d, w, Ts)
19
20 f = w/(2*pi);
21
22 %% dealing with timestep
23 if(Ts ~= 0)
24 % Compute discrete freq:
25 f_disc = f*Ts;
26 % Compute z = e^jw (with w discrete)
27 z = exp(1i*2*pi*f_disc);
28 end
29
30 %% getting matrices properly ordered
31
32 %% stop warnings
33 s = warning;
34 warning('off','MATLAB:nearlySingularMatrix'); % turn all warnings off
35
36 %% looping over SISO systems
37 Ni = size(b,2);
38 No = size(c,1);
39 G = zeros(No,Ni,length(f));
40 G2 = zeros(No,Ni,length(f));
41 %% compute responses
42 % To conserve accuracy, the system is transformed to Hessenberg form [1]
43 % Laub, A.J., "Efficient Multivariable Frequency Response Computations",
44 % IEEE Transactions on Automatic Control, AC-26 (1981), pp. 407-408
45
46 % This transformation leads to a faster resolution and at the
47 % same time to a more accurate and precise answer of C*(jw - A)^-1 *B + D
48
49 reduce_a = a;
50 reduce_c = c;
51
52 %% loop on inputs
53 if (Ts == 0)
54 for ii=1:Ni
55 reduce_b = b(:,ii);
56 reduce_d = d(:,ii);
57 [num,den] = ss2tf(reduce_a,reduce_b,reduce_c,reduce_d);
58 denr = polyval(den, 2*pi*1i*f);
59 for jj = 1:No
60 numr = polyval(num(jj,:), 2*pi*1i*f);
61 G(jj,ii,:) = numr./denr;
62 end
63 end
64 else
65 [T,H]=hess(reduce_a);
66 % Step 1:
67 P = reduce_c*T;
68 reduce_b = b(:,1:Ni);
69 reduce_d = d(1:No,1:Ni);
70 Q = T'*reduce_b;
71 I = eye(size(reduce_a));
72
73 %% time discrete
74 for ff = 1:length(f)
75 G(1:No,1:Ni,ff) = (P/(z(ff)*I - H))*(Q) + reduce_d;
76 end
77 end
78
79 %% execute warnings
80 [msg, msgid] = lastwarn;
81 if(strcmp(msgid,'MATLAB:nearlySingularMatrix') == 1)
82 %warning('This frequency response may be unaccurate because the system matrix is close to singular or badly scaled.') %#ok<WNTAG>
83 end
84 warning(s) % restore the warning state
85
86 %% output
87 varargout = {G};
88 end
89
90