Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/vcfit.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 % VCFIT Fits continuous models to frequency responses
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % DESCRIPTION:
4 %
5 % Fits a continuous model to a frequency response using relaxed
6 % s-domain vector fitting algorithm [1, 2]. The function is able to fit more
7 % than one frequency response per time. In case that more than one
8 % frequency response is passed as input, they are fitted with a set of
9 % common poles. Model functions are expanded in partial fractions:
10 %
11 % r1 rN
12 % f(s) = ------- + ... + ------- + d
13 % s - p1 s - pN
14 %
15 % CALL:
16 %
17 % [res,poles,dterm,mresp,rdl] = vcfit(y,f,poles,weight,fitin)
18 %
19 % INPUTS:
20 %
21 % - y: Is a vector wuth the frequency response data.
22 % - f: Is the frequency vector in Hz.
23 % - poles: are a set of starting poles.
24 % - weight: are a set of weights used in the fitting procedure.
25 % - fitin: is a struct containing fitting options and parameters. fitin
26 % fields are:
27 % - fitin.stable = 0; fit without forcing poles to be stable.
28 % - fitin.stable = 1; force poles to be stable flipping unstable
29 % poles in the left side of the complex plane. s -> s - 2*conj(s).
30 % - fitin.dterm = 0; fit with d = 0.
31 % - fitin.dterm = 1; fit with d different from 0.
32 % - fitin.polt = 0; fit without plotting results.
33 % - fitin.plot = 1; plot fit results in loglog scale.
34 % - fitin.plot = 2; plot fit results in semilogx scale.
35 % - fitin.plot = 3; plot fit results in semilogy scale.
36 % - fitin.plot = 4; plot fit results in linear xy scale.
37 % - fitin.ploth = #; a plot handle to define the figure target for
38 % plotting. Default: [1]
39 %
40 % OUTPUT:
41 %
42 % - res: vector or residues.
43 % - poles: vector of poles.
44 % - dterm: direct term d.
45 % - mresp: frequency response of the fitted model
46 % - rdl: residuals y - mresp
47 % - mse: normalized mean squared error
48 %
49 % EXAMPLES:
50 %
51 % - Fit on a single transfer function:
52 %
53 % INPUT
54 % y is a (Nx1) or (1xN) vector
55 % f is a (Nx1) or (1xN) vector
56 % poles is a (Npx1) or (1xNp) vector
57 % weight is a (Nx1) or (1xN) vector
58 %
59 % OUTPUT
60 % res is a (Npx1) vector
61 % poles is a (Npx1) vector
62 % dterm is a constant
63 % mresp is a (Nx1) vector
64 % rdl is a (Nx1) vector
65 % mse is a constant
66 %
67 % - Fit Nt transfer function at the same time:
68 %
69 % INPUT
70 % y is a (NxNt) or (NtxN) vector
71 % f is a (Nx1) or (1xN) vector
72 % poles is a (Npx1) or (1xNp) vector
73 % weight is a (NxNt) or (NtxN) vector
74 %
75 % OUTPUT
76 % res is a (NpxNt) vector
77 % poles is a (Npx1) vector
78 % dterm is a (1xNt) vector
79 % mresp is a (NxNt) vector
80 % rdl is a (NxNt) vector
81 % mse is a (1xNt) vector
82 %
83 % REFERENCES:
84 %
85 % [1] B. Gustavsen and A. Semlyen, "Rational approximation of frequency
86 % domain responses by Vector Fitting", IEEE Trans. Power Delivery
87 % vol. 14, no. 3, pp. 1052-1061, July 1999.
88 % [2] B. Gustavsen, "Improving the Pole Relocating Properties of Vector
89 % Fitting", IEEE Trans. Power Delivery vol. 21, no. 3, pp.
90 % 1587-1592, July 2006.
91 %
92 % NOTE:
93 %
94 %
95 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
96 % HISTORY: 12-09-2008 L Ferraioli
97 % Creation
98 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
99 % VERSION: '$Id: vcfit.m,v 1.11 2010/02/08 19:05:59 luigi Exp $';
100 %
101 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
102 function [res,poles,dterm,mresp,rdl,mse] = vcfit(y,f,poles,weight,fitin)
103
104 warning off all
105 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
106 % Collecting inputs
107
108 % Default input struct
109 defaultparams = struct('stable',0, 'dterm',0, 'plot',0, 'ploth', 1);
110
111 names = {'stable','dterm','plot','ploth'};
112
113 % collecting input and default params
114 if ~isempty(fitin)
115 for jj=1:length(names)
116 if isfield(fitin, names(jj)) && ~isempty(fitin.(names{1,jj}))
117 defaultparams.(names{1,jj}) = fitin.(names{1,jj});
118 end
119 end
120 end
121
122 stab = defaultparams.stable; % Enforce pole stability is is 1
123 dt = defaultparams.dterm; % 1 to fit with direct term
124 plotting = defaultparams.plot; % set to 1 if plotting is required
125 plth = defaultparams.ploth; % set the figure target
126
127 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
128 % Inputs in row vectors
129
130 [a,b] = size(y);
131 if a > b % shifting to row
132 y = y.';
133 end
134
135 [a,b] = size(f);
136 if a > b % shifting to row
137 f = f.';
138 end
139
140 [a,b] = size(poles);
141 if a > b % shifting to row
142 poles = poles.';
143 end
144
145 clear w
146 w = weight;
147 [a,b] = size(w);
148 if a > b % shifting to row
149 w = w.';
150 end
151
152 N = length(poles); % Model order
153
154 if dt
155 dl = 1; % Fit with direct term
156 else
157 dl = 0; % Fit without direct term
158 end
159
160 % definition of s
161 s = 1i.*2.*pi.*f;
162
163 Nz = length(s);
164 [Nc,Ny] = size(y);
165 if Ny ~= Nz
166 error(' Number of data points different from number of frequency points! ')
167 end
168
169 %Tolerances used by relaxed version of vector fitting
170 TOLlow=1e-8;
171 TOLhigh=1e8;
172
173 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
174 % Marking complex and real poles
175
176 % cindex = 1; pole is complex, next conjugate pole is marked with cindex
177 % = 2. cindex = 0; pole is real
178 cindex=zeros(1,N);
179 for m=1:N
180 if imag(poles(m))~=0
181 if m==1
182 cindex(m)=1;
183 else
184 if cindex(m-1)==0 || cindex(m-1)==2
185 cindex(m)=1; cindex(m+1)=2;
186 else
187 cindex(m)=2;
188 end
189 end
190 end
191 end
192
193 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
194 % Augmented problem
195 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
196
197 % Matrix initialinzation
198 BA = zeros(Nc*Nz+1,1);
199 Ak=zeros(Nz,N+1);
200 AA=zeros(Nz*Nc+1,(N+dl)*Nc+N+1);
201 nf = zeros(1,Nc*(N+dl)+N+1); % Normalization factor
202
203 % Defining Ak
204 for jj = 1:N
205 if cindex(jj) == 1 % conjugate complex couple of poles
206 Ak(:,jj) = 1./(s-poles(jj)) + 1./(s-conj(poles(jj)));
207 Ak(:,jj+1) = 1i./(s-poles(jj)) - 1i./(s-conj(poles(jj)));
208 elseif cindex(jj) == 0 % real pole
209 Ak(:,jj) = 1./(s-poles(jj));
210 end
211 end
212
213 Ak(1:Nz,N+1) = 1;
214
215 % Scaling factor
216 sc = 0;
217 for mm = 1:Nc
218 sc = sc + (norm(w(mm,:).*y(mm,:)))^2;
219 end
220 sc=sqrt(sc)/Nz;
221
222 for nn = 1:Nc
223
224 wg = w(nn,:).'; % Weights
225
226 ida=(nn-1)*Nz+1;
227 idb=nn*Nz;
228 idc=(nn-1)*(N+dl)+1;
229
230 for mm =1:N+dl % Diagonal blocks
231 AA(ida:idb,idc-1+mm) = wg.*Ak(1:Nz,mm);
232 end
233 for mm =1:N+1 % Last right blocks
234 AA(ida:idb,Nc*(N+dl)+mm) = -wg.*(Ak(1:Nz,mm).*y(nn,1:Nz).');
235 end
236
237 end
238
239 % setting the last row of AA and BA for the relaxion condition
240 for qq = 1:N+1
241 AA(Nc*Nz+1,Nc*(N+dl)+qq) = real(sc*sum(Ak(:,qq)));
242 end
243
244 AA = [real(AA);imag(AA)];
245
246 % Last element of the solution vector
247 BA(Nc*Nz+1) = Nz*sc;
248
249 xBA = real(BA);
250 xxBA = imag(BA);
251
252 Nrow = Nz*Nc+1;
253
254 BA = zeros(2*Nrow,1);
255
256 BA(1:Nrow,1) = xBA;
257 BA(Nrow+1:2*Nrow,1) = xxBA;
258
259 % Normalization factor
260 % nf = zeros(2*N+dl+1,1);
261 for pp = 1:length(AA(1,:))
262 nf(pp) = norm(AA(:,pp),2); % Euclidean norm
263 AA(:,pp) = AA(:,pp)./nf(pp); % Normalization
264 end
265
266 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
267 % Solving augmented problem
268
269 % XA = pinv(AA)*BA;
270 % XA = inv((AA.')*AA)*(AA.')*BA;
271
272 % XA = AA.'*AA\AA.'*BA;
273
274 XA = AA\BA;
275
276 XA = XA./nf.'; % renormalization
277
278 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
279 % Checking the tolerance
280
281 if abs(XA(end))<TOLlow || abs(XA(end))>TOLhigh
282
283 if XA(end)==0
284 Dnew=1;
285 elseif abs(XA(end))<TOLlow
286 Dnew=sign(XA(end))*TOLlow;
287 elseif abs(XA(end))>TOLhigh
288 Dnew=sign(XA(end))*TOLhigh;
289 end
290
291 for pp = 1:length(AA(1,:))
292 AA(:,pp) = AA(:,pp).*nf(pp); %removing previous scaling
293 end
294
295 ind=length(AA(:,1))/2; %index to additional row related to relaxation
296
297 AA(ind,:)=[]; % removing relaxation term
298
299 BA=-Dnew*AA(:,end); %new right side
300
301 AA(:,end)=[];
302
303 nf(end)=[];
304
305 for pp = 1:length(AA(1,:))
306 nf(pp) = norm(AA(:,pp),2); % Euclidean norm
307 AA(:,pp) = AA(:,pp)./nf(pp); % Normalization
308 end
309
310 % XA=(AA.'*AA)\(AA.'*BA); % using normal equation
311
312 XA=AA\BA;
313
314 XA = XA./nf.'; % renormalization
315
316 XA=[XA;Dnew];
317
318 end
319
320 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
321 % Finding zeros of sigma
322
323 lsr = XA((N+dl)*Nc+1:(N+dl)*Nc+N,1); % collect the least square results
324
325 D = XA(end); % direct term of sigma
326
327 CPOLES = diag(poles);
328 B = ones(N,1);
329 C = lsr.';
330
331 for kk = 1:N
332 if cindex(kk) == 1
333 CPOLES(kk,kk)=real(poles(kk));
334 CPOLES(kk,kk+1)=imag(poles(kk));
335 CPOLES(kk+1,kk)=-1*imag(poles(kk));
336 CPOLES(kk+1,kk+1)=real(poles(kk));
337 B(kk,1) = 2;
338 B(kk+1,1) = 0;
339 end
340 end
341
342 H = CPOLES-B*C/D;
343
344 szeros=eig(H);
345
346 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
347 % Ruling out unstable poles
348
349 % This option force the poles of f to stay inside the left side of the
350 % complex plane
351
352 if stab
353 unst = real(szeros)>0;
354 szeros(unst) = szeros(unst)-2*real(szeros(unst)); % Mirroring respect to the complex axes
355 end
356 szeros = sort(szeros);
357 N = length(szeros);
358
359 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
360 % Separating complex poles from real poles and ordering
361
362 rnpoles = [];
363 inpoles = [];
364 for tt = 1:N
365 if imag(szeros(tt)) == 0
366 % collecting real poles
367 rnpoles = [rnpoles; szeros(tt)];
368 else
369 % collecting complex poles
370 inpoles = [inpoles; szeros(tt)];
371 end
372 end
373
374 % Sorting complex poles in order to have them in the expected order a+jb
375 % and a-jb a>0 b>0
376 inpoles = sort(inpoles);
377 npoles = [rnpoles;inpoles];
378 npoles = npoles - 2.*1i.*imag(npoles);
379
380 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
381 % Marking complex and real poles
382
383 cindex=zeros(N,1);
384 for m=1:N
385 if imag(npoles(m))~=0
386 if m==1
387 cindex(m)=1;
388 else
389 if cindex(m-1)==0 || cindex(m-1)==2
390 cindex(m)=1; cindex(m+1)=2;
391 else
392 cindex(m)=2;
393 end
394 end
395 end
396 end
397
398 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
399 % Direct problem
400 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
401
402 % Matrix initialinzation
403 nB(1:Nz,1:Nc) = real(w.*y).';
404 nB(Nz+1:2*Nz,1:Nc) = imag(w.*y).';
405
406 B = zeros(2*Nz,1);
407 nAD = zeros(Nz,N+dl);
408 AD = zeros(2*Nz,N+dl);
409 Ak = zeros(Nz,N+dl);
410
411 % Defining Ak
412 for jj = 1:N
413 if cindex(jj) == 1 % conjugate complex couple of poles
414 Ak(:,jj) = 1./(s-npoles(jj)) + 1./(s-npoles(jj+1));
415 Ak(:,jj+1) = 1i./(s-npoles(jj)) - 1i./(s-npoles(jj+1));
416 elseif cindex(jj) == 0 % real pole
417 Ak(:,jj) = 1./(s-npoles(jj));
418 end
419 end
420
421 if dt
422 Ak(1:Nz,N+dl) = ones(Nz,1); % considering the direct term
423 end
424
425 XX = zeros(Nc,N+dl);
426 for nn = 1:Nc
427
428 % Defining AD
429 for m=1:N+dl
430 nAD(1:Nz,m) = w(nn,:).'.*Ak(1:Nz,m);
431 end
432
433 B(1:2*Nz,1) = nB(1:2*Nz,nn);
434
435 AD(1:Nz,:) = real(nAD);
436 AD(Nz+1:2*Nz,:) = imag(nAD);
437
438 % Normalization factor
439 nf = zeros(N+dl,1);
440 for pp = 1:N+dl
441 nf(pp,1) = norm(AD(:,pp),2); % Euclidean norm
442 AD(:,pp) = AD(:,pp)./nf(pp,1); % Normalization
443 end
444
445 % Solving direct problem
446
447 % XD = inv((AD.')*AD)*(AD.')*B;
448 % XD = AD.'*AD\AD.'*B;
449 % XD = pinv(AD)*B;
450 XD = AD\B;
451
452 XD = XD./nf; % Renormalization
453
454 XX(nn,1:N) = XD(1:N).';
455
456 end
457
458 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
459 % Final residues and poles of f
460
461 lsr = XX(:,1:N);
462
463 res = zeros(N,Nc);
464 % Real poles have real residues, complex poles have comples residues
465 for nn = 1:Nc
466 for tt = 1:N
467 if cindex(tt) == 1 % conjugate complex couple of poles
468 res(tt,nn) = lsr(nn,tt)+1i*lsr(nn,tt+1);
469 res(tt+1,nn) = lsr(nn,tt)-1i*lsr(nn,tt+1);
470 elseif cindex(tt) == 0 % real pole
471 res(tt,nn) = lsr(nn,tt);
472 end
473 end
474 end
475
476 poles = npoles;
477
478 if dt
479 dterm = XX(:,N+dl).';
480 else
481 dterm = zeros(1,Nc);
482 end
483
484 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
485 % Calculating responses and residuals
486
487 mresp = zeros(Nz,Nc);
488 rdl = zeros(Nz,Nc);
489 mse = zeros(1,Nc);
490
491 for nn = 1:Nc
492 % freq resp of the fit model
493 r = res(:,nn);
494 p = poles;
495 d = dterm(:,nn);
496
497 Nf = length(f);
498 N = length(p);
499
500 rsp = zeros(Nf,1);
501 indx = (0:length(d)-1).';
502 for ii = 1:Nf
503 for jj = 1:N
504 rsptemp = r(jj)/(1i*2*pi*f(ii)-p(jj));
505 rsp(ii) = rsp(ii) + rsptemp;
506 end
507 % Direct terms response
508 rsp(ii) = rsp(ii) + sum((((1i*2*pi*f(ii))*ones(length(d),1)).^indx).*d);
509 end
510
511 % Model response
512 mresp(:,nn) = rsp;
513
514 % Residual
515 rdl(:,nn) = y(nn,:).' - rsp;
516
517 % RMS error
518 % mse(:,nn) = sqrt(sum((abs(rdl(:,nn)./(y(nn,:).')).^2))/(Nf-N));
519 % Chi Square or mean squared error
520 % Note that this error is normalized to the input data in order to
521 % comparable between different sets of data
522 mse(:,nn) = sum((rdl(:,nn)./y(nn,:).').*conj((rdl(:,nn)./y(nn,:).')))/(Nf-N);
523
524 end
525
526 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
527 % Plotting response
528
529 switch plotting
530
531 case 0
532 % No plot
533
534 case 1
535 % LogLog plot for absolute value
536 figure(plth)
537 subplot(2,1,1);
538 p1 = loglog(f,abs(y.'),'k');
539 hold on
540 p2 = loglog(f,abs(mresp),'r');
541 p3 = loglog(f,abs(rdl),'b');
542 xlabel('Frequency [Hz]')
543 ylabel('Amplitude')
544 legend([p1(1) p2(1) p3(1)],'Original', 'VCFIT','Residual')
545 hold off
546
547 subplot(2,1,2);
548 p4 = semilogx(f,(180/pi).*unwrap(angle(y.')),'k');
549 hold on
550 p5 = semilogx(f,(180/pi).*unwrap(angle(mresp)),'r');
551 xlabel('Frequency [Hz]')
552 ylabel('Phase [Deg]')
553 legend([p4(1) p5(1)],'Original', 'VCFIT')
554 hold off
555
556 case 2
557 % Semilogx plot for absolute value
558 figure(plth)
559 subplot(2,1,1);
560 p1 = semilogx(f,abs(y.'),'k');
561 hold on
562 p2 = semilogx(f,abs(mresp),'r');
563 p3 = semilogx(f,abs(rdl),'b');
564 xlabel('Frequency [Hz]')
565 ylabel('Amplitude')
566 legend([p1(1) p2(1) p3(1)],'Original', 'VCFIT','Residual')
567 hold off
568
569 subplot(2,1,2);
570 p4 = semilogx(f,(180/pi).*unwrap(angle(y.')),'k');
571 hold on
572 p5 = semilogx(f,(180/pi).*unwrap(angle(mresp)),'r');
573 xlabel('Frequency [Hz]')
574 ylabel('Phase [Deg]')
575 legend([p4(1) p5(1)],'Original', 'VCFIT')
576 hold off
577
578 case 3
579 % Semilogy plot for absolute value
580 figure(plth)
581 subplot(2,1,1);
582 p1 = semilogy(f,abs(y.'),'k');
583 hold on
584 p2 = semilogy(f,abs(mresp),'r');
585 p3 = semilogy(f,abs(rdl),'b');
586 xlabel('Frequency [Hz]')
587 ylabel('Amplitude')
588 legend([p1(1) p2(1) p3(1)],'Original', 'VCFIT','Residual')
589 hold off
590
591 subplot(2,1,2);
592 p4 = semilogy(f,(180/pi).*unwrap(angle(y.')),'k');
593 hold on
594 p5 = semilogy(f,(180/pi).*unwrap(angle(mresp)),'r');
595 xlabel('Frequency [Hz]')
596 ylabel('Phase [Deg]')
597 legend([p4(1) p5(1)],'Original', 'VCFIT')
598 hold off
599
600 case 4
601 % Linear plot for absolute value
602 figure(plth)
603 subplot(2,1,1);
604 p1 = plot(f,abs(y.'),'k');
605 hold on
606 p2 = plot(f,abs(mresp),'r');
607 p3 = plot(f,abs(rdl),'b');
608 xlabel('Frequency [Hz]')
609 ylabel('Amplitude')
610 legend([p1(1) p2(1) p3(1)],'Original', 'VCFIT','Residual')
611 hold off
612
613 subplot(2,1,2);
614 p4 = plot(f,(180/pi).*unwrap(angle(y.')),'k');
615 hold on
616 p5 = plot(f,(180/pi).*unwrap(angle(mresp)),'r');
617 xlabel('Frequency [Hz]')
618 ylabel('Phase [Deg]')
619 legend([p4(1) p5(1)],'Original', 'VCFIT')
620 hold off
621
622 end