Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/firwhiten.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % FIRWHITEN whitens the input time-series by building an FIR whitening filter.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: FIRWHITEN whitens the input time-series by building an FIR
5 % whitening filter. The algorithm ultimately uses fir2() to
6 % build the whitening filter.
7 %
8 % ALGORITHM:
9 % 1) Make ASD of time-series
10 % 2) Perform running median to get noise-floor estimate (ao/smoother)
11 % 3) Invert noise-floor estimate
12 % 4) Call mfir() on noise-floor estimate to produce whitening filter
13 % 5) Filter data
14 %
15 % CALL: b = firwhiten(a, pl) % returns whitened time-series AOs
16 % [b, filts] = firwhiten(a, pl) % returns the mfir filters used
17 % [b, filts, nfs] = firwhiten(a, pl) % returns the noise-floor
18 % % estimates as fsdata AOs
19 %
20 %
21 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'firwhiten')">Parameters Description</a>
22 %
23 % VERSION: $Id: firwhiten.m,v 1.29 2011/11/11 15:21:19 luigi Exp $
24 %
25 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26
27 function varargout = firwhiten(varargin)
28
29 callerIsMethod = utils.helper.callerIsMethod;
30
31 % Check if this is a call for parameters
32 if utils.helper.isinfocall(varargin{:})
33 varargout{1} = getInfo(varargin{3});
34 return
35 end
36
37 import utils.const.*
38 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
39
40 % Collect input variable names
41 in_names = cell(size(varargin));
42 for ii = 1:nargin,in_names{ii} = inputname(ii);end
43
44 % Collect all AOs and plists
45 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
46 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
47
48 % Decide on a deep copy or a modify
49 bs = copy(as, nargout);
50 inhists = copy([as.hist],1);
51
52 % combine plists
53 pl = parse(pl, getDefaultPlist());
54
55 % Extract necessary parameters
56 iNfft = find(pl, 'Nfft');
57 bw = find(pl, 'bw');
58 hc = find(pl, 'hc');
59 swin = find(pl, 'win');
60 order = find(pl, 'order');
61 fwin = find(pl, 'FIRwin');
62 Ntaps = find(pl, 'Ntaps');
63
64 % Loop over input AOs
65 filts = [];
66 nfs = [];
67
68 for j=1:numel(bs)
69 if ~isa(bs(j).data, 'tsdata')
70 warning('!!! %s expects ao/tsdata objects. Skipping AO %s', mfilename, ao_invars{j});
71 bs(j) = [];
72 else
73 % get Nfft
74 if iNfft < 0 || isempty(iNfft)
75 Nfft = length(bs(j).data.y);
76 else
77 Nfft = iNfft;
78 end
79 utils.helper.msg(msg.PROC1, 'building spectrum');
80 % Make spectrum
81 axx = psd(bs(j), plist('Nfft', Nfft, 'Win', swin, 'Order', order, 'Scale', 'ASD'));
82 % make noise floor estimate
83 utils.helper.msg(msg.PROC1, 'estimating noise-floor');
84 nxx = smoother(axx, plist('width', bw, 'hc', hc));
85 % collect noise-floor estimates for output
86 nfs = [nfs nxx];
87 % invert and make weights
88 w = 1./nxx;
89 % Make mfir object
90 utils.helper.msg(msg.PROC1, 'building filter');
91 ff = mfir(w, plist('Win', fwin, 'N', Ntaps));
92 % collect filters for output
93 filts = [filts ff];
94 % Filter data
95 utils.helper.msg(msg.PROC1, 'filter data');
96 filter(bs(j), ff);
97 % Set name
98 bs(j).name = sprintf('firwhiten(%s)', ao_invars{j});
99 % add history
100 if ~callerIsMethod
101 bs(j).addHistory(getInfo('None'), pl, ao_invars(j), inhists(j));
102 end
103 % clear errors
104 bs(j).clearErrors;
105
106 end
107 end
108
109 % Any errors are meaningless after this process, so clear them on both
110 % axes.
111 bs.clearErrors(plist('axis', 'xy'));
112
113 % Set outputs
114 if nargout > 0
115 varargout{1} = bs;
116 end
117 if nargout > 1
118 varargout{2} = filts;
119 end
120 if nargout > 2
121 varargout{3} = nfs;
122 end
123 end
124
125 %--------------------------------------------------------------------------
126 % Get Info Object
127 %--------------------------------------------------------------------------
128 function ii = getInfo(varargin)
129 if nargin == 1 && strcmpi(varargin{1}, 'None')
130 sets = {};
131 pl = [];
132 else
133 sets = {'Default'};
134 pl = getDefaultPlist;
135 end
136 % Build info object
137 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: firwhiten.m,v 1.29 2011/11/11 15:21:19 luigi Exp $', sets, pl);
138 end
139
140 %--------------------------------------------------------------------------
141 % Get Default Plist
142 %--------------------------------------------------------------------------
143 function plout = getDefaultPlist()
144 persistent pl;
145 if exist('pl', 'var')==0 || isempty(pl)
146 pl = buildplist();
147 end
148 plout = pl;
149 end
150
151 function pl = buildplist()
152
153 pl = plist();
154
155 % Nfft
156 p = param({'Nfft', ['The number of points in the FFT used to estimate<br>'...
157 'the power spectrum. If unspecified, this is calculated as Ndata/4.']}, paramValue.DOUBLE_VALUE(-1));
158 pl.append(p);
159
160 % BW
161 p = param({'bw', ['The bandwidth of the running median filter used to<br>'...
162 'estimate the noise-floor.']}, {1, {20}, paramValue.OPTIONAL});
163 pl.append(p);
164
165 % HC
166 p = param({'hc', 'The cutoff used to reject outliers (0-1).'}, {1, {0.8}, paramValue.OPTIONAL});
167 pl.append(p);
168
169 % Win
170 p = param({'Win', 'Spectral window used in spectral estimation.'}, paramValue.WINDOW);
171 pl.append(p);
172
173 % Order
174 p = param({'Order',['The order of segment detrending:<ul>', ...
175 '<li>-1 - no detrending</li>', ...
176 '<li>0 - subtract mean</li>', ...
177 '<li>1 - subtract linear fit</li>', ...
178 '<li>N - subtract fit of polynomial, order N</li></ul>']}, paramValue.DETREND_ORDER);
179 pl.append(p);
180
181 % FIR win
182 p = param({'FIRwin', 'The window to use in the filter design.'}, paramValue.WINDOW);
183 pl.append(p);
184
185 % Ntaps
186 p = param({'Ntaps', 'The length of the FIR filter to build.'}, {1, {256}, paramValue.OPTIONAL});
187 pl.append(p);
188
189 end
190
191 % PARAMETERS:
192 %
193 % 'Ntaps' - the length of the FIR filter to build [default: 256].
194 % 'FIRwin' - the window to use in the filter design. Pass a
195 % specwin object of the desired type and of any length.
196 % [default: Hanning]
197 %
198 % parameters passed to ltpda_pwelch()
199 %
200 % 'Nfft' - The number of points in the FFT used to estimate
201 % the power spectrum.
202 % [default: Ndata/4]
203 % 'Win' - Spectral window used in spectral estimation.
204 % [default: Kaiser -150dB]
205 % 'Order' - order of segment detrending:
206 % -1 - no detrending
207 % 0 - subtract mean [default]
208 % 1 - subtract linear fit
209 % N - subtract fit of polynomial, order N
210 %
211 % (Segment overlap is taken from the window function.)
212 %
213 % parameters passed to ltpda_nfest()
214 %
215 % 'bw' - The bandwidth of the running median filter used to
216 % estimate the noise-floor.
217 % [default: 20 samples]
218 % 'hc' - The cutoff used to reject outliers (0-1)
219 % [default: 0.8]
220