text:
function rchy = realized_hayashi_yoshida(priceA,timeA,priceB,timeB,timeType,samplingType,samplingInterval,overlap)
% Computed the Hayashi-Yoshida estimator of quadratic covariation, and allows for
% the empirical-performance motivated K-lead-and-lag version similar to Drost and Nijman (1997).
%
% USAGE:
%   [RCHY] = realized_hayashi_yoshida_(PRICEA,TIMEA,PRICEB,TIMEB,TIMETYPE,SAMPLINGTYPE,SAMPLINGINTERVAL,OVERLAP)
%
% INPUTS:
%   PRICEA            - mA by 1 vector of high frequency prices
%   TIMEA             - mA by 1 vector of times where TIMEB(i) corresponds to PRICEA(i)
%   PRICEA            - mB by 1 vector of high frequency prices
%   TIMEA             - mB by 1 vector of times where TIMEB(i) corresponds to PRICEB(i)
%   TIMETYPE          - String describing the way times are measured
%                        'wall'    24-hour clock of the form HHMMSS, e.g. 101543 or 153217
%                        'seconds' Time measures in seconds past midnight.
%                          TIME must satisfy 0<=TIME<86400
%                        'unit'  Unit normalized date format, e.g. .1, .234, .9
%                          Unit normalized times are more general than the
%                          other types and can be applied to data from more
%                          than one calendar day
%   SAMPLINGTYPE      - String describing the type of sampling to use when
%                         filtering PRICE
%                         'CalendarTime' - Sample in calendar time using
%                           observations separated by SAMPLINGINTERVAL
%                           seconds. If TIMETYPE is 'unit',
%                           SAMPLINGINTERVAL must be between 0 and 1 and
%                           represents the fraction of the sample to skip
%                           when sampling.
%                         'CalendarUniform' - Sample in calendar time using
%                           SAMPLINGINTERVAL observations spread uniformly
%                           between TIME(1) and TIME(m)
%                         'BusinessTime' - Sample in business (tick) time
%                           using observation separated by SAMPLINGINTERVAL
%                           ticks
%                         'BusinessUniform' - Sample in business (tick)
%                           time using observations uniformly spaced in
%                           business time.
%                         'Fixed' - Sample at specific points in time. When
%                           using fixed, SAMPLINGINTERVAL must be a n by 1 vector
%                           of times with the same TIMETYPE as TIME (i.e.
%                           seconds if TIME is in seconds)
%   SAMPLINGINTERVAL   - Scalar integer or n by 1 vector whose meaning depends on the
%                          selected SAMPLINGTYPE
%   OVERLAP            - Number of ticks to overlap when computing the HY estimator.  The original HY
%                          estimator uses 0, which corresponds to the maximum likelihood estimator for
%                          price processes whose observation times are driven by independent Poisson
%                          processes.  Empirically this estimator performs poorly because prices are
%                          not a vector semi-martingale and using a larger number of lags can alleviate
%                          this problem.
%
% OUTPUTS:
%   RCHY             - The K-lead-and-lag Hayashi-Yoshida covariance estimator
%
% COMMENTS:
%   Filtering the price in calendar time allow the creation of a Hayashi-Yoshida corrected
%   calendar-time sampled realized covariance.  The value in SAMPLINGTYPE is applied to price A and
%   then the realized HY estimator is computed using the filtered price of A and the filtered times
%   of A, and all observations of B.  Price filtering is done using realized_price_filter
%
% EXAMPLES:
%   % Standard use with all prices with wall time prices
%   RCHY = realized_hayashi_yoshida(PRICEA,TIMEA,PRICEB,TIMEB,'wall','BusinessTime',1,0)
%
%   % 10 lead and lag RCHY with all prices with wall time prices
%   RCHY = realized_hayashi_yoshida(PRICEA,TIMEA,PRICEB,TIMEB,'wall','BusinessTime',1,10)
%
%   % 1-minute realized covariance with a HY correction
%   RCHY = realized_hayashi_yoshida(PRICEA,TIMEA,PRICEB,TIMEB,'wall','CalendarTime',60,0)
%
%  See also REALIZED_MULTIVARIATE_KERNEL, REALIZED_COVARIANCE, REALIZED_KERNEL, REALIZED_VARIANCE,
%  REALIZED_RANGE, REALIZED_QUANTILE_VARIANCE

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 5/1/2008

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<7 || nargin>8
    error('Seven or eight inputs required.')
end
if size(priceA,2)>size(priceA,1)
    priceA=priceA';
end
if size(priceA,2)>1
    error('PRICEA must be a m by 1 vector.')
end
if size(timeA,2)>size(timeA,1)
    timeA=timeA';
end
if any(diff(timeA)<0)
    error('TIMEA must be sorted and increasing')
elseif any(diff(timeA)==0)
    warning('oxfordRealized:realizedPriceFilter','TIMEA contains multiple entries with the same value. This creates an ambiguity and FILTEREDPRICE will contain the last price if TIMEA does not only unique elements.')
end
if size(timeA,2)>1 || length(timeA)~=length(priceA)
    error('TIMEA must be a m by 1 vector.')
end
% Inserted to protect against inputing integer times
timeA = double(timeA);
if size(priceB,2)>size(priceB,1)
    priceB=priceB';
end
if size(priceB,2)>1
    error('PRICEB must be a m by 1 vector.')
end
if size(timeB,2)>size(timeB,1)
    timeB=timeB';
end
if any(diff(timeB)<0)
    error('TIMEB must be sorted and increasing')
elseif any(diff(timeB)==0)
    warning('oxfordRealized:realizedPriceFilter','TIMEB contains multiple entries with the same value. This creates an ambiguity and FILTEREDPRICE will contain the last price if TIMEB does not only unique elements.')
end
if size(timeB,2)>1 || length(timeB)~=length(priceB)
    error('TIMEB must be a m by 1 vector.')
end
% Inserted to protect against inputing integer times
timeB = double(timeB);

% make sure the intersection is non-empty
if ~(min(timeB)<max(timeA) || min(timeA)<max(timeB))
    warning('oxfordRealized:realizedHayashiYoshida','The intersection of TIMESA and TIMESB is empty.  The RCHY will necessarily be 0.');
    intersectionIsEmpty = true;
else
    intersectionIsEmpty = false;
end

timeType=lower(timeType);
if ~ismember(timeType,{'wall','seconds','unit'})
    error('TIMETYPE must be one of ''wall'', ''seconds'' or ''unit''.');
end

samplingType=lower(samplingType);
if ~ismember(samplingType,{'calendartime','calendaruniform','businesstime','businessuniform','fixed'})
    error('SAMPLINGTYPE must be one of ''CalendarTime'', ''CalendarUniform'', ''BusinessTime'', ''BusinessUniform'' or ''Fixed''.');
end

m=size(priceA,1);
t0Original=timeA(1);
tTOriginal=timeA(m);
if ismember(samplingType,{'calendartime','calendaruniform','businesstime','businessuniform'})
    % Must be a scalar integer, unless using unit times
    if (~isscalar(samplingInterval) || floor(samplingInterval)~=samplingInterval || samplingInterval<1) && ~strcmp(timeType,'unit')
        error('SAMPLINGINTERVAL must be a positive integer for the SAMPLINGTYPE selected.')
    end
else
    if size(samplingInterval,2)>size(samplingInterval,1)
        samplingInterval=samplingInterval';
    end
    if ~(any(samplingInterval>=t0Original) && any(samplingInterval<=tTOriginal))
        error('At least one sampling interval must be between min(TIME) and max(TIME) when using ''Fixed'' as SAMPLINGTYPE.')
    end
    if any(diff(samplingInterval)<=0)
        error('When using ''Fixed'' as SAMPLINGTYPE the vector of sampling times in SAMPLINGINTERVAL must be sorted and strictly increasing.')
    end
end

if strcmp(timeType,'unit') && strcmp(samplingType,'calendartime')
    % samplingInterval must be between 0 and 1
    if samplingInterval>1
        error('When TIMETYPE is ''unit'' and SAMPLINGTYPE is ''CalendarTime'', SAMPLINGINTERVAL must also be in ''unit'' terms, and so must be between 0 and 1')
    end
end

if nargin==7
    overlap = 0;
elseif overlap<0 || floor(overlap)~=overlap || max(size(overlap))>1
    error('OVERLAP must be a non-negative integer.')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Price A is the base price, price B is the other price.  First sample priceA according to
% samplingType and samplingInterval, then use the actual times of these to estimate the HY
% respecting the value in overlap

if ~intersectionIsEmpty
    % First filter the price
    [filteredPriceA,filteredTimeA,actualTimeA] = realized_price_filter(priceA,timeA,timeType,samplingType,samplingInterval);

    % Then call the core routine
    rchy = realized_hayashi_yoshida_core(filteredPriceA,actualTimeA,priceB,timeB,overlap);
else
    rchy = 0;
end

function [rchy,times] = realized_hayashi_yoshida_core(priceA,timeA,priceB,timeB,overlap)
% Core routine for computing the Hayashi-Yoshida estimator of quadratic covariation, and allows for
% the empirical-performance motivated K-lead-and-lag version similar to Drost and Nijman (1997).
%
% USAGE:
%   [RCHY,TIMES] = realized_hayashi_yoshida_core(PRICEA,TIMEA,PRICEB,TIMEB,OVERLAP)
%
% INPUTS:
%   PRICEA           - mA by 1 vector of high frequency prices
%   TIMEA            - mA by 1 vector of times where TIMEB(i) corresponds to PRICEA(i)
%   PRICEA           - mB by 1 vector of high frequency prices
%   TIMEA            - mB by 1 vector of times where TIMEB(i) corresponds to PRICEB(i)
%   OVERLAP          - [OPTIONA] Number of ticks to overlap when computing the HY estimator.  The
%                        original HY estimator uses 0, which corresponds to the maximum likelihood
%                        estimator for price processes whose observation times are driven by
%                        independent Poisson processes.  Empirically this estimator performs poorly
%                        because prices are not a vector semi-martingale and using a larger number
%                        of lags can alleviate this problem. If omitted OVERLAP = 0.
%
% OUTPUTS:
%   RCHY             - The K-lead-and-lag Hayashi-Yoshida covariance estimator
%   TIMES            - A mA by 1 matrix of time stamps contains the times where the prices were
%                        sampled for computing the returns in the HY estimator.  This is mostly for
%                        diagnostic purposes.
%
% COMMENTS:
%   This is a helper function of realized_hayashi_yoshida and does no input checking.  In general
%   it should not be directly called.
%
%  See also REALIZED_MULTIVARIATE_KERNEL, REALIZED_COVARIANCE, REALIZED_KERNEL, REALIZED_VARIANCE,
%  REALIZED_RANGE, REALIZED_QUANTILE_VARIANCE

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 5/1/2008

priceA = log(priceA);
priceB = log(priceB);
timeA = double(timeA);
timeB = double(timeB);
% Price A is the base price, price B is the one which will move

% First find the first time that B is available before timeA(1).  If it is empty then find the first
% timeA(1) which is weakly after the timeB(1)
nA = size(priceA,1);
nB = size(priceB,1);

if timeA(1)<timeB(1)
    % A starts before B, so the easy solution is to project backward the price of the first B to the
    % time of the first A.  This will generate a 0 return but makes the algorithm easier.
    timeB=[timeA(1);timeB];
    priceB=[priceB(1);priceB];
    nB = size(priceB,1);
end

if timeA(nA)>timeB(nB)
    % A ends after B, so the easy solution is to project forward the last price B to the
    % time of the last A.  This will generate a 0 return but makes the algorithm easier.
    timeB=[timeB;timeA(nA)];
    priceB=[priceB;priceB(nB)];
    nB = size(priceB,1);
end

% Initialize the indices for A and B
indexA = 1;
indexB = find(timeB<=timeA(1), 1,'last');
if overlap>0
    indexB = max(indexB-overlap,1); % Make sure that indexB is >= 1
end

% Initialize rchy and the times
rchy = 0;
times = zeros(nA-1,4);
while indexA<nA
    % Get the two prices at indexA and indexA+1
    pA1 = priceA(indexA);
    pA2 = priceA(indexA+1);
    % get the times
    times(indexA-1,1) = timeA(indexA);
    times(indexA-1,2) = timeA(indexA+1);
    % Increment the index to A
    indexA = indexA+1;

    % Compute the index, including overlap if needed, make sure it doesn't end up before the first
    % index
    indexB_minus_overlap = max(indexB - overlap,1);
    % Get the first price of B
    pB1 = priceB(indexB_minus_overlap);
    % Get the time
    times(indexA-1,3) = timeB(indexB_minus_overlap);
    while timeB(indexB)<timeA(indexA)
        % The time stamp of B is strictly less than the time stamp of the second price in the return
        % to asset A so increment the counter.  Since we do the pre- and post- pend trick above,
        % there is no need to worry about uneven ending.
        indexB = indexB + 1;
    end
    % Compute the index plus the overlap, make sure it doesn't go outsie the number of obs
    indexB_plus_overlap = min(indexB + overlap,nB);
    % Get the second price
    pB2 = priceB(indexB_plus_overlap);
    % Get the final time
    times(indexA-1,4) = timeB(indexB_plus_overlap);
    % Add to rchy
    rchy = rchy + (pA2-pA1)*(pB2-pB1);

    % Rewind indexB by 1 since it is now after A
    if timeB(indexB)>timeA(indexA)
        indexB = indexB - 1;
    end
end

author: bashtage

repo: mfe-toolbox


text:
function [M, TF] = padcat(varargin)
% PADCAT - concatenate vectors with different lengths by padding with NaN
%
%   M = PADCAT(V1, V2, V3, ..., VN) concatenates the vectors V1 through VN
%   into one large matrix. All vectors should have the same orientation,
%   that is, they are all row or column vectors. The vectors do not need to
%   have the same lengths, and shorter vectors are padded with NaNs.
%   The size of M is determined by the length of the longest vector. For
%   row vectors, M will be a N-by-MaxL matrix and for column vectors, M
%   will be a MaxL-by-N matrix, where MaxL is the length of the longest
%   vector.
%
%   Examples:
%      a = 1:5 ; b = 1:3 ; c = [] ; d = 1:4 ;
%      padcat(a,b,c,d) % row vectors
%         % ->   1     2     3     4     5
%         %      1     2     3   NaN   NaN
%         %    NaN   NaN   NaN   NaN   NaN
%         %      1     2     3     4   NaN
%      CC = {d.' a.' c.' b.' d.'} ;
%      padcat(CC{:}) % column vectors
%         %      1     1   NaN     1     1
%         %      2     2   NaN     2     2
%         %      3     3   NaN     3     3
%         %      4     4   NaN   NaN     4
%         %    NaN     5   NaN   NaN   NaN
%
%   [M, TF] = PADCAT(..) will also return a logical matrix TF with the same
%   size as R having true values for those positions that originate from an
%   input vector. This may be useful if any of the vectors contain NaNs.
%
%   Example:
%       a = 1:3 ; b = [] ; c = [1 NaN] ;
%       [M,tf] = padcat(a,b,c)
%       % find the original NaN
%       [Vev,Pos] = find(tf & isnan(M))
%       % -> Vec = 3 , Pos = 2
%
%   This second output can also be used to change the padding value into
%   something else than NaN.
%
%       [M, tf] = padcat(1:3,1,1:4)
%       M(~tf) = 99 % change the padding value into 99
%
%   Scalars will be concatenated into a single column vector.
%
%   See also CAT, RESHAPE, STRVCAT, CHAR, HORZCAT, VERTCAT, ISEMPTY
%            NONES, GROUP2CELL (Matlab File Exchange)

% for Matlab 2008 and up (tested in R2011a)
% version 1.2 (oct 2011)
% (c) Jos van der Geest
% email: jos@jasen.nl

% History
% 1.0 (feb 2009) created
% 1.1 (feb 2011) improved comments
% 1.2 (oct 2011) added help on changing the padding value into something
%     else than NaN

% Acknowledgements:
% Inspired by padadd.m (feb 2000) Fex ID 209 by Dave Johnson

error(nargchk(1,Inf,nargin)) ;

% check the inputs
SZ = cellfun(@size,varargin,'UniformOutput',false) ; % sizes
Ndim = cellfun(@ndims,varargin) ; %

if ~all(Ndim==2)
    error([mfilename ':WrongInputDimension'], ...
        'Input should be vectors.') ;
end

TF = [] ; % default second output so we do not have to check all the time

% for 2D matrices (including vectors) the size is a 1-by-2 vector
SZ = cat(1,SZ{:}) ;
maxSZ = max(SZ) ;    % probable size of the longest vector
% maxSZ equals :
%  - [1 1] for all scalars input
%  - [X 1] for column vectors
%  - [1 X] for all row vectors
%  - [X Y] otherwise (so padcat will not work!)

if ~any(maxSZ == 1),  % hmm, not all elements are 1-by-N or N-by-1
    % 2 options ...
    if any(maxSZ==0),
        % 1) all inputs are empty
        M  = [] ;
        return
    else
        % 2) wrong input
        % Either not all vectors have the same orientation (row and column
        % vectors are being mixed) or an input is a matrix.
        error([mfilename ':WrongInputSize'], ...
            'Inputs should be all row vectors or all column vectors.') ;
    end
end

if nargin == 1,
    % single input, nothing to concatenate ..
    M = varargin{1} ;
