LINEDETECT find spectral lines in the ao/fsdata objects. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DESCRIPTION: LINEDETECT find spectral lines in the ao/fsdata objects. CALL: b = linedetect(a, pl) PARAMETERS: N - max number of lines to return [default: 10] fsearch - freqeuncy search interval [default: all] thresh - a threshold to test normalised amplitude spectrum against [default: 2] bw - bandwidth over which to compute median [default: 20 samples] hc - percent of outliers to exclude from median estimation (0-1) [default: 0.8] M-FILE INFO: Get information about this methods by calling >> ao.getInfo('linedetect') Get information about a specified set-plist by calling: >> ao.getInfo('linedetect', 'None') VERSION: $Id: linedetect.m,v 1.3 2008/09/05 11:05:29 ingo Exp $ HISTORY: 14-05-2007 M Hewitson Creation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0001 % LINEDETECT find spectral lines in the ao/fsdata objects. 0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0003 % 0004 % DESCRIPTION: LINEDETECT find spectral lines in the ao/fsdata objects. 0005 % 0006 % CALL: b = linedetect(a, pl) 0007 % 0008 % PARAMETERS: N - max number of lines to return [default: 10] 0009 % fsearch - freqeuncy search interval [default: all] 0010 % thresh - a threshold to test normalised amplitude spectrum against 0011 % [default: 2] 0012 % bw - bandwidth over which to compute median [default: 20 samples] 0013 % hc - percent of outliers to exclude from median estimation (0-1) 0014 % [default: 0.8] 0015 % 0016 % M-FILE INFO: Get information about this methods by calling 0017 % >> ao.getInfo('linedetect') 0018 % 0019 % Get information about a specified set-plist by calling: 0020 % >> ao.getInfo('linedetect', 'None') 0021 % 0022 % VERSION: $Id: linedetect.m,v 1.3 2008/09/05 11:05:29 ingo Exp $ 0023 % 0024 % HISTORY: 14-05-2007 M Hewitson 0025 % Creation 0026 % 0027 % 0028 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0029 0030 function varargout = linedetect(varargin) 0031 0032 % Check if this is a call for parameters 0033 if utils.helper.isinfocall(varargin{:}) 0034 varargout{1} = getInfo(varargin{3}); 0035 return 0036 end 0037 0038 if nargout == 0 0039 error('### cat cannot be used as a modifier. Please give an output variable.'); 0040 end 0041 0042 % Collect input variable names 0043 in_names = cell(size(varargin)); 0044 for ii = 1:nargin,in_names{ii} = inputname(ii);end 0045 0046 % Collect all AOs 0047 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names); 0048 [pli, pl_invars] = utils.helper.collect_objects(varargin(:), 'plist', in_names); 0049 0050 % Decide on a deep copy or a modify 0051 bs = copy(as, nargout); 0052 0053 Na = numel(bs); 0054 if isempty(bs) 0055 error('### Please input at least one AO.'); 0056 end 0057 0058 % Combine plists 0059 if ~isempty(pli) 0060 pl = combine(pli, getDefaultPlist()); 0061 else 0062 pl = getDefaultPlist(); 0063 end 0064 0065 % Get parameters from plist 0066 N = find(pl, 'N'); 0067 fsearch = find(pl, 'fsearch'); 0068 thresh = find(pl, 'thresh'); 0069 0070 % Loop over input AOs 0071 for j=1:Na 0072 if isa(bs(j).data, 'fsdata') 0073 0074 % Make noise-floor estimate 0075 nf = smoother(bs(j), pl); 0076 0077 % Make ratio 0078 r = bs(j)./nf; 0079 0080 % find lines 0081 lines = findLines(bs(j).data.getY, r.data.getX, r.data.getY, thresh, N, fsearch); 0082 0083 f = [lines(:).f]; 0084 y = [lines(:).a]; 0085 0086 % Keep the data shpare of the input AO 0087 if size(bs(j).data.y, 2) == 1 0088 f = f.'; 0089 y = y.'; 0090 end 0091 0092 % Make output data 0093 fs = fsdata(f, y, bs(j).data.fs); 0094 0095 else 0096 error('### I can only find lines in frequency-series AOs.'); 0097 end 0098 0099 %------- Make output AO 0100 0101 % make output analysis object 0102 bs(j).data = fs; 0103 0104 bs(j).setName(sprintf('lines(%s)', ao_invars{1})); 0105 bs(j).addHistory(getInfo('None'), pl, ao_invars(j), bs(j).hist); 0106 end 0107 0108 % return output 0109 varargout{1} = bs; 0110 end 0111 0112 %-------------------------------------------------------------------------- 0113 % find spectral lines 0114 function lines = findLines(ay, f, nay, thresh, N, fsearch) 0115 0116 % look for spectral lines 0117 l = 0; 0118 pmax = 0; 0119 pmaxi = 0; 0120 line = []; 0121 idx = find( f>=fsearch(1) & f<=fsearch(2) ); 0122 for j=1:length(idx) 0123 v = nay(idx(j)); 0124 if v > thresh 0125 if v > pmax 0126 pmax = v; 0127 pmaxi = idx(j); 0128 end 0129 else 0130 if pmax > 0 0131 % find index when we have pmax 0132 fidx = pmaxi; %(find(nay(1:idx(j))==pmax)); 0133 l = l+1; 0134 line(l).idx = fidx; 0135 line(l).f = f(fidx); 0136 line(l).a = ay(fidx); 0137 end 0138 pmax = 0; 0139 end 0140 end 0141 0142 % Select largest peaks 0143 lines = []; 0144 if ~isempty(line) 0145 [bl, lidx] = sort([line.a], 'descend'); 0146 lidxs = lidx(1:min([N length(lidx)])); 0147 lines = line(lidxs); 0148 disp(sprintf(' + found %d lines.', length([lines.f]))); 0149 else 0150 disp(' + found 0 lines.'); 0151 end 0152 end 0153 0154 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0155 % Local Functions % 0156 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0157 0158 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0159 % 0160 % FUNCTION: getInfo 0161 % 0162 % DESCRIPTION: Get Info Object 0163 % 0164 % HISTORY: 11-07-07 M Hewitson 0165 % Creation. 0166 % 0167 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0168 0169 function ii = getInfo(varargin) 0170 if nargin == 1 && strcmpi(varargin{1}, 'None') 0171 sets = {}; 0172 pl = []; 0173 else 0174 sets = {'Default'}; 0175 pl = getDefaultPlist; 0176 end 0177 % Build info object 0178 ii = minfo(mfilename, 'm', '', utils.const.categories.sigproc, '$Id: linedetect.m,v 1.3 2008/09/05 11:05:29 ingo Exp $', sets, pl); 0179 end 0180 0181 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0182 % 0183 % FUNCTION: getDefaultPlist 0184 % 0185 % DESCRIPTION: Get Default Plist 0186 % 0187 % HISTORY: 11-07-07 M Hewitson 0188 % Creation. 0189 % 0190 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0191 0192 function plo = getDefaultPlist() 0193 plo = plist('N', 10, 'fsearch', [0 1e10], 'thresh', 2, 'bw', 20, 'hc', 0.8); 0194 end 0195