Home > m > sigproc > frequency_domain > ltpda_linedetect.m

ltpda_linedetect

PURPOSE ^

LTPDA_LINEDETECT find spectral lines in the ao/fsdata objects.

SYNOPSIS ^

function varargout = ltpda_linedetect(varargin)

DESCRIPTION ^

 LTPDA_LINEDETECT find spectral lines in the ao/fsdata objects.
 
 Usage b = ltpda_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] 
 
 % The following call returns a parameter list object that contains the
 default parameter values:
 
 >> pl = ltpda_linedetect('Params')
 
 
 M Hewitson 14-05-07

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function varargout = ltpda_linedetect(varargin)
0002 
0003 % LTPDA_LINEDETECT find spectral lines in the ao/fsdata objects.
0004 %
0005 % Usage b = ltpda_linedetect(a, pl)
0006 %
0007 % Parameters:
0008 %   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) [default: 0.8]
0014 %
0015 % % The following call returns a parameter list object that contains the
0016 % default parameter values:
0017 %
0018 % >> pl = ltpda_linedetect('Params')
0019 %
0020 %
0021 % M Hewitson 14-05-07
0022 %
0023 %
0024 
0025 % Check if this is a call for parameters
0026 if nargin == 1
0027   in = char(varargin{1});
0028   if strcmp(in, 'Params')
0029     varargout{1} = getDefaultPL();
0030     return
0031   end
0032 end
0033 
0034 % capture input variable names
0035 invars = {};
0036 for j=1:nargin
0037   if isa(varargin{j}, 'ao')
0038     invars = [invars cellstr(inputname(j))];
0039   end
0040 end
0041 
0042 ALGONAME = mfilename;
0043 VERSION  = '$Id: ltpda_linedetect.m,v 1.5 2007/06/12 12:32:14 hewitson Exp $';
0044 
0045 % Get inputs
0046 as      = [];
0047 ps      = [];
0048 queries = [];
0049 for j=1:nargin
0050   if isa(varargin{j}, 'ao')
0051     as = [as varargin{j}];
0052   end
0053   if isa(varargin{j}, 'plist')
0054     ps = [ps varargin{j}];
0055   end
0056 end
0057 
0058 Na = length(as);
0059 if isempty(as)
0060   error('### Please input at least one AO.');
0061 end
0062 
0063 % Combine plists
0064 if ~isempty(ps)
0065   pl = combine(ps, getDefaultPL());
0066 else
0067   pl = getDefaultPL();
0068 end
0069 
0070 % Get parameters from plist
0071 N       = find(pl, 'N');
0072 fsearch = find(pl, 'fsearch');
0073 thresh  = find(pl, 'thresh');
0074 bw      = find(pl, 'bw');
0075 hc      = find(pl, 'hc');
0076 
0077 % Loop over input AOs
0078 bo = [];
0079 for j=1:Na
0080 
0081   a = as(j);
0082   d = a.data;
0083   
0084   if isa(d, 'fsdata')
0085 
0086     % Make noise-floor estimate
0087     nf = ltpda_nfest(a, pl);
0088     
0089     % Make ratio
0090     r = a./nf;
0091     
0092     % find lines
0093     lines = findLines(d.xx, get(r.data, 'f'), get(r.data, 'xx'), thresh, N, fsearch);
0094 
0095     % Make output data
0096     fs = fsdata([lines(:).f], [lines(:).a], d.fs);
0097     
0098   else
0099     error('### I can only find lines in frequency-series AOs.');
0100   end
0101   
0102   %------- Make output AO
0103   
0104   % create new output history
0105   h = history(ALGONAME, VERSION, pl, a.hist);
0106   h = set(h, 'invars', invars);
0107 
0108   % make output analysis object
0109   b = ao(fs, h);
0110 
0111   % set name
0112   % name for this object
0113   if isempty(invars{j})
0114     n1 = a.name;
0115   else
0116     n1 = invars{j};
0117   end
0118   b = set(b, 'name', sprintf('lines(%s)', n1));
0119 
0120   % Add to output array
0121   bo = [bo b];
0122   
0123 end
0124 
0125 % return output
0126 varargout{1} = bo;
0127 
0128 %--------------------------------------------------------------------------
0129 % find spectral lines
0130 function lines = findLines(axx, f, naxx, thresh, N, fsearch)
0131 
0132 % look for spectral lines
0133 l       = 0;
0134 pmax    = 0;
0135 pmaxi   = 0;
0136 line    = [];
0137 idx     = find( f>=fsearch(1) & f<=fsearch(2) );
0138 for j=1:length(idx)
0139   v = naxx(idx(j));
0140   if v > thresh
0141     if v > pmax
0142       pmax  = v;
0143       pmaxi = idx(j);
0144     end
0145   else
0146     if pmax > 0
0147       % find index when we have pmax
0148       fidx = pmaxi; %(find(naxx(1:idx(j))==pmax));
0149       l = l+1;
0150       line(l).idx = fidx;
0151       line(l).f   = f(fidx);
0152       line(l).a   = axx(fidx);
0153     end
0154     pmax = 0;
0155   end
0156 end
0157 
0158 % Select largest peaks
0159 lines = [];
0160 if ~isempty(line)
0161   [bl, lidx] = sort([line.a], 'descend');
0162   lidxs = lidx(1:min([N length(lidx)]));
0163   lines = line(lidxs);
0164   disp(sprintf('   + found %d lines.', length([lines.f])));
0165 else
0166   disp('   + found 0 lines.');
0167 end
0168 
0169 
0170 
0171 %--------------------------------------------------------------------------
0172 % get default parameter list
0173 function pl = getDefaultPL()
0174 
0175 pl = plist();
0176 pl = append(pl, param('N', 10));
0177 pl = append(pl, param('fsearch', [0 1e10]));
0178 pl = append(pl, param('thresh', 2));
0179 pl = append(pl, param('bw', 20));
0180 pl = append(pl, param('hc', 0.8));
0181 
0182 
0183 
0184 
0185 % END

Generated on Mon 02-Jul-2007 12:19:41 by m2html © 2003