comparison m-toolbox/classes/+utils/@math/ctfit.m @ 0:f0afece42f48

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 % CTFIT fits a continuous model to a frequency response.
2 %
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 % DESCRIPTION:
5 %
6 % Fits a continuous model to a frequency response using relaxed
7 % s-domain vector fitting algorithm [1, 2]. Model function is expanded
8 % in partial fractions:
9 %
10 % r1 rN
11 % f(s) = ------- + ... + ------- + d
12 % s - p1 s - pN
13 %
14 % CALL:
15 %
16 % [res,poles,dterm,mresp,rdl] = ctfit(y,f,poles,weight,fitin)
17 %
18 % INPUTS:
19 %
20 % - y: Is a vector wuth the frequency response data.
21 % - f: Is the frequency vector in Hz.
22 % - poles: are a set of starting poles.
23 % - weight: are a set of weights used in the fitting procedure.
24 % - fitin: is a struct containing fitting options and parameters. fitin
25 % fields are:
26 % - fitin.stable = 0; fit without forcing poles to be stable.
27 % - fitin.stable = 1; force poles to be stable flipping unstable
28 % poles in the left side of the complex plane. s -> s - 2*conj(s).
29 % - fitin.dterm = 0; fit with d = 0.
30 % - fitin.dterm = 1; fit with d different from 0.
31 % - fitin.polt = 0; fit without plotting results.
32 % - fitin.plot = 1; plot fit results.
33 %
34 % OUTPUT:
35 %
36 % - res: vector or residues.
37 % - poles: vector of poles.
38 % - dterm: direct term d.
39 % - mresp: frequency response of the fitted model
40 % - rdl: residuals y - mresp
41 %
42 % REFERENCES:
43 %
44 % [1] B. Gustavsen and A. Semlyen, "Rational approximation of frequency
45 % domain responses by Vector Fitting", IEEE Trans. Power Delivery
46 % vol. 14, no. 3, pp. 1052-1061, July 1999.
47 % [2] B. Gustavsen, "Improving the Pole Relocating Properties of Vector
48 % Fitting", IEEE Trans. Power Delivery vol. 21, no. 3, pp.
49 % 1587-1592, July 2006.
50 %
51 % NOTE:
52 %
53 % This function cannot handle more than one frequency response per time
54 %
55 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
56 % VERSION: $Id: ctfit.m,v 1.2 2008/10/24 06:19:23 hewitson Exp $
57 %
58 % HISTORY: 12-09-2008 L Ferraioli
59 % Creation
60 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
61
62 function [res,poles,dterm,mresp,rdl] = ctfit(y,f,poles,weight,fitin)
63
64 %% Collecting inputs
65
66 % Default input struct
67 defaultparams = struct('stable',0, 'dterm',0, 'plot',0);
68
69 names = {'stable','dterm','plot'};
70
71 % collecting input and default params
72 if ~isempty(fitin)
73 for jj=1:length(names)
74 if isfield(fitin, names(jj))
75 defaultparams.(names{1,jj}) = fitin.(names{1,jj});
76 end
77 end
78 end
79
80 stab = defaultparams.stable; % Enforce pole stability is is 1
81 dt = defaultparams.dterm; % 1 to fit with direct term
82 plotting = defaultparams.plot; % set to 1 if plotting is required
83
84 %% Inputs in column vectors
85
86 [a,b] = size(y);
87 if a < b % shifting to column
88 y = y.';
89 end
90
91 [a,b] = size(f);
92 if a < b % shifting to column
93 f = f.';
94 end
95
96 [a,b] = size(poles);
97 if a < b % shifting to column
98 poles = poles.';
99 end
100
101 clear w
102 w = weight;
103 [a,b] = size(w);
104 if a < b % shifting to column
105 w = w.';
106 end
107
108 N = length(poles); % Model order
109
110 if dt
111 dl = 1; % Fit with direct term
112 else
113 dl = 0; % Fit without direct term
114 end
115
116 % definition of s
117 s = 1i.*2.*pi.*f;
118
119 Nz = length(s);
120
121 %% Marking complex and real poles
122
123 % cindex = 1; pole is complex, next conjugate pole is marked with cindex
124 % = 2. cindex = 0; pole is real
125 cindex=zeros(N,1);
126 for m=1:N
127 if imag(poles(m))~=0
128 if m==1
129 cindex(m)=1;
130 else
131 if cindex(m-1)==0 || cindex(m-1)==2
132 cindex(m)=1; cindex(m+1)=2;
133 else
134 cindex(m)=2;
135 end
136 end
137 end
138 end
139
140 %% Initializing the augmented problem matrices
141
142
143
144 % Matrix initialinzation
145 BA = zeros(Nz+1,1);
146 AA = zeros(Nz+1,2*N+dl+1);
147 Ak=zeros(Nz,N+1);
148
149 % Defining Ak
150 for jj = 1:N
151 if cindex(jj) == 1 % conjugate complex couple of poles
152 Ak(:,jj) = (1./(s-poles(jj)))+(1./(s-poles(jj+1)));
153 Ak(:,jj+1) = (1i./(s-poles(jj)))-(1i./(s-poles(jj+1)));
154 elseif cindex(jj) == 0 % real pole
155 Ak(:,jj) = 1./(s-poles(jj));
156 end
157 end
158
159
160 Ak(1:Nz,N+1) = ones(Nz,1);
161
162
163 for m=1:N+dl % left columns
164 AA(1:Nz,m)=w.*Ak(1:Nz,m);
165 end
166 for m=1:N+1 %Rightmost blocks
167 AA(1:Nz,N+dt+m)=-w.*(Ak(1:Nz,m).*y);
168 end
169
170 % Scaling factor
171 clear sc
172 sc = norm(w.*y)/Nz;
173
174 % setting the last row of AA and BA for the relaxion condition
175 for qq = 1:N+1
176 AA(Nz+1,N+dl+qq) = real(sc*sum(Ak(:,qq)));
177 end
178
179 AA = [real(AA);imag(AA)];
180
181 % AAstr1 = AA; % storing
182
183 % Last element of the solution vector
184 BA(Nz+1) = Nz*sc;
185
186 % solving for real and imaginary part of the solution vector
187 nBA = [real(BA);imag(BA)];
188
189 % Normalization factor
190 nf = zeros(2*N+dl+1,1);
191 for pp = 1:2*N+dl+1
192 nf(pp,1) = norm(AA(:,pp),2); % Euclidean norm
193 AA(:,pp) = AA(:,pp)./nf(pp,1); % Normalization
194 end
195
196
197 %% Solving augmented problem
198
199 % XA = pinv(AA)*nBA;
200 % XA = inv((AA.')*AA)*(AA.')*nBA;
201
202 % XA = AA.'*AA\AA.'*nBA;
203
204 XA = AA\nBA;
205
206 XA = XA./nf; % renormalization
207
208 %% Finding zeros of sigma
209
210 lsr = XA(N+dl+1:2*N+dl,1); % collect the least square results
211
212 Ds = XA(end); % direct term of sigma
213
214 % Real poles have real residues, complex poles have comples residues
215 rs = zeros(N,1);
216 for tt = 1:N
217 if cindex(tt) == 1 % conjugate complex couple of poles
218 rs(tt,1) = lsr(tt)+1i*lsr(tt+1);
219 rs(tt+1,1) = lsr(tt)-1i*lsr(tt+1);
220 elseif cindex(tt) == 0 % real pole
221 rs(tt,1) = lsr(tt);
222 end
223 end
224
225 % [snum, sden] = residue(rs,poles,Ds);
226 %
227 % % ceking for numerical calculation errors
228 % for jj = 1:length(snum)
229 % if ~isequal(imag(snum(jj)),0)
230 % snum(jj)=real(snum(jj));
231 % end
232 % end
233 %
234 % % Zeros of sigma are poles of F
235 % szeros = roots(snum);
236
237 DPOLES = diag(poles);
238 B = ones(N,1);
239 C = rs.';
240 for kk = 1:N
241 if cindex(kk) == 1
242 DPOLES(kk,kk)=real(DPOLES(kk,kk));
243 DPOLES(kk,kk+1)=imag(DPOLES(kk,kk));
244 DPOLES(kk+1,kk)=-1*imag(DPOLES(kk,kk));
245 DPOLES(kk+1,kk+1)=real(DPOLES(kk,kk));
246 B(kk,1) = 2;
247 B(kk+1,1) = 0;
248 C(1,kk) = real(C(1,kk));
249 C(1,kk+1) = imag(C(1,kk));
250 end
251 end
252
253 H = DPOLES-B*C/Ds;
254 szeros = eig(H);
255
256 %% Ruling out unstable poles
257
258 % This option force the poles of f to stay inside the left side of the
259 % complex plane
260
261 if stab
262 unst = real(szeros)>0;
263 szeros(unst) = szeros(unst)-2*real(szeros(unst)); % Mirroring respect to the complex axes
264 end
265 N = length(szeros);
266
267 %% Separating complex poles from real poles and ordering
268
269 rnpoles = [];
270 inpoles = [];
271 for tt = 1:N
272 if imag(szeros(tt)) == 0
273 % collecting real poles
274 rnpoles = [rnpoles; szeros(tt)];
275 else
276 % collecting complex poles
277 inpoles = [inpoles; szeros(tt)];
278 end
279 end
280
281 % Sorting complex poles in order to have them in the expected order a+jb
282 % and a-jb a>0 b>0
283 inpoles = sort(inpoles);
284 npoles = [rnpoles;inpoles];
285 npoles = npoles - 2.*1i.*imag(npoles);
286
287 %% Marking complex and real poles
288
289 cindex=zeros(N,1);
290 for m=1:N
291 if imag(npoles(m))~=0
292 if m==1
293 cindex(m)=1;
294 else
295 if cindex(m-1)==0 || cindex(m-1)==2
296 cindex(m)=1; cindex(m+1)=2;
297 else
298 cindex(m)=2;
299 end
300 end
301 end
302 end
303
304 %% Initializing direct problem
305
306 % Matrix initialinzation
307 B = w.*y;
308 AD = zeros(Nz,N+dl);
309 Ak=zeros(Nz,N+dl);
310
311 % Defining Ak
312 for jj = 1:N
313 if cindex(jj) == 1 % conjugate complex couple of poles
314 Ak(:,jj) = (1./(s-npoles(jj)))+(1./(s-npoles(jj+1)));
315 Ak(:,jj+1) = (1i./(s-npoles(jj)))-(1i./(s-npoles(jj+1)));
316 elseif cindex(jj) == 0 % real pole
317 Ak(:,jj) = 1./(s-npoles(jj));
318 end
319 end
320
321 if dt
322 Ak(1:Nz,N+dl) = ones(Nz,1); % considering the direct term
323 end
324
325 % Defining AD
326 for m=1:N+dl
327 AD(1:Nz,m)=w.*Ak(1:Nz,m);
328 end
329
330
331 AD = [real(AD);imag(AD)];
332 nB = [real(B);imag(B)];
333
334 % Normalization factor
335 nf = zeros(N+dl,1);
336 for pp = 1:N+dl
337 nf(pp,1) = norm(AD(:,pp),2); % Euclidean norm
338 AD(:,pp) = AD(:,pp)./nf(pp,1); % Normalization
339 end
340
341 %% Solving direct problem
342
343 % XD = inv((AD.')*AD)*(AD.')*nB;
344 % XD = AD.'*AD\AD.'*nB;
345 % XD = pinv(AD)*nB;
346 XD = AD\nB;
347
348 XD = XD./nf; % Renormalization
349
350 %% Final residues and poles of f
351
352 if dt
353 lsr = XD(1:end-1); % Fitting with direct term
354 else
355 lsr = XD(1:end); % Fitting without direct term
356 end
357
358 res = zeros(N,1);
359 % Real poles have real residues, complex poles have comples residues
360 for tt = 1:N
361 if cindex(tt) == 1 % conjugate complex couple of poles
362 res(tt) = lsr(tt)+1i*lsr(tt+1);
363 res(tt+1) = lsr(tt)-1i*lsr(tt+1);
364 elseif cindex(tt) == 0 % real pole
365 res(tt) = lsr(tt);
366 end
367 end
368
369 clear poles
370 poles = npoles;
371
372 if dt
373 dterm = XD(end);
374 else
375 dterm = 0;
376 end
377
378 %% freq resp of the fit model
379
380 % f = pfparams.freq;
381 r = res;
382 p = poles;
383 d = dterm;
384
385 Nf = length(f);
386 N = length(p);
387
388 rsp = zeros(Nf,1);
389 indx = (0:length(d)-1).';
390 for ii = 1:Nf
391 for jj = 1:N
392 rsptemp = r(jj)/(1i*2*pi*f(ii)-p(jj));
393 rsp(ii) = rsp(ii) + rsptemp;
394 end
395 % Direct terms response
396 rsp(ii) = rsp(ii) + sum(((1i*2*pi*f(ii))*ones(length(d),1).^indx).*d);
397 end
398
399 mresp = rsp;
400
401 % Residual
402 rdl = y - mresp;
403
404 %% Plotting response
405
406 if plotting
407 figure(1)
408 subplot(2,1,1);
409 loglog(f,abs(y),'k')
410 hold on
411 loglog(f,abs(mresp),'r')
412 loglog(f,abs(rdl),'b')
413 xlabel('Frequency [Hz]')
414 ylabel('Amplitude')
415 legend('Original', 'CTFIT','Residual')
416 hold off
417
418 subplot(2,1,2);
419 semilogx(f,angle(y),'k')
420 hold on
421 semilogx(f,angle(mresp),'r')
422 xlabel('Frequency [Hz]')
423 ylabel('Phase [Rad]')
424 legend('Original', 'CTFIT')
425 end
426 hold off
427 end