Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/pfallps2.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 % PFALLPS2 all pass filtering to stabilize TF poles and zeros.
2 %
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 % DESCRIPTION:
5 %
6 % All pass filtering in order to stabilize transfer function poles and
7 % zeros. It inputs a partial fraction expanded discrete model and
8 % outputs a pole-zero minimum phase system
9 %
10 % CALL:
11 %
12 % [resp,np] = pfallps2(ip,mresp,f)
13 %
14 % INPUTS:
15 %
16 % ip: are poles
17 % f: is the frequancies vector in (Hz)
18 % fs: is the sampling frequency in (Hz)
19 %
20 % OUTPUTS:
21 %
22 % resp: is the functions phase frequency response
23 % np: are the new stable poles
24 %
25 % NOTE:
26 %
27 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
28 % VERSION: $Id: pfallpz.m,v 1.6 2009/06/10 15:47:00 luigi Exp $
29 %
30 % HISTORY: 12-09-2008 L Ferraioli
31 % Creation
32 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33 function varargout = pfallps2(ip,mresp,f)
34
35 [a,b] = size(ip);
36 if a<b
37 ip = ip.'; % reshape as a column vector
38 end
39
40 [a,b] = size(f);
41 if a<b
42 f = f.'; % reshape as a column vector
43 end
44
45 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
46
47 Nb = numel(ip);
48 for nn = 1:Nb
49
50 p = ip(nn).poles;
51
52 % stabilizing poles
53 sp = p;
54 unst = real(p) > 0;
55
56 sp(unst) = conj(sp(unst));
57
58
59
60 pp = p(unst);
61 psp = sp(unst);
62 allpstr = '(1';
63 for jj = 1:numel(sp(unst))
64 allpstr = [allpstr sprintf('.*((s-%0.20d)./(s+%0.20d))',pp(jj),psp(jj))];
65 end
66 allpstr = [allpstr ')'];
67
68 funcell{nn} = allpstr;
69
70 end
71
72
73 s = (1i*2*pi).*f;
74 fullallprsp = 1;
75
76 for nn = 1:Nb
77
78 nterm = eval(funcell{nn});
79 % willing to work with columns
80 if size(nterm,2)>1
81 nterm = nterm.';
82 end
83
84 allprsp(:,nn) = nterm;
85
86 fullallprsp = fullallprsp.*nterm;
87
88 end
89
90 phs = angle(fullallprsp);
91
92 for kk=1:Nb
93 resp(:,kk) = mresp(:,kk).*(cos(phs)+1i.*sin(phs));
94 end
95
96
97 % output
98 if nargout == 1
99 varargout{1} = resp;
100 else
101 error('Too many output arguments!')
102 end
103 end
104
105
106
107