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 b = 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]

 VERSION: $Id: ltpda_nfest.html,v 1.1 2007/06/08 14:15:10 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 b = 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 % VERSION: $Id: ltpda_nfest.html,v 1.1 2007/06/08 14:15:10 hewitson Exp $
0017 %
0018 % HISTORY: 14-05-07 M Hewitson
0019 %             Creation
0020 %
0021 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0022 
0023 
0024 % capture input variable names
0025 invars = {};
0026 for j=1:nargin
0027   if isa(varargin{j}, 'ao')
0028     invars = [invars cellstr(inputname(j))];
0029   end
0030 end
0031 
0032 ALGONAME = mfilename;
0033 VERSION  = '$Id: ltpda_nfest.html,v 1.1 2007/06/08 14:15:10 hewitson Exp $';
0034 
0035 % Get inputs
0036 as      = [];
0037 ps      = [];
0038 queries = [];
0039 for j=1:nargin
0040   if isa(varargin{j}, 'ao')
0041     as = [as varargin{j}];
0042   end
0043   if isa(varargin{j}, 'plist')
0044     ps = [ps varargin{j}];
0045   end
0046 end
0047 
0048 Na = length(as);
0049 if isempty(as)
0050   error('### Please input at least one AO.');
0051 end
0052 
0053 % Combine plists
0054 if ~isempty(ps)
0055   pl = combine(ps, getDefaultPlist());
0056 else
0057   pl = getDefaultPlist();
0058 end
0059 
0060 % Get parameters from plist
0061 bw      = find(pl, 'bw');
0062 hc      = find(pl, 'hc');
0063 
0064 % Loop over input AOs
0065 bo = [];
0066 for j=1:Na
0067 
0068   a = as(j);
0069   d = a.data;
0070 
0071   if isa(d, 'fsdata')
0072 
0073     nf = nfest(d.xx, bw, hc);
0074 
0075   else
0076     error('### I can only find lines in frequency-series AOs.');
0077   end
0078 
0079   %------- Make output AO
0080 
0081   % make fsdata
0082   d = set(d, 'xx', nf);
0083   d = set(d, 'name', sprintf('nfest(%s)', d.name));
0084 
0085   % create new output history
0086   h = history(ALGONAME, VERSION, pl, a.hist);
0087   h = set(h, 'invars', invars);
0088 
0089   % make output analysis object
0090   b = ao(d, h);
0091 
0092   % set name
0093   % name for this object
0094   if isempty(invars{j})
0095     n1 = a.name;
0096   else
0097     n1 = invars{j};
0098   end
0099   b = set(b, 'name', sprintf('nfest(%s)', n1));
0100 
0101   % Add to output array
0102   bo = [bo b];
0103 
0104 end
0105 
0106 %--------------------------------------------------------------------------
0107 % get default parameter list
0108 function pl = getDefaultPlist()
0109 
0110 pl = plist();
0111 pl = append(pl, param('bw', 20));
0112 pl = append(pl, param('hc', 0.8));
0113 
0114 
0115 %--------------------------------------------------------------------------
0116 % compute noise floor estimate
0117 function nf = nfest(pxx, bw, hc)
0118 
0119 N = length(pxx);
0120 nf = [];
0121 
0122 for j=1:N
0123 
0124   % Determine the interval we are looking in
0125   interval = j-bw/2:j+bw/2;
0126   idx = find(interval<=0);
0127   interval(idx)=1;
0128   idx = find(interval>N);
0129   interval(idx)=N;
0130 
0131   % calculate median value of interval
0132   % after throwing away outliers
0133   trial = sort(pxx(interval));
0134   b = round(hc*length(trial));
0135   nf(j) = median(trial(1:b));
0136 
0137 end
0138 
0139 nf = nf(:);
0140 
0141 % END
0142

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