Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/startpoles.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 % STARTPOLES defines starting poles for fitting procedures ctfit, dtfit.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % DESCRIPTION
4 %
5 % Defines the starting poles for the fitting procedure with ctfit or
6 % dtfit. Starting poles definition process is different for s-domain
7 % and z-domain.
8 % s-domain identification:
9 % Starting poles can be real or complex in conjugate couples. Real
10 % poles are chosen on the [-2*pi*f(1),-2*pi*f(end)] intervall. Complex poles can
11 % be defined with the real and imaginary parts logspaced or linespaced
12 % on the intervall [2\pi f(1),2\pi f(end)]. Ratio between real and
13 % imaginary part can be setted by the user.
14 % z-domain identification:
15 % Starting poles can be real or come in complex conjugate couples. Real
16 % poles are chosen on the [-1,1] intervall. Complex poles are
17 % chosen inside the unit circle as:
18 % \alfa e^{j\pi\theta}
19 % where \theta is linespaced inside the intervall [0,2\pi]. In this
20 % case two different methods are used: the first method define angles
21 % as \theta = linspace(0,pi,N/2+1), then take out the first element and
22 % construct the complex conjugate couples. If N is odd the first
23 % element is added as real pole. With this method the last two
24 % conjugates poles have a real part very similar to that of the real
25 % pole. This may generate problems so a second method is implemented in
26 % which the angle are generated as \theta = linspace(0,pi,N/2+2). Then
27 % the first and the last elements of the set are taken out and the
28 % first element is used only for N odd. The last element instead is
29 % never used. This allow to have well distributed poles on the unit
30 % circle. The amplitude parameter \alfa can be set by the user.
31 %
32 % CALL:
33 %
34 % spoles = startpoles(order,f,params)
35 %
36 % INPUT:
37 %
38 % order: is the function order, ie. the number of desired poles
39 % f: is the frequency vector in Hz
40 % params: is a struct with the setting parameters
41 %
42 % params.type = 'CONT' --> Output a set of poles for s-domain
43 % identification
44 % params.type = 'DISC' --> Output a set of poles for z-domain
45 % identification
46 %
47 % params.spolesopt = 1 --> generate linear spaced real poles on
48 % the intervall [-2*pi*f(1),-2*pi*f(end)] for s-domain and [-1,1] for z-domain.
49 % params.spolesopt = 2 --> in case of s-domain generates logspaced
50 % complex conjugates poles. In case of z-domain generates complex
51 % conjugates poles with \theta = linspace(0,pi,N/2+1).
52 % params.spolesopt = 3 --> in case of s-domain generates linespaced
53 % complex conjugates poles. In case of z-domain generates complex
54 % conjugates poles with \theta = linspace(0,pi,N/2+2). We advice to
55 % make use of this option for z-domain identification.
56 %
57 % params.pamp = # --> s-domain: set the ratio \alfa/\beta between
58 % poles real and imaginary parts. Adviced value is 0.01.
59 % params.pamp = # --> z-domain: set the amplitude of the poles.
60 % Adviced value is 0.98.
61 %
62 % OUTPUT:
63 %
64 % spoles: is the set of starting poles
65 %
66 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
67 % VERSION: $Id: startpoles.m,v 1.6 2010/04/29 09:00:00 luigi Exp $
68 % HISTORY: 08-10-2008 L Ferraioli
69 % Creation
70 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
71
72 function spoles = startpoles(order,f,params)
73
74 % Default input struct
75 defaultparams = struct('spolesopt',1, 'type','CONT', 'pamp', 0.98);
76
77 names = {'spolesopt', 'type', 'pamp'};
78
79 % collecting input and default params
80 if ~isempty(params)
81 for jj=1:length(names)
82 if isfield(params, names(jj))
83 defaultparams.(names{1,jj}) = params.(names{1,jj});
84 end
85 end
86 end
87
88 type = defaultparams.type;
89 spolesopt = defaultparams.spolesopt;
90 pamp = defaultparams.pamp;
91
92 N = order;
93
94 % switching between continuous and discrete
95 switch type
96 case 'CONT'
97 switch spolesopt
98
99 case 0
100 disp(' Using external starting poles')
101
102 case 1 % real starting poles
103 spoles = -1.*2.*pi.*linspace(f(1),f(end),N).';
104
105 case 2 % complex logspaced starting poles
106 if f(1)==0
107 bet=2.*pi.*logspace(log10(f(2)),log10(f(end)),N/2);
108 else
109 bet=2.*pi.*logspace(log10(f(1)),log10(f(end)),N/2);
110 end
111 spoles=[];
112 for n=1:length(bet)
113 alf=-bet(n)*pamp;
114 spoles=[spoles;(alf-1i*bet(n));(alf+1i*bet(n))];
115 end
116 if (N-2*floor(N/2)) ~= 0
117 % rpl = rand(1,1);
118 % if rpl > 0
119 % rpl = -1*rpl;
120 % end
121 rpl = -1;
122 spoles = [rpl; spoles];
123 end
124
125 case 3 % complex linspaced starting poles
126 bet=linspace(2*pi*f(1),2*pi*f(end),N/2);
127 spoles=[];
128 for n=1:length(bet)
129 alf=-bet(n)*pamp;
130 spoles=[spoles;(alf-1i*bet(n));(alf+1i*bet(n))];
131 end
132 if (N-2*floor(N/2)) ~= 0
133 % rpl = rand(1,1);
134 % if rpl > 0
135 % rpl = -1*rpl;
136 % end
137 rpl = -1;
138 spoles = [rpl; spoles];
139 end
140
141 end
142 case 'DISC'
143 switch spolesopt
144
145 case 0
146 disp(' Using external starting poles')
147
148 case 1 % real starting poles
149 spoles = linspace(-0.99,0.99,N).';
150
151 case 2 % complex starting poles
152 ang = linspace(0,pi,N/2+1);
153 spoles=[];
154 for nn=2:length(ang)
155 spoles=[spoles; exp(1i*ang(nn)); exp(-1i*ang(nn))]; % Taking complex conjugate pairs on the unit circle
156 end
157 if (N-2*floor(N/2)) ~= 0
158 rpl = exp(1i*ang(1));
159 spoles = [rpl; spoles];
160 end
161 spoles = spoles.*pamp; % shifting starting ?poles a little inside the unit circle
162
163 case 3 % complex starting poles
164 ang = linspace(0,pi,N/2+2);
165 spoles=[];
166 for nn=2:length(ang)-1
167 spoles=[spoles; exp(1i*ang(nn)); exp(-1i*ang(nn))]; % Taking complex conjugate pairs on the unit circle
168 end
169 if (N-2*floor(N/2)) ~= 0
170 rpl = exp(1i*ang(1));
171 spoles = [rpl; spoles];
172 end
173 spoles = spoles.*pamp; % shifting starting ?poles a little inside the unit circle
174
175 end
176
177 end
178
179 end