MITT v1.0

Original version as published in Computers and Geosciences (authors:
Bruce MacVicar, Scott Dilling, and Jay Lacey) along with test files

CalcArrayStats.m

@@ -0,0 +1,21 @@
+function xdat = CalcArrayStats(dat,typeanalysisi)
+% calculates a statistical parameter from a data array
+% called from AFiltArrayGUI
+
+% 
+switch typeanalysisi
+    case 'sum '
+        xdat = sum(dat);
+    case 'mean'
+        xdat = mean(dat);
+    case 'std '
+        xdat = std(dat);
+    case 'skew'
+        xdat = skewness(dat);
+    case 'kurt'
+        xdat = kurtosis(dat);
+    case 'box '
+        xdat = mean(dat);
+end
+
+end

CalcChannelMesh.m

@@ -0,0 +1,28 @@
+function GUIControl = CalcChannelMesh(GUIControl,CSVControl)
+% 3rd level program for organization of calculation of test volume
+% topography and water surface grid
+% Called from OrganizeInput
+% Calls InterpolateUniformChannel, InterpolateNonUniformChannel, or custom subprogram
+
+% check how channel information will be input
+% if it is a uniform channel
+if GUIControl.IsUniform == 1
+    % send to subprogram
+    [oneD,twoD] = InterpUniformChan(GUIControl.Slope,GUIControl.Width,GUIControl.Depth,GUIControl.Length,GUIControl.Sideslope,GUIControl.Widthgrid,GUIControl.Depthgrid,GUIControl.Lengthgrid);
+% if it is non-uniform
+else
+    % get
+    if GUIControl.Channel == 1 % a *.csv file has been specified
+        filename1 = [GUIControl.CalcChannelpathname,GUIControl.CalcChannelfile];
+        [oneD,twoD] = InterpNonUniformChan(filename1);
+    elseif GUIControl.Channel == 2 % a subprogram has been specified
+        eval(['[oneD,twoD] = ',GUIControl.CalcChannelfile(1:end-2),'(GUIControl,CSVControl);']);
+
+    end
+end
+
+% save one and two dimensional meshes to C for safekeeping
+GUIControl.oneD = oneD;
+GUIControl.twoD = twoD;
+
+end

CalcGoodCells.m

