Home > m > sigproc > frequency_domain > ltpda_nfest.m

ltpda_nfest

PURPOSE ^

LTPDA_NFEST Estimate the noise floor of the given frequency-series AO.

SYNOPSIS ^

function varargout = ltpda_nfest(varargin)

DESCRIPTION ^

 LTPDA_NFEST Estimate the noise floor of the given frequency-series AO.

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

 DESCRIPTION: LTPDA_NFEST Estimate the noise floor of the given
              frequency-series AO. This essentially implements a running
              median filter with some outlier selection.

 CALL:        b = ltpda_nfest(a, pl)

 PARAMETERS:
   bw  - the bandwidth of the noise floor estimate [default: 20 samples]
   hc  - a cutoff to throw away outliers (0-1)  [default: 0.8]

 The following call returns a parameter list object that contains the
 default parameter values:
 
 >> pl = ltpda_nfest('Params')
 
 
 VERSION: $Id: ltpda_nfest.m,v 1.5 2007/06/12 12:32:14 hewitson Exp $

 HISTORY: 14-05-07 M Hewitson
             Creation

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function varargout = ltpda_nfest(varargin)
0002 % LTPDA_NFEST Estimate the noise floor of the given frequency-series AO.
0003 %
0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0005 %
0006 % DESCRIPTION: LTPDA_NFEST Estimate the noise floor of the given
0007 %              frequency-series AO. This essentially implements a running
0008 %              median filter with some outlier selection.
0009 %
0010 % CALL:        b = ltpda_nfest(a, pl)
0011 %
0012 % PARAMETERS:
0013 %   bw  - the bandwidth of the noise floor estimate [default: 20 samples]
0014 %   hc  - a cutoff to throw away outliers (0-1)  [default: 0.8]
0015 %
0016 % The following call returns a parameter list object that contains the
0017 % default parameter values:
0018 %
0019 % >> pl = ltpda_nfest('Params')
0020 %
0021 %
0022 % VERSION: $Id: ltpda_nfest.m,v 1.5 2007/06/12 12:32:14 hewitson Exp $
0023 %
0024 % HISTORY: 14-05-07 M Hewitson
0025 %             Creation
0026 %
0027 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0028 
0029 % Check if this is a call for parameters
0030 if nargin == 1
0031   in = char(varargin{1});
0032   if strcmp(in, 'Params')
0033     varargout{1} = getDefaultPL();
0034     return
0035   end
0036 end
0037 
0038 % capture input variable names
0039 invars = {};
0040 for j=1:nargin
0041   if isa(varargin{j}, 'ao')
0042     invars = [invars cellstr(inputname(j))];
0043   end
0044 end
0045 
0046 ALGONAME = mfilename;
0047 VERSION  = '$Id: ltpda_nfest.m,v 1.5 2007/06/12 12:32:14 hewitson Exp $';
0048 
0049 % Get inputs
0050 as      = [];
0051 ps      = [];
0052 queries = [];
0053 for j=1:nargin
0054   if isa(varargin{j}, 'ao')
0055     as = [as varargin{j}];
0056   end
0057   if isa(varargin{j}, 'plist')
0058     ps = [ps varargin{j}];
0059   end
0060 end
0061 
0062 Na = length(as);
0063 if isempty(as)
0064   error('### Please input at least one AO.');
0065 end
0066 
0067 % Combine plists
0068 if ~isempty(ps)
0069   pl = combine(ps, getDefaultPL());
0070 else
0071   pl = getDefaultPL();
0072 end
0073 
0074 % Get parameters from plist
0075 bw      = find(pl, 'bw');
0076 hc      = find(pl, 'hc');
0077 
0078 % Loop over input AOs
0079 bo = [];
0080 for j=1:Na
0081 
0082   a = as(j);
0083   d = a.data;
0084 
0085   if isa(d, 'fsdata')
0086 
0087     nf = nfest(d.xx, bw, hc);
0088 
0089   else
0090     error('### I can only find lines in frequency-series AOs.');
0091   end
0092 
0093   %------- Make output AO
0094 
0095   % make fsdata
0096   d = set(d, 'xx', nf);
0097   d = set(d, 'name', sprintf('nfest(%s)', d.name));
0098 
0099   % create new output history
0100   h = history(ALGONAME, VERSION, pl, a.hist);
0101   h = set(h, 'invars', invars);
0102 
0103   % make output analysis object
0104   b = ao(d, h);
0105 
0106   % set name
0107   % name for this object
0108   if isempty(invars{j})
0109     n1 = a.name;
0110   else
0111     n1 = invars{j};
0112   end
0113   b = set(b, 'name', sprintf('nfest(%s)', n1));
0114 
0115   % Add to output array
0116   bo = [bo b];
0117 
0118 end
0119 
0120 varargout{1} = bo;
0121 
0122 %--------------------------------------------------------------------------
0123 % get default parameter list
0124 function pl = getDefaultPL()
0125 
0126 pl = plist();
0127 pl = append(pl, param('bw', 20));
0128 pl = append(pl, param('hc', 0.8));
0129 
0130 
0131 %--------------------------------------------------------------------------
0132 % compute noise floor estimate
0133 function nf = nfest(pxx, bw, hc)
0134 
0135 N = length(pxx);
0136 nf = [];
0137 
0138 for j=1:N
0139 
0140   % Determine the interval we are looking in
0141   interval = j-bw/2:j+bw/2;
0142   idx = find(interval<=0);
0143   interval(idx)=1;
0144   idx = find(interval>N);
0145   interval(idx)=N;
0146 
0147   % calculate median value of interval
0148   % after throwing away outliers
0149   trial = sort(pxx(interval));
0150   b = round(hc*length(trial));
0151   nf(j) = median(trial(1:b));
0152 
0153 end
0154 
0155 nf = nf(:);
0156 
0157 % END
0158

Generated on Mon 03-Sep-2007 12:12:34 by m2html © 2003