Home > m > sigproc > frequency_domain > ltpda_linedetect.m

ltpda_linedetect

PURPOSE ^

LTPDA_LINEDETECT find spectral lines in the ao/fsdata objects.

SYNOPSIS ^

function b = 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] 
 
 
 M Hewitson 14-05-07

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function b = 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 %
0016 % M Hewitson 14-05-07
0017 %
0018 %
0019 
0020 % capture input variable names
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 % Get inputs
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 % Combine plists
0050 if ~isempty(ps)
0051   pl = combine(ps, getDefaultPlist());
0052 else
0053   pl = getDefaultPlist();
0054 end
0055 
0056 % Get parameters from plist
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 % Loop over input AOs
0064 bo = [];
0065 for j=1:Na
0066 
0067   a = as(j);
0068   d = a.data;
0069   
0070   if isa(d, 'fsdata')
0071 
0072     % Make noise-floor estimate
0073     nf = ltpda_nfest(a, pl);
0074     
0075     % Make ratio
0076     r = a./nf;
0077     
0078     % find lines
0079     lines = findLines(d.xx, get(r.data, 'f'), get(r.data, 'xx'), thresh, N, fsearch);
0080 
0081     % Make output data
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   %------- Make output AO
0089   
0090   % create new output history
0091   h = history(ALGONAME, VERSION, pl, a.hist);
0092   h = set(h, 'invars', invars);
0093 
0094   % make output analysis object
0095   b = ao(fs, h);
0096 
0097   % set name
0098   % name for this object
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   % Add to output array
0107   bo = [bo b];
0108 
0109   
0110   
0111 end
0112 
0113 
0114 %--------------------------------------------------------------------------
0115 % find spectral lines
0116 function lines = findLines(axx, f, naxx, thresh, N, fsearch)
0117 
0118 % look for spectral lines
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       % find index when we have pmax
0134       fidx = pmaxi; %(find(naxx(1:idx(j))==pmax));
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 % Select largest peaks
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 % get default parameter list
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 % END

Generated on Fri 08-Jun-2007 16:09:11 by m2html © 2003