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