Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/linedetect.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 % LINEDETECT find spectral lines in the ao/fsdata objects.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: LINEDETECT find spectral lines in the ao/fsdata objects.
5 %
6 % CALL: b = linedetect(a, pl)
7 %
8 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'linedetect')">Parameters Description</a>
9 %
10 % VERSION: $Id: linedetect.m,v 1.15 2011/04/08 08:56:16 hewitson Exp $
11 %
12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13
14 function varargout = linedetect(varargin)
15
16 % Check if this is a call for parameters
17 if utils.helper.isinfocall(varargin{:})
18 varargout{1} = getInfo(varargin{3});
19 return
20 end
21
22 if nargout == 0
23 error('### cat cannot be used as a modifier. Please give an output variable.');
24 end
25
26 % Collect input variable names
27 in_names = cell(size(varargin));
28 for ii = 1:nargin,in_names{ii} = inputname(ii);end
29
30 % Collect all AOs
31 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
32 [pli, pl_invars] = utils.helper.collect_objects(varargin(:), 'plist', in_names);
33
34 % Decide on a deep copy or a modify
35 bs = copy(as, nargout);
36
37 Na = numel(bs);
38 if isempty(bs)
39 error('### Please input at least one AO.');
40 end
41
42 % Combine plists
43 if ~isempty(pli)
44 pl = parse(pli, getDefaultPlist());
45 else
46 pl = getDefaultPlist();
47 end
48
49 % Get parameters from plist
50 N = find(pl, 'N');
51 fsearch = find(pl, 'fsearch');
52 thresh = find(pl, 'thresh');
53
54 % Loop over input AOs
55 for jj = 1:Na
56 if isa(bs(jj).data, 'fsdata')
57
58 % Make noise-floor estimate
59 nf = smoother(bs(jj), pl);
60
61 % Make ratio
62 r = bs(jj)./nf;
63
64 % find lines
65 lines = findLines(bs(jj).data.getY, r.data.getX, r.data.getY, thresh, N, fsearch);
66
67 f = [lines(:).f];
68 y = [lines(:).a];
69
70 % Keep the data shpare of the input AO
71 if size(bs(jj).data.y, 2) == 1
72 f = f.';
73 y = y.';
74 end
75
76 % Make output data: copy the fsdata object so to inherit all the feautures
77 fs = copy(bs(jj).data, 1);
78
79 % Make output data: set the values
80 fs.setX(f);
81 fs.setY(y);
82
83 else
84 error('### I can only find lines in frequency-series AOs.');
85 end
86
87 %------- Make output AO
88
89 % make output analysis object
90 bs(jj).data = fs;
91
92 bs(jj).name = sprintf('lines(%s)', ao_invars{1});
93 bs(jj).addHistory(getInfo('None'), pl, ao_invars(jj), bs(jj).hist);
94 end
95
96 % Set output
97 if nargout == numel(bs)
98 % List of outputs
99 for ii = 1:numel(bs)
100 varargout{ii} = bs(ii);
101 end
102 else
103 % Single output
104 varargout{1} = bs;
105 end
106 end
107
108 %--------------------------------------------------------------------------
109 % find spectral lines
110 function lines = findLines(ay, f, nay, thresh, N, fsearch)
111
112 % look for spectral lines
113 l = 0;
114 pmax = 0;
115 pmaxi = 0;
116 line = [];
117 idx = find( f>=fsearch(1) & f<=fsearch(2) );
118 for jj = 1:length(idx)
119 v = nay(idx(jj));
120 if v > thresh
121 if v > pmax
122 pmax = v;
123 pmaxi = idx(jj);
124 end
125 else
126 if pmax > 0
127 % find index when we have pmax
128 fidx = pmaxi; %(find(nay(1:idx(jj))==pmax));
129 l = l+1;
130 line(l).idx = fidx;
131 line(l).f = f(fidx);
132 line(l).a = ay(fidx);
133 end
134 pmax = 0;
135 end
136 end
137
138 % Select largest peaks
139 lines = [];
140 if ~isempty(line)
141 [bl, lidx] = sort([line.a], 'descend');
142 lidxs = lidx(1:min([N length(lidx)]));
143 lines = line(lidxs);
144 disp(sprintf(' + found %d lines.', length([lines.f])));
145 else
146 disp(' + found 0 lines.');
147 end
148 end
149
150 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
151 % Local Functions %
152 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
153
154 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
155 %
156 % FUNCTION: getInfo
157 %
158 % DESCRIPTION: Get Info Object
159 %
160 % HISTORY: 11-07-07 M Hewitson
161 % Creation.
162 %
163 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
164
165 function ii = getInfo(varargin)
166 if nargin == 1 && strcmpi(varargin{1}, 'None')
167 sets = {};
168 pl = [];
169 else
170 sets = {'Default'};
171 pl = getDefaultPlist;
172 end
173 % Build info object
174 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: linedetect.m,v 1.15 2011/04/08 08:56:16 hewitson Exp $', sets, pl);
175 ii.setModifier(false);
176 end
177
178 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
179 %
180 % FUNCTION: getDefaultPlist
181 %
182 % DESCRIPTION: Get Default Plist
183 %
184 % HISTORY: 11-07-07 M Hewitson
185 % Creation.
186 %
187 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
188
189 function plout = getDefaultPlist()
190 persistent pl;
191 if exist('pl', 'var')==0 || isempty(pl)
192 pl = buildplist();
193 end
194 plout = pl;
195 end
196
197 function pl = buildplist()
198
199 pl = plist();
200
201 % N
202 p = param({'N', 'The maximum number of lines to return.'}, {1, {10}, paramValue.OPTIONAL});
203 pl.append(p);
204
205 % fsearch
206 p = param({'fsearch', 'The frequency search interval.'}, {1, {[0 1e10]}, paramValue.OPTIONAL});
207 pl.append(p);
208
209 % thresh
210 p = param({'thresh', 'A threshold to test normalised amplitude against. (A sort-of SNR threshold.)'}, {1, {2}, paramValue.OPTIONAL});
211 pl.append(p);
212
213 % BW
214 p = param({'bw', ['The bandwidth of the running median filter used to<br>'...
215 'estimate the noise-floor.']}, {1, {20}, paramValue.OPTIONAL});
216 pl.append(p);
217
218 % HC
219 p = param({'hc', 'The cutoff used to reject outliers (0-1).'}, {1, {0.8}, paramValue.OPTIONAL});
220 pl.append(p);
221
222
223 end
224
225 % PARAMETERS: N - max number of lines to return [default: 10]
226 % fsearch - freqeuncy search interval [default: all]
227 % thresh - a threshold to test normalised amplitude spectrum against
228 % [default: 2]
229 % bw - bandwidth over which to compute median [default: 20 samples]
230 % hc - percent of outliers to exclude from median estimation (0-1)
231 % [default: 0.8]