Mercurial > hg > ltpda

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