Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/cpf.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 % CPF finds the partial fraction expansion of the ratio of two polynomials A(s)/B(s).
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: CPF finds the residues, poles and direct terms of the
5 % partial fraction expansion of the ratio of two polynomials A(s)/B(s).
6 % This function assumes that the input continous filter is written in the
7 % rational form or in poles, zeros and gain factorization:
8 %
9 % A(s) a(1)s^m + a(2)s^{m-1} + ... + a(m+1)
10 % H(s)= ---- = --------------------------------------
11 % B(s) b(1)s^n + b(2)s^{n-1} + ... + b(n+1)
12 %
13 % or
14 %
15 % A(s) (s-z1)...(s-zn)
16 % H(s)= ---- = g ----------------
17 % B(s) (s-p1)...(s-pn)
18 %
19 %
20 % It inputs a plist containing the coefficients vectors and the
21 % repeated-root tolerance.
22 % Eg:
23 % A = [a(1), a(2), ..., a(m+1)]
24 % B = [b(1), b(2), ..., b(m+1)]
25 %
26 % or
27 %
28 % Z = [z(1), z(2), ..., z(m)]
29 % P = [p(1), p(2), ..., p(m)]
30 % G = g (Gain is a scalar)
31 %
32 %
33 % If there are no multiple roots,
34 %
35 % A(s) R(1) R(2) R(n)
36 % ---- = -------- + -------- + ... + -------- + K(s)
37 % B(s) s - P(1) s - P(2) s - P(n)
38 %
39 % The number of poles is n = length(B)-1 = length(R) = length(P).
40 % The direct term coefficient vector is empty if length(A) < length(B),
41 % otherwise length(K) = length(A)-length(B)+1.
42 % K(s) is returned in the form:
43 %
44 % K(s) = k(1)*s^(m-n) + ... + k(m-n)*s + k(m-n+1)
45 %
46 % so that the output vector of direct terms is:
47 % K = [k(1), ..., k(m-n), k(m-n+1)]
48 %
49 % If P(j) = ... = P(j+m-1) is a pole of multplicity m, then the
50 % expansion includes terms of the form
51 % R(j) R(j+1) R(j+m-1)
52 % -------- + ------------ + ... + ------------
53 % s - P(j) (s - P(j))^2 (s - P(j))^m
54 %
55 % The function is also capable to convert a partial fraction expanded
56 % function to its rational form by setting the 'PARFRACT' input option.
57 % In this case the output is composed by a plist containing the vactors
58 % of numerator and denominator polynomial coefficients.
59 %
60 %
61 %
62 % CALL: varargout = cpf(varargin)
63 %
64 %
65 %
66 % INPUTS:
67 % Input options are:
68 % 'INOPT' define the input function type
69 % 'RAT' input the continous function in rational form.
70 % then you have to input the vector of coefficients:
71 % 'NUM' is the vector with numerator coefficients.
72 % 'DEN' is the vector of denominator coefficienets.
73 % 'PZ' input the continuous function in poles and
74 % zeros form. Then you have to input the vectors with poles
75 % and zeros:
76 % 'POLES' the vector with poles
77 % 'ZEROS' the vector with zeros
78 % 'GAIN' the value of the gain
79 % 'PF' input the coefficients of a partial fraction
80 % expansion of the transfer function. When this option is
81 % setted the function performs the conversion from partial
82 % fraction to rational transfer function. You have to input the
83 % vectors containing the residues, poles and direct terms:
84 % 'RES' the vector with residues
85 % 'POLES' the vector with poles
86 % 'DTERMS' the vector with direct terms
87 % 'MODE' Is the used mode for the calculation of the roots of a
88 % polynomial. It is an useful option only with rational functions
89 % at the input. Admitted values are:
90 % 'SYM' uses symbolic roots calculation (you need symbolic
91 % math toolbox to use this option)
92 % 'DBL' uses the standard numerical matlab style roots
93 % claculation (double precision)
94 % 'RRTOL' the repeated-root tolerance default value is
95 % 1e-15. If two roots differs less than rrtolerance value, they
96 % are reported as multiple roots
97 %
98 %
99 %
100 % OUTPUTS:
101 %
102 % When 'INOPT' is set to 'RAT' or 'PZ', outputs
103 % are:
104 % RES vector of residues coefficients
105 % POLES vector of poles coefficients
106 % DTERMS vector of direct terms coefficients
107 % PMul vector of poles multiplicity
108 %
109 % When 'INOPT' is setted to 'PF', outputs are:
110 % NUM the vector with the numerator polynomial
111 % coefficients
112 % DEN the vector with the denominator polynomial
113 % coefficents
114 %
115 %
116 %
117 % NOTE:
118 % - 'SYM' option for 'MODE' requires the Symblic Math Toolbox. It
119 % is used only for rational function input
120 %
121 %
122 %
123 % EXAMPLES:
124 % - Input a function in rational form and output the partial
125 % fraction expansion
126 % [Res, Poles, DTerms, PMul] = cpf('INOPT', 'RAT',
127 % 'NUM', [], 'DEN', [], 'MODE','SYM', 'RRTOL', 1e-15)
128 % - Input a function in poles and zeros and output the partial
129 % fraction expansion
130 % [Res, Poles, DTerms, PMul] = cpf('INOPT', 'PZ',
131 % 'POLES', [], 'ZEROS', [], 'GAIN', #, 'RRTOL', 1e-15)
132 % - Input a function in partial fractions and output the rational
133 % expression
134 % [Num, Den] = cpf('INOPT', 'PF', 'POLES', [], 'RES',
135 % [], 'DTERMS', [], 'RRTOL', 1e-15)
136 %
137 %
138 %
139 % REFERENCES:
140 % [1] Alan V. Oppenheim, Allan S. Willsky and Ian T. Young, Signals
141 % and Systems, Prentice-Hall Signal Processing Series, Prentice
142 % Hall (June 1982), ISBN-10: 0138097313. Pages 767 - 776.
143 %
144 %
145 %
146 % VERSION: $Id: cpf.m,v 1.4 2008/11/19 11:37:51 luigi Exp $
147 %
148 %
149 % HISTORY: 16-04-2008 L Ferraioli
150 % Creation
151 %
152 %
153 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
154
155 function varargout = cpf(varargin)
156
157 %% VERSION
158
159 VERSION = '$Id: cpf.m,v 1.4 2008/11/19 11:37:51 luigi Exp $';
160
161 %% Extracting parameters
162
163 % default parameters
164 inopt = 'RAT';
165 mode = 'DBL';
166 tol = 1e-15;
167
168 % Finding input parameters
169 if ~isempty(varargin)
170 for j=1:length(varargin)
171 if strcmp(varargin{j},'INOPT')
172 inopt = varargin{j+1};
173 end
174 if strcmp(varargin{j},'MODE')
175 mode = varargin{j+1};
176 end
177 if strcmp(varargin{j},'RRTOL')
178 tol = varargin{j+1};
179 end
180 end
181 end
182
183 % Switching between input options and setup inputs for proper calculation
184 switch inopt;
185 case 'RAT'
186 % etracting numerator and denominator vectors
187 for jj=1:length(varargin)
188 if strcmp(varargin{jj},'NUM')
189 u = varargin{jj+1};
190 end
191 if strcmp(varargin{jj},'DEN')
192 v = varargin{jj+1};
193 end
194 end
195
196 % For the conversion we need the denominator factored in poles and the
197 % numerator in polynomial form
198 switch mode
199 case 'DBL'
200 % adopt numerical calculation
201 poles_vect = roots(v);
202 case 'SYM'
203 % adopt symbolic calculation
204 % syms s
205 % % Construct the symbolic polynomial
206 % numel = length(v);
207 % PowerVector = [];
208 % for ii=1:numel
209 % PowerVector = [PowerVector v(ii)*s^(numel-ii)];
210 % end
211 % PowerMatrix = diag(PowerVector);
212 % Polyv = trace(PowerMatrix);
213 % % Solve the polynomial in order to find the roots
214 % Sp = solve(Polyv,s);
215 % % output of the poles vector in Matlab double format
216 % numpoles = length(Sp);
217 % poles_vect = zeros(1,numpoles);
218 % for jj=1:numpoles
219 % poles_vect(jj) = double(Sp(jj));
220 % end
221
222 n = length(v);
223 cN = -1.*v(2:end)./v(1);
224 A = sym(diag(ones(1,n-2),-1));
225 A(1,:) = cN;
226 sol = eig(A);
227 poles_vect = double(sol);
228 end
229 % setting the output option to residues and poles
230 outopt = 'RP';
231
232 case 'PZ'
233 % Extracting zeros, poles and gain from inputs
234 for jj=1:length(varargin)
235 if strcmp(varargin{jj},'ZEROS')
236 zeros_vect = varargin{jj+1};
237 end
238 if strcmp(varargin{jj},'POLES')
239 poles_vect = varargin{jj+1};
240 end
241 if strcmp(varargin{jj},'GAIN')
242 gain = varargin{jj+1};
243 end
244 end
245
246 u = poly(zeros_vect).*gain;
247 v = poly(poles_vect);
248 if ((~isempty(v)) && (v(1)~=1))
249 u = u ./ v(1); v = v ./ v(1); % Normalize.
250 end
251 % setting the output option to residues and poles
252 outopt = 'RP';
253
254 case 'PF'
255 % Calculate numerator and denominator of a transfer function expanded
256 % in partial fractions
257 % etracting residues, poles and direct terms
258 for jj=1:length(varargin)
259 if strcmp(varargin{jj},'RES')
260 u = varargin{jj+1};
261 end
262 if strcmp(varargin{jj},'POLES')
263 v = varargin{jj+1};
264 end
265 if strcmp(varargin{jj},'DTERMS')
266 k = varargin{jj+1};
267 end
268 end
269 % setting the output option to Transfer Function
270 outopt = 'TF';
271 end
272
273
274
275 %% Partial fractions expansion
276
277 % Switching between output cases
278 % Note: rational input and poles zeros are equivalent from this point on,
279 % so PF expansion is calculated in the same way
280 switch outopt
281
282 case 'RP'
283 % Direct terms calculation
284 if length(u) >= length(v)
285 [dterms,new_u]=deconv(u,v);
286 else
287 dterms = 0;
288 new_u = u;
289 end
290
291 % identification of multiple poles
292 poles_vect = sort(poles_vect); % sort the poles in ascending order
293 mul = mpoles(poles_vect,tol,0); % find the multiplicity
294
295 mmul = mul;
296 for kk=1:length(mmul)
297 if mmul(kk)>1
298 for hh=1:mmul(kk)
299 mmul(kk-hh+1)=mmul(kk);
300 end
301 end
302 end
303
304 % finding the residues
305 resids = zeros(length(poles_vect),1);
306
307 for ii=1:length(poles_vect)
308
309 den = v;
310 p = [1 -poles_vect(ii)];
311 for hh=1:mmul(ii)
312 den = deconv(den,p);
313 end
314
315 dnum = new_u;
316 dden = den;
317
318 c = 1;
319 if mmul(ii)>mul(ii)
320 c = prod(1:(mmul(ii)-mul(ii)));
321
322 for jj=1:(mmul(ii)-mul(ii))
323 [dnum,dden] = polyder(dnum,dden);
324 end
325
326 end
327
328 resids(ii)=(polyval(dnum,poles_vect(ii))./polyval(dden,poles_vect(ii)))./c;
329 end
330
331 % Converting from partial fractions to rational function
332 case 'TF'
333 % This code is directly taken from matlab 'residue' function
334 [mults,i]=mpoles(v,tol,0);
335 p=v(i); r=u(i);
336 n = length(p);
337 q = [p(:).' ; mults(:).']; % Poles and multiplicities.
338 v = poly(p); u = zeros(1,n,class(u));
339 for indx = 1:n
340 ptemp = q(1,:);
341 i = indx;
342 for j = 1:q(2,indx), ptemp(i) = nan; i = i-1; end
343 ptemp = ptemp(find(~isnan(ptemp))); temp = poly(ptemp);
344 j = length(temp);
345 if j < n, temp = [zeros(1,n-j) temp]; end
346 u = u + (r(indx) .* temp);
347 end
348 if ~isempty(k)
349 if any(k ~= 0)
350 u = [zeros(1,length(k)) u];
351 k = k(:).';
352 temp = conv(k,v);
353 u = u + temp;
354 end
355 end
356 num = u; den = v; % Rename.
357 end
358
359 %% Output data
360
361 switch outopt
362 case 'RP'
363 if nargout == 1
364 varargout{1} = [resids poles_vect dterms mul];
365 elseif nargout == 2
366 varargout{1} = resids;
367 varargout{2} = poles_vect;
368 elseif nargout == 3
369 varargout{1} = resids;
370 varargout{2} = poles_vect;
371 varargout{3} = dterms;
372 elseif nargout == 4
373 varargout{1} = resids;
374 varargout{2} = poles_vect;
375 varargout{3} = dterms;
376 varargout{4} = mul;
377 else
378 error('Unespected number of outputs! Set 1, 2, 3 or 4')
379 end
380 % plout = plist('RESIDUES', resids, 'POLES', poles_vect, 'PMul', mul, 'DIRECT_TERMS', dterms);
381 % % plout = combine(plout, pl);
382 case 'TF'
383 if nargout == 1
384 varargout{1} = [num den];
385 elseif nargout == 2
386 varargout{1} = num;
387 varargout{2} = den;
388 else
389 error('Unespected number of outputs! Set 1 or 2')
390 end
391 % plout = plist('NUMERATOR', num, 'DENOMINATOR', den);
392 end
393
394 end
395 % END
396
397
398
399