0001 function b = ltpda_linedetect(varargin)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021 invars = {};
0022 for j=1:nargin
0023 if isa(varargin{j}, 'ao')
0024 invars = [invars cellstr(inputname(j))];
0025 end
0026 end
0027
0028 ALGONAME = mfilename;
0029 VERSION = '$Id: ltpda_linedetect.html,v 1.1 2007/06/08 14:15:10 hewitson Exp $';
0030
0031
0032 as = [];
0033 ps = [];
0034 queries = [];
0035 for j=1:nargin
0036 if isa(varargin{j}, 'ao')
0037 as = [as varargin{j}];
0038 end
0039 if isa(varargin{j}, 'plist')
0040 ps = [ps varargin{j}];
0041 end
0042 end
0043
0044 Na = length(as);
0045 if isempty(as)
0046 error('### Please input at least one AO.');
0047 end
0048
0049
0050 if ~isempty(ps)
0051 pl = combine(ps, getDefaultPlist());
0052 else
0053 pl = getDefaultPlist();
0054 end
0055
0056
0057 N = find(pl, 'N');
0058 fsearch = find(pl, 'fsearch');
0059 thresh = find(pl, 'thresh');
0060 bw = find(pl, 'bw');
0061 hc = find(pl, 'hc');
0062
0063
0064 bo = [];
0065 for j=1:Na
0066
0067 a = as(j);
0068 d = a.data;
0069
0070 if isa(d, 'fsdata')
0071
0072
0073 nf = ltpda_nfest(a, pl);
0074
0075
0076 r = a./nf;
0077
0078
0079 lines = findLines(d.xx, get(r.data, 'f'), get(r.data, 'xx'), thresh, N, fsearch);
0080
0081
0082 fs = fsdata([lines(:).f], [lines(:).a], d.fs);
0083
0084 else
0085 error('### I can only find lines in frequency-series AOs.');
0086 end
0087
0088
0089
0090
0091 h = history(ALGONAME, VERSION, pl, a.hist);
0092 h = set(h, 'invars', invars);
0093
0094
0095 b = ao(fs, h);
0096
0097
0098
0099 if isempty(invars{j})
0100 n1 = a.name;
0101 else
0102 n1 = invars{j};
0103 end
0104 b = set(b, 'name', sprintf('lines(%s)', n1));
0105
0106
0107 bo = [bo b];
0108
0109
0110
0111 end
0112
0113
0114
0115
0116 function lines = findLines(axx, f, naxx, thresh, N, fsearch)
0117
0118
0119 l = 0;
0120 pmax = 0;
0121 pmaxi = 0;
0122 line = [];
0123 idx = find( f>=fsearch(1) & f<=fsearch(2) );
0124 for j=1:length(idx)
0125 v = naxx(idx(j));
0126 if v > thresh
0127 if v > pmax
0128 pmax = v;
0129 pmaxi = idx(j);
0130 end
0131 else
0132 if pmax > 0
0133
0134 fidx = pmaxi;
0135 l = l+1;
0136 line(l).idx = fidx;
0137 line(l).f = f(fidx);
0138 line(l).a = axx(fidx);
0139 end
0140 pmax = 0;
0141 end
0142 end
0143
0144
0145 lines = [];
0146 if ~isempty(line)
0147 [bl, lidx] = sort([line.a], 'descend');
0148 lidxs = lidx(1:min([N length(lidx)]));
0149 lines = line(lidxs);
0150 disp(sprintf(' + found %d lines.', length([lines.f])));
0151 else
0152 disp(' + found 0 lines.');
0153 end
0154
0155
0156
0157
0158
0159 function pl = getDefaultPlist()
0160
0161 pl = plist();
0162 pl = append(pl, param('N', 10));
0163 pl = append(pl, param('fsearch', [0 1e10]));
0164 pl = append(pl, param('thresh', 2));
0165 pl = append(pl, param('bw', 20));
0166 pl = append(pl, param('hc', 0.8));
0167
0168
0169
0170
0171