Dspike.m
Jump to navigation
Jump to search
function dspike()
%created by Erik Stassinos
%3/16/10
%%the infile will will be the whole lcd file
%despike function with hard coded trend of .01 and flag of 15
%first difference methode is used for despike processing
%if difference between data value "a" and "b" is greater than .01 than the
%new field in line a will contain the average of 15 pts surrounding "a"
%as same as old processing, running dspike twice will be hard coded as well
%as derived parameters d-a650, d-d-a650 which is in bbop scripts
%list_sdespike
%cd(pb)
id = fopen('list_ac9');
if id == -1
id = fopen('list_hs6')
if id ~= -1
hs6 =1;
else
disp('Cannot find list file. Exiting')
return
end
end
h1=0;
while feof(id) ~=1
tu = fgets(id);
h1 = h1+1; %to get # loop iterations
end
frewind(id);
for i = 1:h1
%clear h1;
r = fgets(id);
length(r);
t = textscan(r,'%s');
p = char(t{1,:}); %p is the name of the lcd file being processed from list_ac9
[in,mer,dep] = index1(p); %has locations in data for processing
%-now need to get to locations
%%------------------------main despike-----------------------------------
linecount = 0;
file = fopen(p);
trend = .01; %<=preset from original processing
flag = 15; %<=preset from original processing
% % --------------------splits upcasts and downcasts %hard coded 4 casts
fflag =0;
endflag =0;
if strncmp(p,'x',1) ==1
endflag =1;
end
ogfile = fopen(p);
file = strrep(p,'z','dz');
nfile1 = fopen(strrep(file,'.lcd','.lcd.1'),'r+');
if nfile1 == -1 %if this is the first time dspike is run it will set the flag to 1 and create a new dz.....lcd file
fflag =1; %if does not exist
nfile1 = fopen(strrep(file,'.lcd','.lcd.1'),'w+');
end
nfile2 = fopen(strrep(file,'.lcd','.lcd.2'),'r+');
if nfile2 == -1
nfile2 = fopen(strrep(file,'.lcd','.lcd.2'),'w+');
end
%=============first time
if fflag ==1
m = strrep(p,'z','');
m = strrep(m,'.lcd','.zcal');
mat = load(m,'.ASC');
%format long
%ptsfordncastatrix = in{2} - in{1}
%in2 = in{2}
%in1 = in{1}
%important these are clear every loop
clear ij
clear downcastmatrix
clear upcastmatrix
clear diffdn
clear diffup
width = length(mat(1,:));
for ij = in{1}:1:in{2};
downcastmatrix(ij-(in{1}-1),:) = mat(ij,1:width); %starts from 1 and goes to in2
end
for r = in{3}:in{4};
upcastmatrix(r-(in{3}-1),:) = mat(r,1:width);
end
end
%============
%============if done before
if fflag ~=1 %gets matrix name to read
m1 = strcat('desp1d',p);
mat1 = load(m1,'.ASC');
m2 = strcat('desp2d',p);
mat2 = load(m2,'.ASC');
% for i = in{1}:1:in{2};
% downcastmatrix(i-(in{1}-1),:) = mat1(i,:); %starts from 1 and goes to in2
downcastmatrix = mat1;
upcastmatrix = mat2;
sizemat1 = size(mat1);
% disp('line 110')
sizeoldup = size(upcastmatrix);
sizeolddn = size(downcastmatrix);
end
%=========== endif done before
%--------------------end splits
%------------------------------------------------------------------------
%------------------------------------------------------------------------
%--------------------computes difference in points for up and down-------
%------------------------------------------------------------------------
clear totalptsdn
clear totalptsup
%disp('line 126')
totalptsdn = in{2} - in{1};
totalptsup = in{4} - in{3};
%-------------if first time processed
%downcast matrix consists only wavelengths which will be processed with dspike
rd=0;
if fflag ==1
for cdn = 2:19 %shifted over by 1 and lost one point due to 1mer removal
%while cdn+1 ~=41
cdn ;
for rd = 1:length(downcastmatrix(:,1))%totalptsdn
if rd == 1
diffdn(rd,cdn) = abs(downcastmatrix(rd,cdn) - downcastmatrix((rd+1),cdn));
else
diffdn(rd,cdn)=abs(downcastmatrix(rd,cdn)-downcastmatrix((rd-1),cdn)); %commented
end
end
clear rd
for rd = 1:length(downcastmatrix(:,1))
if diffdn(rd,cdn) > trend;%all lines are the same after column 28 since they have been inserted from every 10th row and coppied therefore there is no diffdn(720,28+)
diffdn(rd,cdn);
downcastmatrix(rd,cdn);
movavgdn(rd,cdn) = moavg(flag,diffdn,rd,cdn,downcastmatrix,trend,totalptsdn);
movavgdn(rd,cdn);
end
end
end
for cu = 2:19
%while cu+1 ~=41
for ru = 1:length(upcastmatrix(:,1))%totalptsup%length(upcastmatrix(:,1))%totalptsup
if ru ==1
%diffup(ru,cu) = abs(upcastmatrix(ru,cu) - upcastmatrix((ru+1),(cu)));
diffup(ru,cu) = abs(upcastmatrix(ru,cu) - upcastmatrix((ru+1),(cu)));
else
diffup(ru,cu) = abs(upcastmatrix(ru,cu) - upcastmatrix((ru-1),(cu)));
end
end
clear ru
for ru = 1:length(upcastmatrix(:,1))
if diffup(ru,cu) > trend
movavgup(ru,cu) = moavg(flag,diffup,ru,cu,upcastmatrix,trend,totalptsup);
end
end
end
end
%-------------end if first time
%disp('line 188 after it calculates movavg and diffup')
%iffirstsizediffdn = size(diffdn);
%iffirstsizediffup = size(diffup);
%iffirstsizemovavgdn = size(movavgdn);
%lengthdowncastmatrix_forloop = length(downcastmatrix(:,1));
%=================if not done before
%disp('berfore here?')
if fflag ~=1
clear diffdn
clear movavgdn
clear diffup
clear movavgup
clear cdn
clear rd
for cdn = 42:59
for rd = 1:length(downcastmatrix(:,1));%totalptsdn
%format long
if rd ==1
diffdn(rd,cdn) = abs(downcastmatrix(rd,cdn) - downcastmatrix((rd+1),(cdn))); %diffdn or up will be empty until the correct column location in the data matrix
else
diffdn(rd,cdn)= abs(downcastmatrix(rd,cdn) - downcastmatrix((rd-1),cdn));
end
end
clear rd
for rd = 1:length(downcastmatrix(:,1));
if diffdn(rd,cdn) > trend ; %all lines are the same after column 28 since they have been inserted from every 10th row and coppied therefore there is no diffdn(720,28+)
movavgdn(rd,cdn) = moavg(flag,diffdn,rd,cdn,downcastmatrix,trend,totalptsdn);
end
end
end
clear ru
clear cu
for cu = 42:59
for ru = 1:length(upcastmatrix)%totalptsup
if ru ==1
diffup(ru,cu) = abs(upcastmatrix(ru,cu) - upcastmatrix((ru+1),(cu)));
else
diffup(ru,cu) = abs(upcastmatrix(ru,cu) - upcastmatrix((ru-1),(cu)));
end
end
clear ru
for ru = 1:length(upcastmatrix)
if diffup(ru,cu) > trend;
movavgup(ru,cu) = moavg(flag,diffup,ru,cu,upcastmatrix,trend,totalptsup);
end
end
end
sizediffdn = size(diffdn);
sizediffup = size(diffup);
diffdn59 = diffdn(1,59);
end
%==============end if done before
%------------------------------------------------------------------------
%------------end computes difference---------------------------
%------------------------------------------------------------------------
%------------------------------------------------------------------------
%------------writes down to derived parameters
ogline = 0; %trip into loop
if fflag ==1
while strncmp(ogline,'<derived_parameters>',20) ~=1
ogline = fgets(ogfile);
fwrite(nfile1,ogline);
fwrite(nfile2,ogline);
derivedloc = ftell(nfile1);
end
end
%------------end writes
%-------check to see if already despiked if so write 'd-d-'
check =0;
u= 0;
clear stuff
stuff = 0;
%=====================for getting oplist if files already dspiked
if fflag ~=1
% loca = ftell(nfile1)
%frewind(nfile1)
while strncmp(stuff,'<derived_parameters>',20) ~=1
% %stuff = fgets(ogfile);
stuff = fgets(nfile1);
if strncmp(stuff,'a650ref',7);
%
oplist(1,:) = (strcat('d-',stuff(:,1:4))');
num = 2;
for tr = num:18
stuff = fgets(nfile1);
% %oplist(2,:) = strcat('d-',stuff(:,1:4))
oplist(tr,:) = strcat('d-',stuff(:,1:4));
end
end
dploc = ftell(nfile1);
end
%=====================end for getting
%=====================rest of list commands for already processed file
while strncmp(check,'<data>',6) ~= 1
u = u+1;
check = fgets(nfile1);
if strncmp(check,'d-',2) ==1
auglist(u,:) = strcat('d-',check(:,1:6));
end
end
end
%====================end of rest for already processed
%-------end check
%--------------------creating new oplist if not already processed
clear u
stuff = 0; %trip into loop that makes list of derived
if fflag == 1;
%frewind(nfile1)
frewind(ogfile)
while strncmp(stuff,'<derived_parameters>',20) ~=1
stuff = fgets(ogfile);
if strncmp(stuff,'a650ref',7);
oplist(1,:) = (strcat('d-',stuff(:,1:4))');
num = 2;
for yt = num:18
stuff = fgets(ogfile);
%oplist(2,:) = strcat('d-',stuff(:,1:4))
oplist(yt,:) = strcat('d-',stuff(:,1:4));
end
end
end
end
%---------------------end of creating new file
%%-------------------------make output file-----------------------------
%downcast 1 upcast 2
%-------------------------make out matricies loop----------------------
%----------------------------------------------------------------------
%----------------------------------------------------------------------
%
% %allocate space for entire matrix information + new
%oplist is the list of wavelengths
if fflag ==1
wholedn = zeros(totalptsdn,(41+length(oplist)));
wholeup = zeros(totalptsup,(41+length(oplist)));
%disp('writes to wholedn')
%sizewholedn = size(wholedn)
op = length(oplist);
end
if fflag ~= 1
wholedn = zeros(totalptsdn,(41+length(oplist) + length(auglist)));
wholeup = zeros(totalptsup,(41+length(oplist) + length(auglist)));
aug = length(auglist);
op = length(oplist);
sizeofwholedn = size(wholedn);
end
%at this point it's needed to write the entire matrix to the new lcd up and
%down files
%to do this, the original data will be inserted in the whole-up/dn matrix
%in its proper location, then depending on weather diff > trend, either the
%movavg points will be inserted or the original points will be
%starts here
%if this is the second time dspike is used the original dspiked matrix will
%be augmented again with new d-d- columns
%~~~~~~~~~~~~~~~~~~~~~~~~~~
%for down segments
if fflag ==1 %first time
for c =1:length(downcastmatrix(1,:)); %for all 58 columns
for k= 1:length(downcastmatrix(:,1))%for all the rows in downcastmatrix
wholedn(k,c) = downcastmatrix(k,c); %puts origina stuff in for first 41 columns
end
end
clear r
clear c
%size(downcastmatrix)
%size(diffdn)
%diffdn(720,2)
count =0;
for c =2:19 ;%length(diffdn(1,:))
% count = 0;
count = count+1;
for r = 1:length(diffdn(:,1))%totalptsdn %will skip last line because difference is one line less; it needs two lines
if diffdn(r,c) > trend;
wholedn(r,41+count) = movavgdn(r,c);
% nowassn = wholedn(r,41+count)
%break
%return%----------------------------------------------
else
wholedn(r,41+count) = downcastmatrix(r,c);
end
end
end
clear c
%for up segments
for c =1:length(upcastmatrix(1,:)); %for all 58 columns
for k= 1:length(upcastmatrix(:,1)) %for all the rows in downcastmatrix
wholeup(k,c) = upcastmatrix(k,c); %puts origina stuff in for first 41 columns
end
end
clear r
count = 0;
for c =2:19; %length(diffdn(1,:))
count = count+1;
for r = 1:length(diffup(:,1));%totalptsup;
if diffup(r,c) > trend;
wholeup(r,41+count) = movavgup(r,c);
else
wholeup(r,41+count) = upcastmatrix(r,c);
end
end
end
clear count
nname = strrep(p,'za','dza');
justbeforeitwrites = size(wholedn);
dlmwrite(strcat('desp1',nname),wholedn,'delimiter',' ','precision',6);
dlmwrite(strcat('desp2',nname),wholeup,'delimiter',' ','precision',6);
end
%~~~~~~~~~~~~~~~~~~~~~~~~~
%adds data to wholedn if already done before
if fflag ~=1
for c =1:length(downcastmatrix(1,:)); %for all 58 columns
for k= 1:length(downcastmatrix(:,1))%totalptsdn ;%for all the rows in downcastmatrix
wholedn(k,c) = downcastmatrix(k,c); %puts origina stuff in for first 41 columns
end
end
clear r
clear c
for r = 1:length(diffdn(:,1));%totalptsdn;
count = 0;
for c =42:59 ;%length(diffdn(1,:))
count = count+1;
if diffdn(r,c) > trend;
wholedn(r,59+count) = movavgdn(r,c);
else
wholedn(r,59+count) = downcastmatrix(r,c);
end
end
clear count
end
clear c
%for up segments
for c =1:length(upcastmatrix(1,:)); %for all 58 columns
for k= 1:length(upcastmatrix(:,1))%totalptsup; %for all the rows in downcastmatrix
wholeup(k,c) = upcastmatrix(k,c); %puts origina stuff in for first 41 columns
end
end
clear r
for r = 1:length(diffup(:,1));%totalptsup;
count = 0;
for c =42:59; %length(diffdn(1,:))
count = count+1;
if diffup(r,c) > trend;
wholeup(r,59+count) = movavgup(r,c);
else
wholeup(r,59+count) = upcastmatrix(r,c);
end
end
clear count
end
nname = strrep(p,'za','dza')%=====================6666666666666666666666666666
%nname = strcat('d',p)
dlmwrite(strcat('desp1',nname),wholedn,'delimiter',' ','precision',6);
dlmwrite(strcat('desp2',nname),wholeup,'delimiter',' ','precision',6);
end
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%--------------------------end make out loop-------------------------------
%------------loop to print oplist or aug after derived parameters-----------------
if fflag == 1; %if creating new file
for j= 1:length(oplist);
fprintf(nfile1,'%s\n',oplist(j,:));
fprintf(nfile2,'%s\n',oplist(j,:));
end
fprintf(nfile1,'<data>\n');
fprintf(nfile2,'<data>\n');
end
%===================
if fflag ~=1 %if working with old file
fseek(nfile1,dploc,-1);
fseek(nfile2,dploc,-1);
for j = 1:length(oplist)
fprintf(nfile1,'%s\n',oplist(j,:));
fprintf(nfile2,'%s\n',oplist(j,:));
end
for k = 1:length(auglist)
fprintf(nfile1,'%s\n',auglist(k,:));
fprintf(nfile2,'%s\n',auglist(k,:));
end
fprintf(nfile1,'<data>\n');
fprintf(nfile2,'<data>\n');
ndatloc = ftell(nfile1);
end
%==================
%------end loop for oplist or aug-----------------------------------------------
%==============part that actually writes matrix to lcd file=============
mindn = fopen(strcat('desp1',nname));
minup = fopen(strcat('desp2',nname));
%working with new file
if fflag ==1
for o= 1:length(diffdn(:,1))%totalptsdn
ln = fgets(mindn);
fprintf(nfile1,ln);
end
for o= 1:length(diffup(:,1))%totalptsup
ln = fgets(minup);
fprintf(nfile2,ln);
end
end
%if working with already created file
if fflag ~=1
fseek(nfile1,ndatloc,-1)
fseek(nfile2,ndatloc,-1)
for o= 1:length(diffdn(:,1));%totalptsdn
ln = fgets(mindn);
fprintf(nfile1,ln);
end
for o= 1:length(diffup(:,1));%totalptsup
ln = fgets(minup);
fprintf(nfile2,ln);
end
end
fclose(mindn);
fclose(minup);
%==============end write matrix to lcd ===================================
%-----------notes filters that have already been used
frewind(ogfile);
%frewind(nfile1)
while feof(ogfile)~=1
loc1 = ftell(ogfile);
lnprev = fgets(ogfile);
if strncmp(lnprev,'<filters_used>',14) == 1 %tests for filters used location
fprintf(nfile1,'<filters_used>\n');
fprintf(nfile2,'<filters_used>\n');
% fseek(ogfile,loc1,-1)
fil = 0;% trip into loop
count =0;
while feof(ogfile) ~=1 %makes variable to hold each filter
count = count +1;
fil = fgets(ogfile);
fwrite(nfile1,fil);
fwrite(nfile2,fil);
eval(['filter',num2str(count) '= fil']);
end
end
end
%-----------end notes
%-----------add to filters the dspiked filters
for b = 1:length(oplist);
fprintf(nfile1,'%s', 'dspike ');
fprintf(nfile1,'%s',oplist(b,:));
fprintf(nfile1,'%d %d\n', trend,flag);
fprintf(nfile2,'%s', 'dspike ');
fprintf(nfile2,'%s',oplist(b,:));
fprintf(nfile2,'%d %d\n', trend,flag);
end
if fflag ~=1
%fseek(nfile1,loc1,-1);
fprintf(nfile1,'%s', 'dspike ');
fprintf(nfile1,'%s',oplist(b,:));
fprintf(nfile1,'%d %d\n', trend,flag);
fprintf(nfile2,'%s', 'dspike ');
fprintf(nfile2,'%s',oplist(b,:));
fprintf(nfile2,'%d %d\n', trend,flag);
for b = 1:length(auglist);
fprintf(nfile1,'%s', 'dspike ');
fprintf(nfile1,'%s',auglist(b,:));
fprintf(nfile1,'%d %d\n', trend,flag);
end
clear b
%fseek(nfile2,loc2,-1)
for b = 1:length(auglist);
fprintf(nfile2,'%s', 'dspike ');
fprintf(nfile2,'%s',auglist(b,:));
fprintf(nfile2,'%d %d\n', trend,flag);
end
end
%-------------update list_ac9--------------
mark = ftell(id);
if fflag ==1
frewind(id);
h=0;
r =1;
for po=1:h1 %-1
h = h+2;
r = h-1;
ac9 = fgets(id);
t = textscan(ac9,'%s');
newlst = char(t{1,:});
newlst = strrep(newlst,'z','dz');
newl(r,:) = strrep(newlst,'.lcd','.lcd.1');
newl(h,:) = strrep(newlst,'.lcd','.lcd.2');
end
%if strncmp(p,'z',1) ==1
dlmwrite('list_ac9dz',newl,'');
%end
end
fseek(id,mark,-1);
fclose(ogfile);
if endflag==1
break
end
end
fclose all;