@@ -0,0 +1,182 @@
+function Config = CalcGoodCells(Config,Data,GUIControl)
+% assesses quality of sampled cells based on various flexible criteria
+% called from ClassifyArrayGUI, ClassifyArrayAuto
+% calls various ClassifyCor, Classifyxcorr, Classifyw1w2xcorr, ClassifyNoise, ClassifySpike, ClassifyPoly
+% NoiseHurtherLemmin, ConvMulti2Struct 
+% Note: criteria can be controlled from launch window for automatic
+% analysis or from interactive array plot
+
+faQCname = fieldnames(Config.faQC);
+xvar = {'Vel','Despiked','Filtered'};
+yvar = {'zZ'};
+% if it is empty, no filter has yet been run, or if you want to reset
+if isempty(faQCname)||GUIControl.resetFilter
+    % set faQC from C structure array
+    if isfield(GUIControl,'faQC')
+        faQC = GUIControl.faQC;
+    else
+        faQC = DefaultfaQC;
+    end
+
+    ncomptot = length(Config.comp);
+    goodCellsi = ones(Config.nCells,ncomptot);
+    goodCellsii = ones(Config.nCells,ncomptot);
+
+    % create matrix to store
+    nQtot = ncomptot+1+faQC.Range+faQC.w1w2xcorr+ncomptot*(faQC.Correlation+faQC.xcorr+faQC.Spikes+faQC.InertialSlope+faQC.PolyFilt);
+    Qdat = zeros(Config.nCells,nQtot);
+    Qcnames = cell(1,nQtot);
+    Qi = ncomptot+1;
+
+    % filter by position in the water column
+    outofwater = Config.zZ<=0 | Config.zZ>=1;
+    goodCellsi(outofwater,:) = 0;
+    Qcnames{Qi} = 'z/Z';
+    Qdat(:,Qi) = Config.zZ;
+    Qi = Qi+1;
+
+    % filter by lateral position
+%     yY = Config.ypos/Config.Y;
+%     outofwater = yY<0 | yY>1;
+%     goodCellsi(outofwater,:) = 0;
+
+    %% filtering subprograms
+    % filter by range
+    if faQC.Range
+        if ~isempty(faQC.nigood)||faQC.nigood>1;
+            if faQC.nigood > Config.nCells
+                goodCellsi(:,:) = 0;
+            else 
+                goodCellsi(1:faQC.nigood-1,:)=0;
+            end
+        end
+        if ~isempty(faQC.negood)||faQC.negood<Config.nCells
+            goodCellsi(faQC.negood+1:end,:)=0;
+        end
+    end
+    % filter by signal correlation
+    if faQC.Correlation && isfield(Data,'Cor')
+        for ncomp = 1:ncomptot
+            compi = char(Config.comp(ncomp));
+            [goodCellsii(:,ncomp),QCor] = ClassifyCor(compi,Data.Cor,faQC.Corthreshold);
+            Qcnames{Qi} = ['CORBeam',num2str(ncomp)];
+            Qdat(:,Qi) = QCor(ncomp);
+            Qi = Qi+1;
+
+        end
+        goodCellsi = goodCellsi.*goodCellsii;
+
+    end
+    % filter by correlation between adjacent cells
+    if faQC.xcorr
+        % for each component
+        for ncomp = 1:ncomptot
+            compi = char(Config.comp(ncomp));
+            %here
+            xdata = Data.(GUIControl.X.var).(compi);
+            [goodCellsii(:,ncomp),Qxcorr] = Classifyxcorr(xdata,faQC.xcorrthreshold);
+            Qcnames{Qi} = ['X corr ',compi,' (%)'];
+            Qdat(:,Qi) = Qxcorr;
+            Qi = Qi+1;
+        end
+        goodCellsi = goodCellsi.*goodCellsii;
+
+    end
+    % filter by Hurther Lemmin w1 w2 
+    if faQC.w1w2xcorr
+        if isfield(Data.Vel,'w2')
+            xdata = Data.(GUIControl.X.var);
+            [goodCellsiii,QNRW] = ClassifyNoiseRatio(xdata,faQC.w1w2xcorrthreshold,Config.transformationMatrix);
+            Qdat(:,Qi) = QNRW;
+            for ncomp = 1:ncomptot
+                goodCellsi(:,ncomp) = goodCellsi(:,ncomp).*goodCellsiii;
+            end
+        else
+            Qdat(:,Qi) = NaN;
+
+        end
+        Qcnames{Qi} = 'Noise Ratio (%)';
+        Qi = Qi+1;
+
+    end
+    % filter by spike numbers
+    if faQC.Spikes 
+        if isfield(Data,'SpikeY')
+            % for each component
+            for ncomp = 1:ncomptot
+                compi = char(Config.comp(ncomp));
+                SpikeYi = Data.SpikeY.(compi);
+                [goodCellsii(:,ncomp),QSpike] = ClassifySpike(SpikeYi,faQC.SpikeThreshold);
+                Qcnames{Qi} = ['Spikes ',compi,'(%)'];
+                Qdat(:,Qi) = QSpike;
+                Qi = Qi+1;
+            end
+            goodCellsi = goodCellsi.*goodCellsii;
+        else
+            Qdat(:,Qi) = NaN;
+
+        end
+
+    end
+    % filter by -5/3 slope in inertial range
+    if faQC.InertialSlope
+        % don't use filtered because it changes the slope
+        if isfield(Data,'Despike')
+            xdata = Data.Despike;
+        else
+            xdata = Data.Vel;
+        end
+        for ncomp = 1:ncomptot
+            compi = char(Config.comp(ncomp));
+            %xdata = Data.(C.X.var).(compi);
+            xdatai = xdata.(compi);
+            [goodCellsii(:,ncomp),QSis] = ClassifyNoiseFloor(xdatai,faQC.InertialSlopeThreshold,Config.Hz,Config.zpos,Config.samplingVolume);
+            Qcnames{Qi} = ['S Inertial ',compi];
+            Qdat(:,Qi) = QSis;
+            Qi = Qi+1;
+        end
+
+        goodCellsi = goodCellsi.*goodCellsii;
+    end
+    % filter by a 3rd order polynomial to the mean, std, and skewness of profile
+    if faQC.PolyFilt
+        ydata = Config.(GUIControl.Y.var);
+
+        % for each component
+        for ncomp = 1:ncomptot
+            compi = char(Config.comp(ncomp));
+            xdata = Data.(GUIControl.X.var).(compi);
+            goodCellsii(:,ncomp) = ClassifyPoly(ydata,xdata,goodCellsi(:,ncomp),faQC.zscore);
+            Qcnames{Qi} = ['Poly fit ',compi];
+            Qdat(:,Qi) = goodCellsii(:,ncomp);
+            Qi = Qi+1;
+        end
+        goodCellsi = goodCellsi.*goodCellsii;
+    end
+    for ncomp=1:ncomptot
+        compi = char(Config.comp(ncomp));
+        eval(['Config.goodCells.',compi,' = goodCellsi(:,ncomp)'';']);% ' added Sept 10 so that goodCells is a 1xnCells matrix, with each value in a column (why?  so annoying!)
+    end
+
+    %% output to table
+    % add in goodCells
+    for ncomp=1:ncomptot
+        compi = char(Config.comp(ncomp));
+        Qcnames{ncomp} = ['goodCells ',compi];
+        Qdat(:,ncomp) = goodCellsi(:,ncomp);
+    end
+
+    %% save data
+    Config.faQC = faQC;
+    Config.Qcnames = Qcnames;
+    Config.Qdat = Qdat;
+end
+
+end
+
+%%%%%
+
+
+
+
+%%%%%

