Mercurial > hg > ltpda

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