text:
function r8lib_test121 ( )

%*****************************************************************************80
%
%% R8LIB_TEST121 tests R8VEC_FRAC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8LIB_TEST121\n' );
  fprintf ( 1, '  R8VEC_FRAC: K-th smallest DVEC entry;\n' );

  seed = 123456789;

  [ a, seed ] = r8vec_uniform_01 ( n, seed );

  r8vec_print ( n, a, '  Array to search:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Fractile  Value\n' );
  fprintf ( 1, '\n' );

  for k = 1 : floor ( n / 2 ) : n

    afrac = r8vec_frac ( n, a, k );

    fprintf ( 1, '  %8d  %14f\n', k, afrac );

  end

  return
end


text:
function test_suite = test_comparisonMessage
%test_comparisonMessage Unit test for comparisonMessage.

%   Steven L. Eddins
%   Copyright 2009 The MathWorks, Inc.

initTestSuite;

function test_happyCase
s = xunit.utils.comparisonMessage('user message', 'assertion message', ...
    [1 2 3], 'hello');
c = xunit.utils.stringToCellArray(s);

expected_output = { 'user message'
    'assertion message'
    ''
    'First input:'
    '     1     2     3'
    ''
    'Second input:'
    'hello'};

assertEqual(c, expected_output);

author: ricedsp


text:
function [xUpdate,SUpdate,innov,Szz,W]=sqrtEKFUpdate(xPred,SPred,z,SR,h,HJacob,innovTrans,stateTrans)
%SQRTEKFUPDATE Perform the measurement update step in a square-root version
%              of the first-order Extended Kalman Filter (EKF).
%
%INPUTS: xPred The xDim X 1 predicted target state.
%        SPred The xDim X xDim lower-triangular square root predicted state
%              covariance matrix.
%            z The zDim X 1 measurement vector.
%           SR The zDim X zDim lower-triangular square root of the
%              measurement covariance matrix in the native coordinate
%              system of the measurement.
%            h A function handle for the measurement function that takes
%              the state as its argument.
%       HJacob A function handle for the measurement Jacobian matrix that
%              takes the target state as a parameter. If not supplied or an
%              empty matrix is passed, then HJacob will be found using
%              numerical differentiation via the numDiff function with
%              default parameters.
%   innovTrans An optional function handle that computes and optionally
%              transforms the value of the difference between the
%              observation and any predicted points. This is called as
%              innovTrans(a,b) and the default if omitted or an empty
%              matrix is passed is @(a,b)bsxfun(@minus,a,b). This must be
%              able to handle sets of values. For a zDimX1 measurement,
%              either of the inputs could be zDimXN in size while one of
%              the inputs could be zDimX1 in size.  This only needs to be
%              supplied when a measurement difference must be restricted
%              to a certain range. For example, the innovation between two
%              angles will be 2*pi if one angle is zero and the other
%              2*pi, even though they are the same direction. In such an
%              instance, a function handle to the
%              wrapRange(bsxfun(@minus,a,b),-pi,pi) function with the
%              appropriate parameters should be passed for innovTrans.
%   stateTrans An optional function that takes a state estimate and
%              transforms it. This is useful if one wishes the elements of
%              the state to be bound to a certain domain. For example, if
%              an element of the state is an angle, one should generally
%              want to bind it to the region +/-pi.
%
%OUTPUTS: xUpdate The xDim X 1 updated state vector.
%         SUpdate The updated xDim X xDim lower-triangular square root
%                 state covariance matrix.
%      innov, Szz The zDimX1 innovation and the zDimXzDim square root
%                 innovation covariance matrix are returned in case one
%                 wishes to analyze the consistency of the estimator or use
%                 those values in gating or likelihood evaluation.
%               W The xDimXzDim gain used in the update. This can be
%                 useful when gating and using the function
%                 calcMissedGateCov.
%
%The first-order EKF is summarized in Figure 10.3.3-1 in Chapter 10.3.3 of
%[1]. The Joseph-form covariance update given in Chapter 5 of the same book
%is used for improved numerical stability. The mathematics behind the
%specific square root implementation used here are those used in the
%standard Kalman filter, which are described in [2] with a flow chart given
%in Appendix G.
%
%The optional parameter innovTrans is not described in the above
%reference, but allows for possible modifications to the filter as
%described in [3].
%
%The parameters have been added to allow the filter to be used with
%angular quantities. For example, if the measurement consisted of
%range and angle, z=[r;theta], then
%innovTrans=@(a,b)[bsxfun(@minus,a(1,:),b(1,:));
%                  wrapRange(bsxfun(@minus,a(2,:),b(2,:)),-pi,pi)];
%should be used to approximately deal with the circular nature of the
%measurements.
%
%REFERENCES:
%[1] Y. Bar-Shalom, X. R. Li, and T. Kirubarajan, Estimation with
%    Applications to Tracking and Navigation. New York: John Wiley and
%    Sons, Inc, 2001.
%[2] David F. Crouse , "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems
%    Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%[3] David F. Crouse, "Cubature/ unscented/ sigma point Kalman filtering
%    with angular measurement models," in Proceedings of the 18th
%    International Conference on Information Fusion, Washington, D.C.,
%    6-9 Jul. 2015.
%
%March 2015 David Karnick, Naval Research Laboratory, Washington D.C.
%Updated with additional inputs, August 2015 David F. Crouse,
%                               Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

zDim=size(z,1);

if(nargin<8||isempty(stateTrans))
    stateTrans=@(x)x;
end

if(nargin<7||isempty(innovTrans))
    innovTrans=@(a,b)bsxfun(@minus,a,b);
end

if(nargin<6||isempty(HJacob))
    HJacob=@(x)numDiff(x,h,zDim);
end

zPred=h(xPred);

H=HJacob(xPred);
Pxz=SPred*SPred'*H';