CalcIllinoisFlumePos.m

@@ -0,0 +1,57 @@
+function [oneD,twoD] = CalcIllinoisFlumePos(GUIControl,CSVControl)
+% calculate 1 and 2D surfaces for plotting and calculation
+
+filename = [GUIControl.CSVControlpathname,CSVControl(1).ffname,'.csv'];
+
+% get flume positions from data file
+[xchannel,ychannel,zchannel,beamh,wsurf]=textread(filename,'%f %f %f %f %f',...
+    'delimiter',',','commentstyle','matlab');
+
+%% 1D way
+% isolate xpositions and interpolate to a finer grid
+xu = unique(xchannel);
+nxutot = length(xu);
+bedi = zeros(nxutot,1);
+wsurfi = zeros(nxutot,1);
+beamhi = zeros(nxutot,1);
+
+% find average bed position
+for nx = 1:nxutot
+    xi = xchannel==xu(nx);
+    bedi(nx) = mean(zchannel(xi));
+    wsurfi(nx) = mean(wsurf(xi));
+    beamhi(nx) = mean(beamh(xi));
+end
+
+% interpolate on regular grid and
+% put data into C structure
+xmem = xu(1):diff(xu)/200:xu(end);
+oneD.xchannel = xmem;
+oneD.bedElevation = interp1(xu,bedi,xmem);
+oneD.waterElevation = interp1(xu,wsurfi,xmem);
+oneD.beamh = interp1(xu,beamhi,xmem);
+oneD.Y = CSVControl(1).Y*ones(1,length(xmem));
+
+% isolate xpositions and interpolate to a finer grid
+ymem = 0:CSVControl(1).Y/20:CSVControl(1).Y;
+oneD.ymem = ymem;
+
+%% 2D way 
+
+% create 2D grid for interpolation
+[xgrid,ygrid]=ndgrid(xmem,ymem); % grid for 2d xy plane
+
+% create 2D interpolants for bed and water surface
+Fbed = scatteredInterpolant(xchannel,ychannel,zchannel,'linear','linear');
+Fwater = scatteredInterpolant(xchannel,ychannel,wsurf,'linear','linear');
+
+
+twoD.xchannel = xgrid;
+twoD.ychannel = ygrid;
+twoD.bedElevation = Fbed(xgrid,ygrid);
+twoD.waterElevation = Fwater(xgrid,ygrid);
+% can get the bed and water surface position of any x and y point
+twoD.Fbed = Fbed; 
+twoD.Fwater = Fwater;
+
+end

CalcXYZMoras.m

@@ -0,0 +1,47 @@
+function Config = CalcXYZMoras(GUIControl,Config)
+% rotate data from local coordinate system to global system
+
+% user parameters
+dmultiplier = 0.01;
+
+% load bridgefile for rotation
+bridgeDate = Config.bridgeDate;
+if length(bridgeDate)<6
+    bridgeDate = ['0',bridgeDate]; % to compensate for zeros that disappear when you go in and out of excel with a text file
+end
+bridgefile = [GUIControl.CSVControlpathname,'bridge',bridgeDate,'.mat'];
+load(bridgefile); 
+
+% choose the benchmark to which the local data points will be translated and
+% around which they will be rotated
+pont = eval(['bridge',bridgeDate]);
+origg = find(pont(:,1) == 12);
+% choose a second benchmark for calculation of rotation angle
+rotg = find(pont(:,1) == 18);
+% identify origin and rotation points
+xgo = pont(origg,2);
+ygo = pont(origg,3);
+xgr = pont(rotg,2);
+ygr = pont(rotg,3);
+% transform rotation point into polar coordinates
+[thetarot,rhorot] = cart2pol(xgr-xgo,ygr-ygo);
+
+% calculate translation distances to global origin
+% translate to origin (using point x=2,y=0)
+xlt = Config.xlocal-0;
+ylt = Config.ylocal-0;
+%convert to polar
+[thetal,rhol] = cart2pol(xlt,ylt);
+% remove local angle (equal to 0), add global
+thetalr = thetal-0+thetarot;
+% transform back into cartesian
+[xltr,yltr] = pol2cart(thetalr,rhol);
+% add constants to obtain new coordinates
+Config.xpos=xltr+xgo;
+Config.ypos=yltr+ygo;
+
+% calculate zpos data
+Config.zpos = [Config.z1 Config.z2 Config.z3 Config.z4]*dmultiplier;
+Config.waterDepth = Config.waterDepth*dmultiplier;
+Config.bedElevation = Config.bridgeElevation-Config.bridgeDepth*dmultiplier; % bed surface 
+end