Mercurial > hg > ltpda

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