0001 function varargout = ltpda_nfest(varargin)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
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
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.html,v 1.14 2008/03/31 10:27:39 hewitson Exp $';
0048
0049
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
0068 if ~isempty(ps)
0069 pl = combine(ps, getDefaultPL());
0070 else
0071 pl = getDefaultPL();
0072 end
0073
0074
0075 bw = find(pl, 'bw');
0076 hc = find(pl, 'hc');
0077
0078
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.y, bw, hc);
0088
0089 else
0090 error('### I can only find lines in frequency-series AOs.');
0091 end
0092
0093
0094
0095
0096 d = set(d, 'y', nf);
0097 d = set(d, 'name', sprintf('nfest(%s)', d.name));
0098
0099
0100 h = history(ALGONAME, VERSION, pl, a.hist);
0101 h = set(h, 'invars', invars);
0102
0103
0104 b = ao(d, h);
0105
0106
0107
0108 if isempty(invars{j})
0109 n1 = a.name;
0110 else
0111 n1 = invars{j};
0112 end
0113 b = setnh(b, 'name', sprintf('nfest(%s)', n1));
0114
0115
0116 bo = [bo b];
0117
0118 end
0119
0120 varargout{1} = bo;
0121
0122
0123
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
0133 function nf = nfest(pxx, bw, hc)
0134
0135 N = length(pxx);
0136 nf = [];
0137
0138 for j=1:N
0139
0140
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
0148
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
0158