Mercurial > hg > ltpda
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 |