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.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 DESCRIPTION: LTPDA_LINEDETECT find spectral lines in the ao/fsdata objects.

 CALL:        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]

 VERSION:    $Id: ltpda_linedetect.m,v 1.8 2007/11/08 10:24:54 hewitson Exp $

 HISTORY:    14-05-2007 M Hewitson
                Creation

 The following call returns a parameter list object that contains the
 default parameter values:

 >> pl = ltpda_linedetect('Params')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Mon 31-Mar-2008 13:54:54 by m2html © 2003