else
    % Concatenate row vectors in a row, and column vectors in a column.
    dim = (maxSZ(1)==1) + 1 ;      % Find out the dimension to work on
    X = cat(dim, varargin{:}) ;    % make one big list

    % we will use linear indexing, which operates along columns. We apply a
    % transpose at the end if the input were row vectors.

    if maxSZ(dim) == 1,
        % if all inputs are scalars, ...
        M = X ;   % copy the list
    elseif all(SZ(:,dim)==SZ(1,dim)),
        % all vectors have the same length
        M = reshape(X,SZ(1,dim),[]) ;% copy the list and reshape
    else
        % We do have vectors of different lengths.
        % Pre-allocate the final output array as a column oriented array. We
        % make it one larger to accommodate the largest vector as well.
        M = zeros([maxSZ(dim)+1 nargin]) ;
        % where do the fillers begin in each column
        M(sub2ind(size(M), SZ(:,dim).'+1, 1:nargin)) = 1 ;
        % Fillers should be put in after that position as well, so applying
        % cumsum on the columns
        % Note that we remove the last row; the largest vector will fill an
        % entire column.
        M = cumsum(M(1:end-1,:),1) ; % remove last row

        % If we need to return position of the non-fillers we will get them
        % now. We cannot do it afterwards, since NaNs may be present in the
        % inputs.
        if nargout>1,
            TF = ~M ;
            % and make use of this logical array
            M(~TF) = NaN ; % put the fillers in
            M(TF)  = X ;   % put the values in
        else
            M(M==1) = NaN ; % put the fillers in
            M(M==0) = X ;   % put the values in
        end
    end

    if dim == 2,
        % the inputs were row vectors, so transpose
        M = M.' ;
        TF = TF.' ; % was initialized as empty if not requested
    end
end % nargin == 1

if nargout > 1 && isempty(TF),
    % in this case, the inputs were all empty, all scalars, or all had the
    % same size.
    TF = true(size(M)) ;
end

author: Sable

repo: mcbench-benchmarks


text:
function rs = scalogram(s, scalemin, scalemax, scalestep, mlen)

%tstoolbox/@signal/scalogram
%   Syntax:
%     * rs = scalogram(s) => scalemin=0.1
%     * rs = scalogram(s, scalemin) => scalemax=1
%     * rs = scalogram(s, scalemin, scalemax) => scalestep=0.1
%     * rs = scalogram(s, scalemin, scalemax, scalestep) => mlen=10
%     * rs = scalogram(s, scalemin, scalemax, scalestep, mlen)
%
%   Scalogram of signal s using morlet wavelet. See also: spec2.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(1,5);

if nargin<2
	scalemin = 0.1;
end
if nargin<3
	scalemax = 1;
end
if nargin<4
	scalestep = 0.1;
end
if nargin<6
	mlen = 10;
end

c = scalogram(s.core, scalemin, scalemax, scalestep, mlen);
rs = signal(c, s);	% special constructor calling syntax for working routines
a = achse(unit, scalemin, scalestep);
rs = setaxis(rs, 2, a);
rs = setplothint(rs, 'spectrogram');
rs = addhistory(rs, ['Calculated scalogram']);
rs = addcommandlines(rs, 's = scalogram(s', scalemin, scalemax, scalestep);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
%{
    Copyright (C) Cambridge Electronic Design Limited 2014
    Author: James Thompson
    Web: www.ced.co.uk email: james@ced.co.uk, softhelp@ced.co.uk

    This file is part of CEDS64ML, a MATLAB interface to the SON64 library.

    CEDS64ML 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.

    CEDS64ML 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 CEDS64ML.  If not, see <http://www.gnu.org/licenses/>.
%}

function [ i64MaxTime ] = CEDS64MaxTime( fhand )
%CEDS64MAXTIME Returns the maximum time in ticks of all channels for a SON file
%   [ i64MaxTime ] = CEDS64MaxTime( fhand )
%   Inputs
%   fhand - Integer file handle
%   Outputs
%   i64MaxTime - the maximum time in ticks of all channels in this file,
%   otherwise a negative error code.

if (nargin ~= 1)
    i64MaxTime = -22;
    return;
end

i64MaxTime = calllib('ceds64int', 'S64MaxTime', fhand);
end

author: brainstorm-tools


text:
function determ = fiedler_determinant ( n, x )

%*****************************************************************************80
%
%% FIEDLER_DETERMINANT returns the determinant of the FIEDLER matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real X(N), the values that define A.
%
%    Output, real DETERM, the determinant.
%
  determ = 2.0 ^ ( n - 2 );

  if ( mod ( n, 2 ) == 1 )
    determ = - determ;
  end

  for i = 1 : n - 1
    for j = i + 1 : n
      if ( x(j) < x(i) )
        t    = x(j);
        x(j) = x(i);
        x(i) = t;
        determ = - determ;
      end
    end
  end

  determ = determ * ( x(n) - x(1) );

  for i = 2 : n
    determ = determ * ( x(i) - x(i-1) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = logmvn(X, Mu, Sigma, options)
%LOGMVN Evaluate the logarithms of normal density functions.
%The function evaluates the logarithm of Gaussian density function given
%data at X(n,:) with center at Mu(n,:) and covariance matrix Sigma(:,:,n).
%The code can handle different means and covariance matrices for each
%sample. It is much faster than the native function mvnpdf.
%
% Input:
%   X - N by D data. N is the number of samples. D is the number of
%       dimensions.
%   Mu - N by D (multiple) or 1 by D (single) mean.
%   Sigma - D by D by N (multiple) or D by D (single) covariance matrice.
%   Options - structure containing indices for transforming the Sigma to a
%       sparse matrix. Usually no need to set.
% Output:
%   y - N by 1 results that contain the logarithm of the Gaussians.
%
% Example:
%
% N = 1000;
% D = 2;
% Mu = [1 -1];
% Mu = repmat(Mu, 1000, 1) + randn(N,D);
% SN = repmat(reshape(abs(randn(N,1)), [1,1,N]), D, D) .* repmat(eye(D), [1,1,N]);
% Sigma = [.9 .4; .4 .3];
% Sigma = repmat(Sigma, [1,1,N]) + SN;
% X = mvnrnd(Mu,Sigma,1000);
%
% tic
% y = logmvn(X, Mu, Sigma);
% y = exp(y);
% toc
%
% tic
% y1 = mvnpdf(X, Mu, Sigma);
% toc
%
% norm(y-y1)
%
% Author: Yuan Zhou (zhouyuanzxcv@gmail.com)

% Copyright (C) 2015
if nargin < 4
    options = [];
end

if size(Mu,1) == 1 && ndims(Sigma) == 2 % single
    y = loggausspdf(X, Mu, Sigma);
elseif size(Mu,1) == size(X,1) && size(Sigma,3) == size(X,1) % multiple means
    y = logmvn_multiple(X, Mu, Sigma, options);
end

end

function y = logmvn_multiple(X, Mu, Sigma, options)
[N,B] = size(X);

% Transform the sigma matrix to a sparse block diagonal matrix
if ~isempty(options) && isfield(options,'sigma_update_in_logmvn')
    sparse_update_inplace(options.sigma_update_in_logmvn, Sigma(:));
    Sigma1 = options.sigma_update_in_logmvn;
else
    if ~isempty(options) && isfield(options,'Is') && isfield(options,'Js')
        Is = options.Is;
        Js = options.Js;
    else
        Is = (1:N*B);
        Is = repmat(reshape(Is, [B,1,N]), 1, B);
        Is = Is(:);

        Js = (1:N*B);
        Js = repmat(reshape(Js, [1,B,N]), B, 1);
        Js = Js(:);
    end

    Sigma1 = sparse(Is,Js,reshape(Sigma,N*B*B,1),N*B,N*B);
end

% Compute the Cholesky decomposition
[R,err] = chol(Sigma1);
if err ~= 0
    error('A Covariance matrix is not positive definite.');
end

% Final calculation
Y1 = X - Mu;

x1 = reshape(Y1', 1, N*B) / R;
x1 = sum(reshape(x1, B, N).^2, 1)';
x2 = sum(reshape(log(diag(R)),B,N), 1)';

y = -0.5 * x1 - x2 - B * log(2*pi) / 2;
end

function y = loggausspdf(X, mu, Sigma)
X = X';
mu = mu';
d = size(X,1);
X = bsxfun(@minus,X,mu);
[U,p]= cholcov(Sigma);
if p ~= 0
    error('ERROR: Sigma is not PD.');
end
Q = U'\X;
q = dot(Q,Q,1);  % quadratic term (M distance)
c = d*log(2*pi)+2*sum(log(diag(U)));   % normalization constant
y = -(c+q)/2;

y = y';
end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function pass = test_harmonic( )
% simple tests for cylindrical harmonic command.

tol = 3e2*chebfunpref().cheb2Prefs.chebfun2eps;

%example: build via jbessel w/known bessel root
 %L = 5, m = 3;
 F = diskfun(@(t,r) sqrt(2)/ (sqrt(pi)*abs(besselj(5+1, 15.7001740797116)))...
     *besselj(5, 15.7001740797116*r).*cos(5*t), 'polar');
 G = diskfun.harmonic(5,3);
 pass(1) = norm( F-G ) < tol;

 F = diskfun(@(t,r) sqrt(2)/ (sqrt(pi)*abs(besselj(5+1, 15.7001740797116)))...
     *besselj(5, 15.7001740797116*r).*sin(5*t), 'polar');
 G = diskfun.harmonic(-5,3);
 pass(2) = norm( F-G ) < tol;

%L = 0, m = 5;
 F = diskfun(@(t,r) sqrt(2)/ (sqrt(2*pi)*abs(besselj(0+1, 14.9309177084877)))...
     *besselj(0, 14.9309177084877*r), 'polar');
 G = diskfun.harmonic(0,5);
 pass(3) = norm(F-G) < tol;

%test orthonormality

A=diskfun.harmonic(29,10);
B=diskfun.harmonic(6,23);
C = diskfun.harmonic(-4, 3);
pass(4) = abs(sum2(A.*A)-1) < tol;
pass(5) = abs(sum2(B.*B)-1) < tol;
pass(6) = abs(sum2(A.*B)) < tol;
pass(7) = abs(sum2(A.*C)) < tol;
pass(8) = abs(sum2(C.*C)-1) < tol;

%neumann conditions, test orthonormal
A = diskfun.harmonic(10,8, 'neumann');
B = diskfun.harmonic(-4, 3, 'neumann');
pass(8) = abs(sum2(A.*A)-1) < tol;
pass(9) = abs(sum2(B.*B)-1) < tol;
pass(10) = abs(sum2(A.*B)) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef MockDatasetDiskCache < handle
    % MockDatasetDiskCache. Part of the PTK test framework
    %
    % This class is used in tests in place of a PTKDatasetDiskCache. It allows
    % expected calls to be verified, while maintaining some of the expected
    % behaviour of a PTKDatasetDiskCache object.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties (Access = private)
        StoredData
    end

    methods
        function obj = MockDatasetDiskCache()
            obj.StoredData = containers.Map;
        end

        function AddEntry(obj, name, value, context)
            obj.StoredData([name char(context)]) = value;
        end

        % Test methods

        function [value, cache_info] = LoadPluginResult(obj, plugin_name, context, reporting)
        end

        % Stores a plugin result in the disk cache and updates cached dependency
        % information
        function SavePluginResult(obj, plugin_name, result, cache_info, context, disk_cache_policy, memory_cache_policy, reporting)
        end

        % Saves additional data associated with this dataset to the cache
        function SaveData(obj, data_filename, value, reporting)
            obj.StoredData(data_filename) = value;
        end

        % Loads additional data associated with this dataset from the cache
        function value = LoadData(obj, data_filename, reporting)
        end

        function cache_path = GetCachePath(obj, ~)
        end

        function RemoveAllCachedFiles(obj, remove_framework_files, reporting)
        end

        function exists = Exists(obj, name, context, reporting)
            exists = obj.StoredData.isKey([name char(context)]);
        end

        function valid = CheckDependencyValid(obj, next_dependency, reporting)
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
%THRESHOLD  Applies a fixed-level threshold to each array element
%
%     dst = cv.threshold(src, thresh)
%     dst = cv.threshold(src, thresh, 'OptionName',optionValue, ...)
%     [dst, thresh] = cv.threshold(src, 'auto', ...)
%
% ## Input
% * __src__ Input array (single or multiple channel, 8-bit, 16-bit, or
%   floating-point).
% * __thresh__ Threshold value. Scalar numeric value or one of the strings:
%   * __Otsu__ use Otsu algorithm to choose the optimal threshold value
%   * __Triangle__ use Triangle algorithm to choose the optimal threshold value
%
% ## Output
% * __dst__ Output array of the same size and type as `src`.
% * __thresh__ the computed threshold value if Otsu's or Triangle methods are
%   used.
%
% ## Options
% * __MaxValue__ Maximum value to use with the 'Binary' and 'BinaryInv'
%   thresholding types. default 255
% * __Type__ Thresholding type, default 'Binary'. One of:
%   * __Binary__    `dst(x,y) = (src(x,y) > thresh) ? maxVal : 0`
%   * __BinaryInv__ `dst(x,y) = (src(x,y) > thresh) ? 0 : maxVal`
%   * __Trunc__     `dst(x,y) = (src(x,y) > thresh) ? thresh : src(x,y)`
%   * __ToZero__    `dst(x,y) = (src(x,y) > thresh) ? src(x,y) : 0`
%   * __ToZeroInv__ `dst(x,y) = (src(x,y) > thresh) ? 0 : src(x,y)`
%
% The function applies fixed-level thresholding to a multiple-channel array.
% The function is typically used to get a bi-level (binary) image out of a
% grayscale image (cv.compare could be also used for this purpose) or for
% removing a noise, that is, filtering out pixels with too small or too large
% values. There are several types of thresholding supported by the function.
% They are determined by `Type` parameter.
%
% When `thresh` is set 'Otsu' or 'Triangle', the function determines the
% optimal threshold value using the Otsu's or Triangle algorithm.
%
% Note: Currently, the Otsu's and Triangle methods are implemented only for
% 8-bit single-channel images.
%
% See also: cv.adaptiveThreshold, cv.findContours, cv.compare,
%  im2bw, graythresh, multithresh, imbinarize, otsuthresh, grayslice
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function pass = test_multipleOutputs(~)
% TEST_MULTIPLEOUTPUTS   Check that syntax [u,v,w] = N\0 works.

%% BVP scalar case
dom  = [0, 2];
N = chebop(@(x,u) diff(u,2) + 0.01*sin(u), dom);
N.lbc = 0;
N.rbc = 2;
u = N\0;
pass(1) = isa(u, 'chebfun');

%% BVP scalar case with info
dom  = [0, 2];
N = chebop(@(x,u) diff(u,2) + 0.01*sin(u), dom);
N.lbc = 0;
N.rbc = 2;
[u, info] = N\0;
pass(2) = isa(u, 'chebfun') && isstruct(info);

%% IVP scalar case
dom  = [0, 2];
N = chebop(@(x,u) diff(u,2) + sin(u), dom);
N.lbc = @(u) [u; diff(u)-1];
u = N\0;
pass(3) = isa(u, 'chebfun');

%% IVP scalar case with info
dom  = [0, 2];
N = chebop(@(x,u) diff(u,2) + sin(u), dom);
N.lbc = @(u) [u; diff(u)-1];
[u, info] = N\0;
pass(4) = isa(u, 'chebfun') && isstruct(info);

%% BVP two var case, one output
dom  = [0, 2];
N = chebop(@(x,u,v) [diff(u,2) + 0.1*sin(u) + v; diff(v,2) + 0.1*cos(v) + u], dom);
N.lbc = @(u,v) [u - 1 ; diff(v)];
N.rbc = @(u,v) [diff(v); v + 1];
u = N\0;
pass(5) = isa(u, 'chebmatrix');

%% BVP two var case, two outputs
dom  = [0, 2];
N = chebop(@(x,u,v) [diff(u,2) + 0.1*sin(u) + v; diff(v,2) + 0.1*cos(v) + u], dom);
N.lbc = @(u,v) [u - 1 ; diff(v)];
N.rbc = @(u,v) [diff(v); v + 1];
[u, v] = N\0;
pass(6) = isa(u, 'chebfun') && isa(v, 'chebfun');

%% BVP two var case with info
dom  = [0, 2];
N = chebop(@(x,u,v) [diff(u,2) + 0.1*sin(u) + v; diff(v,2) + 0.1*cos(v) + u], dom);
N.lbc = @(u,v) [u - 1 ; diff(v)];
N.rbc = @(u,v) [diff(v); v + 1];
[u, v, info] = N\0;
pass(7) = isa(u, 'chebfun') && isa(v, 'chebfun') && isstruct(info);

%% IVP two var case
dom  = [0, 2];
N = chebop(@(x,u,v) [diff(u,2) + 0.1*sin(u) - v; diff(v,2) + 0.1*cos(v) - u], dom);
N.lbc = @(u,v) [u - 1 ; v + 1; diff(u); diff(v)];
u = N\0;
pass(8) = isa(u, 'chebmatrix');

%% IVP two var case, two outputs
dom  = [0, 2];
N = chebop(@(x,u,v) [diff(u,2) + 0.1*sin(u) - v; diff(v,2) + 0.1*cos(v) - u], dom);
N.lbc = @(u,v) [u - 1 ; v + 1; diff(u); diff(v)];
[u, v] = N\0;
pass(9) = isa(u, 'chebfun') && isa(v, 'chebfun');

%% IVP two var case with info
dom  = [0, 2];
N = chebop(@(x,u,v) [diff(u,2) + 0.1*sin(u) - v; diff(v,2) + 0.1*cos(v) - u], dom);
N.lbc = @(u,v) [u - 1 ; v + 1; diff(u); diff(v)];
[u, v, info] = N\0;
pass(10) = isa(u, 'chebfun') && isa(v, 'chebfun') && isstruct(info);

%% BVP four var case
dom  = [0, 2];
N = chebop(@(x,u,v,w,y) [ ...
    diff(u) + 0.1*sin(u) + w; ...
    diff(v) + 0.1*sin(w) + y;
    diff(w) + 0.1*sin(y) - u;
    diff(y) + 0.1*sin(u) - v], dom);
N.lbc = @(u,v,w,y) [u - 1 ; v + 1];
N.rbc = @(u,v,w,y) [w + .5; y - .5];
u = N\0;
pass(11) = isa(u, 'chebmatrix');

%% BVP four var case, multiple outputs
dom  = [0, 2];
N = chebop(@(x,u,v,w,y) [ ...
    diff(u) + 0.1*sin(u) + w; ...
    diff(v) + 0.1*sin(w) + y;
    diff(w) + 0.1*sin(y) - u;
    diff(y) + 0.1*sin(u) - v], dom);
N.lbc = @(u,v,w,y) [u - 1 ; v + 1];
N.rbc = @(u,v,w,y) [w + .5; y - .5];
[u,v,w,y] = N\0;
pass(12) = isa(u, 'chebfun') && isa(v, 'chebfun') && ...
    isa(w, 'chebfun') && isa(y, 'chebfun');

%% BVP four var case with info
dom  = [0, 2];
N = chebop(@(x,u,v,w,y) [ ...
    diff(u) + 0.1*sin(u) + w; ...
    diff(v) + 0.1*sin(w) + y;
    diff(w) + 0.1*sin(y) - u;
    diff(y) + 0.1*sin(u) - v], dom);
N.lbc = @(u,v,w,y) [u - 1 ; v + 1];
N.rbc = @(u,v,w,y) [w + .5; y - .5];
[u,v,w,y,info] = N\0;
pass(13) = isa(u, 'chebfun') && isa(v, 'chebfun') && ...
    isa(w, 'chebfun') && isa(y, 'chebfun') && isstruct(info);

%% IVP four var case
dom  = [0, 2];
N = chebop(@(x,u,v,w,y) [ ...
    diff(u) + 0.1*sin(u) + w; ...
    diff(v) + 0.1*sin(w) + y;
    diff(w) + 0.1*sin(y) - u;
    diff(y) + 0.1*sin(u) - v], dom);
N.lbc = @(u,v,w,y) [u - 1 ; v + 1; w-.5; y+.5];
u = N\0;
pass(14) = isa(u, 'chebmatrix');

%% IVP four var case, multiple outputs
dom  = [0, 2];
N = chebop(@(x,u,v,w,y) [ ...
    diff(u) + 0.1*sin(u) + w; ...
    diff(v) + 0.1*sin(w) + y;
    diff(w) + 0.1*sin(y) - u;
    diff(y) + 0.1*sin(u) - v], dom);
N.lbc = @(u,v,w,y) [u - 1 ; v + 1; w-.5; y+.5];
[u,v,w,y] = N\0;
pass(15) = isa(u, 'chebfun') && isa(v, 'chebfun') && ...
    isa(w, 'chebfun') && isa(y, 'chebfun');

%% IVP four var case with info
dom  = [0, 2];
N = chebop(@(x,u,v,w,y) [ ...
    diff(u) + 0.1*sin(u) + w; ...
    diff(v) + 0.1*sin(w) + y;
    diff(w) + 0.1*sin(y) - u;
    diff(y) + 0.1*sin(u) - v], dom);
N.lbc = @(u,v,w,y) [u - 1 ; v + 1; w-.5; y+.5];
[u,v,w,y] = N\0;
pass(16) = isa(u, 'chebfun') && isa(v, 'chebfun') && ...
    isa(w, 'chebfun') && isa(y, 'chebfun') && isstruct(info);

end


text:
function sdpvarExpr = parseLMI(X)
%PARSELMI Internal function for dirty parsing of lmi string

% Assume Error

% TypeofConstraint = 1;

% Check for obsolete notation .<, .>, =
single_equality = strfind(X,'=');
if isempty(strfind(X,'>=')) & isempty(strfind(X,'<=')) & (rem(length(single_equality),2) == 1) % There is a single =
    error('Obsolete constraint =, use == in equalities.');
    %X = strrep(X,'=','==');
    %X = strrep(X,'====','=='); % Whoops, == replaced with =====!
end

if ~isempty(strfind(X,'.<'))
    disp(' ');
    disp('Warning: Obsolete constraint .<, use < in equalities.');
    disp('If you have a Hermitian matrix which you want to')
    disp('constrain to have negative values, use ,e.g.,  P(:)<0')
    disp('The constraint has been changed to <')
    disp(' ');
    X = strrep(X,'.<','<');
end

if ~isempty(strfind(X,'.>'))
    disp(' ');
    disp('Warning: Obsolete constraint .>, use > in equalities.');
    disp('If you have a Hermitian matrix which you want to')
    disp('constrain to have negative values, use e.g.  P(:)>0')
    disp('The constraint has been changed to >')
    disp(' ');
    X = strrep(X,'.>','>');
end

% Any norm? If not, we're done!
if isempty(strfind(X,'||'))
    sdpvarExpr = X;
    return
end

% The only actual parsing needed is for the ||Axplusb||<cx+d notation
delimiters = {'.<','.>','<.','>.','<','>','==',''};
delindex = 0;indtodel = [];
while isempty(indtodel) & (delindex<=7)
    delindex = delindex+1;
    indtodel = strfind(X,delimiters{delindex});
end

switch delindex
    case 5 %<
        TypeofConstraint = 1;
        ind_start =  indtodel-1;
        ind_end   =  indtodel+1;
        REVERSE   = 1;
    case 6 %>
        TypeofConstraint = 1;
        ind_start =  indtodel-1;
        ind_end   =  indtodel+1;
        REVERSE   = 0;
    case 7%=
        TypeofConstraint = 3;
        ind_start = indtodel-1;
        ind_end   = indtodel+2;
        REVERSE   = 0;
    case {3,4}
        error('Incorrect argument. Perhaps you mean .> or .<');
    otherwise
        error('Incorrect argument. Could not find <=>.>.<')
end

LeftHand   = X(1:ind_start);
RightHand  = X(ind_end:end);

if REVERSE
    temp = LeftHand;
    LeftHand = RightHand;
    RightHand = temp;
end

% Search for a norm expression
ind_norm_Right = strfind(RightHand,'||');
ind_norm_Left  = strfind(LeftHand,'||');

% Any norm at all, if not, we're done!
if isempty(ind_norm_Right) & isempty(ind_norm_Left)
    sdpvarExpr = [LeftHand '-(' RightHand ')'];
    return
end

% Equality constrained norm?
if (TypeofConstraint == 3) & (ind_norm_Right | ind_norm_Left)
    error('Equality constraints cannot be used with ||..|| operator');
end

% Convex normconstraint?
if ~isempty(ind_norm_Left)
    error('Norm ||..|| must look like || column vector ||< scalar, or scalar>|| vector ||')
end

% Everthing seem ok in norm
TypeofConstraint = 4;
Axplusb = RightHand(2+ind_norm_Right(1):ind_norm_Right(2)-1);
WithoutNorm = strrep(RightHand,RightHand(ind_norm_Right(1):ind_norm_Right(2)+1),'0');
if length(WithoutNorm)==1
    cxplusd = LeftHand;
else
    cxplusd = [LeftHand '-(' WithoutNorm ')'];
end

sdpvarExpr = ['cone( ' Axplusb ',' cxplusd ')'];

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function r8col_sorted_tol_undex_test ( )

%*****************************************************************************80
%
%% R8COL_SORTED_TOL_UNDEX_TEST tests R8COL_SORTED_TOL_UNDEX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 July 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8COL_SORTED_TOL_UNDEX_TEST\n' );
  fprintf ( 1, '  R8COL_SORTED_TOL_UNDEX produces index vectors which \n' );
  fprintf ( 1, '  create a sorted list of the tolerably unique columns of \n' );
  fprintf ( 1, '  a sorted R8COL and a map from the original R8COL to the \n' );
  fprintf ( 1, '  (implicit) R8COL of sorted tolerably unique elements.\n' );

  m = 3;
  n = 22;

  a = [ ...
    1.9,  0.0, 10.0; ...
    2.0,  6.0, 10.0; ...
    4.0,  8.0, 12.0; ...
    1.0,  5.0,  9.0; ...
    3.0,  7.0, 11.0; ...
    2.0,  6.0,  0.0; ...
    2.0,  0.0, 10.1; ...
    2.0,  0.1, 10.0; ...
    3.0,  4.0, 18.0; ...
    1.9,  8.0, 10.0; ...
    0.0,  0.0,  0.0; ...
    0.0,  6.0, 10.0; ...
    2.1,  0.0, 10.0; ...
    2.0,  6.0, 10.0; ...
    3.0,  7.0, 11.0; ...
    2.0,  0.0, 10.0; ...
    2.0,  0.0, 10.0; ...
    2.0,  6.0, 10.0; ...
    1.0,  5.0,  9.0; ...
    2.0,  0.0, 10.1; ...
    1.0,  5.0,  9.1; ...
    1.0,  5.1,  9.0 ]';

  r8mat_transpose_print ( m, n, a, '  The unsorted R8COL (transposed):' );

  a = r8col_sort_heap_a ( m, n, a );

  r8mat_transpose_print ( m, n, a, '  The sorted R8COL (transposed):' );

  tol = 0.25;
  n_unique = r8col_sorted_tol_unique_count ( m, n, a, tol );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Tolerance for equality is %e\n', tol );
  fprintf ( 1, '  Number of unique entries in X is %d\n', n_unique );

  [ undx, xdnu ] = r8col_sorted_tol_undex ( m, n, a, n_unique, tol );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  XDNU points to the representative for each item.\n' );
  fprintf ( 1, '  UNDX selects the representatives..\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n_unique
    fprintf ( 1, '  %4d  %4d  %4d\n', i, xdnu(i), undx(i) );
  end
  for i = n_unique + 1 : n
    fprintf ( 1, '  %4d  %4d\n', i, xdnu(i) );
  end

  for j = 1 : n_unique
    au(1:m,j) = a(1:m,undx(j));
  end

  r8mat_transpose_print ( m, n_unique, au, ...
    '  The tolerably unique R8COL (transposed):' );

  return
end


text:
function [ image ] = DeleteWindowsImage(image,binaryImage)
%UNTITLED3 Summary of this function goes here
%   Detailed explanation goes here
discHeight = size(binaryImage,1);
discWidth = size(binaryImage,2);
Height = size(image,1);
s = Height/discHeight;

for i=1:discHeight
    for j=1:discWidth
        if (binaryImage(i,j)==1)
            image((i-1)*s+1:i*s,(j-1)*s+1:j*s,:) = 0;
        end
    end
end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function loadlayerprofile(varargin)
%*************************************************
%DATE: 10/8/2007 (created); 	10/8/2007 (modified)
%FUNCTION: 	load layer profiles
%INPUTS: 	NULL
%OUTPUTS:	NULL
%muwn.gu@gmail.com
%*************************************************
global Layerprofile
Layerprofile.color.silicon =  	[53, 53, 53]/255;
Layerprofile.color.sidioxi=  	[85, 85, 85]/255;
Layerprofile.color.metal =  	[146, 146, 146]/255;
Layerprofile.color.sini=  		[190, 190, 190]/255;
Layerprofile.color.posi=  		[220, 220, 200]/255;
Layerprofile.color.phore =  	[0.1, 0.1, 0.4];

Layerprofile.origin.x = 20;
Layerprofile.origin.y = 10;

Layerprofile.stepSpace = 6;

Layerprofile.currentY = Layerprofile.origin.y;

Layerprofile.diameter = 54;

Layerprofile.lineWidth = 0.3;

global layerMatrix;
global layerColor;


text:
function [xn] = normalize(x_kk,fc,cc,kc,alpha_c)

%normalize
%
%[xn] = normalize(x_kk,fc,cc,kc,alpha_c)
%
%Computes the normalized coordinates xn given the pixel coordinates x_kk
%and the intrinsic camera parameters fc, cc and kc.
%
%INPUT: x_kk: Feature locations on the images
%       fc: Camera focal length
%       cc: Principal point coordinates
%       kc: Distortion coefficients
%       alpha_c: Skew coefficient
%
%OUTPUT: xn: Normalized feature locations on the image plane (a 2XN matrix)
%
%Important functions called within that program:
%
%comp_distortion_oulu: undistort pixel coordinates.

if nargin < 5,
   alpha_c = 0;
   if nargin < 4;
      kc = [0;0;0;0;0];
      if nargin < 3;
         cc = [0;0];
         if nargin < 2,
            fc = [1;1];
         end;
      end;
   end;
end;

% First: Subtract principal point, and divide by the focal length:
x_distort = [(x_kk(1,:) - cc(1))/fc(1);(x_kk(2,:) - cc(2))/fc(2)];

% Second: undo skew
x_distort(1,:) = x_distort(1,:) - alpha_c * x_distort(2,:);

if norm(kc) ~= 0,
	% Third: Compensate for lens distortion:
	% xn = comp_distortion_oulu(x_distort,kc);
  xn = undistort(x_distort, kc);
else
   xn = x_distort;
end;

author: jianxiongxiao

repo: ProfXkit


text:
function out = isempty(f)
% overloads isempty

out = isempty(f.h);

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function Y = imag2real(Y)
%IMAG2REAL

if isreal(Y)
    y = Y;
    return
end

realBase = real(Y.basis);
imagBase = imag(Y.basis);

lmi_variables = getvariables(Y);
nv = length(lmi_variables);

% [re im;-im re] = kron(I,re) + kron([0 1;-1 0],im)
sparse_X1 = [1 0;0 1];
sparse_X2 = [0 1;-1 0];
temp = kron(sparse_X1,reshape(realBase(:,1),Y.dim(1),Y.dim(1)))+ kron(sparse_X2,reshape(imagBase(:,1),Y.dim(1),Y.dim(1)));
temp = temp(:);

Y.basis = temp(:);
for i = 1:nv
    temp1 = kron(sparse_X1,reshape(realBase(:,i+1),Y.dim(1),Y.dim(1)));
    temp2 = kron(sparse_X2,reshape(imagBase(:,i+1),Y.dim(1),Y.dim(1)));
    Y.basis(:,i+1) = temp1(:) + temp2(:);
end;

Y.dim(1) = size(temp1,1);
Y.dim(2) = size(temp1,2);
Y = clean(Y);
% Reset info about conic terms
if isa(Y,'sdpvar')
    Y.conicinfo = [0 0];
end


text:
%  meandpth  finds the average depth for a predefined running
%  window in terms of number of events, and a selected step
%  and plots the results.
%                                                     R.Z. 6/94
%                          Operates on newcat
%
report_this_filefun(mfilename('fullpath'));

% Find out of figure already exists
%
[existFlag,figNumber]=figure_exists('Mean Depth2',1);
newDep2WindowFlag=~existFlag;

% Set up the Seismicity Map window Enviroment
%
if newDep2WindowFlag

    figure_w_normalized_uicontrolunits(...
        'Name','Mean Depth2',...
        'visible','off',...
        'NumberTitle','off', ...
        'MenuBar','none', ...
        'Color',[ 1 1 1], ...
        'NextPlot','new', ...
        'Units','Pixel',  'Position',[wex wey 550 400'])
    depfg2 = gcf;
    hold on
    axis off
    matdraw

    uicontrol('Style','Pushbutton',...
        'Position',[.9 .80 .10 .05],...
        'Units','normalized',...
        'Callback','sta = ''ast'';medispas1','String','AS');
    uicontrol('Style','Pushbutton',...
        'Position',[.9 .70 .10 .05],...
        'Units','normalized',...
        'Callback','sta = ''lta'';medispas1','String','LTA');
    uicontrol('Style','Pushbutton',...
        'Position',[.9 .90 .10 .05],...
        'Units','normalized',...
        'Callback','dispma4','String','Com');

    new = uicontrol('style','edit','value',iwl,...
        'string',num2str(iwl), 'background','y',...
        'Callback','iwl=str2num(get(new,''String''));''String'',num2str(iwl);medispas1',...
        'units','norm','pos',[.90 .30 .10 .06],'min',0.1,'max',100);

    newlabel = uicontrol('style','text','units','norm','pos',[.85 .30 .05 .06]);
    set(newlabel,'string','iwl:','background',color_fbg);

    uicontrol('Units','normal',...
        'Position',[.90 .25 .08 .06],'String','Go',...
         'Callback','medispas1')

end  % if figure exist

figure_w_normalized_uicontrolunits(depfg2)
delete(gca);delete(gca);delete(gca);
set(gca,'visible','off');

%orient tall
set(gcf,'Units','centimeter','PaperPosition',[1 1 8 7])
rect = [0.15, 0.15, 0.65, 0.30];
axes('position',rect)
p5 = gca;

% plot errbar
%errorbar(xt2,meand,er)
%plot(xt2,meand,'co')
%for i = 1:length(xt2)
%boxutil(-me(:,i),1,xt2(i),0.5, 'r.',1,1.5);
%hold on
%end

pl = plot(xt2,meand,'-k')
set(pl,'LineWidth',1.0)
hold on
pl = plot(xt2,meand,'ok')
set(pl,'LineWidth',1.0,'MarkerSize',6)
if isempty(maepi) == 0
    pl =   plot(maepi(:,3),-maepi(:,7),'xk');
    set(pl,'LineWidth',2.0)
end

axis([min(xt2) max(xt2+0.5) min([meand*1.1 ])  max([meand*0.9 ])])
v = axis;
grid
xlabel('Time (years)','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m,'Color','k')
ylabel('Mean Depth (km)','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m,'Color','k')
stro = [' ' file1 '; ' num2str(iwln) ' / ' num2str(step)];
title2(stro,'FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m,'Color','k')

set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)

hold off

set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)

set(gca,'visible','on');
set(gcf,'visible','on');

ic = 1;

author: CelsoReyes


text:
close all; clear all; clc;
yf = spx.data.image.EhsanYaleFaces();
[faces, sizes, labels] = yf.all_faces();
% normalize the faces
faces_normalized = spx.norm.normalize_l2(faces);
[M, S]  = size(faces_normalized);

% only between 2 subjects
% interesting cases are (10, 11), (28, 29), (29, 30)
% (35, 36), (36, 37)
s1 = 35;
s2 = s1+1;
columns_s1 = (1:64) + (s1-1)*64;
columns_s2 = (1:64) + (s2-1)*64;
all_columns = [columns_s1 columns_s2];
fprintf('\n\n\n Statistics for the first pair of subjects:\n\n');
faces_normalized = faces_normalized(:, all_columns);
[M, S]  = size(faces_normalized);
sizes = [64 64];
angle_result = spx.cluster.subspace.nearest_same_subspace_neighbors_by_inner_product(faces_normalized, sizes);

fprintf('Within Neighbor Counts:\n %s\n', spx.stats.format_descriptive_statistics(angle_result.within_neighbor_counts));
fprintf('Nearest Within Neighbor Indices:\n %s\n', spx.stats.format_descriptive_statistics(angle_result.nearst_within_neighbor_indices));
nearst_within_neighbor_indices = angle_result.nearst_within_neighbor_indices;
% club all those cases where nearest within neighbor is far away.
nearst_within_neighbor_indices(nearst_within_neighbor_indices > 5) = -1;
tabulate(nearst_within_neighbor_indices);
fprintf('First In Out Angle Spreads:\n %s\n', spx.stats.format_descriptive_statistics(angle_result.first_in_out_angle_spreads));
fioas = angle_result.first_in_out_angle_spreads;
%fioas = round(fioas * 10) / 10;
fioas = round(fioas/4)*4;
fioas(fioas > 8) = 100;
fioas(fioas < -8) = 100;
angle_gt_zero_flags = angle_result.first_in_out_angle_spreads >= 0;
angle_gte_zero = angle_result.first_in_out_angle_spreads(angle_gt_zero_flags);
angle_lte_two_flags =  angle_gte_zero <= 2;
angle_lte_four_flags =  angle_gte_zero <= 4;
angle_lte_eight_flags =  angle_gte_zero <= 8;
angle_lt_zero_flags = angle_result.first_in_out_angle_spreads < 0;
angle_lt_zero =  angle_result.first_in_out_angle_spreads(angle_lt_zero_flags);
angle_gte_minus_two_flags =  angle_lt_zero >= -2;
angle_gte_minus_four_flags =  angle_lt_zero >= -4;
abs_angle_lte_four_flags = abs(angle_result.first_in_out_angle_spreads) <=4;
fprintf('angle greater than equal to 0: %.2f %%\n', sum(angle_gt_zero_flags) * 100 / S);
fprintf('angle less than 2: %.2f %%\n', sum(angle_lte_two_flags) * 100 / S);
fprintf('angle less than 4: %.2f %%\n', sum(angle_lte_four_flags) * 100 / S);
fprintf('angle less than 8: %.2f %%\n', sum(angle_lte_four_flags) * 100 / S);
fprintf('angle less than 0: %.2f %%\n', sum(angle_lt_zero_flags) * 100 / S);
fprintf('angle greater than -2: %.2f %%\n', sum(angle_gte_minus_two_flags) * 100 / S);
fprintf('angle greater than -4: %.2f %%\n', sum(angle_gte_minus_four_flags) * 100 / S);
fprintf('abs angle less than 4: %.2f %%\n', sum(abs_angle_lte_four_flags) * 100 / S);

author: indigits

repo: sparse-plex


text:
function varargout = openNEV(varargin)

% openNEV
%
% Opens an .nev file for reading, returns all file information in a NEV
% structure. Works with File Spec 2.1 & 2.2 & 2.3.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use OUTPUT = openNEV(fname, 'noread', 'report', 'noparse', 'nowarning',
%                             'nosave', 'nomat', 'uV', 'overwrite').
%
% NOTE: All input arguments are optional. Input arguments may be in any order.
%
%   fname:        Name of the file to be opened. If the fname is omitted
%                 the user will be prompted to select a file using an open
%                 file user interface.
%                 DEFAULT: Will open Open File UI.
%
%   'noread':     Will not read the spike waveforms if user passes this argument.
%                 DEFAULT: will read spike waveform.
%
%   'report':     Will show a summary report if user passes this argument.
%                 DEFAULT: will not show report.
%
%   'parse':    The code will not parse the experimental parameters in digital I/O.
%                 See below for guidelines on how to format your parameters.
%                 DEFAULT: will not parse the parameters.
%
%   'nowarning':  The code will not give a warning if there is an error in
%                 parsing.
%                 DEFAULT: will give warning message.
%
%   'nosave':     The code will not save a copy of the NEV structure as a
%                 MAT file. By default the code will save a copy in the same
%                 folder as the NEV file for easy future access.
%                 DEFAULT: will save the MAT file.
%
%   'nomat':      Will not look for a MAT file. This option will force
%                 openNEV to open a NEV file instead of any available MAT
%                 files.
%                 DEFAULT: will load the MAT file if available.
%
%   'uV':         Will read the spike waveforms in unit of uV instead of
%                 raw values. Note that this conversion may lead to loss of
%                 information (e.g. 15/4 = 4) since the waveforms type will
%                 stay in int16. It's recommended to read raw spike
%                 waveforms and then perform the conversion at a later
%                 time.
%                 DEFAULT: will read waveform information in raw.
%
%   '8bits':      Indicates that 8 bits on the digital IO port was used
%                 instead of 16 bits.
%                 DEFAULT: will assumes that 16 bits of digital IO were used.
%
%   't:':         Indicats the time window of the NEV file to be read. For
%                 example, if t: is set to 2 (i.e. 't:0.6')
%                 then only the first 2 seconds of the file is to be read. If set
%                 to 2-50 (i.e. 't:2:50) then the time between 2 seconds
%                 and 50 seconds will be read.
%                 DEFAULT: the entire file will be read if 't:xx:xx' is not
%                 passed to the function.
%
%   'overwrite':  If MATLAB loads a NEV file using 'nomat' and a MAT file
%                 already exists, by default it will prompt the user to
%                 allow for overwriting the old MAT. Passing the
%                 'overwrite' flag will automatically overwrite the newly
%                 opened NEV file ont the old MAT file.
%                 DEFAULT: will ask the user whether to overwrite the old
%                 MAT.
%
%   OUTPUT:       Contains the NEV structure.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   USAGE EXAMPLE:
%
%   openNEV('report','read');
%
%   In the example above, the file dialogue will prompt for a file. A
%   report of the file contents will be shown. The digital data will not be
%   parsed. The data needs to be in the proper format (refer below). The
%   spike waveforms are in raw units and not in uV.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   DIGITAL PARAMETERS/MARKERS FORMAT:
%
%   In order for this function to parse your experimental parameters they
%   need to be in the following format:
%
%   *ParamLabel:Parameter1=value1;Parameter2=value2;Parameter3=value3;#
%
%   TWO EXAMPLES:
%   *ExpParameter:Intensity=1.02;Duration=400;Trials=1;PageSegment=14;#
%
%   *Stimulation:StimCount=5;Duration=10;#
%
%   In the first example, the parameter is of type "ExpParameter". The
%   parameters are, "Intensity, Duration, Trials, and PageSement." The
%   values of those parameters are, "1.02, 400, 1, and 14," respectively.
%   The second example is of type "Stimulation". The name of the parameters
%   are "StimCount" and "Duration" and the values are "5" and "10"
%   respectively.
%   -----------------------------------------------------------------------
%   It can also read single value markers that follow the following format.
%
%   *MarkerName=Value;#
%
%   EXAMPLES:  *WaitSeconds=10;# OR  *JuiceStatus=ON;#
%
%   The above line is a "Marker". The marker value is 10 in the first
%   and it's ON in the second example.
%   -----------------------------------------------------------------------
%   Moreover, the marker could be a single value:
%
%   *MarkerValue#
%
%   EXAMPLES: *JuiceOff#  OR  *HandsOnSwitches#
%   -----------------------------------------------------------------------
%   The label, parameter name, and values are flexible and can be anything.
%   The only required formatting is that the user needs to have a label
%   followed by a colon ':', followed by a field name 'MarkerVal', followed
%   by an equal sign '=', followed by the parameter value '10', and end
%   with a semi-colon ';'.
%
%   NOTE:
%   Every line requires a pound-sign '#' at the very end.
%   Every line requires a star sign '*' at the very beginning. If you
%   use LabVIEW SendtoCerebus.vi by Kian Torab then there is no need for
%   a '*' in the beginning.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Kian Torab
%   kian@blackrockmicro.com
%   Blackrock Microsystems
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Version History
%
% 4.4.0.0:
%   - Major performance boost in reading NEV files when tracking data is
%   stored in the file.
%
% 4.4.0.2:
%   - Updated documentation.
%
% 4.4.0.3: 5 January 2014
%   - Fixed the way DayOfWeek is read in MetaTags.
%   - Fixed 'noread' argument, so when passed, openNEV will not read the
%   spike waveforms.
%
% 4.4.1.0: 25 January 2014
%   - Fixed a bug that resulted from passing 'read' to openNEV.
%
% 4.4.2.0: 28 February 2014
%   - Fixed bug related to loading data with t:XX:XX argument.
%
% 4.4.3.0: 12 June 2014
%   - Fixed a typo in the help.
%
% 4.4.3.1: 13 June 2014
%   - Updated the version numbers in the help and in the function itself.
%
% 4.4.3.2: 21 June 2014
%   - Fixed a bug where Application name wasn't being read properly.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Defining structures
NEV = struct('MetaTags',[], 'ElectrodesInfo', [], 'Data', []);
NEV.MetaTags = struct('Subject', [], 'Experimenter', [], 'DateTime', [],...
    'SampleRes',[],'Comment',[],'FileTypeID',[],'Flags',[], 'openNEVver', [], ...
    'DateTimeRaw', [], 'FileSpec', [], 'PacketBytes', [], 'HeaderOffset', [], ...
    'DataDuration', [], 'DataDurationSec', [], 'PacketCount', [], ...
    'TimeRes', [], 'Application', [], 'Filename', [], 'FilePath', []);
NEV.Data = struct('SerialDigitalIO', [], 'Spikes', [], 'Comments', [], 'VideoSync', [], ...
    'Tracking', [], 'TrackingEvents', [], 'PatientTrigger', [], 'Reconfig', []);
NEV.Data.Spikes = struct('TimeStamp', [],'Electrode', [],...
    'Unit', [],'Waveform', [], 'WaveformUnit', []);
NEV.Data.SerialDigitalIO = struct('InputType', [], 'TimeStamp', [],...
    'TimeStampSec', [], 'Type', [], 'Value', [], 'InsertionReason', [], 'UnparsedData', []);
NEV.Data.VideoSync = struct('TimeStamp', [], 'FileNumber', [], 'FrameNumber', [], 'ElapsedTime', [], 'SourceID', []);
NEV.Data.Comments = struct('TimeStamp', [], 'TimeStampSec', [], 'CharSet', [], 'Color', [], 'Text', []);
NEV.Data.Tracking = [];
NEV.Data.TrackingEvents = struct('TimeStamp', [], 'TimeStampSec', [], 'Text', []);
NEV.Data.PatientTrigger = struct('TimeStamp', [], 'TriggerType', []);
NEV.Data.Reconfig = struct('TimeStamp', [], 'ChangeType', [], 'CompName', [], 'ConfigChanged', []);
Flags = struct;

NEV.MetaTags.openNEVver = '4.4.3.1';

%% Check for multiple versions of openNEV in path
if size(which('openNEV', '-ALL'),1) > 1
    disp('WARNING: There are multiple openNEV functions in the path. Use which openNEV -ALL for more information.');
end

%% Validating input arguments
for i=1:length(varargin)
    switch lower(varargin{i})
        case 'report'
            Flags.Report = varargin{i};
        case 'noread'
            Flags.ReadData = varargin{i};
        case 'read'
            Flags.ReadData = varargin{i};
        case 'nosave'
            Flags.SaveFile = varargin{i};
        case 'nomat'
            Flags.NoMAT = varargin{i};
        case 'nowarning'
            Flags.WarningStat = varargin{i};
        case 'parse'
            Flags.ParseData = 'parse';
        case 'uv'
            Flags.waveformUnits = 'uV';
        case '8bits'
            Flags.digIOBits = '8bits';
        case 'overwrite'
            Flags.Overwrite = 'overwrite';
        otherwise
            temp = varargin{i};
            if length(temp)>3 && ...
                    (strcmpi(temp(3),'\') || ...
                     strcmpi(temp(1),'/') || ...
                     strcmpi(temp(2),'/'))
                fileFullPath = varargin{i};
                if exist(fileFullPath, 'file') ~= 2
                    disp('The file does not exist.');
                    varargout{1} = [];
                    return;
                end
            elseif length(temp)>3 && strcmpi(temp(1:2),'t:') && ~strcmpi(temp(3), '\') && ~strcmpi(temp(3), '/')
                temp(1:2) = [];
                temp = str2num(temp);
                if length(temp) == 1
                    fprintf('Only one timepoint (%0.0f) was passed to the function.\n', temp);
                    fprintf('The initial timepoint is set to 0, so data between 0 and %0.0f will be read.\n', temp);
                    temp(2) = temp;
                    temp(1) = 0;
                end
                readTime = [temp(1), temp(end)];
                Flags.SaveFile = 'nosave';
                Flags.NoMAT = 'nomat';
            else
                if ~isnumeric(varargin{i})
                    disp(['Invalid argument ''' varargin{i} ''' .']);
                else
                    disp(['Invalid argument ''' num2str(varargin{i}) ''' .']);
                end
                clear variables;
                if nargout
                    varargout{1} = [];
                end
                return;
            end
            clear temp;
    end
end; clear i;

%% Defining and validating variables
if ~exist('fileFullPath', 'var')
    if exist('getFile.m', 'file') == 2
        [fileName pathName] = getFile('*.nev', 'Choose a NEV file...');
    else
        [fileName pathName] = uigetfile;
    end
    fileFullPath = [pathName fileName];
    if fileFullPath==0;
        clear variables;
        if nargout
            varargout{1} = [];
        end
        disp('No file was selected.');
        return
    end
end

if ~isfield(Flags, 'Report');        Flags.Report = 'noreport'; end
if ~isfield(Flags, 'WarningStat');   Flags.WarningStat = 'warning'; end;
if ~isfield(Flags, 'ReadData');      Flags.ReadData = 'read'; end
if ~isfield(Flags, 'ParseData');     Flags.ParseData = 'noparse'; end
if ~isfield(Flags, 'SaveFile');      Flags.SaveFile = 'save'; end;
if ~isfield(Flags, 'NoMAT');         Flags.NoMAT = 'yesmat'; end;
if ~isfield(Flags, 'waveformUnits'); Flags.waveformUnits = 'raw'; end;
if ~isfield(Flags, 'digIOBits');     Flags.digIOBits = '16bits'; end;
if ~isfield(Flags, 'Overwrite');     Flags.Overwrite = 'noOverwrite'; end;
if strcmpi(Flags.Report, 'report')
    disp(['openNEV ' NEV.MetaTags.openNEVver]);
end
%%  Validating existance of parseCommand
if strcmpi(Flags.ParseData, 'parse')
    if exist('parseCommand.m', 'file') ~= 2
        disp('This version of openNEV requires function parseCommand.m to be placed in path.');
        clear variables;
        if nargout
            varargout{1} = [];
        end
        return;
    end
end

tic;
matPath = [fileFullPath(1:end-4) '.mat'];

%% Check for a MAT file and load that instead of NEV
if exist(matPath, 'file') == 2 && strcmpi(Flags.NoMAT, 'yesmat')
    disp('MAT file corresponding to selected NEV file already exists. Loading MAT instead...');
    load(matPath);
    if isempty(NEV.Data.Spikes.Waveform) && strcmpi(Flags.ReadData, 'read')
        disp('The MAT file does not contain waveforms. Loading NEV instead...');
    else
        if ~nargout
            assignin('base', 'NEV', NEV);
            clear variables;
        else
            varargout{1} = NEV;
        end
        return;
    end
end

%% Reading BasicHeader information from file
FID                       = fopen(fileFullPath, 'r', 'ieee-le');
BasicHeader               = fread(FID, 336, '*uint8');
NEV.MetaTags.FileTypeID   = char(BasicHeader(1:8)');
NEV.MetaTags.FileSpec     = [num2str(double(BasicHeader(9))) '.' num2str(double(BasicHeader(10)))];
NEV.MetaTags.Flags        = dec2bin(double(typecast(BasicHeader(11:12), 'uint16')),16);
Trackers.fExtendedHeader  = double(typecast(BasicHeader(13:16), 'uint32'));
NEV.MetaTags.HeaderOffset = Trackers.fExtendedHeader;
Trackers.countPacketBytes = double(typecast(BasicHeader(17:20), 'uint32'));
NEV.MetaTags.PacketBytes  = Trackers.countPacketBytes;
NEV.MetaTags.TimeRes      = double(typecast(BasicHeader(21:24), 'uint32'));
NEV.MetaTags.SampleRes    = typecast(BasicHeader(25:28), 'uint32');
t                         = double(typecast(BasicHeader(29:44), 'uint16'));
tempApp                   = BasicHeader(45:76)';
tempApp(find(tempApp == 0):end) = [];
NEV.MetaTags.Application  = char(tempApp); clear tempApp;
NEV.MetaTags.Comment      = char(BasicHeader(77:332)');
[NEV.MetaTags.FilePath, NEV.MetaTags.Filename] = fileparts(fileFullPath);
Trackers.countExtHeader   = typecast(BasicHeader(333:336), 'uint32');
clear BasicHeader;

if strcmpi(NEV.MetaTags.FileTypeID, 'NEURALEV')
    if exist([fileFullPath(1:end-8) '.sif'], 'file') == 2
        METATAGS = textread([fileFullPath(1:end-8) '.sif'], '%s');
        NEV.MetaTags.Subject      = METATAGS{3}(5:end-5);
        NEV.MetaTags.Experimenter = [METATAGS{5}(8:end-8) ' ' METATAGS{6}(7:end-7)];
    end
end
if ~any(strcmpi(NEV.MetaTags.FileSpec, {'2.1', '2.2', '2.3'}))
    disp('Unknown filespec. Cannot open file.');
    fclose FID;
    clear variables;
    if nargout
        varargout{1} = [];
    end
    return;
end
clear fileFullPath;

%% Parsing and validating FileSpec and DateTime variables
NEV.MetaTags.DateTimeRaw = t.';
NEV.MetaTags.DateTime = [num2str(t(2)) '/'  num2str(t(4)+2) '/' num2str(t(1))...
    ' ' datestr(t(3), 'dddd') ' ' num2str(t(5)) ':'  ...
    num2str(t(6)) ':'  num2str(t(7)) '.' num2str(t(8))] ;
clear t;

%% Removing extra garbage characters from the Comment field.
NEV.MetaTags.Comment(find(NEV.MetaTags.Comment==0,1):end) = 0;

%% Recording after BasicHeader file position
Trackers.fBasicHeader = ftell(FID); %#ok<NASGU>

% Calculating the length of the data
currentLocation = ftell(FID);
fseek(FID, -Trackers.countPacketBytes, 'eof');
NEV.MetaTags.DataDuration = fread(FID, 1, 'uint32=>double');
NEV.MetaTags.DataDurationSec = double(NEV.MetaTags.DataDuration) / double(NEV.MetaTags.SampleRes);
fseek(FID, currentLocation, 'bof');

%% Reading ExtendedHeader information
for ii=1:Trackers.countExtHeader
    ExtendedHeader = fread(FID, 32, '*uint8');
    PacketID = char(ExtendedHeader(1:8)');
    switch PacketID
        case 'ARRAYNME'
            NEV.ArrayInfo.ElectrodeName    = char(ExtendedHeader(9:end));
        case 'ECOMMENT'
            NEV.ArrayInfo.ArrayComment     = char(ExtendedHeader(9:end));
        case 'CCOMMENT'
            NEV.ArrayInfo.ArrayCommentCont = char(ExtendedHeader(9:end));
        case 'MAPFILE'
            NEV.ArrayInfo.MapFile          = char(ExtendedHeader(9:end));
        case 'NEUEVWAV'
            ElectrodeID                       = typecast(ExtendedHeader(9:10), 'uint16');
            NEV.ElectrodesInfo(ElectrodeID).ElectrodeID     = ElectrodeID;
            NEV.ElectrodesInfo(ElectrodeID).ConnectorBank   = char(ExtendedHeader(11)+64);
            NEV.ElectrodesInfo(ElectrodeID).ConnectorPin    = ExtendedHeader(12);
            df   = typecast(ExtendedHeader(13:14),'int16');
            % This is a workaround for the DigitalFactor overflow in NEV
            % files. Remove once Central is updated
            if df == 21516
                NEV.ElectrodesInfo(ElectrodeID).DigitalFactor = 152592.547;
            else
                NEV.ElectrodesInfo(ElectrodeID).DigitalFactor = df;
            end
            % End of workaround
            NEV.ElectrodesInfo(ElectrodeID).EnergyThreshold = typecast(ExtendedHeader(15:16),'uint16');
            NEV.ElectrodesInfo(ElectrodeID).HighThreshold   = typecast(ExtendedHeader(17:18),'int16');
            NEV.ElectrodesInfo(ElectrodeID).LowThreshold    = typecast(ExtendedHeader(19:20),'int16');
            NEV.ElectrodesInfo(ElectrodeID).Units           = ExtendedHeader(21);
            NEV.ElectrodesInfo(ElectrodeID).WaveformBytes   = ExtendedHeader(22);
            clear ElectrodeID;
        case 'NEUEVLBL'
            ElectrodeID                       = typecast(ExtendedHeader(9:10), 'uint16');
            NEV.ElectrodesInfo(ElectrodeID).ElectrodeLabel = char(ExtendedHeader(11:26));
            clear ElectrodeID;
        case 'NEUEVFLT'
            ElectrodeID                       = typecast(ExtendedHeader(9:10), 'uint16');
            NEV.ElectrodesInfo(ElectrodeID).HighFreqCorner = typecast(ExtendedHeader(11:14),'uint32');
            NEV.ElectrodesInfo(ElectrodeID).HighFreqOrder  = typecast(ExtendedHeader(15:18),'uint32');
            NEV.ElectrodesInfo(ElectrodeID).HighFilterType = typecast(ExtendedHeader(19:20),'uint16');
            NEV.ElectrodesInfo(ElectrodeID).LowFreqCorner  = typecast(ExtendedHeader(21:24),'uint32');
            NEV.ElectrodesInfo(ElectrodeID).LowFreqOrder   = typecast(ExtendedHeader(25:28),'uint32');
            NEV.ElectrodesInfo(ElectrodeID).LowFilterType  = typecast(ExtendedHeader(29:30),'uint16');
            clear ElectrodeID;
        case 'DIGLABEL'
            Mode                    = ExtendedHeader(25);
            NEV.IOLabels{Mode+1}    = char(ExtendedHeader(9:24).');
            clear Mode;
        case 'NSASEXEV' %% Not implemented in the Cerebus firmware.
                        %% Needs to be updated once implemented into the
                        %% firmware by Blackrock Microsystems.
            NEV.NSAS.Freq          = typecast(ExtendedHeader(9:10),'uint16');
            NEV.NSAS.DigInputConf  = char(ExtendedHeader(11));
            NEV.NSAS.AnalCh1Conf   = char(ExtendedHeader(12));
            NEV.NSAS.AnalCh1Detect = typecast(ExtendedHeader(13:14),'uint16');
            NEV.NSAS.AnalCh2Conf   = char(ExtendedHeader(15));
            NEV.NSAS.AnalCh2Detect = typecast(ExtendedHeader(16:17),'uint16');
            NEV.NSAS.AnalCh3Conf   = char(ExtendedHeader(18));
            NEV.NSAS.AnalCh3Detect = typecast(ExtendedHeader(19:20),'uint16');
            NEV.NSAS.AnalCh4Conf   = char(ExtendedHeader(21));
            NEV.NSAS.AnalCh4Detect = typecast(ExtendedHeader(22:23),'uint16');
            NEV.NSAS.AnalCh5Conf   = char(ExtendedHeader(24));
            NEV.NSAS.AnalCh5Detect = typecast(ExtendedHeader(25:26),'uint16');
        case 'VIDEOSYN'
            cnt = 1;
            if (isfield(NEV, 'VideoSyncInfo'))
                cnt = size(NEV.VideoSyncInfo, 2) + 1;
            end
            NEV.VideoSyncInfo(cnt).SourceID     = typecast(ExtendedHeader(9:10),'uint16');
            NEV.VideoSyncInfo(cnt).SourceName   = char(ExtendedHeader(11:26))';
            NEV.VideoSyncInfo(cnt).FrameRateFPS = typecast(ExtendedHeader(27:30),'single')';
            clear cnt;
        case 'TRACKOBJ'
            cnt = 1;
            if (isfield(NEV, 'ObjTrackInfo'))
                cnt = size(NEV.ObjTrackInfo, 2) + 1;
            end
            NEV.ObjTrackInfo(cnt).TrackableType = typecast(ExtendedHeader(9:10),'uint16');
            NEV.ObjTrackInfo(cnt).TrackableID   = typecast(ExtendedHeader(11:14), 'uint32');
            NEV.ObjTrackInfo(cnt).TrackableName = char(ExtendedHeader(15:30))';
            clear cnt;
        otherwise
            disp(['PacketID ' PacketID ' is invalid.']);
            disp('Please make sure this version of openNEV is compatible with your current NSP firmware.')
            fclose(FID);
            clear variables;
            if nargout
                varargout{1} = [];
            end
            return;
    end
end
NEV.MetaTags.ChannelID = [NEV.ElectrodesInfo.ElectrodeID];
clear ExtendedHeader PacketID ii;

%% Recording after ExtendedHeader file position and calculating Data Length
%  and number of data packets
fseek(FID, 0, 'eof');
Trackers.fData = ftell(FID);
Trackers.countDataPacket = (Trackers.fData - Trackers.fExtendedHeader)/Trackers.countPacketBytes;
NEV.MetaTags.PacketCount = Trackers.countDataPacket;

%%
Flags.UnparsedDigitalData = 0;

%% Reading packet headers and digital values
Timestamp = [];
PacketIDs = [];
tempClassOrReason = [];
tempDigiVals = [];
if NEV.MetaTags.PacketCount ~= 0
    fseek(FID, Trackers.fExtendedHeader, 'bof');
    tRawData  = fread(FID, [10 Trackers.countDataPacket], '10*uint8=>uint8', Trackers.countPacketBytes - 10);
    Timestamp = tRawData(1:4,:);
    Timestamp = typecast(Timestamp(:), 'uint32').';

    %% Calculate the number of packets that need to be read based on the time input parameters
    if ~exist('readTime', 'var')
        Trackers.readPackets = [1, length(Timestamp)];
    else
        [tmp,tempReadPackets] = find(Timestamp > readTime(1)*NEV.MetaTags.SampleRes,1,'first');
        if ~isempty(tempReadPackets)
            Trackers.readPackets(1) = tempReadPackets;
        else
            Trackers.readPackets(1) = NaN;
        end
        if isnan(Trackers.readPackets(1))
            fprintf('The file contains %0.2f seconds of data. The requested begining timestamp of %0.2f seconds is longer than the duration.\n', ...
                double(Timestamp(end))/double(NEV.MetaTags.SampleRes), ...
                readTime(1));
            clear variables;
            if nargout
                varargout{1} = [];
            end
            return;
        end
        [tmp,tempReadPackets] = find(Timestamp < readTime(2)*NEV.MetaTags.SampleRes,1,'last');
        if ~isempty(tempReadPackets)
            if readTime(2)*NEV.MetaTags.SampleRes > Timestamp(end)
                fprintf('The file contains %0.2f seconds of data. The requested end duration of %0.2f seconds will be adjusted to %0.2f seconds.\n', ...
                    double(Timestamp(end))/double(NEV.MetaTags.SampleRes), ...
                    readTime(2),...
                    double(Timestamp(end))/double(NEV.MetaTags.SampleRes));
            end
            Trackers.readPackets(2) = tempReadPackets;
        else
            Trackers.readPackets(2) = NaN;
        end
        clear tmp, tempReadPackets;
    end

    PacketIDs = tRawData(5:6,Trackers.readPackets(1):Trackers.readPackets(2));
    PacketIDs = typecast(PacketIDs(:), 'uint16').';
    tempClassOrReason = uint8(tRawData(7,Trackers.readPackets(1):Trackers.readPackets(2)));
    if strcmpi(Flags.digIOBits, '16bits')
        tempDigiVals      = tRawData(9:10,Trackers.readPackets(1):Trackers.readPackets(2));
        tempDigiVals      = typecast(tempDigiVals(:), 'uint16');
    else
        tempDigiVals      = uint16(tRawData(9,Trackers.readPackets(1):Trackers.readPackets(2)));
    end
    clear tRawData;
else
    Trackers.readPackets = zeros(1,2);
end

%% Defining PacketID constants
digserPacketID = 0;
neuralIndicesPacketIDBounds = [1, 16384];
commentPacketID = 65535;
videoSyncPacketID = 65534;
trackingPacketID = 65533;
patientTrigPacketID = 65532;
reconfigPacketID = 65531;

%% Parse read digital data. Please refer to help to learn about the proper
% formatting if the data.
digserIndices              = find(PacketIDs == digserPacketID);
neuralIndices              = find(neuralIndicesPacketIDBounds(2) >= PacketIDs & PacketIDs >= neuralIndicesPacketIDBounds(1));
commentIndices             = find(PacketIDs == commentPacketID);
videoSyncPacketIDIndices   = find(PacketIDs == videoSyncPacketID);
trackingPacketIDIndices    = find(PacketIDs == trackingPacketID);
patientTrigPacketIDIndices = find(PacketIDs == patientTrigPacketID);
reconfigPacketIDIndices    = find(PacketIDs == reconfigPacketID);
clear digserPacketID neuralIndicesPacketIDBounds commentPacketID ...
      videoSyncPacketID trackingPacketID patientTrigPacketID reconfigPacketID;
digserTimestamp            = Timestamp(digserIndices);
NEV.Data.Spikes.TimeStamp  = Timestamp(neuralIndices);
NEV.Data.Spikes.Electrode  = PacketIDs(neuralIndices);
clear PacketIDs;
NEV.Data.Spikes.Unit       = tempClassOrReason(neuralIndices);
%clear neuralIndices;
NEV.Data.SerialDigitalIO.InsertionReason   = tempClassOrReason(digserIndices);
clear tempClassOrReason;
DigiValues                 = tempDigiVals(digserIndices);
clear tempDigiVals;

%% Reads the waveforms if 'read' is passed to the function
if strcmpi(Flags.ReadData, 'read')
    allExtraDataPacketIndices  = [commentIndices, ...
                                  videoSyncPacketIDIndices, ...
                                  trackingPacketIDIndices, ...
                                  patientTrigPacketIDIndices, ...
                                  reconfigPacketIDIndices];

    if ~isempty(allExtraDataPacketIndices) % if there is any extra packets
        fseek(FID, Trackers.fExtendedHeader, 'bof');
        fseek(FID, (Trackers.readPackets(1)-1) * Trackers.countPacketBytes, 'cof');
        tRawData  = fread(FID, [Trackers.countPacketBytes Trackers.readPackets(2)], ...
            [num2str(Trackers.countPacketBytes) '*uint8=>uint8'], 0);
        if ~isempty(commentIndices)
            NEV.Data.Comments.TimeStamp = Timestamp(commentIndices);
            NEV.Data.Comments.TimeStampSec = double(NEV.Data.Comments.TimeStamp)/double(NEV.MetaTags.TimeRes);
            NEV.Data.Comments.CharSet = tRawData(7, commentIndices);
            NEV.Data.Comments.Color = tRawData(9:12, commentIndices);
            NEV.Data.Comments.Color = typecast(NEV.Data.Comments.Color(:), 'uint32').';
            NEV.Data.Comments.Text  = char(tRawData(13:Trackers.countPacketBytes, commentIndices).');

            % Transferring NeuroMotive Events to its own structure
            neuroMotiveEvents = find(NEV.Data.Comments.CharSet == 255);
            NEV.Data.TrackingEvents.TimeStamp = NEV.Data.Comments.TimeStamp(neuroMotiveEvents);
            NEV.Data.TrackingEvents.TimeStampSec = double(NEV.Data.TrackingEvents.TimeStamp)/double(NEV.MetaTags.TimeRes);
            NEV.Data.TrackingEvents.Text = NEV.Data.Comments.Text(neuroMotiveEvents,:);
            NEV.Data.Comments.TimeStamp(neuroMotiveEvents) = [];
            NEV.Data.Comments.TimeStampSec(neuroMotiveEvents) = [];
            NEV.Data.Comments.CharSet(neuroMotiveEvents) = [];
            NEV.Data.Comments.Color(neuroMotiveEvents) = [];
            NEV.Data.Comments.Text(neuroMotiveEvents) = [];

            clear commentIndices;
        end
        if ~isempty(videoSyncPacketIDIndices)
            NEV.Data.VideoSync.TimeStamp       = Timestamp(videoSyncPacketIDIndices);
            NEV.Data.VideoSync.FileNumber      = tRawData(7:8, videoSyncPacketIDIndices);
            NEV.Data.VideoSync.FileNumber      = typecast(NEV.Data.VideoSync.FileNumber(:), 'uint16').';
            NEV.Data.VideoSync.FrameNumber     = tRawData(9:12, videoSyncPacketIDIndices);
            NEV.Data.VideoSync.FrameNumber     = typecast(NEV.Data.VideoSync.FrameNumber(:), 'uint32').';
            NEV.Data.VideoSync.ElapsedTime     = tRawData(13:16, videoSyncPacketIDIndices);
            NEV.Data.VideoSync.ElapsedTime     = typecast(NEV.Data.VideoSync.ElapsedTime(:), 'uint32').';
            NEV.Data.VideoSync.SourceID        = tRawData(17:20, videoSyncPacketIDIndices);
            NEV.Data.VideoSync.SourceID        = typecast(NEV.Data.VideoSync.SourceID(:), 'uint32').';
            clear videoSyncPacketIDIndices;
        end
        if ~isempty(trackingPacketIDIndices)
            tmp.TimeStamp     = Timestamp(trackingPacketIDIndices);
            tmp.TimeStampSec  = double(tmp.TimeStamp)/30000;
            % This portion is commented out because it does not contain any
            % information as of yet.
            tmp.ParentID      = tRawData(7:8, trackingPacketIDIndices);
            tmp.ParentID      = typecast(tmp.ParentID(:), 'uint16').';
            tmp.NodeID        = tRawData(9:10, trackingPacketIDIndices);
            tmp.NodeID        = typecast(tmp.NodeID(:), 'uint16').';
            tmp.NodeCount     = tRawData(11:12, trackingPacketIDIndices);
            tmp.NodeCount     = typecast(tmp.NodeCount(:), 'uint16').';
            tmp.MarkerCount   = tRawData(13:14, trackingPacketIDIndices);
            tmp.MarkerCount   = typecast(tmp.MarkerCount(:), 'uint16').';

            tmp.rigidBodyPoints = tRawData(15:102, trackingPacketIDIndices);
            tmp.rigidBodyPoints = reshape(typecast(tmp.rigidBodyPoints(:), 'uint16'), 44, length(tmp.ParentID));

            objectIndex = [0 '1' '2' '3' '4' '1' '2' '3' '4'];
            if (isfield(NEV, 'ObjTrackInfo'))
                for IDX = 1:size(NEV.ObjTrackInfo,2)
                    emptyChar = find(NEV.ObjTrackInfo(IDX).TrackableName == 0, 1);
                    NEV.ObjTrackInfo(IDX).TrackableName(emptyChar) = objectIndex(IDX);
                    indicesOfEvent = find(tmp.NodeID == IDX-1);
                    if ~isempty(indicesOfEvent)
                        NEV.Data.Tracking.(NEV.ObjTrackInfo(IDX).TrackableName).TimeStamp = tmp.TimeStamp(indicesOfEvent);
                        NEV.Data.Tracking.(NEV.ObjTrackInfo(IDX).TrackableName).TimeStampSec = tmp.TimeStampSec(indicesOfEvent);
                        NEV.Data.Tracking.(NEV.ObjTrackInfo(IDX).TrackableName).ParentID =    tmp.ParentID(indicesOfEvent);
                        NEV.Data.Tracking.(NEV.ObjTrackInfo(IDX).TrackableName).NodeCount =   tmp.NodeCount(indicesOfEvent);
                        NEV.Data.Tracking.(NEV.ObjTrackInfo(IDX).TrackableName).MarkerCount = tmp.MarkerCount(indicesOfEvent);
                        NEV.Data.Tracking.(NEV.ObjTrackInfo(IDX).TrackableName).MarkerCoordinates(size(NEV.Data.Tracking.(NEV.ObjTrackInfo(IDX).TrackableName).TimeStamp,2)).X = [];
                        NEV.Data.Tracking.(NEV.ObjTrackInfo(IDX).TrackableName).MarkerCoordinates(size(NEV.Data.Tracking.(NEV.ObjTrackInfo(IDX).TrackableName).TimeStamp,2)).X = [];
                        for xyIDX = 1:size(NEV.Data.Tracking.(NEV.ObjTrackInfo(IDX).TrackableName).TimeStamp,2)
                            NEV.Data.Tracking.(NEV.ObjTrackInfo(IDX).TrackableName).MarkerCoordinates(xyIDX).X = ...
                                tmp.rigidBodyPoints(1:2:NEV.Data.Tracking.(NEV.ObjTrackInfo(IDX).TrackableName).MarkerCount(xyIDX)*2, indicesOfEvent(xyIDX));
                            NEV.Data.Tracking.(NEV.ObjTrackInfo(IDX).TrackableName).MarkerCoordinates(xyIDX).Y = ...
                                tmp.rigidBodyPoints(2:2:NEV.Data.Tracking.(NEV.ObjTrackInfo(IDX).TrackableName).MarkerCount(xyIDX)*2, indicesOfEvent(xyIDX));
                        end
                    end
                end
            end
            clear trackingPacketIDIndices tmp;
        end
        if ~isempty(patientTrigPacketIDIndices)
            NEV.Data.PatientTrigger.TimeStamp    = Timestamp(patientTrigPacketIDIndices);
            NEV.Data.PatientTrigger.TriggerType  = tRawData(7:8, patientTrigPacketIDIndices);
            NEV.Data.PatientTrigger.TriggerType  = typecast(NEV.Data.PatientTrigger.TriggerType(:), 'uint16').';
            clear patientTrigPacketIDIndices;
        end
        if ~isempty(reconfigPacketIDIndices)
            NEV.Data.Reconfig.TimeStamp     = Timestamp(reconfigPacketIDIndices);
            NEV.Data.Reconfig.ChangeType    = tRawData(7:8, reconfigPacketIDIndices);
            NEV.Data.Reconfig.ChangeType    = typecast(NEV.Data.Reconfig.ChangeType(:), 'uint16').';
            NEV.Data.Reconfig.CompName      = char(tRawData(9:24, reconfigPacketIDIndices));
            NEV.Data.Reconfig.ConfigChanged = char(tRawData(25:Trackers.countPacketBytes, reconfigPacketIDIndices));
            clear reconfigPacketIDIndices;
        end
    end % end if ~isempty(allExtraDataPacketIndices)

    clear Timestamp tRawData count idx;

   % now read waveform
    fseek(FID, Trackers.fExtendedHeader + 8, 'bof'); % Seek to location of spikes
    fseek(FID, (Trackers.readPackets(1)-1) * Trackers.countPacketBytes, 'cof');
    NEV.Data.Spikes.WaveformUnit = Flags.waveformUnits;
    NEV.Data.Spikes.Waveform = fread(FID, [(Trackers.countPacketBytes-8)/2 Trackers.readPackets(2)], ...
        [num2str((Trackers.countPacketBytes-8)/2) '*int16=>int16'], 8);
    NEV.Data.Spikes.Waveform(:, [digserIndices allExtraDataPacketIndices]) = [];

    clear allExtraDataPacketIndices;
    if strcmpi(Flags.waveformUnits, 'uv')
        elecDigiFactors = double(1000./[NEV.ElectrodesInfo(NEV.Data.Spikes.Electrode).DigitalFactor]);
        NEV.Data.Spikes.Waveform = bsxfun(@rdivide, double(NEV.Data.Spikes.Waveform), elecDigiFactors);
        if strcmpi(Flags.WarningStat, 'warning')
            fprintf(1,'\nThe spike waveforms are in unit of uV.\n');
            fprintf(2,'WARNING: This conversion may lead to loss of information.');
            fprintf(1,'\nRefer to help for more information.\n');
        end
    end
end
clear digserIndices;

%% Parse digital data if requested
if ~isempty(DigiValues)
    if strcmpi(Flags.ParseData, 'parse')
        try
            DigiValues = char(DigiValues);
            Inputs                     = {'Digital'; 'AnCh1'; 'AnCh2'; 'AnCh3'; 'AnCh4'; 'AnCh5'; 'PerSamp'; 'Serial'};
            AsteriskIndices   = find(DigiValues == '*');
            DataBegTimestamp = digserTimestamp(AsteriskIndices);
            splitDigiValues = regexp(DigiValues(2:end), '*', 'split')';
            for idx = 1:length(splitDigiValues)
                try
                    if isempty(find(splitDigiValues{idx} == ':', 1))
                        splitDigiValues{idx}(find(splitDigiValues{idx} == '#')) = [];
                        NEV.Data.SerialDigitalIO(idx).Value = splitDigiValues{idx};
                        NEV.Data.SerialDigitalIO(idx).Type = 'Marker';
                    else
                        [tempParsedCommand error] = parseCommand(splitDigiValues{idx});
                        if ~error
                            pcFields = fields(tempParsedCommand);
                            NEV.Data.SerialDigitalIO(idx).Value = splitDigiValues{idx};
                            for fidx = 1:length(pcFields)
                                NEV.Data.SerialDigitalIO(idx).(pcFields{fidx}) = tempParsedCommand.(pcFields{fidx});
                            end
                        else
                            NEV.Data.SerialDigitalIO(idx).Value = splitDigiValues{idx};
                            NEV.Data.SerialDigitalIO(idx).Type = 'UnparsedData';
                            Flags.UnparsedDigitalData = 1;
                        end
                    end
                catch
                    disp(['Error parsing: ' splitDigiValues{idx}]);
                    disp('Please refer to the help for more information on how to properly format the digital data for parsing.');
                end
            end
            % Populate the NEV structure with Timestamp and inputtypes for the
            % digital data
            if ~isempty(DataBegTimestamp)
                c = num2cell(DataBegTimestamp); [NEV.Data.SerialDigitalIO(1:length(NEV.Data.SerialDigitalIO)).TimeStamp] = deal(c{1:end});
                c = num2cell(DataBegTimestamp/NEV.MetaTags.SampleRes); [NEV.Data.SerialDigitalIO.TimeStampSec] = deal(c{1:end});
                c = {Inputs{NEV.Data.SerialDigitalIO.InsertionReason(AsteriskIndices)}}; [NEV.Data.SerialDigitalIO.InputType] = deal(c{1:end});
            end
            clear Inputs DigiValues digserTimestamp;
        catch
            disp(lasterr);
            disp('An error occured during reading digital data. This is due to a problem with formatting digital data.');
            disp('Refer to help ''help openNEV'' for more information on how to properly format the digital data.');
            disp('Try using openNEV with ''noparse'', i.e. openNEV(''noparse'').');
        end
    else
        NEV.Data.SerialDigitalIO.TimeStamp = digserTimestamp;
        NEV.Data.SerialDigitalIO.TimeStampSec = double(digserTimestamp)/30000;
        clear digserTimestamp;
        NEV.Data.SerialDigitalIO.UnparsedData = DigiValues;
        clear DigiValues;
    end
else
    if strcmpi(Flags.ReadData, 'read')
        if strcmpi(Flags.Report, 'report')
            disp('No digital data to read.');
        end
    end
end

if strcmpi(Flags.ParseData, 'parse')
    if Flags.UnparsedDigitalData && strcmpi(Flags.WarningStat, 'warning')
        fprintf(2, 'WARNING: The NEV file contains unparsed digital data.\n');
    end
end

%% Show a report if 'report' is passed as an argument
if strcmpi(Flags.Report, 'report')
    % Displaying report
    disp( '*** FILE INFO **************************');
    disp(['File Name           = ' NEV.MetaTags.Filename]);
    disp(['Filespec            = ' NEV.MetaTags.FileSpec]);
    disp(['Data Duration (min) = ' num2str(round(NEV.MetaTags.DataDuration/NEV.MetaTags.SampleRes/60))]);
    disp(['Packet Counts       = ' num2str(Trackers.countDataPacket)]);
    disp(' ');
    disp( '*** BASIC HEADER ***********************');
    disp(['Sample Resolution   = ' num2str(NEV.MetaTags.SampleRes)]);
    disp(['Date and Time       = '         NEV.MetaTags.DateTime]);
    disp(['Comment             = '         NEV.MetaTags.Comment(1:64)   ]);
    disp(['                      '         NEV.MetaTags.Comment(65:128) ]);
    disp(['                      '         NEV.MetaTags.Comment(129:192)]);
    disp(['                      '         NEV.MetaTags.Comment(193:256)]);
    disp(['The load time was for NEV file was ' num2str(toc, '%0.1f') ' seconds.']);
end

%% Saving the NEV structure as a MAT file for easy access
if strcmpi(Flags.SaveFile, 'save')
    if exist(matPath, 'file') == 2 && strcmpi(Flags.Overwrite, 'nooverwrite')
        disp(['File ' matPath ' already exists.']);
        overWrite = input('Would you like to overwrite (Y/N)? ', 's');
        if strcmpi(overWrite, 'y')
            disp('Saving MAT file. This may take a few seconds...');
            save(matPath, 'NEV', '-v7.3');
        else
            disp('File was not overwritten.');
        end
    elseif exist(matPath, 'file') == 2 && strcmpi(Flags.Overwrite, 'overwrite')
        disp(['File ' matPath ' already exists.']);
        disp('Overwriting the old MAT file. This may take a few seconds...');
        save(matPath, 'NEV', '-v7.3');
    else
        disp('Saving MAT file. This may take a few seconds...');
        save(matPath, 'NEV', '-v7.3');
    end
    clear overWrite;
end

if ~nargout
    assignin('base', 'NEV', NEV);
else
    varargout{1} = NEV;
end

fclose(FID);
clear Flags Trackers FID matPath;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function geometry_test036 ( )

%*****************************************************************************80
%
%% TEST036 tests SEGMENT_CONTAINS_POINT_1D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  ntest = 4;

  ptest =  [ 3.0,   7.5, 20.0,  5.0 ];
  p1test = [ 2.0,  10.0,  8.0, 88.0 ];
  p2test = [ 6.0, -10.0, 10.0, 88.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST036\n' );
  fprintf ( 1, '  SEGMENT_CONTAINS_POINT_1D determines if a point\n' );
  fprintf ( 1, '    lies within a line segment in 1D.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       P1     P       T\n' );
  fprintf ( 1, '\n' );

  for itest = 1 : ntest

    p1 = p1test(itest);
    p2 = p2test(itest);
    p = ptest(itest);

    t = segment_contains_point_1d ( p1, p2, p );
    fprintf ( 1, '  %10f  %10f  %10f  %12f\n', p1, p2, p, t );

  end

  return
end


text:
function stroud_test425 ( )

%*****************************************************************************80
%
%% TEST425 tests TRIANGLE_UNIT_SET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  order_max = 64;
  rule_max = 20;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST425\n' );
  fprintf ( 1, '  TRIANGLE_UNIT_SET sets up a quadrature \n' );
  fprintf ( 1, '    in the unit triangle,\n' )
  fprintf ( 1, '\n' );

  for a = 0 : 10

    for b = 0 : 10 - a

      coef = ( a + b + 2 ) * ( a + b + 1 );
      for i = 1 : b
        coef = coef * ( a + i ) / i;
      end

      fprintf ( 1, '\n' );
      fprintf ( 1, '  A = %f  B = %f\n', a, b );
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Rule       QUAD           ERROR\n' );
      fprintf ( 1, '\n' );

      for rule = 1 : rule_max

        order = triangle_unit_size ( rule );

        [ xtab, ytab, weight ] = triangle_unit_set ( rule, order );

        quad = 0.0;

        for i = 1 : order

          x = xtab(i);
          y = ytab(i);

          if ( a == 0 & b == 0 )
            value = coef;
          elseif ( a == 0 & b~= 0 )
            value = coef             * ytab(i)^b;
          elseif ( a ~= 0 & b == 0 )
            value = coef * xtab(i)^a;
          elseif ( a ~= 0 & b ~= 0 )
            value = coef * xtab(i)^a * ytab(i)^b;
          end

          quad = quad + 0.5 * weight(i) * value;

        end

        exact = 1.0;
        err = abs ( exact - quad );

        fprintf ( 1, '  %4d  %14f  %14f\n', rule, quad, err );

      end

    end

  end

  return
end


text:
function [ x, w ] = hermite_ek_compute ( n )

%*****************************************************************************80
%
%% HERMITE_EK_COMPUTE computes a Gauss-Hermite quadrature rule.
%
%  Discussion:
%
%    The code uses an algorithm by Elhay and Kautsky.
%
%    The abscissas are the zeros of the N-th order Hermite polynomial.
%
%    The integral:
%
%      integral ( -oo < x < +oo ) exp ( - x * x ) * f(x) dx
%
%    The quadrature rule:
%
%      sum ( 1 <= i <= n ) w(i) * f ( x(i) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2012
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer N, the number of abscissas.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights.
%

%
%  Define the zero-th moment.
%
  zemu = gamma ( 0.5 );
%
%  Define the Jacobi matrix.
%
  bj = zeros ( n, 1 );
  for i = 1 : n
    bj(i) = i / 2.0;
  end
  bj(1:n) = sqrt ( bj(1:n) );

  x = zeros ( n, 1 );

  w = zeros ( n, 1 );
  w(1) = sqrt ( zemu );
%
%  Diagonalize the Jacobi matrix.
%
  [ x, w ] = imtqlx ( n, x, bj, w );
%
%  If N is odd, force the center X to be 0.
%
  if ( mod ( n, 2 ) == 1 )
    x((n+1)/2) = 0.0;
  end

  w(1:n) = w(1:n).^2;

  return
end