Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/sineParams.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 % SINEPARAMS estimates parameters of sinusoids
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: SINEPARAMS estimates the parameters of the sinusoids in the
5 % time series. Number of sinusoids needs to be specified.
6 %
7 % CALL: b = sineParams(a,pl)
8 %
9 % INPUTS: a - input AOs
10 % pl - parameter list (see below)
11 %
12 % OUTPUTs: b - pest object
13 %
14 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'sineParams')">Parameters Description</a>
15 %
16 % VERSION: $Id: sineParams.m,v 1.9 2011/04/08 08:56:13 hewitson Exp $
17 %
18 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
19
20 function varargout = sineParams(varargin)
21
22 % Check if this is a call for parameters
23 if utils.helper.isinfocall(varargin{:})
24 varargout{1} = getInfo(varargin{3});
25 return
26 end
27
28 import utils.const.*
29 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
30
31 % Collect input variable names
32 in_names = cell(size(varargin));
33 for ii = 1:nargin,in_names{ii} = inputname(ii);end
34
35 % Collect all AOs and plists
36 [aos, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
37 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
38
39 %%% Decide on a deep copy or a modify
40 bs = copy(aos, nargout);
41
42 % combine plists
43 pl = parse(pl, getDefaultPlist());
44
45 % check (current version with only one AO)
46 if numel(aos) >1
47 error('### Current version of ''sineParams'' works only with one AO')
48 end
49
50 % get parameters
51 Nsine = find(pl, 'N');
52 realSine = find(pl, 'real');
53 method = find(pl,'method');
54 nonlin = find(pl,'non-linear');
55
56 % check number of sinusoids
57 if isempty(Nsine)
58 error('### You need to set the number of sinusoids ''N'' ')
59 end
60
61 if realSine
62 Nexp = 2*Nsine;
63 else
64 Nexp = Nsine;
65 end
66 f =[];A=[];
67 % for i = 1:numel(bs)
68
69 switch method
70 case 'MUSIC'
71 % MUSIC algorithm
72 [w,pow,w_mse,pow_mse] = utils.math.rootmusic(bs.y,Nexp,bs.fs);
73
74 if realSine
75 f = w(1:2:end);
76 df = w_mse(1:2:end);
77 % A = pow;
78 % dA = pow_mse;
79 A = sqrt(2*pow); % from the expression for power in a sinusoid
80 dA = 1/sqrt(2*pow)*pow_mse;
81 else
82 f = w;
83 df = w_mse;
84 A = sqrt(2*pow); % from the expression for power in a sinusoid
85 dA = A*sqrt(2/pow)*pow_mse;
86 end
87
88 case 'ESPRIT'
89
90 % p = 1; % num. sinusoid
91 % N = 50;
92 %
93 % m = xcov(c.y,N-1);
94 %
95 % clear cmat
96 % for i =1:N
97 % cmat(i,:) = m(N-(i-1):end-(i-1));
98 % end
99 %
100 % [U,S,UT] = svd(cmat);
101 %
102 % D = diag(S);
103 %
104 % Ls = D(1:2*p);
105 % Us = U(:,1:2*p);
106 %
107 % A1 = Us(1:end-1,:);
108 % A2 = Us(2:end,:);
109 %
110 % M = A1\A2;
111 %
112 % freq = imag(-log(eig(M))/2/pi*fs)
113 % error = real(-log(eig(M))/2/pi*fs)
114
115 case 'IFFT'
116
117 % % fft
118 % p = psd(phi(i),plist('win','Hanning','scale','AS','Navs',1));
119 %
120 % [m,index] = max(p.y);
121 % ratio between max. and next value
122 % alpha = p.y(index+1)/p.y(index);
123 %
124 % xm = (2*alpha-1)/(alpha+1);
125 % frequency
126 % f(i) = (i+xm)/phi(i).nsecs;
127 % amplitude
128 % A(i) = abs(2*pi*xm*(1-xm)/sin(pi*xm)*exp(-pi*1i*xm)*(1*xm)*p.y(i));
129 %
130 %
131 %
132 %
133 end
134 % end
135
136 % create output pest
137 mdl = [];
138 p = pest();
139 for i = 1:Nsine
140 Cname = sprintf('C%d',i);
141 fname = sprintf('f%d',i);
142 Aname = sprintf('A%d',i);
143 pname = sprintf('phi%d',i);
144 % setY
145 p.setYforParameter(Cname,0);
146 p.setYforParameter(Aname,A(i));
147 p.setYforParameter(fname,f(i));
148 p.setYforParameter(pname,0);
149 % errors for single sinusoid only, error is the sqrt(mse)
150 if Nsine == 1
151 p.setDyForParameter(Cname,inf);
152 p.setDyForParameter(Aname,sqrt(dA(i)));
153 p.setDyForParameter(fname,sqrt(df(i)));
154 p.setDyForParameter(pname,inf);
155 end
156 % smodel for each sinusoid
157 m = smodel(sprintf('%s + %s*sin(2*pi*%s*t+%s)',Cname,Aname,fname,pname));
158 m.setXunits('s');
159 m.setXvar('t');
160 m.setXvals(bs.x);
161 m.setParams({Cname,Aname,fname,pname},[0 A(i) f(i) 0]);
162 % optional non-linear
163 if nonlin
164 pl = plist('Function',m);
165 pnl = xfit(bs,pl);
166 p.setY(pnl.y);
167 p.setDy(pnl.dy);
168 p.setCorr(pnl.corr);
169 p.setCov(pnl.cov);
170 p.setDof(pnl.dof);
171 p.setChi2(pnl.chi2);
172 % set values in model as well
173 m.setValues(pnl.y);
174 end
175 mdl = [mdl m];
176 end
177 % set name
178 p.name = sprintf('sineParams(%s)', ao_invars{:});
179 % set models
180 p.setModels(mdl);
181 % Add history
182 p.addHistory(getInfo('None'), pl, ao_invars(:), bs(:).hist);
183
184 % Set outputs
185 if nargout > 0
186 varargout{1} = p;
187 end
188 end
189
190
191 %--------------------------------------------------------------------------
192 % Get Info Object
193 %--------------------------------------------------------------------------
194 function ii = getInfo(varargin)
195 if nargin == 1 && strcmpi(varargin{1}, 'None')
196 sets = {};
197 pl = [];
198 else
199 sets = {'Default'};
200 pl = getDefaultPlist;
201 end
202 % Build info object
203 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: sineParams.m,v 1.9 2011/04/08 08:56:13 hewitson Exp $', sets, pl);
204 end
205
206 %--------------------------------------------------------------------------
207 % Get Default Plist
208 %--------------------------------------------------------------------------
209 function plout = getDefaultPlist()
210 persistent pl;
211 if exist('pl', 'var')==0 || isempty(pl)
212 pl = buildplist();
213 end
214 plout = pl;
215 end
216
217 function pl = buildplist()
218
219 pl = plist();
220
221 % number of sinusoids
222 p = param({'N', 'Number of sinusoids/complex exp. in the time series'},paramValue.EMPTY_DOUBLE);
223 pl.append(p);
224
225 % number of sinusoids
226 p = param({'real', 'Set to true if working with real sinusoids'},paramValue.TRUE_FALSE);
227 pl.append(p);
228
229 % number of sinusoids
230 p = param({'method', ['Choose one of the following methods:<ul>', ...
231 '<li>MUSIC - MUltiple SIgnal Classification algorithm </li>']}, {1, {'MUSIC'}, paramValue.SINGLE});
232 pl.append(p);
233
234 % non-linear
235 p = param({'non-linear','Set to true to perform non-linear fit'},...
236 paramValue.FALSE_TRUE);
237 pl.append(p);
238
239 end
240 % END
241