Szz=tria([H*SPred,SR]);
W=(Pxz/Szz')/Szz;

innov=innovTrans(z,zPred);
xUpdate=stateTrans(xPred+W*innov);

temp=W*H;
SUpdate=tria([(eye(size(temp))-temp)*SPred,W*SR]);

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% The COBRAToolbox: testSelMatrix.m
%
% Purpose:
%     - testSelMatrix tests the functionality of testSelMatrix()
%
% Authors:
%     - Lemmer El Assal February 2017
%

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testSelMatrix'));
cd(fileDir);

% load reference data
load('refData_selMatrix.mat');

selMat1 = selMatrix(selVec1);
selMat2 = selMatrix(selVec2);
assert(isequal(selMat1, ref_selMat1))
assert(isequal(selMat2, ref_selMat2))

%return to original directory
cd(currentDir);

author: opencobra

repo: cobratoolbox


text:
function J=normVecJacob(x)
%%NORVECJACOB Consider a unit vector created as u=x/norm(x). This function
%         returns the partial derivatives of u with respect to the elements
%         of x. This partial derivative
%
%INPUTS: x A real NX1 vector.
%
%OUTPUTS: J The NXN partial derivative matrix The columns are the element
%           of x with respect to which the partial derivative is taken and
%           the rows are the elements of u.
%
%EXAMPLE:
%Here, we verify that the output of this function agrees with the numerical
%differentiation result.
% x=[4;-18;12];
% u=@(x)x/norm(x);
% JNumDiff=numDiff(x,u,3);
% J=normVecJacob(x);
% RelErr=max(abs((J(:)-JNumDiff(:))./J(:)))
%The relative error will be about 3.8488e-10, indicating a good level of
%agreement between the numerical differentiation result and the analytic
%result.
%
%November 2020 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

xDim=length(x);
normX=norm(x);

J=-x*x'/normX^3+eye(xDim,xDim)/normX;

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
pairs = [0.1 2.2527126517342059598697; ...
			0.6 .39823385806923489961685; ...
			0.7 .26086724653166651438573; ...
			1.0 0; ...
			2.0 0; ...
			3.4 1.0923280598027415674947; ...
			4.0 1.791759469228055000812477; ...
			8.0 8.525161361065414300165531; ...
			64.0 201.00931639928152667928; ...
			256.0 1161.71210111840065079];
err = [];
for i = 1:rows(pairs)
	err(i) = abs(gammaln(pairs(i,1)) - pairs(i,2))/pairs(i,1);
end
if max(err) > 1e-10
	err
	error('maximum err > 1e-10')
end
err = [];
for i = 1:rows(pairs)
	err(i) = abs(gammaln(pairs(i,1),1) - pairs(i,2))/pairs(i,1);
end
if max(err) > 1e-10
	err
	error('maximum err > 1e-10')
end
err = abs(gammaln(1.1,2) - 0.920726359734123);
if err > 1e-10
	error('gammaln(1.1,2) != 0.920726359734123');
end
if gammaln(0) ~= Inf
  error('gammaln(0) should be Inf');
end
if ~isnan(gammaln(-1))
  %error('gammaln(-1) should be NaN');
end
if gammaln(Inf) ~= Inf
  error('gammaln(Inf) should be Inf');
end
% should be NaN?
%gammaln(-Inf)
if ~isnan(gammaln(NaN))
  error('gammaln(NaN) should be NaN');
end

author: andrewssobral

repo: lrslibrary


text:
run = true;
last_mrk = '';

state_acquire= bbci_acquire_lsl('init');
cnt = zeros(64, 100);
mrk = zeros(1, 100);

step = 1;

while state_acquire.running
    [cnt_new, cntTime, mrkTime, mrkDesc, state_acquire]= bbci_acquire_lsl(state_acquire);

    if not(isempty(mrkDesc))
        fprintf('RECEIVED %s\n',  mrkDesc{1});
        fprintf('TIME %f \n',  mrkTime);
    end

    %   fprintf('Data %s at %f \n', cnt_new, cntTime);
%     fprintf('step %f \n', step)
%     last_mrk = mrkDesc;
%     last_mrkTime = mrkTime;
%
%     cnt(:,step) = cnt_new;
    step = step + 1;
end
display('Streams broke off')

author: bbci


text:
function launcher(alignmentType)
% clc; clear; close all
addpath(genpath('thirdparty'));

%%
if ~exist('alignmentType','var')
    alignmentType = 1; %0 for nowarp, 1 for optical flow, 2 for homography, 3 for similarity
end
path2dataset = '../dataset/qualitative_datasets';

%%
datasets = dir(path2dataset);
dirFlags = [datasets.isdir];
datasets(~dirFlags) = [];
datasets = {datasets(3:end).name};

% permute and divide into training and testing sets
num_datasets = numel(datasets);
seed = 12;
rng(seed);
ds_idx = randperm(num_datasets);

alignment = '';
if alignmentType == 0
    alignment = '_nowarp';
elseif alignmentType == 1
    alignment = '_OF';
elseif alignmentType == 2
    alignment = '_homography';
end

path2output = @(stage,d) sprintf('../data/testing_real_all_nostab%s/%s/image_%s',alignment,stage,d);

fn_in = @(dsi,type,fri) [path2dataset '/' datasets{ds_idx(dsi)} '/' type '/' sprintf('%05d.jpg',fri)];
fn_out = @(stage,fri,frid) sprintf('%s/%05d.jpg',path2output(stage,frid),fri);

ds_range = 1:num_datasets;

for l = -2:2
    for i = 1:num_datasets
        checkDir(path2output(datasets{i},num2str(l)));
    end
end

%%
for ii = 1:length(ds_range)
    fr_cnt = 0;
    i = ds_range(ii);
    % get the frame range
    datasets{ds_idx(i)}
    files = dir([path2dataset '/' datasets{ds_idx(i)} '/input/*.jpg']);
    if isempty(files)
        files = dir([path2dataset '/' datasets{ds_idx(i)} '/input/*.png']);
    end
    if ~isempty(files)
        [~,ststr,~] = fileparts(files(1).name);
        [~,enstr,~] = fileparts(files(end).name);
        start_frame = str2num(ststr);
        end_frame = str2num(enstr);
        frame_range = start_frame:min(start_frame+99,end_frame);
        num_frame = numel(frame_range);

        fr_idx = floor(linspace(frame_range(1),frame_range(end),num_frame));
        for j = 1:num_frame
            fr_cnt = fr_cnt+1;
            % save image_1 to image_5
            v0 = im2double(imread(fn_in(i,'input',fr_idx(j)+0)));
            v0g = single(rgb2gray(v0));
            [h,w,~] = size(v0);

            for l = -2:2
                if l ~= 0
                    vi = im2double(imread(fn_in(i,'input',max(min(fr_idx(j)+l,frame_range(end)),frame_range(1)))));
                    vig = single(rgb2gray(vi));
                    if alignmentType == 0
                        v_i0 = vi;
                    elseif alignmentType == 1
                        flo_i0 = genFlow(v0g, vig);
                        [v_i0, ~] = warpToRef(v0, vi, flo_i0);
                    elseif alignmentType == 2
                        v_i0 = homographyAlignment(v0,vi,0);
                    elseif alignmentType == 3
                        v_i0 = similarityAlignment(v0,vi,0);
                    end
                else
                    v_i0 = v0;
                end
                imwrite(v_i0, fn_out(datasets{ds_idx(i)},fr_cnt,num2str(l)));
            end
        end
    end
end

author: shuochsu

repo: DeepVideoDeblurring


text:
filename='Bridge_Quadrilateral_Bilinear_Structured';
ptype = 'MACRO';
method = 'SIMP_P3'; % !! Instead of proportional to material density !!
materialType = 'ISOTROPIC';
initial_case = 'holes';
cost = {'compliance'};
weights = [1];
constraint = {'volumeConstraint'};
optimizer = 'HAMILTON-JACOBI'; incrementFactor = 1;
filterType = 'P1';
constraint_case = 'INEQUALITY';

HJiter0 = 1;
e2 = 30;
N_holes = [5 6];
R_holes = 0.7;
phase_holes = [0 pi/2];

nsteps = 1;
Vfrac_final = 0.5;
Perimeter_target=3.5;
optimality_final =1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 0.3;
TOL.nu_minus = 0.3;

plotting = 1;
printing = 0;
monitoring = 1;
monitoring_interval = 1;

maxiter = 50;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: compact version of matrix-based regularization
%
%  S(y) = alpha/2 * norm(B*(y-yRef))^2,
%
% where
%  alpha regularization parameter, weights regularization versus
%        distance in the joint objective function, alpha = 1 here
%  yRef  is a reference configuration, e.g. yRef = x or a
%        pre-registration result
%  B     a discrete partial differential operator either in explicit
%        matrix form or as a structure containing the necessary
%        parameters to compute B*y
% see also regularizer E8_regularization_MF
%==============================================================================

clear, close all, help(mfilename);

% (c) Jan Modersitzki 2009/03/25, see FAIR.2 and FAIRcopyright.m.
% note: '%' is used so hide comments in the book
% Example for usage of regularization
% (c) Jan Modersitzki 2009/04/02, see FAIR.2 and FAIRcopyright.m.
% illustrates the usage of L2-norm based regularization
%
%

% initialize the regularization and create a starting  point
regularizer('reset','regularizer','mbElastic','alpha',1,'mu',1,'lambda',0);
y0 = @(omega,m) randn(size(getStaggeredGrid(omega,m)));

% 2D example, initialize physical domain and number of discretization points
omega = [0,1,0,1]; m = [16,12];   %

% test derivative of 2D implementation
fctn = @(yc) regularizer(yc,omega,m);  checkDerivative(fctn,y0(omega,m));

% 3D example, initialize physical domain and number of discretization points
omega = [0,1,0,1,0,1]; m  = [16,12,8];

% test derivative of 3D implementation
fctn = @(yc) regularizer(yc,omega,m);
checkDerivative(fctn,y0(omega,m));
%==============================================================================


text:
function [] = anim8_DIC_image_3Dmeasure_faces_2n(ImSet,DIC2DpairResults,DIC3DpairResults,faceMeasureString,varargin)
%% function for plotting 2D-DIC results imported from Ncorr in step 2
% called inside plotNcorrPairResults
% plotting the images chosen for stereo DIC (2 views) with the
% correlated points results plotted on top, colored as their correlation
% coefficient.
% on the left side the images from the reference camera (reference image and current images), and on the right side the
% images from the deformed camera
% requirements: GIBBON toolbox
%
%
% INPUT:
% * IMset - a 2nX1 cell array containing 2n grayscale images. The first n
% images are from camera A (the "reference" camera), and the last n images
% are from camera B (the "deformed" camera). The first image in the set is
% considered as the reference image, on which the reference grid of points
% is defined, and all the correlated points and consequent displacements
% and strains, are relative to this image.
% * DIC_2Dpair_results - containig the correlated points, correlation
% coefficients, faces..
% * optional: CorCoeffCutOff - - maximal correlation coefficient to plot
% points
% * optional: CorCoeffDispMax - maximal correlation coefficient in colorbar

%%
DIC2DpairResultsL=DIC2DpairResults{1};
DIC2DpairResultsR=DIC2DpairResults{2};
DIC3DpairResultsL=DIC3DpairResults{1};
DIC3DpairResultsR=DIC3DpairResults{2};

nImages=DIC2DpairResultsL.nImages;
nCam=DIC2DpairResultsL.nCamDef;

PointsL=DIC2DpairResultsL.Points(nImages+1:end);
PointsR=DIC2DpairResultsR.Points(1:nImages);
FL=DIC2DpairResultsL.Faces;
FR=DIC2DpairResultsR.Faces;
FaceCorrL=DIC3DpairResultsL.FaceCorrComb;
FaceCorrR=DIC3DpairResultsR.FaceCorrComb;

switch nargin
    case 4 % in case no results were entered
        optStruct=struct;
    case 5
        optStruct=varargin{1};
    otherwise
        error('wrong number of input arguments');
end

%% cut out point with large correlation coefficient
if ~isfield(optStruct,'CorCoeffCutOff')
    CorCoeffCutOff=max([cell2mat(FaceCorrL); cell2mat(FaceCorrR)]);
else
    CorCoeffCutOff=optStruct.CorCoeffCutOff;
end

for ii=1:nImages
    FaceCorrL{ii}(FaceCorrL{ii}>CorCoeffCutOff)=NaN;
    FaceCorrR{ii}(FaceCorrR{ii}>CorCoeffCutOff)=NaN;
end

%%
switch faceMeasureString
    case {'J','Lamda1','Lamda2'}
        FCL=DIC3DpairResultsL.Deform.(faceMeasureString);
        FCR=DIC3DpairResultsR.Deform.(faceMeasureString);
        cMap=coldwarm;
        if ~isfield(optStruct,'FClimits')
            FCmax=0;
            for ii=1:nImages
                FCmax=max([max(abs(FCL{ii}(~isnan(FaceCorrL{ii}))-1)) max(abs(FCR{ii}(~isnan(FaceCorrR{ii}))-1)) FCmax]);
            end
            FClimits=[1-FCmax 1+FCmax];
        else
            FClimits=optStruct.FClimits;
        end
    case {'Emgn','emgn'}
        FCL=DIC3DpairResultsL.Deform.(faceMeasureString);
        FCR=DIC3DpairResultsR.Deform.(faceMeasureString);        cMap='parula';
        if ~isfield(optStruct,'FClimits')
            FCmax=0;
            for ii=1:nImages
                FCmax=max([max(FCL{ii}(~isnan(FaceCorrL{ii}))) max(FCR{ii}(~isnan(FaceCorrR{ii}))) FCmax]);
            end
            FClimits=[0 FCmax];
        else
            FClimits=optStruct.FClimits;
        end
    case {'Epc1','Epc2','epc1','epc2'}
        FCL=DIC3DpairResultsL.Deform.(faceMeasureString);
        FCR=DIC3DpairResultsR.Deform.(faceMeasureString);        cMap=coldwarm;
        if ~isfield(optStruct,'FClimits')
            FCmax=0;
            for ii=1:nImages
                FCmax=max([max(abs(FCL{ii}(~isnan(FaceCorrL{ii})))) max(abs(FCR{ii}(~isnan(FaceCorrR{ii})))) FCmax]);
            end
            FClimits=[-FCmax FCmax];
        else
            FClimits=optStruct.FClimits;
        end
    otherwise
        error('unexpected face measure string. plots not created');
end

for ii=1:nImages
    FCL{ii}(isnan(FaceCorrL{ii}))=NaN;
    FCR{ii}(isnan(FaceCorrR{ii}))=NaN;
end
%%
hf=figure; hold all;
hf.Units='normalized'; hf.OuterPosition=[.05 .05 .9 .9]; hf.Units='pixels';

ii=1;
hp1=imagesc(repmat(ImSet{ii},1,1,3)); hold on; axis ij
hp2=gpatch(FL,PointsL{1},FCL{ii},'none',0.5);
hp3=gpatch(FR,PointsR{1},FCR{ii},'none',0.5);
pbaspect([size(ImSet{ii},2) size(ImSet{ii},1) 1])
hs1=title(['Cam ' num2str(nCam) ' frame ' num2str(1)]);
colormap(cMap);
hc1=colorbar;
caxis(FClimits)
title(hc1, faceMeasureString);
hc1.FontSize=16;
axis off
drawnow

%Create the time vector
animStruct.Time=linspace(0,1,nImages);

for ii=1:nImages
    PnowL=PointsL{ii};
    PnowR=PointsR{ii};

    cNowL=FCL{ii};
    cNowR=FCR{ii};

    TitleNow=['Cam ' num2str(nCam) ' frame ' num2str(ii)];

    %Set entries in animation structure
    animStruct.Handles{ii}=[hp1,hp2,hp2,hp3,hp3,hs1]; %Handles of objects to animate
    animStruct.Props{ii}={'CData','CData','Vertices','CData','Vertices','String'}; %Properties of objects to animate
    animStruct.Set{ii}={repmat(ImSet{ii},1,1,3),cNowL,PnowL,cNowR,PnowR,TitleNow}; %Property values for to set in order to animate

end

anim8(hf,animStruct);

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function     [position,janelas,messages]= twavef_detector(heasig,samp,time,position,w,intervalo,messages)
% script which identifies T wave, as well as its onset and offset
% implementation of changes to improve FN and reduce mean error utilizando a escala 5
%
%
%Input Parameters:
%   heasig: struct vector with header information
%   samp: samples included in the current excerpt (borders excluded)
%   time:  QRS times in inedexes refering the interval included in the current excerpt (borders excluded)
%   position: struct vector with the detected points
%   w: matrix with WT scales 1 to 5
%   intervalo: numeration of the beats processed in this segment
%
%Output Parameters:
%   janelas: T wave seach windows
%   actualized parameters: position
%
% Rute Almeida
% based on twave.m by Juan Pablo Mart�nez Cort�s
% Last update: Rute Almeida  07FEB2012
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
%
% Designed for MATLAB Version R12; tested with MATLAB Version R13
%
%%%%%% Constants and Thresholds  !!!!!!!!!!!!!!!!!!!!!!!!!
if ~isfield(messages.setup.wavedet,'Kton')
    messages.setup.wavedet.Kton = 4;      % 2 4!
end
if ~isfield(messages.setup.wavedet,'Ktoff')
    messages.setup.wavedet.Ktoff = 2.5;% 3.5 3!
end
if ~isfield(messages.setup.wavedet,'umbraldetT')
    messages.setup.wavedet.umbraldetT = 0.25;  % We use umbraldet*sqrt(mean(w(time(i):time(i+1),4).^2))
end
if ~isfield(messages.setup.wavedet,'umbralsig')
    messages.setup.wavedet.umbralsig  =  1/8;    % To decide if there are really 2, 1 or no peak
    %threshold to decide if there is or not a significative
    % T wave.3!
end
if ~isfield(messages.setup.wavedet,'rrant_mintol')
    messages.setup.wavedet.rrant_mintol  =  0.5;    % minimum time in sec admited for current RR to be used in T wave delineation or in Exponentially averaged RR
end
if ~isfield(messages.setup.wavedet,'rrant_maxtol')
    messages.setup.wavedet.rrant_maxtol  =  1.5;    % max time in sec admited for current RR to be used in Exponentially averaged RR
end
if ~isfield(messages.setup.wavedet,'rrant_average')
    messages.setup.wavedet.rrant_average  =  [0.8 0.2];    % Exponentially averaged RR weights
end
if ~isfield(messages.setup.wavedet,'inivent_tol')
    messages.setup.wavedet.inivent_tol  =  0.1;
end
if ~isfield(messages.setup.wavedet,'inivent_tol_S')
    messages.setup.wavedet.inivent_tol_S=0.05; % sec
end
if ~isfield(messages.setup.wavedet,'finvent_tol')
    messages.setup.wavedet.finvent_tol=0.240;% sec % it would be nioce to make it depend on round(rrmed/2)
end
if ~isfield(messages.setup.wavedet,'finvent_max')
    messages.setup.wavedet.finvent_max=0.6;% sec
end
if ~isfield(messages.setup.wavedet,'scale')
    messages.setup.wavedet.scale=4;
end
if ~isfield(messages.setup.wavedet,'scale2')
    messages.setup.wavedet.scale2=5;
end
if ~isfield(messages.setup.wavedet,'scalezerocros')
    messages.setup.wavedet.scalezerocros=3;
end
if ~isfield(messages.setup.wavedet,'scalezerocros2') %for twavetask5
    messages.setup.wavedet.scalezerocros2=4;
end

if ~isfield(messages.setup.wavedet,'min_vent')
    messages.setup.wavedet.min_vent=0.1;%sec
end
if ~isfield(messages.setup.wavedet,'Tmax_Tmin_time_min')
    messages.setup.wavedet.Tmax_Tmin_time_min=0.2;%0.15;%sec %10 ENE2012: 150 ms es poco
end
if ~isfield(messages.setup.wavedet,'Tmax_Tmin_bifasic')
    messages.setup.wavedet.Tmax_Tmin_bifasic=2.5;
end
if ~isfield(messages.setup.wavedet,'T_bound_tol')
    messages.setup.wavedet.T_bound_tol=0.12;%sec
end
Kton=messages.setup.wavedet.Kton;
Ktoff=messages.setup.wavedet.Ktoff;
umbraldetT = messages.setup.wavedet.umbraldetT;
umbralsig = messages.setup.wavedet.umbralsig;
rrant_mintol=messages.setup.wavedet.rrant_mintol;
rrant_maxtol=messages.setup.wavedet.rrant_maxtol;
rrant_average=messages.setup.wavedet.rrant_average;
inivent_tol= messages.setup.wavedet.inivent_tol;
inivent_tol_S=messages.setup.wavedet.inivent_tol_S;
finvent_tol=messages.setup.wavedet.finvent_tol;
finvent_max= messages.setup.wavedet.finvent_max;
scale=messages.setup.wavedet.scale;
scale2=messages.setup.wavedet.scale2;
scalezerocros=messages.setup.wavedet.scalezerocros;
scalezerocros2=messages.setup.wavedet.scalezerocros2; %#ok<NASGU>
min_vent=messages.setup.wavedet.min_vent;%sec
Tmax_Tmin_time_min= messages.setup.wavedet.Tmax_Tmin_time_min;%sec
Tmax_Tmin_bifasic=messages.setup.wavedet.Tmax_Tmin_bifasic;% extra criteria
T_bound_tol=messages.setup.wavedet.T_bound_tol;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isempty(heasig)

end
% janelas=[];
% endwinT=[];

% Initialization of auxiliary variables
if length(time)>1
    rrant = time(2)-time(1);
else
    aux=find(position.qrs~=0);
    if length(aux)>1
        rrant= diff(position.qrs(aux(end-1):aux(end))); %18MAR09
    else
        rrant=messages.setup.wavedet.freq;
    end
end

%18MAR09 % for the case in which the first RR is very wrong AGO2011
if rrant<rrant_mintol*messages.setup.wavedet.freq || rrant>rrant_maxtol*messages.setup.wavedet.freq;
    messages.warnings=[messages.warnings {['RR lower than' rrant_mintol*messages.setup.wavedet.freq ' ms or higher than' rrant_maxtol*messages.setup.wavedet.freq ' found in twavef']}];
    rrant=rrant_mintol*messages.setup.wavedet.freq;%Jul2011
end

T = []; Tprima=[]; picon=[]; picoff=[]; Ton=[]; Toff=[]; tipoT=[];
% minapos = []; minppos=[]; maxapos=[]; maxppos=[]; mina=[]; minp=[]; maxa=[];
% maxp=[];
picon_keep=[];%multileadchange
picoff_keep=[];%multileadchange
lead_keep=[];%multileadchange

%figure
%%%plot(messages.lixo)
% hold on
janelas=NaN*ones(length(time),3);
endwinT=NaN*ones(length(time),1);
for i = 1:length(time)   % For each beat

    if (i>1),
        if (rrant_maxtol*rrant>(time(i)-time(i-1)))&&(time(i)-time(i-1)>rrant_mintol*rrant),
            rrmed = rrant_average(1)*rrant + rrant_average(2)*(time(i)-time(i-1));% Exponentially averaged RR new value
        else
            rrmed = rrant;  % Exponentially averaged RR old value
        end
    else                  % Only for the first in each segment of the ECG
        rrmed = rrant;
    end
    rrant = rrmed;        % For next segment
    inivent = round(inivent_tol*messages.setup.wavedet.freq);   % Begining of window

    %Rute multilead 03.Dec.04

    %if ~isempty(pos.S(i)),              % If there is an S wave
    %         inivent = max(inivent, pos.S(i)-pos.qrs(i)+round(0.05*messages.setup.wavedet.freq));
    %     end

    if ~isempty(position.S(i+intervalo(1)-1));
        inivent = max(inivent, position.S(i+intervalo(1)-1)-samp(1)+1-time(i)+round(inivent_tol_S*messages.setup.wavedet.freq));
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% changed 13/06/02 Rute
    %if rrmed >= messages.setup.wavedet.freq,  %
    %if rrmed <= messages.setup.wavedet.freq,  %   21AGO09 % for rrmed < 1 seg finiven=finvent_max seg
    if rrmed <= messages.setup.wavedet.freq && i~=length(time) && (rrant_maxtol*rrant>(time(i+1)-time(i)))&&(time(i+1)-time(i)>rrant_mintol*rrant) %02FEB2011
        finvent = round(finvent_max*messages.setup.wavedet.freq);     % End of window
    else
        finvent = round(rrmed*finvent_max);
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %% changed 13/06/02 Rute
    if i~=length(time),                 % For last beat in the segment
        %finvent = min(finvent,time(i+1)-time(i)-round(finvent_tol*messages.setup.wavedet.freq));
        %finvent = min(finvent,time(i+1)-time(i)-min(round((time(i+1)-time(i))*1/3),round(finvent_tol*messages.setup.wavedet.freq)));
        finvent = min(finvent,time(i+1)-time(i)-min(round((time(i+1)-time(i))*0.4),round(finvent_tol*messages.setup.wavedet.freq)));
    elseif i~=1%Rute 18MAR09
        finvent = min(finvent,time(i)-time(i-1)-round(finvent_tol*messages.setup.wavedet.freq));
    elseif length(aux)>1
        finvent = min(finvent,diff(position.qrs(aux(end-1):aux(end)))-round(finvent_tol*messages.setup.wavedet.freq));%Rute 18MAR09
    end

    % We work at scale 4 (in general).
    if isempty(position.QRSoff(i+intervalo(1)-1)),   % Should never happen, but...
        begwin = min(inivent + time(i),length(w)); %Rute 4Jul2011
    else
        begwin = min(max(inivent + time(i), position.QRSoff(i+intervalo(1)-1)-samp(1)+1+1),length(w));%Rute 4Jul2011
        % if QRSoff is bad anoated it can miss T wave...
    end

    endwin = max(1,min(finvent + time(i),length(w))); %% Rute 20/05/02
    janelas(i,:)=[i begwin endwin];
    endwinT(i)=endwin;
    %%%%plot([begwin endwin ], messages.lixo([begwin endwin]),'*g')

    % Positive maxima and negative minima in the window
    maxpos = begwin + modmax(w(begwin+1:endwin,scale),2,0,+1);
    minpos = begwin + modmax(w(begwin+1:endwin,scale),2,0,-1);
    [maxim ind] = max(w(maxpos,scale));   % The biggest of the positive
    maxpos = maxpos(ind);
    [minim ind] = min(w(minpos,scale));   % The biggest of the negative
    minpos = minpos(ind);

    if isempty(maxpos),               % If no local positive maximum
        % The maximum will be the first
        % or the last sample
        if (w(begwin,scale)>=w(endwin,scale)) && w(begwin,scale)>0,
            maxpos = begwin; maxim = w(maxpos,scale);
        elseif (w(endwin,scale)>=w(begwin,scale)) && w(endwin,scale)>0,
            maxpos = endwin; maxim = w(maxpos,scale);
        end
    end
    if isempty(minpos),               % if no local negative minimum
        % the minimum will be the first
        % or the last sample
        if (w(begwin,scale)<=w(endwin,scale)) && w(begwin,scale)<0,
            minpos = begwin; minim = w(minpos,scale);
        elseif (w(endwin,scale)<=w(begwin,scale)) && w(endwin,scale)<0,
            minpos = endwin; minim = w(minpos,scale);
        end
    end

    absmax = abs(maxim);
    absmin = abs(minim);

    if i<length(time),
        veficaz = sqrt(mean(w(time(i):time(i+1),scale).^2));
    else                        % if last beat of the segment
        veficaz = sqrt(mean(w(time(i):end,scale).^2));
    end

    hay_onda = ((absmax>umbraldetT*veficaz)|(absmin>umbraldetT*veficaz));
    if isempty(hay_onda)
        hay_onda = 0;
    end
    % Rute 18/06/02

    if endwin-begwin<(min_vent*messages.setup.wavedet.freq) || isnan(hay_onda)% se a janela tem amplitude menor do que min_vent seg enato nao existe onda T
        hay_onda =0;
    end

    % Is there a wave?
    if hay_onda,
        if absmax >= absmin,    % the greatest modulus maximum is the maximum
            % Now we search the two minima nearest to maxpos, one before and one after
            minapos = max(modmax(w(begwin+1:maxpos-1,scale),2,0,-1));
            minapos = begwin + minapos; % Position of the negative minimum before the maximum
            if isempty(minapos) && (maxpos ~= begwin) && (w(begwin,scale)<0),
                minapos = begwin;% If no local minimum before the maximum, take the first sample
            end
            minppos = min(modmax(w(maxpos+1:endwin,scale),2,0,-1));
            minppos = maxpos + minppos; % Position of the positive maximum after the minimum
            if isempty(minppos) && (maxpos ~= endwin) && (w(endwin,scale)<0),
                minppos = endwin;        % If no local minimum after the maximum, take the last sample
            end

            mina = abs(w(minapos,scale));     % Amplitude of minimum before maximum
            minp = abs(w(minppos,scale));     % Amplitude of minimum after maximum

            if (mina < umbralsig*absmax), % If mina is not big enough
                mina =[];                  % forget it

            elseif (maxpos-minapos>Tmax_Tmin_time_min*messages.setup.wavedet.freq), %or if ther are more than Tmax_Tmin_time_min sec to maxpos
                mina = [];
            end
            if (minp < umbralsig*absmax), % If minp is not big enough
                minp =[];                  % forget it
            elseif (minppos-maxpos>Tmax_Tmin_time_min*messages.setup.wavedet.freq), % and also if there are more than Tmax_Tmin_time_min sec to maxpos
                minp =[];
            end

            if ~isnan(mina)&~isnan(minp), %#ok<AND2>     %%% NUEVO JP
                if (mina >= minp)&&(minp < umbralsig*absmax*Tmax_Tmin_bifasic),
                    minp = [];
                elseif (minp> mina)&&(mina < umbralsig*absmax*Tmax_Tmin_bifasic),
                    mina = [];
                end
            end

            % Test which modulus maxima are significative and find zero crossings
            if isempty(mina),
                if isempty(minp),
                    tipoT = 2;      % only upwards T wave
                    if maxpos - minapos > 2,         % if not !!!!!!!!!!!
                        ind = zerocros(flipud(w(minapos:maxpos,scalezerocros)));   %Scale scalezerocros !!!
                        T = maxpos - ind +1;                           % Zero crossing = T wave position
                        picoff = maxpos;                               %wavelet  peak to detect offset
                    elseif isempty(minapos);  % If there were no minimum, there is no zero crossing
                        T = picant (w(begwin:maxpos,scale),maxpos);       % Take the minimum at scale scale
                        if ~isempty(T) %%%%%%%%%%%% 14/06/02 Rute
                            picoff = maxpos;
                        else
                            picoff = []; % if did not exist a peak in scale scale there is no T
                        end %%%%%%%%%%%% 14/06/02 Rute
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave' }];
                    end
                else      % minp exists (is significative) but mina not
                    tipoT = 0;  		%normal T wave
                    if minppos -maxpos >2,       % if not!!!!!!???
                        ind = zerocros(w(maxpos:minppos,scalezerocros));  %% 07/06/02 Rute
                        if isempty(ind) %%%%%%%%%%%%%%Rute 03/09/02
                            ind = zerocros(w(maxpos:minppos,scale));
                        end %%%%%%%%%%%%%%Rute 03/09/02
                        T = maxpos + ind -1;
                        picon = maxpos;		% For determining onset and offset
                        picoff = minppos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave'}];
                    end
                end
            else
                if isempty(minp),   %  mina exists (is significative) but minp not
                    tipoT = 1;  	%inverted T wave
                    if maxpos -minapos >2,  % if not !!!!!!!!!
                        ind = zerocros(w(minapos:maxpos,scalezerocros));  % wavelet zero crossing is T wave peak  %% 07/06/02 Rute
                        if isempty(ind)%%%%%%%%%%%%%%Rute 03/09/02
                            ind = zerocros(w(minapos:maxpos,scale));  % wavelet zero crossing is T wave peak  %% 03/09/02 Rute
                        end %%%%%%%%%%%%%%Rute 03/09/02
                        T = minapos + ind -1;
                        picon = minapos;
                        picoff = maxpos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave'}];
                    end
                else    % both mina and minp are significative.  Biphasic wave.
                    tipoT = 5;	% biphasic neg-pos T wave
                    if maxpos - minapos > 2,    %!!!!!!!!!!!!
                        ind = zerocros(flipud(w(minapos:maxpos,scalezerocros))); %% 07/06/02 Rute
                        if isempty(ind) %%%%%%%%%%%%%%Rute 03/09/02
                            ind = zerocros(flipud(w(minapos:maxpos,scale))); %% 03/09/02 Rute

                        end%%%%%%%%%%%%%%Rute 03/09/02
                        T = maxpos - ind +1;
                        picon = minapos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave'}];
                    end
                    if minppos - maxpos > 2,    %!!!!!!!!!!!!
                        ind = zerocros(flipud(w(maxpos:minppos,scalezerocros))); %% 07/06/02 Rute
                        if isempty(ind)%%%%%%%%%%%%%%Rute 03/09/02
                            %ind = zerocros(flipud(w(minapos:maxpos,scale))); %% 03/09/02 Rute %18.11.05 DUVIDA
                            %DUVIDA 18Nov05
                            ind = zerocros(flipud(w(maxpos:minppos,scale))); %% 03/09/02 Rute
                        end %%%%%%%%%%%%%%Rute 03/09/02
                        %Tprima = maxpos + ind -1;
                        %DUVIDA 18Nov05
                        Tprima = minppos- ind +1;
                        picoff = minppos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave'}];
                    end
                end
            end
        else        % If the greatest modulus maximum is the minimum
            % Search two maxima, one before and one after the minimum
            maxapos = max(modmax(w(begwin+1:minpos-1,scale),2,0,1));
            maxapos = begwin + maxapos;
            if isempty(maxapos) && (minpos ~= begwin) && (w(begwin,scale)>0),
                maxapos = begwin;
            end
            maxppos = min(modmax(w(minpos+1:endwin,scale),2,0,1));
            maxppos = minpos + maxppos;
            if isempty(maxppos) && (minpos ~= endwin) && (w(endwin,scale)>0),
                maxppos = endwin;
            end
            maxa = abs(w(maxapos,scale));
            maxp = abs(w(maxppos,scale)) ;                                   % See if they are significative
            if (maxa < umbralsig*absmin)
                maxa =[];
            elseif (minpos-maxapos>Tmax_Tmin_time_min*messages.setup.wavedet.freq),
                maxa = [];
            end
            if (maxp < umbralsig*absmin),
                maxp =[];
            elseif (maxppos-minpos>Tmax_Tmin_time_min*messages.setup.wavedet.freq),
                maxp = [];
            end
            if ~isnan(maxa)&~isnan(maxp), %#ok<AND2>     %%% NUEVO JP
                if (maxa >= maxp)&&(maxp < umbralsig*absmin*Tmax_Tmin_bifasic),
                    maxp = [];
                elseif (maxp> maxa)&&(maxa < umbralsig*absmin*Tmax_Tmin_bifasic),
                    maxa = [];
                end
            end
            % Test which modulus maxima are significative and find zero crossings
            if isempty(maxa),
                if isempty(maxp),
                    tipoT = 3;      % only downwards T wave
                    if minpos - maxapos > 2, %!!!!!!!!
                        ind = zerocros(flipud(w(maxapos:minpos,scalezerocros)));   %Scale scalezerocros
                        if isempty(ind) %%%%%%%%%%%%%%Rute 03/09/02
                            ind = zerocros(flipud(w(maxapos:minpos,scale)));
                        end %%%%%%%%%%%%%%Rute 03/09/02
                        T = minpos - ind +1;
                        picoff = minpos;
                    elseif isempty(maxapos);  % If there were no maximum, there is no zero crossing.
                        T = picant (w(begwin:minpos,scale),minpos);  % minimo en escala scale.
                        if ~isempty(T) %%%%%%%%%%%% 14/06/02 Rute
                            picoff = minpos;
                        else
                            picoff = []; % if did not exist a peak in scale scale there is no T
                        end %%%%%%%%%%%% 14/06/02 Rute
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave'}];
                    end
                else      % maxp is signficative, but not maxa
                    tipoT = 1;  %inverted T wave
                    if maxppos -minpos >2,  % !!!!!!
                        ind = zerocros(w(minpos:maxppos,scalezerocros));  %% 07/06/02 Rute
                        if isempty(ind) %%%%%%%%%%%%%%Rute 03/09/02
                            ind = zerocros(w(minpos:maxppos,scale));  %% 03/09/02 Rute
                        end %%%%%%%%%%%%%%Rute 03/09/02
                        T = minpos + ind -1;
                        picon = minpos;		% For calculating onset and offset
                        picoff = maxppos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave'}];
                    end
                end
            else
                if isempty(maxp),   % maxa is significative, but not maxp
                    tipoT = 0;       %normal T wave
                    if minpos -maxapos >2,
                        ind = zerocros(w(maxapos:minpos,scalezerocros)); %% 07/06/02 Rute
                        if isempty(ind) %%%%%%%%%%%%%%%%%%%% 03/09/02 Rute
                            ind = zerocros(w(maxapos:minpos,scale)); %% 03/09/02 Rute
                        end %%%%%%%%%%%%%%%%%%%% 03/09/02 Rute
                        T = maxapos + ind -1;
                        picon = maxapos;
                        picoff = minpos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave'}];
                    end
                else    % both maxa and maxp are significative.  Biphasic wave.
                    tipoT = 4;	% biphasic pos-neg T wave
                    if minpos - maxapos > 2,  %!!!!!!!!!!!
                        ind = zerocros(flipud(w(maxapos:minpos,scalezerocros))); %% 07/06/02 Rute
                        if isempty(ind) %%%%%%%%%%%%%%%%%%%% 03/09/02 Rute
                            ind = zerocros(flipud(w(maxapos:minpos,scale))); %% 03/09/02 Rute
                        end %%%%%%%%%%%%%%%%%%%% 03/09/02 Rute
                        T = minpos - ind +1;
                        picon = maxapos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave'}];
                    end
                    if maxppos - minpos > 2,    %!!!!!!!!!!!!
                        ind = zerocros(flipud(w(minpos:maxppos,scalezerocros))); %% 07/06/02 Rute
                        if isempty(ind) %%%%%%%%%%%%%%%%%%%% 03/09/02 Rute
                            %ind = zerocros(flipud(w(minapos:maxpos,scale))); %% 03/09/02 Rute %18.11.05 DUVIDA
                            %DUVIDA 18Nov05
                            ind = zerocros(flipud(w(minpos:maxppos,scale))); %% 03/09/02 Rute
                        end %%%%%%%%%%%%%%%%%%%% 03/09/02 Rute
                        %Tprima = minpos + ind -1; %18.11.05 DUVIDA
                        %DUVIDA 18Nov05
                        Tprima = maxppos- ind +1;
                        picoff = maxppos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave'}];
                    end
                end
            end
        end

        % T wave onset and offset detection
        %if isempty(T), picon=[]; picoff=[]; end

        picon_keep=[picon_keep picon]; %#ok<AGROW> %multileadchange
        picoff_keep=[picoff_keep picoff];  %#ok<AGROW> %multileadchange
        lead_keep=[lead_keep scale];  %#ok<AGROW> %multileadchange

        if ~isempty(picon),
            Ton = searchon (picon, w(max(begwin,picon-round(T_bound_tol*messages.setup.wavedet.freq)):picon,scale), Kton);
        end
        if ~isempty(picoff),
            Toff=searchoff(picoff, w(picoff:min([size(w,1) picoff+round(T_bound_tol*messages.setup.wavedet.freq) ]) ,scale) , Ktoff);
            if (Toff > endwin),
                Toff = endwin;
            end
        end
    else
        %% using scale 5 07/06/02 Rute %%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% Neste momento usa-se a mesma janela e os mesmos limiares que para a escala scale

        maxpos = begwin + modmax(w(begwin+1:endwin,scale2),2,0,+1);
        minpos = begwin + modmax(w(begwin+1:endwin,scale2),2,0,-1);
        [maxim ind] = max(w(maxpos,scale2));   % The biggest of the positive
        maxpos = maxpos(ind);
        [minim ind] = min(w(minpos,scale2));   % The biggest of the negative
        minpos = minpos(ind);
        if isempty(maxpos),               % If no local positive maximum
            % The maximum will be the first
            % or the last sample
            if (w(begwin,scale2)>=w(endwin,scale2)) && w(begwin,scale2)>0,
                maxpos = begwin; maxim = w(maxpos,scale2);
            elseif (w(endwin,scale2)>=w(begwin,scale2)) && w(endwin,scale2)>0,
                maxpos = endwin; maxim = w(maxpos,scale2);
            end
        end
        if isempty(minpos),               % if no local negative minimum
            % the minimum will be the first
            % or the last sample
            if (w(begwin,scale2)<=w(endwin,scale2)) && w(begwin,scale2)<0,
                minpos = begwin; minim = w(minpos,scale2);
            elseif (w(endwin,scale2)<=w(begwin,scale2)) && w(endwin,scale2)<0,
                minpos = endwin; minim = w(minpos,scale2);
            end
        end

        absmax = abs(maxim);
        absmin = abs(minim);
        if i<length(time),
            veficaz = sqrt(mean(w(time(i):time(i+1),scale2).^2));
        else                        % if last beat of the segment
            veficaz = sqrt(mean(w(time(i):end,scale2).^2));
        end
        hay_ondascale2 = ((absmax>umbraldetT*veficaz)|(absmin>umbraldetT*veficaz));

        if endwin-begwin<(min_vent*messages.setup.wavedet.freq) % se a janela tem amplitude menor do que min_vent seg enato nao existe onda T
            hay_ondascale2 =0;
        end
        if hay_ondascale2,
            t5; % procurar na escala scale2 10/07/02 Rute
            picon_keep=[picon_keep picon];%#ok<AGROW> %multileadchange
            picoff_keep=[picoff_keep picoff];%#ok<AGROW>%multileadchange
            lead_keep=[lead_keep scale2];%#ok<AGROW>%multileadchange
        end
    end
    % Is there a wave
    %% utilizar a escala scale2 07/06/02 Rute %%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Filling the structure with positions
    if isempty(Ton), Ton = NaN; end;
    if isempty(Toff), Toff = NaN;  end;
    if isempty(T), T=NaN; end;
    if isnan(T),
        picon_keep=[picon_keep NaN]; %#ok<AGROW>%multileadchange
        picoff_keep=[picoff_keep NaN]; %#ok<AGROW>%multileadchange
        lead_keep=[lead_keep NaN]; %#ok<AGROW>%multileadchange
    end
    if isempty(Tprima), Tprima=NaN; end;
    if isempty(tipoT), tipoT=NaN; end;
    pos.Tscale(i)=scale;
    pos.Ton(i) = Ton;
    pos.Toff(i) = Toff;
    pos.T(i) = T;
    pos.Tprima(i) = Tprima;
    pos.Ttipo(i) = tipoT;
    %%%plot([Toff(~isnan(Toff))], messages.lixo([Toff(~isnan(Toff))]),'*r')
    T = []; Tprima=[]; picon=[]; picoff=[]; Ton=[]; Toff=[]; tipoT=[];
%     minapos = []; minppos=[]; maxapos=[]; maxppos=[]; mina=[]; minp=[]; maxa=[];
%     maxp=[];

end
position.Tscale(intervalo(1):intervalo(2)) = pos.Tscale;
position.Ton(intervalo(1):intervalo(2)) = pos.Ton+samp(1)-1;
position.Toff(intervalo(1):intervalo(2))= pos.Toff+samp(1)-1;
position.T(intervalo(1):intervalo(2)) = pos.T+samp(1)-1;
position.Tprima(intervalo(1):intervalo(2)) = pos.Tprima+samp(1)-1;
position.Ttipo(intervalo(1):intervalo(2)) = pos.Ttipo;

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function mbasis = basis_matrix_hermite ( )

%*****************************************************************************80
%
%% BASIS_MATRIX_HERMITE sets up the Hermite spline basis matrix.
%
%  Discussion:
%
%    This basis matrix assumes that the data points are stored as
%    ( P1, P2, P1', P2' ), with P1 and P1' being the data value and
%    the derivative dP/dT at T = 0, while P2 and P2' apply at T = 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Foley, van Dam, Feiner, Hughes,
%    Computer Graphics: Principles and Practice,
%    page 484.
%
%  Parameters:
%
%    Output, real MBASIS(4,4), the basis matrix.
%
  mbasis(1,1) =  2.0;
  mbasis(1,2) = -2.0;
  mbasis(1,3) =  1.0;
  mbasis(1,4) =  1.0;

  mbasis(2,1) = -3.0;
  mbasis(2,2) =  3.0;
  mbasis(2,3) = -2.0;
  mbasis(2,4) = -1.0;

  mbasis(3,1) =  0.0;
  mbasis(3,2) =  0.0;
  mbasis(3,3) =  1.0;
  mbasis(3,4) =  0.0;

  mbasis(4,1) =  1.0;
  mbasis(4,2) =  0.0;
  mbasis(4,3) =  0.0;
  mbasis(4,4) =  0.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function events = LoadEvents(filename)

%LoadEvents - Read events from file.
%
%  USAGE
%
%    events = LoadEvents(filename)
%
%    filename            event file name

% Copyright (C) 2004-2011 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

events.time = [];
events.description = [];
if ~exist(filename),
	error(['File ''' filename ''' not found.']);
end

file = fopen(filename,'r');
if file == -1,
	error(['Cannot read ' filename ' (insufficient access rights?).']);
end

while ~feof(file),
	time = fscanf(file,'%f',1);
	if isempty(time),
		if feof(file), break; end
		error(['Failed to read events from ' filename ' (possibly an empty file).']);
	end
	events.time(end+1,1) = time;
	line = fgetl(file);
	start = regexp(line,'[^\s]','once');
	events.description{end+1,1} = sscanf(line(start:end),'%c');
end
fclose(file);

% Convert to seconds
if ~isempty(events.time), events.time(:,1) = events.time(:,1) / 1000; end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function y = previewsong(frequencies, durations, gap, Fs)
% PREVIEWSONG plays the output of song() on the host PC.
%
%   PREVIEWSONG converts frequency domain data in the
%   form of a frequency vector and a duration vector
%   into a single time domain vector and plays it over
%   the PC speaker using the sound command.
%
%   frequencies - Vector of frequencies.
%   durations   - Vector of times.
%   gap         - Pause between notes. Default 0.01.
%   Fs          - Sample rate. Default 8192.
%
%   Example:
%      [sng.tune, sng.rhythm, sng.message, sng.stats] = song('random');
%      previewsong(sng.tune, sng.rhythm);

   if nargin < 4, Fs          = 8192;     end
   if nargin < 3, gap         = 0.01;     end
   if nargin < 2, durations   = 1;        end
   if nargin < 1, frequencies = 261.6256; end

   if length(durations) < length(frequencies)
      for i=length(durations)+1:length(frequencies)
         durations(i) = 0;
      end
   end

   if gap
      paddedFrequencies = [];
      paddedDurations = [];
      for i=1:length(frequencies)
         paddedFrequencies = [paddedFrequencies frequencies(i) 0]; %#ok<AGROW>
         paddedDurations = [paddedDurations durations(i) gap]; %#ok<AGROW>
      end
   else
      paddedFrequencies = frequencies;
      paddedDurations = durations;
   end

   deltaT = 1.0/Fs;

   startTime = 0;
   y = [];
   for i=1:length(paddedFrequencies)
%     timeIntervals = 0 : deltaT : paddedDurations(i);
      timeIntervals = startTime : deltaT : startTime+paddedDurations(i);
      startTime = timeIntervals(1)+deltaT;
      yCurrent = sin(2*pi*paddedFrequencies(i)*timeIntervals);
      y = [y yCurrent]; %#ok<AGROW>
   end

   sound(y, Fs);

return   % end of previewsong()

author: Sable

repo: mcbench-benchmarks


text:
function [ dist ] = euclidean_dist( X_gallery, X_probe, model_para )
    dist=pdist2(X_probe,X_gallery,'euclidean');
end

author: wuancong


text:
% Calculate adaptive smoothed rate map
%
% Calculates an adaptive smoothed rate map as described in "Skaggs et al 1996 -
% Theta Phase Precession in Hippocampal Neuronal Population and the Compression of Temporal Sequences"
%
%  USAGE
%   [map, posPdf] = analyses.ratemapAdaptiveSmoothing(pos, spkPos, <options>)
%   pos         Position samples (either Nx3 or Nx2 matrix). X and Y coordinates are used.
%   spkPos      Spike position samples (either Nx2 or Nx2 matrix). X and Y coordinates are used.
%   <options>   optional list of property-value pairs (see table below)
%
%   =========================================================================
%    Properties     Values
%   -------------------------------------------------------------------------
%    'binWidth'     width (and height) of bins in firing rate map. Value units
%                   are [cm]. binWidth has priority over parameter nBins. See
%                   description of nBins.
%
%    'nBins'        number of horizontal and vertical bins (default = [50 50]).
%                   If single value is provided, it is used for both horizontal
%                   and vertical number of bins (i.e. nBins = 50 => [50 50]).
%                   Either binWidth or nBins should be provided. If non is provided,
%                   then the default value of binWidth is used. If both are provided,
%                   then the value of binWidth is used.
%
%    'minTime'      minimum time spent in each bin (in s, default = 0).
%    'alphaValue'   scaling parameter. Default value is 10000. See original paper
%                   for the details.
%    'shape'        shape of the arena. Possible values are: 1 for square box,
%                   2 for cylinder. Only square box is currently supported!
%    =========================================================================
%
%  OUTPUT
%
%   map.x       x bins
%   map.y       y bins
%   map.z       adaptively smoothed firing rate map
%   posPdf      position probability density function

function [map, posPdf] = mapAdaptiveSmoothing(pos, spkPos, varargin)
    % Check number of parameters
    if nargin < 2 || mod(length(varargin), 2) ~= 0,
        error('BNT:numArgs', 'Incorrect number of parameters (type ''help analyses.mapAdaptiveSmoothing'' for details).');
    end

    % Check parameter sizes
    if size(pos, 2) < 2
        error('Parameter ''pos'' should have at least 2 columns (type ''help analyses.mapAdaptiveSmoothing'' for details).');
    end
    if size(spkPos, 2) < 2
        error('Parameter ''spkPos'' should have 1 or 2 columns (type ''help analyses.mapAdaptiveSmoothing'' for details).');
    end

    if isempty(pos)
        return;
    end

    % Default values
    alphaValue = 10000;
    binWidth = 5;
    minTime = 0;
    limits = [];
    shape = 1;
    sampleTime = helpers.sampleTimeFromData(pos);
    nBins = [50 50];

    haveBinWidth = false;
    haveNBins = false;

    % Parse parameter list
    i = 1;
    while i < length(varargin)
        if ~ischar(varargin{i}),
            error(['Parameter ' num2str(i+2) ' is not a property (type ''help analyses.mapAdaptiveSmoothing'' for details).']);
        end

        switch(lower(varargin{i})),
            case 'binwidth',
                binWidth = varargin{i+1};
                if ~helpers.isdscalar(binWidth, '>0')
                    error('Incorrect value for property ''nBins'' (type ''help analyses.mapAdaptiveSmoothing'' for details).');
                end
                haveBinWidth = true;
                i = i + 2;

            case 'nbins'
                nBins = varargin{i+1};
                if ~helpers.isivector(nBins, '>0') || length(nBins) > 2,
                    error('Incorrect value for property ''nBins'' (type ''help analyses.mapAdaptiveSmoothing'' for details).');
                end
                haveNBins = true;
                i = i + 2;

            case 'mintime',
                minTime = varargin{i+1};
                if ~helpers.isdscalar(minTime, '>=0'),
                    error('Incorrect value for property ''minTime'' (type ''help analyses.mapAdaptiveSmoothing'' for details).');
                end
                i = i + 2;

            case 'alphavalue'
                alphaValue = varargin{i+1};
                if ~helpers.isiscalar(alphaValue, '>0'),
                    error('Incorrect value for property ''alphaValue'' (type ''help analyses.mapAdaptiveSmoothing'' for details).');
                end
                i = i + 2;

            case 'limits'
                limits = varargin{i+1};
                if ~isvector(limits)
                    error('Incorrect value for property ''limits'' (type ''help analyses.mapAdaptiveSmoothing'' for details).');
                end
                i = i + 2;

            otherwise,
                error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help analyses.mapAdaptiveSmoothing'' for details).']);
        end
    end

    % Some info about x, y and z
    if size(pos, 2) > 2
        posx = pos(:, 2);
        posy = pos(:, 3);
    else
        posx = pos(:, 1);
        posy = pos(:, 2);
    end

    if size(spkPos, 2) > 2
        spkx = spkPos(:, 2);
        spky = spkPos(:, 3);
    else
        spkx = spkPos(:, 1);
        spky = spkPos(:, 2);
    end

    if length(nBins) == 1
        nBins(2) = nBins(1);
    end

    % Calculate the border coordinates
    % xStart        Minimum x-coordinate for the path
    %
    % yStart        Minimum y-coordinate for the path
    %
    % xLength       Length of the arena in the x-direction [cm](for cylinder
    %               this equals the diameter)
    % yLength       Length of the arena in the y-direction [cm] (for cylinder
    %               this equals the diameter)
    maxX = nanmax(posx);
    maxY = nanmax(posy);
    xStart = nanmin(posx);
    yStart = nanmin(posy);
    xLength = maxX - xStart;
    yLength = maxY - yStart;

    if haveNBins && ~haveBinWidth
        binWidth = ceil(xLength / nBins(1));
        numColBins = nBins(1);
        numRowBins = nBins(2);
    else
        % Number of bins in each direction of the map
        numColBins = ceil(xLength / binWidth);
        numRowBins = ceil(yLength / binWidth);
    end

    halfBin = binWidth / 2;

    rowAxis = halfBin:binWidth:(numRowBins*binWidth)-halfBin;
    rowAxis = yStart + rowAxis';
    colAxis = halfBin:binWidth:(numColBins*binWidth)-halfBin;
    colAxis = xStart + colAxis';

    maxBins = max([numColBins, numRowBins]);

    map.x = colAxis;
    map.y = rowAxis;
    map.z = zeros(numRowBins, numColBins);
    posPdf = zeros(numRowBins, numColBins);

    if shape(1) == 1
        % Overall clue:
        %     - grow circle from r=1:maxBins (mult. of binWidth), tracking inside
        %     - stop at smallest rad. such that r >= alpha/samples*sqrt(spikes)
        radsqs = ((1:maxBins) * binWidth) .^ 2; % square radius in cm

        for i = 1:numColBins
            binPosX = colAxis(i);

            dist_sample_xdir = (posx - binPosX) .^ 2;
            dist_spike_xdir = (spkx - binPosX) .^ 2;

            for j = 1:numRowBins
                binPosY = rowAxis(j);

                % Calculate sample and spike distances from bin center
                dist_sample = dist_sample_xdir + (posy-binPosY).^2;
                dist_spike = dist_spike_xdir + (spky-binPosY).^2;

                found = 0;
                % Grow circle in increments of binWidth
                for r = 1:maxBins
                   n = nnz(dist_sample <= radsqs(r));
                   s = nnz(dist_spike <= radsqs(r));

                   if r >= alphaValue/(n*sqrt(s))
                       found = 1;
                       break;
                   end
                end

                % Set the rate for this bin
                map.z(j, i) = found * s/(n*sampleTime);
                posPdf(j, i) = found * n*sampleTime;
            end
        end
    else
        for ii = 1:numColBins
            binPosY = (yStart + binWidth/2);
            for jj = 1:numRowBins
                currentPosition = sqrt(binPosX^2 + binPosY^2);
                if currentPosition > shape(2)/2
                    map.z(numRowBins-jj+1, ii) = NaN;
                    posPdf(numRowBins-jj+1, ii) = NaN;
                else
                    n = 0;
                    s = 0;
                    for r = 1:maxBins
                        % Set the current radius of the circle
                        radius = r * binWidth;
                        % Number of samples inside the circle
                        n = insideCircle(binPosX, binPosY, radius, posx, posy);
                        % Number of spikes inside the circle
                        s = insideCircle(binPosX, binPosY, radius, spkx, spky);

                        if r >= alphaValue/(n*sqrt(s))
                            break;
                        end

                    end
                    % Set the rate for this bin
                    map.z(jj,ii) = s/(n*sampleTime);
                    posPdf(jj,ii) = n*sampleTime;

                end
                binPosY = binPosY + binWidth;
            end

            binPosX = binPosX + binWidth;
        end
    end

    map.z(posPdf < minTime) = NaN;
    posPdf = posPdf / nansum(nansum(posPdf));
end

author: kavli-ntnu


text:
function [ Y_hat, W_hat] = MCWNNM_ADMM2_Estimation( NL_mat, Sigma_arr, CurPat, Par )

Y_hat = zeros(size(CurPat));
W_hat    = zeros(size(CurPat));
for  i      =  1 : length(Par.SelfIndex) % For each keypatch group
    Y    =   CurPat(:, NL_mat(1:Par.nlsp,i)); % Non-local similar patches to the keypatch
    mY  =   repmat(mean( Y, 2 ),1,Par.nlsp);
    Y    =   Y-mY;
    X 	=   MCWNNM_ADMM2( Y, Sigma_arr(:, Par.SelfIndex(i)), Par); % WNNM Estimation
    Y_hat(:,NL_mat(1:Par.nlsp,i))  = Y_hat(:,NL_mat(1:Par.nlsp,i))+X+mY;
    W_hat(:,NL_mat(1:Par.nlsp,i))     = W_hat(:,NL_mat(1:Par.nlsp,i))+ones(Par.ps2ch, Par.nlsp);
end
end

author: csjunxu


text:
% DROSRUNODERANKINGS Runs the baseline quadrature ranking presented in the paper
% FORMAT
% DESC Runs the baseline quadrature rankings presented in the
% paper for comparison
%
% COPYRIGHT : Antti Honkela, 2010

% SHEFFIELDML

drosLoadData
z_scores = drosexp.fitmean ./ sqrt(drosexp.fitvar);
active_genes = (mean(z_scores, 2) >= 1.8);

[r_twi, lls_twi, params_twi] = drosODERank(drosexp, drosTF, 'twi', drosexp.probes(active_genes));
r = r_twi; lls = lls_twi; params = params_twi;
save results/dros_mlode_twi_list_expro3_active_genes_results.mat r lls params

[r_mef2, lls_mef2, params_mef2] = drosODERank(drosexp, drosTF, 'mef2', drosexp.probes(active_genes));
r = r_mef2; lls = lls_mef2; params = params_mef2;
save results/dros_mlode_mef2_list_expro3_active_genes_results.mat r lls params

oderank.twi = drosRemoveDuplicateGenes(drosexp, r_twi);
oderank.mef2 = drosRemoveDuplicateGenes(drosexp, r_mef2);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ Velocity_top_lid ] = CreateLidVelocity( alpha,N )
%CREATEALPHAINPUT creates a parabolic velcoity input on the top lid with
%amplitude alpha

% form the grid if not already done it
persistent Grid

if isempty(Grid)
    [Grid]=CollocationGrid_q(N);    % the computational grid
end

Velocity_top_lid = alpha*(1-Grid.x.^2).^2;

end

author: arbabiha


text:
Rotation_Data = get(logsout, 'Rotation_Data');
Translation_Data = get(logsout, 'Translation_Data');
INS_Out = get(logsout, 'INS_Out');

acc_meas_O = get(logsout, 'acc_meas_O_mPs2');
acc_est_O = get(logsout, 'acc_est_O_mPs2');

mag_enable = get(logsout, 'mag_enable');
acc_meas_valid = get(logsout, 'acc_meas_valid');

%% plot attitude data
phi_deg = INS_Out.Values.phi;
phi_deg.Data = single(phi_deg.Data)*2^-12*180/pi;
theta_deg = INS_Out.Values.theta;
theta_deg.Data = single(theta_deg.Data)*2^-12*180/pi;
psi_deg = INS_Out.Values.psi;
psi_deg.Data = single(psi_deg.Data)*2^-12*180/pi;

acc_O_meas_x = acc_meas_O.Values;
acc_O_meas_x.Data = acc_O_meas_x.Data(1,:);
acc_O_meas_y = acc_meas_O.Values;
acc_O_meas_y.Data = acc_O_meas_y.Data(2,:);
acc_O_meas_z = acc_meas_O.Values;
acc_O_meas_z.Data = acc_O_meas_z.Data(3,:);
acc_O_est_x = acc_est_O.Values;
acc_O_est_x.Data = acc_O_est_x.Data(1,:);
acc_O_est_y = acc_est_O.Values;
acc_O_est_y.Data = acc_O_est_y.Data(2,:);
acc_O_est_z = acc_est_O.Values;
acc_O_est_z.Data = acc_O_est_z.Data(3,:);

figure;
ax1 = subplot(3,2,1);
plot(phi_deg, 'Color', 'r');
hold on;
grid on;
legend('\phi')

ax2 = subplot(3,2,3);
plot(theta_deg, 'Color', 'r');
hold on;
grid on;
legend('\theta')

ax3 = subplot(3,2,5);
plot(psi_deg, 'Color', 'r');
hold on;
grid on;
colorAxes(ax3, findLogicalArea(mag_enable.Values.Time, mag_enable.Values.Data > 0), 'green');
colorAxes(ax3, findLogicalArea(acc_meas_valid.Values.Time, acc_meas_valid.Values.Data > 0), 'blue');
legend('\psi')

ax4 = subplot(3,2,2);
plot(acc_O_est_x, 'Color', 'r');
hold on;
plot(acc_O_meas_x, 'Color', 'b');
grid on;
legend('acc\_est\_x', 'acc\_meas\_x')

ax5 = subplot(3,2,4);
plot(acc_O_est_y, 'Color', 'r');
hold on;
plot(acc_O_meas_y, 'Color', 'b');
grid on;
legend('acc\_est\_y', 'acc\_meas\_y')

ax6 = subplot(3,2,6);
plot(acc_O_est_z, 'Color', 'r');
hold on;
plot(acc_O_meas_z, 'Color', 'b');
grid on;
legend('acc\_est\_z', 'acc\_meas\_z')

linkaxes([ax1,ax2,ax3,ax4,ax5,ax6], 'x');

author: Firmament-Autopilot

repo: FMT-Model


text:
function m = r8mat_shortest_path ( n, m )

%*****************************************************************************80
%
%% R8MAT_SHORTEST_PATH computes the shortest distance between all pairs of points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 March 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Floyd,
%    Algorithm 97, Shortest Path,
%    Communications of the ACM,
%    Volume 5, Number 6, June 1962, page 345.
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input/output, real M(N,N).
%    On input, M(I,J) contains the length of the direct link between
%    nodes I and J, or Inf if there is no direct link.
%    On output, M(I,J) contains the distance between nodes I and J,
%    that is, the length of the shortest path between them.  If there
%    is no such path, then M(I,J) will remain Inf.
%
  for i = 1 : n
    for j = 1 : n
      if ( m(j,i) < Inf )
        for k = 1 : n
          if ( m(i,k) < Inf )
            s = m(j,i) + m(i,k);
            if ( s < m(j,k) )
              m(j,k) = s;
            end
          end
        end
      end
    end
  end

  return
end


text:
 function FM = fmeasure(Image, Measure, ROI)
%This function measures the relative degree of focus of
%an image. It may be invoked as:
%
%   FM = fmeasure(Image, Method, ROI)
%
%Where
%   Image,  is a grayscale image and FM is the computed
%           focus value.
%   Method, is the focus measure algorithm as a string.
%           see 'operators.txt' for a list of focus
%           measure methods.
%   ROI,    Image ROI as a rectangle [xo yo width heigth].
%           if an empty argument is passed, the whole
%           image is processed.
%
%  Said Pertuz
%  Abr/2010

if ~isempty(ROI)
    Image = imcrop(Image, ROI);
end

WSize = 15; % Size of local window (only some operators)

switch upper(Measure)
    case 'ACMO' % Absolute Central Moment (Shirvaikar2004)
        if ~isinteger(Image), Image = im2uint8(Image);
        end
        FM = AcMomentum(Image);

    case 'BREN' % Brenner's (Santos97)
        [M N] = size(Image);
        DH = Image;
        DV = Image;
        DH(1:M-2,:) = diff(Image,2,1);
        DV(:,1:N-2) = diff(Image,2,2);
        FM = max(DH, DV);
        FM = FM.^2;
        FM = mean2(FM);

    case 'CONT' % Image contrast (Nanda2001)
        ImContrast = inline('sum(abs(x(:)-x(5)))');
        FM = nlfilter(Image, [3 3], ImContrast);
        FM = mean2(FM);

    case 'CURV' % Image Curvature (Helmli2001)
        if ~isinteger(Image), Image = im2uint8(Image);
        end
        M1 = [-1 0 1;-1 0 1;-1 0 1];
        M2 = [1 0 1;1 0 1;1 0 1];
        P0 = imfilter(Image, M1, 'replicate', 'conv')/6;
        P1 = imfilter(Image, M1', 'replicate', 'conv')/6;
        P2 = 3*imfilter(Image, M2, 'replicate', 'conv')/10 ...
            -imfilter(Image, M2', 'replicate', 'conv')/5;
        P3 = -imfilter(Image, M2, 'replicate', 'conv')/5 ...
            +3*imfilter(Image, M2, 'replicate', 'conv')/10;
        FM = abs(P0) + abs(P1) + abs(P2) + abs(P3);
        FM = mean2(FM);

    case 'DCTE' % DCT energy ratio (Shen2006)
        FM = nlfilter(Image, [8 8], @DctRatio);
        FM = mean2(FM);

    case 'DCTR' % DCT reduced energy ratio (Lee2009)
        FM = nlfilter(Image, [8 8], @ReRatio);
        FM = mean2(FM);

    case 'GDER' % Gaussian derivative (Geusebroek2000)
        N = floor(WSize/2);
        sig = N/2.5;
        [x,y] = meshgrid(-N:N, -N:N);
        G = exp(-(x.^2+y.^2)/(2*sig^2))/(2*pi*sig);
        Gx = -x.*G/(sig^2);Gx = Gx/sum(Gx(:));
        Gy = -y.*G/(sig^2);Gy = Gy/sum(Gy(:));
        Rx = imfilter(double(Image), Gx, 'conv', 'replicate');
        Ry = imfilter(double(Image), Gy, 'conv', 'replicate');
        FM = Rx.^2+Ry.^2;
        FM = mean2(FM);

    case 'GLVA' % Graylevel variance (Krotkov86)
        FM = std2(Image);

    case 'GLLV' %Graylevel local variance (Pech2000)
        LVar = stdfilt(Image, ones(WSize,WSize)).^2;
        FM = std2(LVar)^2;

    case 'GLVN' % Normalized GLV (Santos97)
        FM = std2(Image)^2/mean2(Image);

    case 'GRAE' % Energy of gradient (Subbarao92a)
        Ix = Image;
        Iy = Image;
        Iy(1:end-1,:) = diff(Image, 1, 1);
        Ix(:,1:end-1) = diff(Image, 1, 2);
        FM = Ix.^2 + Iy.^2;
        FM = mean2(FM);

    case 'GRAT' % Thresholded gradient (Snatos97)
        Th = 0; %Threshold
        Ix = Image;
        Iy = Image;
        Iy(1:end-1,:) = diff(Image, 1, 1);
        Ix(:,1:end-1) = diff(Image, 1, 2);
        FM = max(abs(Ix), abs(Iy));
        FM(FM<Th)=0;
        FM = sum(FM(:))/sum(sum(FM~=0));

    case 'GRAS' % Squared gradient (Eskicioglu95)
        Ix = diff(Image, 1, 2);
        FM = Ix.^2;
        FM = mean2(FM);

    case 'HELM' %Helmli's mean method (Helmli2001)
        MEANF = fspecial('average',[WSize WSize]);
        U = imfilter(Image, MEANF, 'replicate');
        R1 = U./Image;
        R1(Image==0)=1;
        index = (U>Image);
        FM = 1./R1;
        FM(index) = R1(index);
        FM = mean2(FM);

    case 'HISE' % Histogram entropy (Krotkov86)
        FM = entropy(Image);

    case 'HISR' % Histogram range (Firestone91)
        FM = max(Image(:))-min(Image(:));

    case 'LAPE' % Energy of laplacian (Subbarao92a)
        LAP = fspecial('laplacian');
        FM = imfilter(Image, LAP, 'replicate', 'conv');
        FM = mean2(FM.^2);

    case 'LAPM' % Modified Laplacian (Nayar89)
        M = [-1 2 -1];
        Lx = imfilter(Image, M, 'replicate', 'conv');
        Ly = imfilter(Image, M', 'replicate', 'conv');
        FM = abs(Lx) + abs(Ly);
        FM = mean2(FM);

    case 'LAPV' % Variance of laplacian (Pech2000)
        LAP = fspecial('laplacian');
        ILAP = imfilter(Image, LAP, 'replicate', 'conv');
        FM = std2(ILAP)^2;

    case 'LAPD' % Diagonal laplacian (Thelen2009)
        M1 = [-1 2 -1];
        M2 = [0 0 -1;0 2 0;-1 0 0]/sqrt(2);
        M3 = [-1 0 0;0 2 0;0 0 -1]/sqrt(2);
        F1 = imfilter(Image, M1, 'replicate', 'conv');
        F2 = imfilter(Image, M2, 'replicate', 'conv');
        F3 = imfilter(Image, M3, 'replicate', 'conv');
        F4 = imfilter(Image, M1', 'replicate', 'conv');
        FM = abs(F1) + abs(F2) + abs(F3) + abs(F4);
        FM = mean2(FM);

    case 'SFIL' %Steerable filters (Minhas2009)
        % Angles = [0 45 90 135 180 225 270 315];
        N = floor(WSize/2);
        sig = N/2.5;
        [x,y] = meshgrid(-N:N, -N:N);
        G = exp(-(x.^2+y.^2)/(2*sig^2))/(2*pi*sig);
        Gx = -x.*G/(sig^2);Gx = Gx/sum(Gx(:));
        Gy = -y.*G/(sig^2);Gy = Gy/sum(Gy(:));
        R(:,:,1) = imfilter(double(Image), Gx, 'conv', 'replicate');
        R(:,:,2) = imfilter(double(Image), Gy, 'conv', 'replicate');
        R(:,:,3) = cosd(45)*R(:,:,1)+sind(45)*R(:,:,2);
        R(:,:,4) = cosd(135)*R(:,:,1)+sind(135)*R(:,:,2);
        R(:,:,5) = cosd(180)*R(:,:,1)+sind(180)*R(:,:,2);
        R(:,:,6) = cosd(225)*R(:,:,1)+sind(225)*R(:,:,2);
        R(:,:,7) = cosd(270)*R(:,:,1)+sind(270)*R(:,:,2);
        R(:,:,7) = cosd(315)*R(:,:,1)+sind(315)*R(:,:,2);
        FM = max(R,[],3);
        FM = mean2(FM);

    case 'SFRQ' % Spatial frequency (Eskicioglu95)
        Ix = Image;
        Iy = Image;
        Ix(:,1:end-1) = diff(Image, 1, 2);
        Iy(1:end-1,:) = diff(Image, 1, 1);
        FM = mean2(sqrt(double(Iy.^2+Ix.^2)));

    case 'TENG'% Tenengrad (Krotkov86)
        Sx = fspecial('sobel');
        Gx = imfilter(double(Image), Sx, 'replicate', 'conv');
        Gy = imfilter(double(Image), Sx', 'replicate', 'conv');
        FM = Gx.^2 + Gy.^2;
        FM = mean2(FM);

    case 'TENV' % Tenengrad variance (Pech2000)
        Sx = fspecial('sobel');
        Gx = imfilter(double(Image), Sx, 'replicate', 'conv');
        Gy = imfilter(double(Image), Sx', 'replicate', 'conv');
        G = Gx.^2 + Gy.^2;
        FM = std2(G)^2;

    case 'VOLA' % Vollath's correlation (Santos97)
        Image = double(Image);
        I1 = Image; I1(1:end-1,:) = Image(2:end,:);
        I2 = Image; I2(1:end-2,:) = Image(3:end,:);
        Image = Image.*(I1-I2);
        FM = mean2(Image);

    case 'WAVS' %Sum of Wavelet coeffs (Yang2003)
        [C,S] = wavedec2(Image, 1, 'db6');
        H = wrcoef2('h', C, S, 'db6', 1);
        V = wrcoef2('v', C, S, 'db6', 1);
        D = wrcoef2('d', C, S, 'db6', 1);
        FM = abs(H) + abs(V) + abs(D);
        FM = mean2(FM);

    case 'WAVV' %Variance of  Wav...(Yang2003)
        [C,S] = wavedec2(Image, 1, 'db6');
        H = abs(wrcoef2('h', C, S, 'db6', 1));
        V = abs(wrcoef2('v', C, S, 'db6', 1));
        D = abs(wrcoef2('d', C, S, 'db6', 1));
        FM = std2(H)^2+std2(V)+std2(D);

    case 'WAVR'
        [C,S] = wavedec2(Image, 3, 'db6');
        H = abs(wrcoef2('h', C, S, 'db6', 1));
        V = abs(wrcoef2('v', C, S, 'db6', 1));
        D = abs(wrcoef2('d', C, S, 'db6', 1));
        A1 = abs(wrcoef2('a', C, S, 'db6', 1));
        A2 = abs(wrcoef2('a', C, S, 'db6', 2));
        A3 = abs(wrcoef2('a', C, S, 'db6', 3));
        A = A1 + A2 + A3;
        WH = H.^2 + V.^2 + D.^2;
        WH = mean2(WH);
        WL = mean2(A);
        FM = WH/WL;
    otherwise
        error('Unknown measure %s',upper(Measure))
end
 end
%************************************************************************
function fm = AcMomentum(Image)
[M N] = size(Image);
Hist = imhist(Image)/(M*N);
Hist = abs((0:255)-255*mean2(Image))'.*Hist;
fm = sum(Hist);
end

%******************************************************************
function fm = DctRatio(M)
MT = dct2(M).^2;
fm = (sum(MT(:))-MT(1,1))/MT(1,1);
end

%************************************************************************
function fm = ReRatio(M)
M = dct2(M);
fm = (M(1,2)^2+M(1,3)^2+M(2,1)^2+M(2,2)^2+M(3,1)^2)/(M(1,1)^2);
end
%******************************************************************


text:
function [ intersect, p ] = plane_imp_line_par_int_3d ( a, b, c, d, x0, y0, z0, ...
  f, g, h )

%*****************************************************************************80
%
%% PLANE_IMP_LINE_PAR_INT_3D: intersection ( implicit plane, parametric line ) in 3D.
%
%  Discussion:
%
%    The implicit form of a plane in 3D is:
%
%      A * X + B * Y + C * Z + D = 0
%
%    The parametric form of a line in 3D is:
%
%      X = X0 + F * T
%      Y = Y0 + G * T
%      Z = Z0 + H * T
%
%    We normalize by choosing F*F+G*G+H*H=1 and 0 <= F.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Adrian Bowyer and John Woodwark,
%    A Programmer's Geometry,
%    Butterworths, 1983, page 111.
%
%  Parameters:
%
%    Input, real A, B, C, D, the implicit plane parameters.
%
%    Input, real X0, Y0, Z0, F, G, H, parameters that define the
%    parametric line.
%
%    Output, logical INTERSECT, is TRUE if the line and the plane
%    intersect.
%
%    Output, real P(3), is a point of intersection of the line
%    and the plane, if INTERSECT is TRUE.
%
  dim_num = 3;
  tol = 0.00001;
%
%  Check.
%
  norm1 = sqrt ( a * a + b * b + c * c );

  if ( norm1 == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PLANE_IMP_LINE_PAR_INT_3D - Fatal error!\n' );
    fprintf ( 1, '  The plane normal vector is null.\n' );
    error ( 'PLANE_IMP_LINE_PAR_INT_3D - Fatal error!' )
  end

  norm2 = sqrt ( f * f + g * g + h * h );

  if ( norm2 == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PLANE_IMP_LINE_PAR_INT_3D - Fatal error!\n' );
    fprintf ( 1, '  The line direction vector is null.\n' );
    error ( 'PLANE_IMP_LINE_PAR_INT_3D - Fatal error!' )
  end

  denom = a * f + b * g + c * h;
%
%  The line and the plane may be parallel.
%
  if ( abs ( denom ) < tol * norm1 * norm2 )

    if ( a * x0 + b * y0 + c * z0 + d == 0.0 )
      intersect = 1;
      p(1) = x0;
      p(2) = y0;
      p(3) = z0;
    else
      intersect = 0;
      p(1:dim_num) = 0.0;
    end
%
%  If they are not parallel, they must intersect.
%
  else

    intersect = 1;
    t = - ( a * x0 + b * y0 + c * z0 + d ) / denom;
    p(1) = x0 + t * f;
    p(2) = y0 + t * g;
    p(3) = z0 + t * h;

  end

  return
end


text:
function [noise_samp, n_noise_samp, noise_seg, D, Dsmth, snre_vad, e]=snre_highenergy(dfdata, nfr10, flen, fsh10, ENERGYFLOOR, pv01)

% Ref:
%   Z-H Tan and B Lindberg
%   Low-Complexity Variable Frame Rate Analysis for Speech Recognition and Voice Activity Detection
%   IEEE Journal of Selected Topics in Signal Processing, 4(5), Oct. 2010.

% square root of a posteriori SNR weighted engergy difference, block based
%
% Modified 01 Mar 2013

Dexpl=18;
Dexpr=18;
segThres = 0.25;
e=zeros(nfr10,1);
for i=1:nfr10
    for j=1:flen
        e(i)=e(i)+dfdata((i-1)*fsh10+j)*dfdata((i-1)*fsh10+j);
    end
    if e(i) <= ENERGYFLOOR
        e(i)=ENERGYFLOOR;
    end
end

emin=ones(nfr10,1);
NESEG = 200;
if nfr10 < NESEG; NESEG=nfr10; end
for i=1:floor(nfr10/NESEG)
	[eY,eI]=sort(e((i-1)*NESEG+1:i*NESEG));
	emin((i-1)*NESEG+1:i*NESEG)=eY(floor(NESEG*0.1));
    if i~=1
        emin((i-1)*NESEG+1:i*NESEG)=0.9*emin((i-1)*NESEG)+0.1*emin((i-1)*NESEG+1);
    end
end
if i*NESEG~=nfr10
    [eY,eI]=sort(e((i-1)*NESEG+1:nfr10));
    emin(i*NESEG+1:nfr10)=eY(floor((nfr10-(i-1)*NESEG)*0.1));
    emin(i*NESEG+1:nfr10)=0.9*emin(i*NESEG)+0.1*emin(i*NESEG+1);
end

D=zeros(nfr10,1);
postsnr=zeros(nfr10,1);
for i=2:nfr10
    postsnr(i) =log10(e(i))-log10(emin(i));
    if postsnr(i)<0
        postsnr(i)=0;
    end
    D(i)=sqrt(abs(e(i)-e(i-1))*postsnr(i));
end
D(1)=D(2);

Dexp = vertcat(ones(Dexpl,1)*D(1), D, ones(Dexpr,1)*D(nfr10));
Dsmth = zeros(nfr10,1);
for i=1:nfr10
    Dsmth(i)=sum(Dexp(i:i+Dexpl+Dexpr));
end

for i=1:floor(nfr10/NESEG)
    Dsmth_max((i-1)*NESEG+1:i*NESEG)=max(e((i-1)*NESEG+1:i*NESEG));
end
if i*NESEG~=nfr10
    Dsmth_max(i*NESEG+1:nfr10)=max(e((i-1)*NESEG+1:nfr10));
end

snre_vad = zeros(nfr10,1);
for i=1:nfr10
   if Dsmth(i)>Dsmth_max(i)*segThres; snre_vad(i)=1; end
end

% block based processing to remove noise part by using snre_vad1.
sign_vad = 0;
noise_seg=zeros(floor(nfr10/1.6),1);
noise_samp=zeros(nfr10,2);
n_noise_samp=0;
for i=1:nfr10
    if snre_vad(i) == 1 && sign_vad == 0 % start of a segment
        sign_vad = 1;
        nstart=i;
    elseif (snre_vad(i) ==0 || i==nfr10) && sign_vad == 1 % end of a segment
        sign_vad = 0;
        nstop=i-1;
        if sum(pv01(nstart:nstop))==0
            noise_seg(round(nstart/1.6):floor(nstop/1.6)) = 1;
            n_noise_samp=n_noise_samp+1;
            noise_samp(n_noise_samp,:)=[(nstart-1)*fsh10+1 nstop*fsh10];
        end
    end
end
noise_samp(n_noise_samp+1:nfr10,:)=[];

author: zhenghuatan

repo: rVAD

save_path: github-repos/MATLAB/zhenghuatan-rVAD


text:
function k = kernDiagGradX(kern, x)

% KERNDIAGGRADX Compute the gradient of the  kernel wrt X.
% FORMAT
% DESC computes the gradient of the diagonal of the kernel matrix
% with respect to the elements of the design matrix given in X.
% ARG kern : the kernel structure for which gradients are being computed.
% ARG X : the input data in the form of a design matrix.
% RETURN gX : the gradients of the diagonal with respect to each
% element of X. The returned matrix has the same dimensions as X.
%
% SEEALSO : kernDiagGradX, kernGradX

% KERN

fhandle = str2func([kern.type 'KernDiagGradX']);
k = fhandle(kern, x);


text:
function value = j_double_product_integral ( i, j, a, b )

%*****************************************************************************80
%
%% J_DOUBLE_PRODUCT_INTEGRAL: integral of J(i,x)*J(j,x)*(1-x)^a*(1+x)^b.
%
%  Discussion:
%
%    VALUE = integral ( -1 <= x <= +1 ) J(i,x)*J(j,x)*(1-x)^a*(1+x)^b dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, J, the polynomial indices.
%
%    Input, real A, B, the parameters.
%    -1 < A, B.
%
%    Output, real VALUE, the value of the integral.
%
  if ( i ~= j )
    value = 0.0;
  else
    value = 2^( a + b + 1.0 ) / ( 2 * i + a + b + 1 ) ...
      * gamma ( i + a + 1 ) * gamma ( i + b + 1 ) ...
      / r8_factorial ( i ) / gamma ( i + a + b + 1 );
  end

  return
end


text:
function  FileC = file2cell(FileNameS, opt)
%Put a file into a cellarray.  Each cell is a
%line converted to a character string.
%If opt is set to 'spreadsheet', each word on each line becomes a cell within a
%cell array for each line, which is in a cell array of all lines.
%If opt is not present, each line becomes a character cell in a cell array.
%copyright (c) 2001, J.O. Deasy and Washington University in St. Louis.
%Use is granted for non-commercial and non-clinical applications.
%No warranty is expressed or implied for any use whatever.
%
%LM:  14 Feb 02, JOD.
%   13 Mar 03, JOD.  Added 'spreadsheet' option.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

spread = 0;
if nargin > 1
  if strcmpi(opt,'spreadsheet')
    spread = 1;
  end
end

try
  fid = fopen(FileNameS,'r');
  FileV = fread(fid);
  fclose(fid);
catch
  FileC = {}
  warning(['Attempt to open ' FileNameS ' failed!'])
  return
end

FileC={};

if ~isempty(FileV)

  if any([FileV == 13])
    Newline = 13; %Define value of newline character; start assuming PC value, but check for UNIX below.
    FileV(FileV == 10) = []; %remove UNIX 'corruption'
  else
    Newline = 10;  %raw UNIX
  end

  %Break into lines.  Find all the newline characters.
  [LocationV] = find(FileV == Newline);

  %Reprocess to a cellarray
  Start = 1;
  for i = 1:length(LocationV)
    Stop = LocationV(i);

    strV = char(FileV(Start:Stop-1))';
    if spread == 0
      FileC{i} = strV;
    elseif spread == 1
      FileC{i} = str2Cell(strV, 'convert');
    end

    Start = Stop + 1;
  end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%% =======================================================================
%  Start Here Script
%  =======================================================================
%
%  This script guides the user through examples of both simulations and
%  Wi-Fi control models for the Parrot ARDrone.
%  Authors:
%       David Escobar Sanabria -> descobar@aem.umn.edu
%       Pieter J. Mosterman -> pieter.mosterman@mathworks.com
%  =======================================================================

%%
%  Cleaning workspace
bdclose all;
clear all;
bdclose all;
clc

%%
disp('This script guides the user through examples of both simulations and');
disp('Wi-Fi control Simulink models for the Parrot ARDrone ');
disp(' ');

disp('Swlect one of the following options:');
disp('    (1) Simulation');
disp('    (2) Wi-Fi control. The computer shoud be connected to the drone');

option = input('');

%%
switch option
    case 1

        disp('Swlect one of the following options for simulation:');
        disp('    (1) Baseline simulation: The ARDrone block with inputs and scopes to visualize outputs');
        disp('    (2) Hover: Vehicle is held at constant position');
        disp('    (3) Waypoint tracking: Vehicle tracks a list of waypoints');
        option2 = input(' ');

        switch option2
            case 1
                cd simulation;
                setupBaseModel;
            case 2
                cd simulation;
                setupHoverSim;
            case 3
                cd simulation;
                setupWPTrackingSim;
            otherwise
                disp('An incorrect option was selected')

        end

    case 2
        disp('Swlect one of the following options for Wi-Fi control:');
        disp('    (1) Hover: Vehicle is held at constant position');
        disp('    (2) Waypoint tracking: Vehicle tracks a list of waypoints');
        option2 = input(' ');

        switch option2
            case 1
                cd wifiControl;
                setupHover;
                % Building model using RTWT. Install rtwt if not installed
                % using 'rtwintgt -setup'
                rtwbuild('ARDroneHover');
            case 2
                cd wifiControl;
                setupWPTracking;
                % Building model using RTWT. Install rtwt if not installed
                % using 'rtwintgt -setup'
                rtwbuild('ARDroneWPTracking');
            otherwise
                disp('An incorrect option was selected')
        end

    otherwise
       disp('An incorrect option was selected')

end


text:
function [success] = writeMatfileToHDF5(inputMatfilePath,varName,outputFilename,varargin)
	% [success] = writeMatfileToHDF5(inputMatfilePath,varName,outputFilename,varargin)
	%
	% Writes a MAT-file dataset to HDF5 by reading the matrix in parts to avoid large overhead.
	%
	% Biafra Ahanonu
	% started: 2022.03.13 [23:49:18]
	%
	% inputs
	% 	inputMatfilePath - path to matfile object containing data.
	% 	varName - Str: name of variable in MAT-file to save. Should be matrix of [x y frames].
	% 	inputFilename - Str: path where HDF5 file should be saved.
	%
	% outputs
	% 	success - Binary: 1 = data saved correctly, 0 = data not saved correctly.

	% changelog
		%
	% TODO
		%

	% ========================
	% Str: name of HDF5 dataset to save data into.
	options.datasetname = '/1';
	% HDF5: append (don't blank HDF5 file) or new (blank HDF5 file)
	options.writeMode = 'new';
	% Int: Number of frames to use for chunking the data to save in parts.
	options.chunkSize = 200;
	% Int: Defines gzip compression level (0-9). 0 = no compression, 9 = most compression.
	options.deflateLevel = 9;
	% Int: chunk size in [x y z] of the dataset, leave empty for auto chunking
	options.dataDimsChunkCopy = [];
	% Struct: structure of information to add. Will create a HDF5 file
	options.addInfo = [];
	% Str: e.g. '/movie/processingSettings'
	options.addInfoName = '';
	% get options
	options = ciapkg.io.getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	% ========================

	try
		success = 0;

		% Create matfile object.
		inputObj = matfile(inputMatfilePath);

		% Get size of the input data, call this way instead of size(inputObj.varName) to avoid loading entire dataset into memory.
		dataDims = size(inputObj,varName);

		if length(dataDims)~=3
			disp('Data is not a 3D matrix of size [:, :, >1], returning...')
			return;
		end

		% Get coordinates to use.
		nFrames = dataDims(3);
		subsetSize = options.chunkSize;
		movieLength = nFrames;
		numSubsets = ceil(movieLength/subsetSize)+1;
		subsetList = round(linspace(1,movieLength,numSubsets));
		nSubsets = (length(subsetList)-1);

		% Write data out in chunks.
		for thisSet = 1:nSubsets
			% subsetStartTime = tic;

			subsetStartIdx = subsetList(thisSet);
			subsetEndIdx = subsetList(thisSet+1);
			disp(repmat('$',1,7))
			if thisSet==nSubsets
				movieSubset = subsetStartIdx:subsetEndIdx;
			else
				movieSubset = subsetStartIdx:(subsetEndIdx-1);
			end
			disp([num2str(movieSubset(1)) '-' num2str(movieSubset(end)) ' ' num2str(thisSet) '/' num2str(nSubsets)])

			% Slice into the desired variable using dynamic field references to avoid loading entire dataset into memory.
			inputDataSlice = inputObj.(varName)(:,:,movieSubset);

			if thisSet==1
				ciapkg.hdf5.createHdf5File(outputFilename, options.datasetname, inputDataSlice,'deflateLevel',options.deflateLevel,'dataDimsChunkCopy',options.dataDimsChunkCopy);
			else
				ciapkg.hdf5.appendDataToHdf5(outputFilename, options.datasetname, inputDataSlice);
			end

			% toc(subsetStartTime)
		end

		% add information about data to HDF5 file
		if strcmp(options.writeMode,'new')
			hdf5write(outputFilename,'/movie/info/dimensions',dataDims,'WriteMode','append');
			currentDateTimeStr = datestr(now,'yyyymmdd_HHMM','local');
			hdf5write(outputFilename,'/movie/info/date',currentDateTimeStr,'WriteMode','append');
			hdf5write(outputFilename,'/movie/info/savePath',outputFilename,'WriteMode','append');
			hdf5write(outputFilename,'/movie/info/Deflate',options.deflateLevel,'WriteMode','append');
		end
		if ~isempty(options.addInfo)
			if ~iscell(options.addInfo)
				options.addInfo = {options.addInfo};
				options.addInfoName = {options.addInfoName};
			end
			addInfoLen = length(options.addInfo);
			for addInfoStructNo = 1:addInfoLen
				thisAddInfo = options.addInfo{addInfoStructNo};
				infoList = fieldnames(thisAddInfo);
				nInfo = length(infoList);
				for fieldNameNo = 1:nInfo
					thisField = infoList{fieldNameNo};
					hdf5write(outputFilename,[options.addInfoName{addInfoStructNo} '/' thisField],thisAddInfo.(thisField),'WriteMode','append');
					% h5write(fileSavePath,[options.addInfoName{addInfoStructNo} '/' thisField],thisAddInfo.(thisField),'WriteMode','append');
				end
			end
		end

		success = 1;
	catch err
		success = 0;
		disp(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		disp(repmat('@',1,7))
	end

	% function [outputs] = nestedfxn_exampleFxn(arg)
	% 	% Always start nested functions with "nestedfxn_" prefix.
	% 	% outputs = ;
	% end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah