0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 % PFALLPSYMZ all pass filtering in order to stabilize TF poles and zeros.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
2 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
4 % DESCRIPTION:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
5 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
6 % All pass filtering in order to stabilize transfer functions poles.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
7 % It inputs a partial fraction expanded discrete model and outpu
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
8 % residues, poles direct terms and frequency response of the stabilized
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
9 % model. Function can handle multiple models with common poles.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
10 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
11 % CALL:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
12 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
13 % [nr,np,nd,resp] = pfallpsymz(r,p,d,f,fs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
14 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
15 % INPUTS:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
16 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
17 % r: are residues. (Npx1) or (NpxM) vector
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
18 % p: are poles. (Npx1) vector
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
19 % d: is direct term (1x1) or (1xM) vector
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
20 % mresp: input model response. (Nx1) or (NxM) vector
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21 % f: is the frequancies vector in (Hz). (Nx1) vector
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22 % fs: is the sampling frequency in (Hz). (1x1)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24 % OUTPUTS:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26 % nr: new residues. (Npx1) or (NpxM) vector
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27 % np: new stable poles. (Npx1) vector
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28 % nd: new direct term. (1x1) or (1xM) vector
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29 % nmresp: new model response. (Nx1) or (NxM) vector
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31 % NOTE:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33 % This function make use of symbolic math toolbox functions
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
35 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
36 % VERSION: $Id: pfallpsymz.m,v 1.5 2008/11/10 15:40:11 luigi Exp $
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
37 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
38 % HISTORY: 12-09-2008 L Ferraioli
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
39 % Creation
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41 function [nr,np,nd,nmresp] = pfallpsymz(r,p,d,mresp,f,fs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43 % Reshaping
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44 [a,b] = size(r);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45 if a<b
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
46 r = r.'; % reshape as a column vector
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49 [a,b] = size(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50 if a<b
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51 p = p.'; % reshape as a column vector
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
52 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
53
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
54 [a,b] = size(f);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
55 if a<b
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
56 f = f.'; % reshape as a column vector
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
57 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
58
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
59 [a,b] = size(f);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
60 if a<b
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
61 f = f.'; % reshape as a column vector
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
62 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
63
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
64 [a,b] = size(d);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
65 if a > b
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
66 d = d.'; % reshape as a row
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
67 d = d(1,:); % taking the first row (the function can handle only simple constant direct terms)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
68 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
69
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
70 if isempty(fs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
71 fs = 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
72 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
73 [a,b] = size(fs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
74 if a ~= b
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
75 disp(' Fs has to be a number. Only first term will be considered! ')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
76 fs = fs(1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
77 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
78
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
79 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
80 f = sym(f,'f');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
81 fs = sym(fs,'f');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
82 syms z
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
83 p = sym(p,'f');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
84
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
85 % stabilize poles
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
86 sp = p;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
87 unst = abs(double(sp)) > 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
88 sp(unst) = 1./conj(sp(unst));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
89 skp = prod(1./conj(p(unst)));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
90
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
91 [Na,Nb] = size(r);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
92
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
93 for nn = 1:Nb
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
94
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
95 % pp = p(unst);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
96 % psp = sp(unst);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
97 % tterm = 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
98 % for ii = 1:length(pp)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
99 % tterm = tterm.*(z-pp(ii))./(z-psp(ii));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
100 % end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
101 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
102 % s = cos((2*pi/fs).*f) + 1i.*sin((2*pi/fs).*f);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
103 % allresp = subs(tterm,z,s);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
104 % phs = angle(allresp);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
105 % nmresp(:,nn) = mresp(:,nn).*(cos(phs)+1i.*sin(phs));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
106 % np = double(sp);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
107
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
108 % digits(100)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
109 % Defining inputs as symbolic variable
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
110 rt = sym(r(:,nn),'f');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
111 % p = sym(p,'f');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
112 dt = sym(d(1,nn),'f');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
113 % f = sym(f,'f');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
114 % fs = sym(fs,'f');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
115 % syms z
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
116
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
117 % Function gain coefficient
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
118 k = sum(rt)+dt;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
119
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
120 Np = length(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
121
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
122 % Defining the symbolic transfer function
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
123 vter = rt.*z./(z-p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
124 vter = [vter; dt];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
125 Mter = diag(vter);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
126 H = trace(Mter);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
127
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
128 % Factorizing the transfer function
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
129 HH = factor(H);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
130
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
131 % Extracting numerator and denominator
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
132 [N,D] = numden(HH);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
133
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
134 % Symbolci calculation of function zeros
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
135 cN = coeffs(expand(N),z);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
136 n = length(cN);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
137 cN2 = -1.*cN(2:end)./cN(1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
138 A = sym(diag(ones(1,n-2),-1));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
139 A(1,:) = cN2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
140 zrs = eig(A);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
141 if double(d) == 0
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
142 zrs = [zrs; sym(0,'f')];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
143 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
144
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
145 % % stabilize zeros
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
146 % szrs = zrs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
147 % unst = abs(double(szrs)) > 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
148 % szrs(unst) = 1./conj(szrs(unst));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
149 % skz = prod(conj(zrs(unst)));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
150
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
151 % % stabilize poles
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
152 % sp = p;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
153 % unst = abs(double(sp)) > 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
154 % sp(unst) = 1./conj(sp(unst));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
155 % skp = prod(1./conj(p(unst)));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
156
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
157 % Correcting for some special cases
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
158 % if isempty(skz)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
159 % skz = sym(1,'f');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
160 % end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
161 if isempty(skp)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
162 skp = sym(1,'f');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
163 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
164 if isempty(k)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
165 k = sym(1,'f');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
166 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
167
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
168 % Calculating new gain
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
169 % sk = real(k*skz*skp);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
170 sk = real(k*skp);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
171
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
172 HHHn = sym(1,'f');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
173
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
174 for jj = 1:Np
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
175 HHHn = HHHn.*(z-zrs(jj));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
176 tsp = sp;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
177 tsp(jj) = [];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
178 tHHHd = sym(1,'f');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
179 for kk = 1:Np-1
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
180 tHHHd = tHHHd.*(z-tsp(kk));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
181 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
182 HHHd(jj,1) = tHHHd;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
183
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
184 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
185
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
186 for jj = 1:Np
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
187 sr(jj,1) = subs(sk*HHHn/(z*HHHd(jj,1)),z,sp(jj));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
188 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
189
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
190 np = double(sp);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
191 for kk = 1:Np
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
192 if imag(np(kk)) == 0
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
193 sr(kk) = real(sr(kk));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
194 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
195 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
196
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
197 nr(:,nn) = double(sr);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
198 nd(:,nn) = double(dt);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
199
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
200 % Model evaluation
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
201 pfparams.type = 'disc';
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
202 pfparams.freq = f;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
203 pfparams.fs = fs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
204 pfparams.res = nr(:,nn);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
205 pfparams.pol = np;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
206 pfparams.dterm = nd(:,nn);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
207 pfr = utils.math.pfresp(pfparams);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
208 resp = pfr.resp;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
209
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
210 ratio = mean(abs(mresp(:,nn))./abs(resp));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
211 resp = resp.*ratio;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
212 nr(:,nn) = nr(:,nn).*ratio;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
213 nd(:,nn) = nd(:,nn).*ratio;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
214 nmresp(:,nn) = resp;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
215
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
216 end
|