0001 function b = 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 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
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
0054 if ~isempty(ps)
0055 pl = combine(ps, getDefaultPlist());
0056 else
0057 pl = getDefaultPlist();
0058 end
0059
0060
0061 bw = find(pl, 'bw');
0062 hc = find(pl, 'hc');
0063
0064
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
0080
0081
0082 d = set(d, 'xx', nf);
0083 d = set(d, 'name', sprintf('nfest(%s)', d.name));
0084
0085
0086 h = history(ALGONAME, VERSION, pl, a.hist);
0087 h = set(h, 'invars', invars);
0088
0089
0090 b = ao(d, h);
0091
0092
0093
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
0102 bo = [bo b];
0103
0104 end
0105
0106
0107
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
0117 function nf = nfest(pxx, bw, hc)
0118
0119 N = length(pxx);
0120 nf = [];
0121
0122 for j=1:N
0123
0124
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
0132
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
0142