Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/lp2z.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 % LP2Z converts a continous TF in to a discrete TF.
2 %
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 %
5 % DESCRIPTION: LP2Z converts a continous (lapalce notation) transfer
6 % function in a discrete (z transfor formalism) partial fractioned transfer
7 % function.
8 %
9 % Input function can be in rational, poles and zeros or partial fractions
10 % form:
11 % The rational case assumes an input function of the type:
12 %
13 % a(1)s^m + a(2)s^{m-1} + ... + a(m+1)
14 % H(s) = --------------------------------------
15 % b(1)s^n + b(2)s^{n-1} + ... + b(n+1)
16 %
17 % The poles and zeros case assumes an impution function of the type:
18 %
19 % (s-Z(1))(s-Z(2))...(s-Z(n))
20 % H(s) = K ---------------------------
21 % (s-P(1))(s-P(2))...(s-P(n))
22 %
23 % The partail fraction case assumes an input function of the type:
24 %
25 % R(1) R(2) R(n)
26 % H(s) = -------- + -------- + ... + -------- + K(s)
27 % s - P(1) s - P(2) s - P(n)
28 %
29 % LP2Z function is also capable to handling input transfer functions with
30 % multiple poles.
31 %
32 %
33 %
34 % CALL:
35 % pfstruct = lp2z(varargin)
36 % [res,poles,dterm] = lp2z(varargin)
37 % [pfstruct,res,poles,dterm] = lp2z(varargin)
38 %
39 %
40 %
41 % INPUTS:
42 % Input parameters are:
43 % 'FS' sampling frequency
44 % 'RRTOL' the repeated-root tolerance, default value is
45 % 1e-15. If two roots differs less than rrtolerance value, they
46 % are reported as multiple roots.
47 % 'INOPT' define the input function type
48 % 'PF' use this option if you want to input a function
49 % expanded in partial fractions. Then you have to input the
50 % vectors containing:
51 % 'RES' the vector containig the residues
52 % 'POLES' The vector containing the poles
53 % 'DTERMS' The vector containing the direct terms
54 % 'RAT' input the continous function in rational form.
55 % then you have to input the vector of coefficients:
56 % 'NUM' is the vector with numerator coefficients.
57 % 'DEN' is the vector of denominator coefficienets.
58 % 'PZ' input the continuous function in poles and
59 % zeros form. Then you have to input the vectors with poles
60 % and zeros:
61 % 'POLES' the vector with poles
62 % 'ZEROS' the vector with zeros
63 % 'GAIN' the value of the gain
64 % 'MODE' Is the mode used for the calculation of the roots of a
65 % polynomial. Admitted values are:
66 % 'SYM' uses symbolic roots calculation (you need Symbolic
67 % Math Toolbox to use this option)
68 % 'DBL' uses the standard numerical matlab style roots
69 % claculation (double precision)
70 %
71 %
72 % OUTPUTS:
73 % Function output is a structure array with two fields: 'num' and
74 % 'den' respectively. Each term of the partial fraction expansion
75 % can be seen as rational transfer function in z domain. These
76 % reduced transfer functions can be used to filter in parallel a
77 % series of data. This representation is particularly useful when
78 % poles with multiplicity higher than one are present. In this
79 % case the corresponding terms of the partial fraction expansion
80 % are really rational transfer functions whose order depends from
81 % the pole multiplicity [1].
82 % As a sencond option residues, poles and direct term are output
83 % as standard vectors (this possibility works only with simple
84 % poles).
85 % The third possibility is to call the function with four output,
86 % 1st is a struct, 2nd are residues, 3rd are ples and 4th is the
87 % direct term.
88 %
89 % NOTE1: The function is capable to convert in partial fractions
90 % only rational functions with numerator of order lower or equal
91 % to the order of the denominator. When the order of the
92 % numerator is higher than the order of the denominator (improper
93 % rational functions) the expansion in partial fractions is not
94 % useful and other methods are preferable.
95 %
96 % NOTE2: 'SYM' calculation mode requires Symbolic Math Toolbox
97 % to work
98 %
99 %
100 %
101 % REFERENCES:
102 % [1] Alan V. Oppenheim, Allan S. Willsky and Ian T. Young, Signals
103 % and Systems, Prentice-Hall Signal Processing Series, Prentice
104 % Hall (June 1982), ISBN-10: 0138097313. Pages 767 - 776.
105 %
106 %
107 % VERSION $Id: lp2z.m,v 1.11 2009/02/10 12:14:29 luigi Exp $
108 %
109 % HISTORY: 16-04-2008 L Ferraioli
110 % Creation
111 %
112 %
113 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
114
115 function varargout = lp2z(varargin)
116
117
118 %% Finding input parameters
119
120 % Default values for the parameters
121 fs = 10;
122 inopt = 'RAT';
123 tol = 1e-15;
124 mode = 'DBL';
125
126 % Finding input parameters
127 if ~isempty(varargin)
128 for j=1:length(varargin)
129 if strcmp(varargin{j},'INOPT')
130 inopt = varargin{j+1};
131 end
132 if strcmp(varargin{j},'FS')
133 fs = varargin{j+1};
134 end
135 if strcmp(varargin{j},'MODE')
136 mode = varargin{j+1};
137 end
138 if strcmp(varargin{j},'RRTOL')
139 tol = varargin{j+1};
140 end
141 end
142 end
143
144
145 %% Conversion to partial fractions and discretization
146
147 % swithing between input options
148 switch inopt
149 case 'RAT'
150 % Finding proper parameters
151 for jj=1:length(varargin)
152 if strcmp(varargin{jj},'NUM')
153 num = varargin{jj+1};
154 end
155 if strcmp(varargin{jj},'DEN')
156 den = varargin{jj+1};
157 end
158 end
159 % conversion in partail fractions cP is the vector of residues, cP is the
160 % vector of poles and cK is the vector of direct terms
161 [cR, cP, cK, mul] = utils.math.cpf('INOPT', 'RAT', 'NUM', num,...
162 'DEN', den, 'MODE', mode, 'RRTOL', tol);
163
164 case 'PZ'
165 % Finding proper parameters
166 for jj=1:length(varargin)
167 if strcmp(varargin{jj},'ZEROS')
168 zer = varargin{jj+1};
169 end
170 if strcmp(varargin{jj},'POLES')
171 pol = varargin{jj+1};
172 end
173 if strcmp(varargin{jj},'GAIN')
174 gain = varargin{jj+1};
175 end
176 end
177 % conversion in partail fractions cP is the vector of residues, cP is the
178 % vector of poles and cK is the vector of direct terms
179 [cR, cP, cK, mul] = utils.math.cpf('INOPT', 'PZ', 'ZEROS', zer,...
180 'POLES', pol, 'GAIN', gain, 'RRTOL', tol);
181
182 case 'PF'
183 % Finding proper parameters
184 for jj=1:length(varargin)
185 if strcmp(varargin{jj},'RES')
186 cR = varargin{jj+1};
187 end
188 if strcmp(varargin{jj},'POLES')
189 cP = varargin{jj+1};
190 end
191 if strcmp(varargin{jj},'DTERMS')
192 cK = varargin{jj+1};
193 end
194 end
195 % finding poles multiplicity
196 mul = mpoles(cP,tol,0);
197 end
198
199 % Checking for not proper functions
200 if length(cK)>1
201 error('Unable to directly discretize not proper continuous functions ')
202 end
203
204 % disctretization
205 T = 1/fs;
206
207 dpls = exp(cP.*T); % poles discretization
208
209 % Construct a struct array in which the elements of the partial
210 % fractions expansion are stored. Each partial fraction can be
211 % considered as rational function with proper numerator and denominator
212 % coefficients
213 pfstruct = struct();
214
215 indxmul = 0;
216 for kk=1:length(mul)
217 if mul(kk)==1
218 pfstruct(kk).num = cR(kk).*T;
219 pfstruct(kk).den = [1 -1*dpls(kk)];
220 else
221 coeffs = PolyLogCoeffs(1-mul(kk),exp(cP(kk)*T));
222 coeffs = flipdim(coeffs,2);
223 pfstruct(kk).num = (cR(kk)*T^mul(kk)/prod(1:mul(kk)-1)).*coeffs;
224 pfstruct(kk).den = BinomialExpansion(mul(kk),exp(cP(kk)*T));
225 indxmul = 1;
226 end
227 end
228
229 % Inserting the coefficients relative to the presence of a direct term
230 if ((length(cK)==1)&&(not(cK(1)==0)))
231 pfstruct(length(mul)+1).num = cK.*T;
232 pfstruct(length(mul)+1).den = [1 0];
233 end
234
235 % Output also the vector of residues, poles and direct term if non multiple
236 % poles are found
237 if indxmul == 0
238 res = cR.*T;
239 poles = dpls;
240 dterm = cK.*T;
241 else
242 res = [];
243 poles = [];
244 dterm = [];
245 end
246
247 % Output empty fields struct in case of no input
248 if isempty(cR)||isempty(cP)
249 pfstruct.num = [];
250 pfstruct.den = [];
251 end
252
253 %% Output data
254
255 if indxmul == 1
256 disp( ' Found poles with multiplicity higher than one. Results are output in struct form only. ')
257 end
258
259 if nargout == 1
260 varargout{1} = pfstruct;
261 elseif nargout == 3
262 varargout{1} = res;
263 varargout{2} = poles;
264 varargout{3} = dterm;
265 elseif nargout == 4
266 varargout{1} = pfstruct;
267 varargout{2} = res;
268 varargout{3} = poles;
269 varargout{4} = dterm;
270 else
271 error('Unespected number of outputs! Output must be 1, 3 or 4')
272 end
273
274
275
276 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
277 %
278 % FUNCTION: Eulerian Number
279 %
280 % DESCRIPTION: Calculates the Eulerian Number
281 %
282 % HISTORY: 12-05-2008 L Ferraioli
283 % Creation
284 %
285 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
286 function enum = eulerian(n,k)
287
288 enum = 0;
289 for jj = 0:k+1;
290 enum = enum + ((-1)^jj)*(prod(1:n+1)/(prod(1:n+1-jj)*prod(1:jj)))*(k-jj+1)^n;
291 end
292
293 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
294 %
295 % FUNCTION: PloyLogCoeffs
296 %
297 % DESCRIPTION: Computes the coefficients of a poly-logarithm expansion for
298 % negative n values
299 %
300 % HISTORY: 12-05-2008 L Ferraioli
301 % Creation
302 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
303 function coeffs = PolyLogCoeffs(n,x)
304
305 if n>=0
306 error('n must be a negative integer!')
307 else
308 % Now we keep the absolute value of n
309 n = -1*n;
310 end
311
312 coeffs = [];
313 for ii=0:n;
314 coeffs = [coeffs eulerian(n,ii)*(x^(n-ii))];
315 end
316
317 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
318 %
319 % FUNCTION: BinomialExpansion
320 %
321 % DESCRIPTION: Performs the binomial expamsion of the expression (1-x)^n
322 %
323 % HISTORY: 12-05-2008 L Ferraioli
324 % Creation
325 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
326 function binexp = BinomialExpansion(n,x)
327
328 binexp = [];
329 for jj=0:n
330 binexp = [binexp ((-1)^jj)*(prod(1:n)/(prod(1:n-jj)*prod(1:jj))*(x^jj))];
331 end