0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 % PFALLPSYMZ2 all pass filtering 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 function poles and
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
7 % zeros. It inputs a partial fraction expanded discrete model and
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
8 % outputs a pole-zero minimum phase system
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
9 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
10 % CALL:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
11 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
12 % resp= pfallpz2(ip,mresp,f,fs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
13 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
14 % INPUTS:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
15 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
16 % ip: is a struct with fields named poles
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
17 % mresp: is a vector with functions response
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
18 % f: is the frequancies vector in (Hz)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
19 % fs: is the sampling frequency in (Hz)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
20 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21 % OUTPUTS:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23 % resp: is the stable functions frequency response
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25 % NOTE:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28 % VERSION: $Id: pfallpz.m,v 1.6 2009/06/10 15:47:00 luigi Exp $
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30 % HISTORY: 12-09-2008 L Ferraioli
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31 % Creation
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33 function varargout = pfallpsymz2(ip,mresp,f,fs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
35 [a,b] = size(ip);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
36 if a<b
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
37 ip = ip.'; % reshape as a column vector
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
38 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
39
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40 [a,b] = size(f);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41 if a<b
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42 f = f.'; % reshape as a column vector
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45 if isempty(fs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
46 fs = 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48 [a,b] = size(fs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49 if a ~= b
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50 disp(' Fs has to be a number. Only first term will be considered! ')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51 fs = fs(1);
|
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
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
55 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
56 Nb = numel(ip);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
57 for nn = 1:Nb
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
58
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
59 p = ip(nn).poles;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
60
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
61 % stabilizing poles
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
62 sp = sym(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
63 unst = abs(p) > 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
64 sp(unst) = conj(sp(unst));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
65
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
66 pp = sym(p(unst));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
67 psp = sp(unst);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
68 syms z
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
69 allpsym = 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
70 for jj=1:numel(psp)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
71 allpsym = allpsym.*((z-pp(jj))./(z*psp(jj)-1));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
72 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
73 funcell{nn} = allpsym;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
74
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
75 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
76
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
77 fullallprsp = 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
78
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
79 for nn = 1:Nb
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
80 symexpr = funcell{nn};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
81 nterm = subs(symexpr,z,cos((2*pi/fs).*f) + 1i.*sin((2*pi/fs).*f));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
82 % willing to work with columns
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
83 if size(nterm,2)>1
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
84 nterm = nterm.';
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
85 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
86
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
87 fullallprsp = fullallprsp.*nterm;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
88
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
89 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
90
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
91 % dballprsp = double(fullallprsp);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
92 % rallp = real(dballprsp);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
93 % iallp = imag(dballprsp);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
94 % ang = angle(dballprsp);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
95
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
96 for kk=1:Nb
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
97 sresp(:,kk) = mresp(:,kk).*fullallprsp;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
98 % resp(:,kk) = mresp(:,kk).*(cos(ang)+1i.*sin(ang));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
99 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
100
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
101 for kk=1:Nb
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
102 resp(:,kk) = double(sresp(:,kk));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
103 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
104
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
105
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
106 % output
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
107 if nargout == 1
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
108 varargout{1} = resp;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
109 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
110 error('Too many output arguments!')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
111 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
112 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
113
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
114
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
115
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
116
|