1 |
function [adcp,cfg,ens]=rdradcp(name,varargin); |
---|
2 |
% RDRADCP Read (raw binary) RDI ADCP files, |
---|
3 |
% ADCP=RDRADCP(NAME) reads the raw binary RDI BB/Workhorse ADCP file NAME and |
---|
4 |
% puts all the relevant configuration and measured data into a data structure |
---|
5 |
% ADCP (which is self-explanatory). This program is designed for handling data |
---|
6 |
% recorded by moored instruments (primarily Workhorse-type but can also read |
---|
7 |
% Broadband) and then downloaded post-deployment. For vessel-mount data I |
---|
8 |
% usually make p-files (which integrate nav info and do coordinate transformations) |
---|
9 |
% and then use RDPADCP. |
---|
10 |
% |
---|
11 |
% This current version does have some handling of VMDAS and WINRIVER output |
---|
12 |
% files, but it is still 'beta'. |
---|
13 |
% |
---|
14 |
% [ADCP,CFG]=RDRADCP(...) returns configuration data in a |
---|
15 |
% separate data structure. |
---|
16 |
% |
---|
17 |
% Various options can be specified on input: |
---|
18 |
% [..]=RDRADCP(NAME,NUMAV) averages NUMAV ensembles together in the result. |
---|
19 |
% [..]=RDRADCP(NAME,NUMAV,NENS) reads only NENS ensembles (-1 for all). |
---|
20 |
% [..]=RDRADCP(NAME,NUMAV,[NFIRST NEND]) reads only the specified range |
---|
21 |
% of ensembles. This is useful if you want to get rid of bad data before/after |
---|
22 |
% the deployment period. |
---|
23 |
% |
---|
24 |
% Note - sometimes the ends of files are filled with garbage. In this case you may |
---|
25 |
% have to rerun things explicitly specifying how many records to read (or the |
---|
26 |
% last record to read). I don't handle bad data very well. |
---|
27 |
% |
---|
28 |
% - I don't read in absolutely every parameter stored in the binaries; |
---|
29 |
% just the ones that are 'most' useful. |
---|
30 |
% |
---|
31 |
% String parameter/option pairs can be added after these initial parameters: |
---|
32 |
% |
---|
33 |
% 'baseyear' : Base century for BB/v8WH firmware (default to 2000). |
---|
34 |
% |
---|
35 |
% 'despike' : [ 'no' | 'yes' | 3-element vector ] |
---|
36 |
% Controls ensemble averaging. With 'no' a simple mean is used |
---|
37 |
% (default). With 'yes' a mean is applied to all values that fall |
---|
38 |
% within a window around the median (giving some outlier rejection). |
---|
39 |
% This is useful for noisy data. Window sizes are [.3 .3 .3] m/s |
---|
40 |
% for [ horiz_vel vert_vel error_vel ] values. If you want to |
---|
41 |
% change these values, set 'despike' to the 3-element vector. |
---|
42 |
% |
---|
43 |
% R. Pawlowicz (rich@ocgy.ubc.ca) - 17/09/99 |
---|
44 |
|
---|
45 |
% R. Pawlowicz - 17/Oct/99 |
---|
46 |
% 5/july/00 - handled byte offsets (and mysterious 'extra" bytes) slightly better, Y2K |
---|
47 |
% 5/Oct/00 - bug fix - size of ens stayed 2 when NUMAV==1 due to initialization, |
---|
48 |
% hopefully this is now fixed. |
---|
49 |
% 10/Mar/02 - #bytes per record changes mysteriously, |
---|
50 |
% tried a more robust workaround. Guess that we have an extra |
---|
51 |
% 2 bytes if the record length is even? |
---|
52 |
% 28/Mar/02 - added more firmware-dependent changes to format; hopefully this |
---|
53 |
% works for everything now (put previous changes on firmer footing?) |
---|
54 |
% 30/Mar/02 - made cfg output more intuitive by decoding things. |
---|
55 |
% - An early version of WAVESMON and PARSE which split out this |
---|
56 |
% data from a wave recorder inserted an extra two bytes per record. |
---|
57 |
% I have removed the code to handle this but if you need it see line 509 |
---|
58 |
% 29/Nov/02 - A change in the bottom-track block for version 4.05 (very old!). |
---|
59 |
% 29/Jan/03 - Status block in v4.25 150khzBB two bytes short? |
---|
60 |
% 14/Oct/03 - Added code to at least 'ignore' WinRiver GPS blocks. |
---|
61 |
% 11/Nov/03 - VMDAS navigation block, added hooks to output |
---|
62 |
% navigation data. |
---|
63 |
|
---|
64 |
num_av=5; % Block filtering and decimation parameter (# ensembles to block together). |
---|
65 |
nens=-1; % Read all ensembles. |
---|
66 |
century=2000; % ADCP clock does not have century prior to firmware 16.05. |
---|
67 |
vels='no'; % Default to simple averaging |
---|
68 |
|
---|
69 |
lv=length(varargin); |
---|
70 |
if lv>=1 & ~isstr(varargin{1}), |
---|
71 |
num_av=varargin{1}; % Block filtering and decimation parameter (# ensembles to block together). |
---|
72 |
varargin(1)=[]; |
---|
73 |
lv=lv-1; |
---|
74 |
if lv>=1 & ~isstr(varargin{1}), |
---|
75 |
nens=varargin{1}; |
---|
76 |
varargin(1)=[]; |
---|
77 |
lv=lv-1; |
---|
78 |
end; |
---|
79 |
end; |
---|
80 |
|
---|
81 |
% Read optional args |
---|
82 |
while length(varargin)>0, |
---|
83 |
switch varargin{1}(1:3), |
---|
84 |
case 'bas', |
---|
85 |
century = varargin{2}; |
---|
86 |
case 'des', |
---|
87 |
if isstr(varargin{2}), |
---|
88 |
if strcmp(varargin{2},'no'), vels='no'; |
---|
89 |
else vels=[.3 .3 .3]; end; |
---|
90 |
else |
---|
91 |
vels=varargin{2}; |
---|
92 |
end; |
---|
93 |
otherwise, |
---|
94 |
error(['Unknown command line option ->' varargin{1}]); |
---|
95 |
end; |
---|
96 |
varargin([1 2])=[]; |
---|
97 |
end; |
---|
98 |
|
---|
99 |
|
---|
100 |
|
---|
101 |
|
---|
102 |
% Check file information first |
---|
103 |
|
---|
104 |
naminfo=dir(name); |
---|
105 |
|
---|
106 |
if isempty(naminfo), |
---|
107 |
fprintf('ERROR******* Can''t find file %s\n',name); |
---|
108 |
return; |
---|
109 |
end; |
---|
110 |
|
---|
111 |
fprintf('\nOpening file %s\n\n',name); |
---|
112 |
fd=fopen(name,'r','ieee-le'); |
---|
113 |
|
---|
114 |
% Read first ensemble to initialize parameters |
---|
115 |
|
---|
116 |
[ens,hdr,cfg]=rd_buffer(fd,-2); % Initialize and read first two records |
---|
117 |
fseek(fd,0,'bof'); % Rewind |
---|
118 |
|
---|
119 |
if (cfg.prog_ver<16.05 & cfg.prog_ver>5.999) | cfg.prog_ver<5.55, |
---|
120 |
fprintf('**************Assuming that the century begins year %d *********** \n\n',century); |
---|
121 |
else |
---|
122 |
century=0; % century included in clock. |
---|
123 |
end; |
---|
124 |
|
---|
125 |
dats=datenum(century+ens.rtc(1,:),ens.rtc(2,:),ens.rtc(3,:),ens.rtc(4,:),ens.rtc(5,:),ens.rtc(6,:)+ens.rtc(7,:)/100); |
---|
126 |
t_int=diff(dats); |
---|
127 |
fprintf('Record begins at %s\n',datestr(dats(1),0)); |
---|
128 |
fprintf('Ping interval appears to be %s\n',datestr(t_int,13)); |
---|
129 |
|
---|
130 |
|
---|
131 |
% Estimate number of records (since I don't feel like handling EOFs correctly, |
---|
132 |
% we just don't read that far!) |
---|
133 |
|
---|
134 |
|
---|
135 |
% Now, this is a puzzle - it appears that this is not necessary in |
---|
136 |
% a firmware v16.12 sent to me, and I can't find any example for |
---|
137 |
% which it *is* necessary so I'm not sure why its there. It could be |
---|
138 |
% a leftoever from dealing with the bad WAVESMON/PARSE problem (now |
---|
139 |
% fixed) that inserted extra bytes. |
---|
140 |
% ...So its out for now. |
---|
141 |
%if cfg.prog_ver>=16.05, extrabytes=2; else extrabytes=0; end; % Extra bytes |
---|
142 |
extrabytes=0; |
---|
143 |
|
---|
144 |
if length(nens)==1, |
---|
145 |
if nens==-1, |
---|
146 |
nens=fix(naminfo.bytes/(hdr.nbyte+2+extrabytes)); |
---|
147 |
fprintf('\nEstimating %d ensembles in this file\nReducing by a factor of %d\n',nens,num_av); |
---|
148 |
else |
---|
149 |
fprintf('\nReading %d ensembles in this file\nReducing by a factor of %d\n',nens,num_av); |
---|
150 |
end; |
---|
151 |
else |
---|
152 |
fprintf('\nReading ensembles %d-%d in this file\nReducing by a factor of %d\n',nens,num_av); |
---|
153 |
fseek(fd,(hdr.nbyte+2+extrabytes)*(nens(1)-1),'bof'); |
---|
154 |
nens=diff(nens)+1; |
---|
155 |
end; |
---|
156 |
|
---|
157 |
if num_av>1, |
---|
158 |
if isstr(vels), |
---|
159 |
fprintf('\n Simple mean used for ensemble averaging\n'); |
---|
160 |
else |
---|
161 |
fprintf('\n Averaging after outlier rejection with parameters [%f %f %f]\n',vels); |
---|
162 |
end; |
---|
163 |
end; |
---|
164 |
|
---|
165 |
% Number of records after averaging. |
---|
166 |
|
---|
167 |
n=fix(nens/num_av); |
---|
168 |
|
---|
169 |
% Structure to hold all ADCP data |
---|
170 |
% Note that I am not storing all the data contained in the raw binary file, merely |
---|
171 |
% things I think are useful. |
---|
172 |
|
---|
173 |
switch cfg.sourceprog, |
---|
174 |
case 'WINRIVER', |
---|
175 |
adcp=struct('name','adcp','config',cfg,'mtime',zeros(1,n),'number',zeros(1,n),'pitch',zeros(1,n),... |
---|
176 |
'roll',zeros(1,n),'heading',zeros(1,n),'pitch_std',zeros(1,n),... |
---|
177 |
'roll_std',zeros(1,n),'heading_std',zeros(1,n),'depth',zeros(1,n),... |
---|
178 |
'temperature',zeros(1,n),'salinity',zeros(1,n),... |
---|
179 |
'pressure',zeros(1,n),'pressure_std',zeros(1,n),... |
---|
180 |
'east_vel',zeros(cfg.n_cells,n),'north_vel',zeros(cfg.n_cells,n),'vert_vel',zeros(cfg.n_cells,n),... |
---|
181 |
'error_vel',zeros(cfg.n_cells,n),'corr',zeros(cfg.n_cells,4,n),... |
---|
182 |
'status',zeros(cfg.n_cells,4,n),'intens',zeros(cfg.n_cells,4,n),... |
---|
183 |
'bt_range',zeros(4,n),'bt_vel',zeros(4,n),... |
---|
184 |
'nav_longitude',zeros(1,n),'nav_latitude',zeros(1,n)); |
---|
185 |
case 'VMDAS', |
---|
186 |
adcp=struct('name','adcp','config',cfg,'mtime',zeros(1,n),'number',zeros(1,n),'pitch',zeros(1,n),... |
---|
187 |
'roll',zeros(1,n),'heading',zeros(1,n),'pitch_std',zeros(1,n),... |
---|
188 |
'roll_std',zeros(1,n),'heading_std',zeros(1,n),'depth',zeros(1,n),... |
---|
189 |
'temperature',zeros(1,n),'salinity',zeros(1,n),... |
---|
190 |
'pressure',zeros(1,n),'pressure_std',zeros(1,n),... |
---|
191 |
'east_vel',zeros(cfg.n_cells,n),'north_vel',zeros(cfg.n_cells,n),'vert_vel',zeros(cfg.n_cells,n),... |
---|
192 |
'error_vel',zeros(cfg.n_cells,n),'corr',zeros(cfg.n_cells,4,n),... |
---|
193 |
'status',zeros(cfg.n_cells,4,n),'intens',zeros(cfg.n_cells,4,n),... |
---|
194 |
'bt_range',zeros(4,n),'bt_vel',zeros(4,n),... |
---|
195 |
'nav_smtime',zeros(1,n),'nav_emtime',zeros(1,n),... |
---|
196 |
'nav_slongitude',zeros(1,n),'nav_elongitude',zeros(1,n),... |
---|
197 |
'nav_slatitude',zeros(1,n),'nav_elatitude',zeros(1,n)); |
---|
198 |
otherwise |
---|
199 |
adcp=struct('name','adcp','config',cfg,'mtime',zeros(1,n),'number',zeros(1,n),'pitch',zeros(1,n),... |
---|
200 |
'roll',zeros(1,n),'heading',zeros(1,n),'pitch_std',zeros(1,n),... |
---|
201 |
'roll_std',zeros(1,n),'heading_std',zeros(1,n),'depth',zeros(1,n),... |
---|
202 |
'temperature',zeros(1,n),'salinity',zeros(1,n),... |
---|
203 |
'pressure',zeros(1,n),'pressure_std',zeros(1,n),... |
---|
204 |
'east_vel',zeros(cfg.n_cells,n),'north_vel',zeros(cfg.n_cells,n),'vert_vel',zeros(cfg.n_cells,n),... |
---|
205 |
'error_vel',zeros(cfg.n_cells,n),'corr',zeros(cfg.n_cells,4,n),... |
---|
206 |
'status',zeros(cfg.n_cells,4,n),'intens',zeros(cfg.n_cells,4,n),... |
---|
207 |
'bt_range',zeros(4,n),'bt_vel',zeros(4,n)); |
---|
208 |
end; |
---|
209 |
|
---|
210 |
|
---|
211 |
% Calibration factors for backscatter data |
---|
212 |
|
---|
213 |
clear global ens |
---|
214 |
% Loop for all records |
---|
215 |
for k=1:n, |
---|
216 |
|
---|
217 |
% Gives display so you know something is going on... |
---|
218 |
|
---|
219 |
if rem(k,50)==0, fprintf('\n%d',k*num_av);end; |
---|
220 |
fprintf('.'); |
---|
221 |
|
---|
222 |
% Read an ensemble |
---|
223 |
|
---|
224 |
ens=rd_buffer(fd,num_av); |
---|
225 |
|
---|
226 |
|
---|
227 |
if ~isstruct(ens), % If aborting... |
---|
228 |
fprintf('Only %d records found..suggest re-running RDRADCP using this parameter\n',(k-1)*num_av); |
---|
229 |
fprintf('(If this message preceded by a POSSIBLE PROGRAM PROBLEM message, re-run using %d)\n',(k-1)*num_av-1); |
---|
230 |
break; |
---|
231 |
end; |
---|
232 |
|
---|
233 |
dats=datenum(century+ens.rtc(1,:),ens.rtc(2,:),ens.rtc(3,:),ens.rtc(4,:),ens.rtc(5,:),ens.rtc(6,:)+ens.rtc(7,:)/100); |
---|
234 |
adcp.mtime(k)=median(dats); |
---|
235 |
adcp.number(k) =ens.number(1); |
---|
236 |
adcp.heading(k) =mean(ens.heading); |
---|
237 |
adcp.pitch(k) =mean(ens.pitch); |
---|
238 |
adcp.roll(k) =mean(ens.roll); |
---|
239 |
adcp.heading_std(k) =mean(ens.heading_std); |
---|
240 |
adcp.pitch_std(k) =mean(ens.pitch_std); |
---|
241 |
adcp.roll_std(k) =mean(ens.roll_std); |
---|
242 |
adcp.depth(k) =mean(ens.depth); |
---|
243 |
adcp.temperature(k) =mean(ens.temperature); |
---|
244 |
adcp.salinity(k) =mean(ens.salinity); |
---|
245 |
adcp.pressure(k) =mean(ens.pressure); |
---|
246 |
adcp.pressure_std(k)=mean(ens.pressure_std); |
---|
247 |
|
---|
248 |
if isstr(vels), |
---|
249 |
adcp.east_vel(:,k) =nmean(ens.east_vel ,2); |
---|
250 |
adcp.north_vel(:,k) =nmean(ens.north_vel,2); |
---|
251 |
adcp.vert_vel(:,k) =nmean(ens.vert_vel ,2); |
---|
252 |
adcp.error_vel(:,k) =nmean(ens.error_vel,2); |
---|
253 |
else |
---|
254 |
adcp.east_vel(:,k) =nmedian(ens.east_vel ,vels(1),2); |
---|
255 |
adcp.north_vel(:,k) =nmedian(ens.north_vel,vels(1),2); |
---|
256 |
adcp.vert_vel(:,k) =nmedian(ens.vert_vel ,vels(2),2); |
---|
257 |
adcp.error_vel(:,k) =nmedian(ens.error_vel,vels(3),2); |
---|
258 |
end; |
---|
259 |
|
---|
260 |
adcp.corr(:,:,k) =nmean(ens.corr,3); % added correlation RKD 9/00 |
---|
261 |
adcp.status(:,:,k) =nmean(ens.status,3); |
---|
262 |
|
---|
263 |
adcp.intens(:,:,k) =nmean(ens.intens,3); |
---|
264 |
|
---|
265 |
adcp.bt_range(:,k) =nmean(ens.bt_range,2); |
---|
266 |
adcp.bt_vel(:,k) =nmean(ens.bt_vel,2); |
---|
267 |
|
---|
268 |
switch cfg.sourceprog, |
---|
269 |
case 'WINRIVER', |
---|
270 |
adcp.nav_longitude(k)=nmean(ens.slongitude); |
---|
271 |
adcp.nav_latitude(k)=nmean(ens.slatitude); |
---|
272 |
case 'VMDAS', |
---|
273 |
adcp.nav_smtime(k) =ens.smtime(1); |
---|
274 |
adcp.nav_emtime(k) =ens.emtime(end); |
---|
275 |
adcp.nav_slatitude(k)=ens.slatitude(1); |
---|
276 |
adcp.nav_elatitude(k)=ens.elatitude(end); |
---|
277 |
adcp.nav_slongitude(k)=ens.slongitude(1); |
---|
278 |
adcp.nav_elongitude(k)=ens.elongitude(end); |
---|
279 |
end; |
---|
280 |
end; |
---|
281 |
|
---|
282 |
fprintf('\n'); |
---|
283 |
fclose(fd); |
---|
284 |
|
---|
285 |
|
---|
286 |
%------------------------------------- |
---|
287 |
function hdr=rd_hdr(fd); |
---|
288 |
% Read config data |
---|
289 |
|
---|
290 |
cfgid=fread(fd,1,'uint16'); |
---|
291 |
if cfgid~=hex2dec('7F7F'), |
---|
292 |
error(['File ID is ' dec2hex(cfgid) ' not 7F7F - data corrupted or not a BB/WH raw file?']); |
---|
293 |
end; |
---|
294 |
|
---|
295 |
hdr=rd_hdrseg(fd); |
---|
296 |
|
---|
297 |
%------------------------------------- |
---|
298 |
function cfg=rd_fix(fd); |
---|
299 |
% Read config data |
---|
300 |
|
---|
301 |
cfgid=fread(fd,1,'uint16'); |
---|
302 |
if cfgid~=hex2dec('0000'), |
---|
303 |
warning(['Fixed header ID ' cfgid 'incorrect - data corrupted or not a BB/WH raw file?']); |
---|
304 |
end; |
---|
305 |
|
---|
306 |
cfg=rd_fixseg(fd); |
---|
307 |
|
---|
308 |
|
---|
309 |
|
---|
310 |
%-------------------------------------- |
---|
311 |
function [hdr,nbyte]=rd_hdrseg(fd); |
---|
312 |
% Reads a Header |
---|
313 |
|
---|
314 |
hdr.nbyte =fread(fd,1,'int16'); |
---|
315 |
fseek(fd,1,'cof'); |
---|
316 |
ndat=fread(fd,1,'int8'); |
---|
317 |
hdr.dat_offsets =fread(fd,ndat,'int16'); |
---|
318 |
nbyte=4+ndat*2; |
---|
319 |
|
---|
320 |
%------------------------------------- |
---|
321 |
function opt=getopt(val,varargin); |
---|
322 |
% Returns one of a list (0=first in varargin, etc.) |
---|
323 |
if val+1>length(varargin), |
---|
324 |
opt='unknown'; |
---|
325 |
else |
---|
326 |
opt=varargin{val+1}; |
---|
327 |
end; |
---|
328 |
|
---|
329 |
% |
---|
330 |
%------------------------------------- |
---|
331 |
function [cfg,nbyte]=rd_fixseg(fd); |
---|
332 |
% Reads the configuration data from the fixed leader |
---|
333 |
|
---|
334 |
%%disp(fread(fd,10,'uint8')) |
---|
335 |
%%fseek(fd,-10,'cof'); |
---|
336 |
|
---|
337 |
cfg.name='wh-adcp'; |
---|
338 |
cfg.sourceprog='instrument'; % default - depending on what data blocks are |
---|
339 |
% around we can modify this later in rd_buffer. |
---|
340 |
cfg.prog_ver =fread(fd,1,'uint8')+fread(fd,1,'uint8')/100; |
---|
341 |
|
---|
342 |
if fix(cfg.prog_ver)==4 | fix(cfg.prog_ver)==5, |
---|
343 |
cfg.name='bb-adcp'; |
---|
344 |
elseif fix(cfg.prog_ver)==8 | fix(cfg.prog_ver)==16, |
---|
345 |
cfg.name='wh-adcp'; |
---|
346 |
else |
---|
347 |
cfg.name='unrecognized firmware version' ; |
---|
348 |
end; |
---|
349 |
|
---|
350 |
config =fread(fd,2,'uint8'); % Coded stuff |
---|
351 |
cfg.config =[dec2base(config(2),2,8) '-' dec2base(config(1),2,8)]; |
---|
352 |
cfg.beam_angle =getopt(bitand(config(2),3),15,20,30); |
---|
353 |
cfg.beam_freq =getopt(bitand(config(1),7),75,150,300,600,1200,2400); |
---|
354 |
cfg.beam_pattern =getopt(bitand(config(1),8)==8,'concave','convex'); % 1=convex,0=concave |
---|
355 |
cfg.orientation =getopt(bitand(config(1),128)==128,'down','up'); % 1=up,0=down |
---|
356 |
cfg.simflag =getopt(fread(fd,1,'uint8'),'real','simulated'); % Flag for simulated data |
---|
357 |
fseek(fd,1,'cof'); |
---|
358 |
cfg.n_beams =fread(fd,1,'uint8'); |
---|
359 |
cfg.n_cells =fread(fd,1,'uint8'); |
---|
360 |
cfg.pings_per_ensemble=fread(fd,1,'uint16'); |
---|
361 |
cfg.cell_size =fread(fd,1,'uint16')*.01; % meters |
---|
362 |
cfg.blank =fread(fd,1,'uint16')*.01; % meters |
---|
363 |
cfg.prof_mode =fread(fd,1,'uint8'); % |
---|
364 |
cfg.corr_threshold =fread(fd,1,'uint8'); |
---|
365 |
cfg.n_codereps =fread(fd,1,'uint8'); |
---|
366 |
cfg.min_pgood =fread(fd,1,'uint8'); |
---|
367 |
cfg.evel_threshold =fread(fd,1,'uint16'); |
---|
368 |
cfg.time_between_ping_groups=sum(fread(fd,3,'uint8').*[60 1 .01]'); % seconds |
---|
369 |
coord_sys =fread(fd,1,'uint8'); % Lots of bit-mapped info |
---|
370 |
cfg.coord=dec2base(coord_sys,2,8); |
---|
371 |
cfg.coord_sys =getopt(bitand(bitshift(coord_sys,-3),3),'beam','instrument','ship','earth'); |
---|
372 |
cfg.use_pitchroll =getopt(bitand(coord_sys,4)==4,'no','yes'); |
---|
373 |
cfg.use_3beam =getopt(bitand(coord_sys,2)==2,'no','yes'); |
---|
374 |
cfg.bin_mapping =getopt(bitand(coord_sys,1)==1,'no','yes'); |
---|
375 |
cfg.xducer_misalign=fread(fd,1,'int16')*.01; % degrees |
---|
376 |
cfg.magnetic_var =fread(fd,1,'int16')*.01; % degrees |
---|
377 |
cfg.sensors_src =dec2base(fread(fd,1,'uint8'),2,8); |
---|
378 |
cfg.sensors_avail =dec2base(fread(fd,1,'uint8'),2,8); |
---|
379 |
cfg.bin1_dist =fread(fd,1,'uint16')*.01; % meters |
---|
380 |
cfg.xmit_pulse =fread(fd,1,'uint16')*.01; % meters |
---|
381 |
cfg.water_ref_cells=fread(fd,2,'uint8'); |
---|
382 |
cfg.fls_target_threshold =fread(fd,1,'uint8'); |
---|
383 |
fseek(fd,1,'cof'); |
---|
384 |
cfg.xmit_lag =fread(fd,1,'uint16')*.01; % meters |
---|
385 |
nbyte=40; |
---|
386 |
|
---|
387 |
if cfg.prog_ver>=8.14, % Added CPU serial number with v8.14 |
---|
388 |
cfg.serialnum =fread(fd,8,'uint8'); |
---|
389 |
nbyte=nbyte+8; |
---|
390 |
end; |
---|
391 |
|
---|
392 |
if cfg.prog_ver>=8.24, % Added 2 more bytes with v8.24 firmware |
---|
393 |
cfg.sysbandwidth =fread(fd,2,'uint8'); |
---|
394 |
nbyte=nbyte+2; |
---|
395 |
end; |
---|
396 |
|
---|
397 |
if cfg.prog_ver>=16.05, % Added 1 more bytes with v16.05 firmware |
---|
398 |
cfg.syspower =fread(fd,1,'uint8'); |
---|
399 |
nbyte=nbyte+1; |
---|
400 |
end; |
---|
401 |
|
---|
402 |
% It is useful to have this precomputed. |
---|
403 |
|
---|
404 |
cfg.ranges=cfg.bin1_dist+[0:cfg.n_cells-1]'*cfg.cell_size; |
---|
405 |
if cfg.orientation==1, cfg.ranges=-cfg.ranges; end |
---|
406 |
|
---|
407 |
|
---|
408 |
%----------------------------- |
---|
409 |
function [ens,hdr,cfg]=rd_buffer(fd,num_av); |
---|
410 |
|
---|
411 |
% To save it being re-initialized every time. |
---|
412 |
global ens hdr |
---|
413 |
|
---|
414 |
% A fudge to try and read files not handled quite right. |
---|
415 |
global FIXOFFSET SOURCE |
---|
416 |
|
---|
417 |
% If num_av<0 we are reading only 1 element and initializing |
---|
418 |
if num_av<0 | isempty(ens), |
---|
419 |
FIXOFFSET=0; |
---|
420 |
SOURCE=0; % 0=instrument, 1=VMDAS, 2=WINRIVER |
---|
421 |
n=abs(num_av); |
---|
422 |
pos=ftell(fd); |
---|
423 |
hdr=rd_hdr(fd); |
---|
424 |
cfg=rd_fix(fd); |
---|
425 |
fseek(fd,pos,'bof'); |
---|
426 |
clear global ens |
---|
427 |
global ens |
---|
428 |
|
---|
429 |
ens=struct('number',zeros(1,n),'rtc',zeros(7,n),'BIT',zeros(1,n),'ssp',zeros(1,n),'depth',zeros(1,n),'pitch',zeros(1,n),... |
---|
430 |
'roll',zeros(1,n),'heading',zeros(1,n),'temperature',zeros(1,n),'salinity',zeros(1,n),... |
---|
431 |
'mpt',zeros(1,n),'heading_std',zeros(1,n),'pitch_std',zeros(1,n),... |
---|
432 |
'roll_std',zeros(1,n),'adc',zeros(8,n),'error_status_wd',zeros(1,n),... |
---|
433 |
'pressure',zeros(1,n),'pressure_std',zeros(1,n),... |
---|
434 |
'east_vel',zeros(cfg.n_cells,n),'north_vel',zeros(cfg.n_cells,n),'vert_vel',zeros(cfg.n_cells,n),... |
---|
435 |
'error_vel',zeros(cfg.n_cells,n),'intens',zeros(cfg.n_cells,4,n),'percent',zeros(cfg.n_cells,4,n),... |
---|
436 |
'corr',zeros(cfg.n_cells,4,n),'status',zeros(cfg.n_cells,4,n),'bt_range',zeros(4,n),'bt_vel',zeros(4,n),... |
---|
437 |
'smtime',zeros(1,n),'emtime',zeros(1,n),'slatitude',zeros(1,n),... |
---|
438 |
'slongitude',zeros(1,n),'elatitude',zeros(1,n),'elongitude',zeros(1,n),... |
---|
439 |
'flags',zeros(1,n)); |
---|
440 |
num_av=abs(num_av); |
---|
441 |
end; |
---|
442 |
|
---|
443 |
k=0; |
---|
444 |
while k<num_av, |
---|
445 |
|
---|
446 |
|
---|
447 |
id1=dec2hex(fread(fd,1,'uint16')); |
---|
448 |
if ~strcmp(id1,'7F7F'), |
---|
449 |
if isempty(id1), % End of file |
---|
450 |
ens=-1; |
---|
451 |
return; |
---|
452 |
end; |
---|
453 |
error(['Not a workhorse/broadband file or bad data encountered: ->' id1]); |
---|
454 |
end; |
---|
455 |
startpos=ftell(fd)-2; % Starting position. |
---|
456 |
|
---|
457 |
|
---|
458 |
% Read the # data types. |
---|
459 |
[hdr,nbyte]=rd_hdrseg(fd); |
---|
460 |
byte_offset=nbyte+2; |
---|
461 |
%%disp(length(hdr.dat_offsets)) |
---|
462 |
% Read all the data types. |
---|
463 |
for n=1:length(hdr.dat_offsets), |
---|
464 |
|
---|
465 |
id=fread(fd,1,'uint16'); |
---|
466 |
%% fprintf('ID=%s\n',dec2hex(id,4)); |
---|
467 |
|
---|
468 |
% handle all the various segments of data. Note that since I read the IDs as a two |
---|
469 |
% byte number in little-endian order the high and low bytes are exchanged compared to |
---|
470 |
% the values given in the manual. |
---|
471 |
% |
---|
472 |
|
---|
473 |
switch dec2hex(id,4), |
---|
474 |
case '0000', % Fixed leader |
---|
475 |
[cfg,nbyte]=rd_fixseg(fd); |
---|
476 |
nbyte=nbyte+2; |
---|
477 |
|
---|
478 |
case '0080' % Variable Leader |
---|
479 |
k=k+1; |
---|
480 |
ens.number(k) =fread(fd,1,'uint16'); |
---|
481 |
ens.rtc(:,k) =fread(fd,7,'uint8'); |
---|
482 |
ens.number(k) =ens.number(k)+65536*fread(fd,1,'uint8'); |
---|
483 |
ens.BIT(k) =fread(fd,1,'uint16'); |
---|
484 |
ens.ssp(k) =fread(fd,1,'uint16'); |
---|
485 |
ens.depth(k) =fread(fd,1,'uint16')*.1; % meters |
---|
486 |
ens.heading(k) =fread(fd,1,'uint16')*.01; % degrees |
---|
487 |
ens.pitch(k) =fread(fd,1,'int16')*.01; % degrees |
---|
488 |
ens.roll(k) =fread(fd,1,'int16')*.01; % degrees |
---|
489 |
ens.salinity(k) =fread(fd,1,'int16'); % PSU |
---|
490 |
ens.temperature(k) =fread(fd,1,'int16')*.01; % Deg C |
---|
491 |
ens.mpt(k) =sum(fread(fd,3,'uint8').*[60 1 .01]'); % seconds |
---|
492 |
ens.heading_std(k) =fread(fd,1,'uint8'); % degrees |
---|
493 |
ens.pitch_std(k) =fread(fd,1,'int8')*.1; % degrees |
---|
494 |
ens.roll_std(k) =fread(fd,1,'int8')*.1; % degrees |
---|
495 |
ens.adc(:,k) =fread(fd,8,'uint8'); |
---|
496 |
nbyte=2+40; |
---|
497 |
|
---|
498 |
if strcmp(cfg.name,'bb-adcp'), |
---|
499 |
|
---|
500 |
if cfg.prog_ver>=5.55, |
---|
501 |
fseek(fd,15,'cof'); % 14 zeros and one byte for number WM4 bytes |
---|
502 |
cent=fread(fd,1,'uint8'); % possibly also for 5.55-5.58 but |
---|
503 |
ens.rtc(:,k)=fread(fd,7,'uint8'); % I have no data to test. |
---|
504 |
ens.rtc(1,k)=ens.rtc(1,k)+cent*100; |
---|
505 |
nbyte=nbyte+15+8; |
---|
506 |
end; |
---|
507 |
|
---|
508 |
elseif strcmp(cfg.name,'wh-adcp'), % for WH versions. |
---|
509 |
|
---|
510 |
ens.error_status_wd(k)=fread(fd,1,'uint32'); |
---|
511 |
nbyte=nbyte+4;; |
---|
512 |
|
---|
513 |
if cfg.prog_ver>=8.13, % Added pressure sensor stuff in 8.13 |
---|
514 |
fseek(fd,2,'cof'); |
---|
515 |
ens.pressure(k) =fread(fd,1,'uint32'); |
---|
516 |
ens.pressure_std(k) =fread(fd,1,'uint32'); |
---|
517 |
nbyte=nbyte+10; |
---|
518 |
end; |
---|
519 |
|
---|
520 |
if cfg.prog_ver>8.24, % Spare byte added 8.24 |
---|
521 |
fseek(fd,1,'cof'); |
---|
522 |
nbyte=nbyte+1; |
---|
523 |
end; |
---|
524 |
|
---|
525 |
if cfg.prog_ver>=16.05, % Added more fields with century in clock 16.05 |
---|
526 |
cent=fread(fd,1,'uint8'); |
---|
527 |
ens.rtc(:,k)=fread(fd,7,'uint8'); |
---|
528 |
ens.rtc(1,k)=ens.rtc(1,k)+cent*100; |
---|
529 |
nbyte=nbyte+8; |
---|
530 |
end; |
---|
531 |
end; |
---|
532 |
|
---|
533 |
case '0100', % Velocities |
---|
534 |
vels=fread(fd,[4 cfg.n_cells],'int16')'*.001; % m/s |
---|
535 |
ens.east_vel(:,k) =vels(:,1); |
---|
536 |
ens.north_vel(:,k)=vels(:,2); |
---|
537 |
ens.vert_vel(:,k) =vels(:,3); |
---|
538 |
ens.error_vel(:,k)=vels(:,4); |
---|
539 |
nbyte=2+4*cfg.n_cells*2; |
---|
540 |
|
---|
541 |
case '0200', % Correlations |
---|
542 |
ens.corr(:,:,k) =fread(fd,[4 cfg.n_cells],'uint8')'; |
---|
543 |
nbyte=2+4*cfg.n_cells; |
---|
544 |
|
---|
545 |
case '0300', % Echo Intensities |
---|
546 |
ens.intens(:,:,k) =fread(fd,[4 cfg.n_cells],'uint8')'; |
---|
547 |
nbyte=2+4*cfg.n_cells; |
---|
548 |
|
---|
549 |
case '0400', % Percent good |
---|
550 |
ens.percent(:,:,k) =fread(fd,[4 cfg.n_cells],'uint8')'; |
---|
551 |
nbyte=2+4*cfg.n_cells; |
---|
552 |
|
---|
553 |
case '0500', % Status |
---|
554 |
% Note in one case with a 4.25 firmware SC-BB, it seems like |
---|
555 |
% this block was actually two bytes short! |
---|
556 |
ens.status(:,:,k) =fread(fd,[4 cfg.n_cells],'uint8')'; |
---|
557 |
nbyte=2+4*cfg.n_cells; |
---|
558 |
|
---|
559 |
case '0600', % Bottom track |
---|
560 |
% In WINRIVER GPS data is tucked into here in odd ways, as long |
---|
561 |
% as GPS is enabled. |
---|
562 |
if SOURCE==2, |
---|
563 |
fseek(fd,2,'cof'); |
---|
564 |
long1=fread(fd,1,'uint16'); |
---|
565 |
fseek(fd,6,'cof'); |
---|
566 |
cfac=180/2^31; |
---|
567 |
ens.slatitude(k) =fread(fd,1,'int32')*cfac; |
---|
568 |
else |
---|
569 |
fseek(fd,14,'cof'); % Skip over a bunch of stuff |
---|
570 |
end; |
---|
571 |
ens.bt_range(:,k)=fread(fd,4,'uint16')*.01; % |
---|
572 |
ens.bt_vel(:,k) =fread(fd,4,'int16'); |
---|
573 |
if SOURCE==2, |
---|
574 |
fseek(fd,12+2,'cof'); |
---|
575 |
ens.slongitude(k)=(long1+65536*fread(fd,1,'uint16'))*cfac; |
---|
576 |
if ens.slongitude(k)>180, ens.slongitude(k)=ens.slongitude(k)-360; end; |
---|
577 |
fseek(fd,71-33-16,'cof'); |
---|
578 |
nbyte=2+68; |
---|
579 |
else |
---|
580 |
fseek(fd,71-33,'cof'); |
---|
581 |
nbyte=2+68; |
---|
582 |
end; |
---|
583 |
if cfg.prog_ver>=5.3, % Version 4.05 firmware seems to be missing these last 11 bytes. |
---|
584 |
fseek(fd,78-71,'cof'); |
---|
585 |
ens.bt_range(:,k)=ens.bt_range(:,k)+fread(fd,4,'uint8')*655.36; |
---|
586 |
nbyte=nbyte+11; |
---|
587 |
if cfg.prog_ver>=16, % RDI documentation claims these extra bytes were added in v 8.17 |
---|
588 |
fseek(fd,4,'cof'); % but they don't appear in my 8.33 data. |
---|
589 |
nbyte=nbyte+4; |
---|
590 |
end; |
---|
591 |
end; |
---|
592 |
|
---|
593 |
% The raw files produced by VMDAS contain a binary navigation data |
---|
594 |
% block. |
---|
595 |
|
---|
596 |
case '2000', % Something from VMDAS. |
---|
597 |
cfg.sourceprog='VMDAS'; |
---|
598 |
SOURCE=1; |
---|
599 |
utim =fread(fd,4,'uint8'); |
---|
600 |
mtime =datenum(utim(3)+utim(4)*256,utim(2),utim(1)); |
---|
601 |
ens.smtime(k) =mtime+fread(fd,1,'uint32')/8640000; |
---|
602 |
fseek(fd,4,'cof'); |
---|
603 |
cfac=180/2^31; |
---|
604 |
ens.slatitude(k) =fread(fd,1,'int32')*cfac; |
---|
605 |
ens.slongitude(k) =fread(fd,1,'int32')*cfac; |
---|
606 |
ens.emtime(k) =mtime+fread(fd,1,'uint32')/8640000; |
---|
607 |
ens.elatitude(k) =fread(fd,1,'int32')*cfac; |
---|
608 |
ens.elongitude(k) =fread(fd,1,'int32')*cfac; |
---|
609 |
fseek(fd,12,'cof'); |
---|
610 |
ens.flags(k) =fread(fd,1,'uint16'); |
---|
611 |
fseek(fd,30,'cof'); |
---|
612 |
nbyte=2+76; |
---|
613 |
|
---|
614 |
% The following blocks come from WINRIVER files, they aparently contain |
---|
615 |
% the raw NMEA data received from a serial port. |
---|
616 |
% |
---|
617 |
% Note that for WINRIVER files somewhat decoded data is also available |
---|
618 |
% tucked into the bottom track block. |
---|
619 |
|
---|
620 |
case '2100', % $xxDBT (Winriver addition) 38 |
---|
621 |
cfg.sourceprog='WINRIVER'; |
---|
622 |
SOURCE=2; |
---|
623 |
str=fread(fd,38,'uchar')'; |
---|
624 |
nbyte=2+38; |
---|
625 |
|
---|
626 |
case '2101', % $xxGGA (Winriver addition) 94 in maanual but 97 seems to work |
---|
627 |
cfg.sourceprog='WINRIVER'; |
---|
628 |
SOURCE=2; |
---|
629 |
str=fread(fd,97,'uchar')'; |
---|
630 |
nbyte=2+97; |
---|
631 |
% disp(setstr(str(1:80))); |
---|
632 |
|
---|
633 |
case '2102', % $xxVTG (Winriver addition) 45 |
---|
634 |
cfg.sourceprog='WINRIVER'; |
---|
635 |
SOURCE=2; |
---|
636 |
str=fread(fd,45,'uchar')'; |
---|
637 |
nbyte=2+45; |
---|
638 |
% disp(setstr(str)); |
---|
639 |
|
---|
640 |
case '2103', % $xxGSA (Winriver addition) 60 |
---|
641 |
cfg.sourceprog='WINRIVER'; |
---|
642 |
SOURCE=2; |
---|
643 |
str=fread(fd,60,'uchar')'; |
---|
644 |
% disp(setstr(str)); |
---|
645 |
nbyte=2+60; |
---|
646 |
|
---|
647 |
case '2104', %xxHDT or HDG (Winriver addition) 38 |
---|
648 |
cfg.sourceprog='WINRIVER'; |
---|
649 |
SOURCE=2; |
---|
650 |
str=fread(fd,38,'uchar')'; |
---|
651 |
% disp(setstr(str)); |
---|
652 |
nbyte=2+38; |
---|
653 |
|
---|
654 |
|
---|
655 |
|
---|
656 |
case '0701', % Number of good pings |
---|
657 |
fseek(fd,4*cfg.n_cells,'cof'); |
---|
658 |
nbyte=2+4*cfg.n_cells; |
---|
659 |
|
---|
660 |
case '0702', % Sum of squared velocities |
---|
661 |
fseek(fd,4*cfg.n_cells,'cof'); |
---|
662 |
nbyte=2+4*cfg.n_cells; |
---|
663 |
|
---|
664 |
case '0703', % Sum of velocities |
---|
665 |
fseek(fd,4*cfg.n_cells,'cof'); |
---|
666 |
nbyte=2+4*cfg.n_cells; |
---|
667 |
|
---|
668 |
% These blocks were implemented for 5-beam systems |
---|
669 |
|
---|
670 |
case '0A00', % Beam 5 velocity (not implemented) |
---|
671 |
fseek(fd,cfg.n_cells,'cof'); |
---|
672 |
nbyte=2+cfg.n_cells; |
---|
673 |
|
---|
674 |
case '0301', % Beam 5 Number of good pings (not implemented) |
---|
675 |
fseek(fd,cfg.n_cells,'cof'); |
---|
676 |
nbyte=2+cfg.n_cells; |
---|
677 |
|
---|
678 |
case '0302', % Beam 5 Sum of squared velocities (not implemented) |
---|
679 |
fseek(fd,cfg.n_cells,'cof'); |
---|
680 |
nbyte=2+cfg.n_cells; |
---|
681 |
|
---|
682 |
case '0303', % Beam 5 Sum of velocities (not implemented) |
---|
683 |
fseek(fd,cfg.n_cells,'cof'); |
---|
684 |
nbyte=2+cfg.n_cells; |
---|
685 |
|
---|
686 |
case '020C', % Ambient sound profile (not implemented) |
---|
687 |
fseek(fd,4,'cof'); |
---|
688 |
nbyte=2+4; |
---|
689 |
|
---|
690 |
otherwise, |
---|
691 |
|
---|
692 |
fprintf('Unrecognized ID code: %s\n',dec2hex(id,4)); |
---|
693 |
nbyte=2; |
---|
694 |
%% ens=-1; |
---|
695 |
%% return; |
---|
696 |
|
---|
697 |
|
---|
698 |
end; |
---|
699 |
|
---|
700 |
% here I adjust the number of bytes so I am sure to begin |
---|
701 |
% reading at the next valid offset. If everything is working right I shouldn't have |
---|
702 |
% to do this but every so often firware changes result in some differences. |
---|
703 |
|
---|
704 |
%%fprintf('#bytes is %d, original offset is %d\n',nbyte,byte_offset); |
---|
705 |
byte_offset=byte_offset+nbyte; |
---|
706 |
|
---|
707 |
if n<length(hdr.dat_offsets), |
---|
708 |
if hdr.dat_offsets(n+1)~=byte_offset, |
---|
709 |
fprintf('%s: Adjust location by %d\n',dec2hex(id,4),hdr.dat_offsets(n+1)-byte_offset); |
---|
710 |
fseek(fd,hdr.dat_offsets(n+1)-byte_offset,'cof'); |
---|
711 |
end; |
---|
712 |
byte_offset=hdr.dat_offsets(n+1); |
---|
713 |
end; |
---|
714 |
end; |
---|
715 |
|
---|
716 |
% Now at the end of the record we have two reserved bytes, followed |
---|
717 |
% by a two-byte checksum = 4 bytes to skip over. |
---|
718 |
|
---|
719 |
readbytes=ftell(fd)-startpos; |
---|
720 |
offset=(hdr.nbyte+2)-byte_offset; % The 2 is for the checksum |
---|
721 |
if offset ~=4 & FIXOFFSET==0, |
---|
722 |
fprintf('\n*****************************************************\n'); |
---|
723 |
fprintf('Adjust location by %d (readbytes=%d, hdr.nbyte=%d)\n',offset,readbytes,hdr.nbyte); |
---|
724 |
fprintf(' NOTE - THIS IS A PROGRAM PROBLEM, POSSIBLY FIXED BY A FUDGE\n'); |
---|
725 |
fprintf(' PLEASE REPORT TO rich@ocgy.ubc.ca WITH DETAILS!!\n'); |
---|
726 |
fprintf(' ATTEMPTING TO RECOVER...\n'); |
---|
727 |
fprintf('******************************************************\n'); |
---|
728 |
FIXOFFSET=offset-4; |
---|
729 |
end; |
---|
730 |
fseek(fd,4+FIXOFFSET,'cof'); |
---|
731 |
|
---|
732 |
% An early version of WAVESMON and PARSE contained a bug which stuck an additional two |
---|
733 |
% bytes in these files, but they really shouldn't be there |
---|
734 |
%if cfg.prog_ver>=16.05, |
---|
735 |
% fseek(fd,2,'cof'); |
---|
736 |
%end; |
---|
737 |
|
---|
738 |
end; |
---|
739 |
|
---|
740 |
% Blank out stuff bigger than error velocity |
---|
741 |
% big_err=abs(ens.error_vel)>.2; |
---|
742 |
big_err=0; |
---|
743 |
|
---|
744 |
% Blank out invalid data |
---|
745 |
ens.east_vel(ens.east_vel==-32.768 | big_err)=NaN; |
---|
746 |
ens.north_vel(ens.north_vel==-32.768 | big_err)=NaN; |
---|
747 |
ens.vert_vel(ens.vert_vel==-32.768 | big_err)=NaN; |
---|
748 |
ens.error_vel(ens.error_vel==-32.768 | big_err)=NaN; |
---|
749 |
|
---|
750 |
|
---|
751 |
|
---|
752 |
|
---|
753 |
%-------------------------------------- |
---|
754 |
function y=nmedian(x,window,dim); |
---|
755 |
% Copied from median but with handling of NaN different. |
---|
756 |
|
---|
757 |
if nargin==2, |
---|
758 |
dim = min(find(size(x)~=1)); |
---|
759 |
if isempty(dim), dim = 1; end |
---|
760 |
end |
---|
761 |
|
---|
762 |
siz = [size(x) ones(1,dim-ndims(x))]; |
---|
763 |
n = size(x,dim); |
---|
764 |
|
---|
765 |
% Permute and reshape so that DIM becomes the row dimension of a 2-D array |
---|
766 |
perm = [dim:max(length(size(x)),dim) 1:dim-1]; |
---|
767 |
x = reshape(permute(x,perm),n,prod(siz)/n); |
---|
768 |
|
---|
769 |
% Sort along first dimension |
---|
770 |
x = sort(x,1); |
---|
771 |
[n1,n2]=size(x); |
---|
772 |
|
---|
773 |
if n1==1, |
---|
774 |
y=x; |
---|
775 |
else |
---|
776 |
if n2==1, |
---|
777 |
kk=sum(finite(x),1); |
---|
778 |
if kk>0, |
---|
779 |
x1=x(max(fix(kk/2),1)); |
---|
780 |
x2=x(max(ceil(kk/2),1)); |
---|
781 |
x(abs(x-(x1+x2)/2)>window)=NaN; |
---|
782 |
end; |
---|
783 |
x = sort(x,1); |
---|
784 |
kk=sum(finite(x),1); |
---|
785 |
x(isnan(x))=0; |
---|
786 |
y=NaN; |
---|
787 |
if kk>0, |
---|
788 |
y=sum(x)/kk; |
---|
789 |
end; |
---|
790 |
else |
---|
791 |
kk=sum(finite(x),1); |
---|
792 |
ll=kk<n1-2; |
---|
793 |
kk(ll)=0;x(:,ll)=NaN; |
---|
794 |
x1=x(max(fix(kk/2),1)+[0:n2-1]*n1); |
---|
795 |
x2=x(max(ceil(kk/2),1)+[0:n2-1]*n1); |
---|
796 |
|
---|
797 |
x(abs(x-ones(n1,1)*(x1+x2)/2)>window)=NaN; |
---|
798 |
x = sort(x,1); |
---|
799 |
kk=sum(finite(x),1); |
---|
800 |
x(isnan(x))=0; |
---|
801 |
y=NaN+ones(1,n2); |
---|
802 |
if any(kk), |
---|
803 |
y(kk>0)=sum(x(:,kk>0))./kk(kk>0); |
---|
804 |
end; |
---|
805 |
end; |
---|
806 |
end; |
---|
807 |
|
---|
808 |
% Permute and reshape back |
---|
809 |
siz(dim) = 1; |
---|
810 |
y = ipermute(reshape(y,siz(perm)),perm); |
---|
811 |
|
---|
812 |
%-------------------------------------- |
---|
813 |
function y=nmean(x,dim); |
---|
814 |
% R_NMEAN Computes the mean of matrix ignoring NaN |
---|
815 |
% values |
---|
816 |
% R_NMEAN(X,DIM) takes the mean along the dimension DIM of X. |
---|
817 |
% |
---|
818 |
|
---|
819 |
kk=finite(x); |
---|
820 |
x(~kk)=0; |
---|
821 |
|
---|
822 |
if nargin==1, |
---|
823 |
% Determine which dimension SUM will use |
---|
824 |
dim = min(find(size(x)~=1)); |
---|
825 |
if isempty(dim), dim = 1; end |
---|
826 |
end; |
---|
827 |
|
---|
828 |
if dim>length(size(x)), |
---|
829 |
y=x; % For matlab 5.0 only!!! Later versions have a fixed 'sum' |
---|
830 |
else |
---|
831 |
ndat=sum(kk,dim); |
---|
832 |
indat=ndat==0; |
---|
833 |
ndat(indat)=1; % If there are no good data then it doesn't matter what |
---|
834 |
% we average by - and this avoid div-by-zero warnings. |
---|
835 |
|
---|
836 |
y = sum(x,dim)./ndat; |
---|
837 |
y(indat)=NaN; |
---|
838 |
end; |
---|
839 |
|
---|
840 |
|
---|
841 |
|
---|
842 |
|
---|
843 |
|
---|
844 |
|
---|
845 |
|
---|
846 |
|
---|
847 |
|
---|
848 |
|
---|
849 |
|
---|
850 |
|
---|
851 |
|
---|
852 |
|
---|
853 |
|
---|
854 |
|
---|
855 |
|
---|
856 |
|
---|
857 |
|
---|
858 |
|
---|
859 |
|
---|
860 |
|
---|
861 |
|
---|
862 |
|
---|
863 |
|
---|
864 |
|
---|