text:
%% plotV
% Below is a demonstration of the features of the |plotV| function

%%
clear; close all; clc;

%% Syntax
% |varargout=plotV(V,varargin);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.


text:
function coef = plotnsdgt(coef,a,varargin)
%PLOTNSDGT Plot non-stationary Gabor coefficients
%   Usage:  plotnsdgt(c,a,fs,dynrange);
%
%   Input parameters:
%         coef     : Cell array of coefficients.
%         a        : Vector of time positions of windows.
%         fs       : signal sample rate in Hz (optional)
%         dynrange : Color scale dynamic range in dB (optional).
%
%   `plotnsdgt(coef,a)` plots coefficients computed using |nsdgt| or
%   |unsdgt|. For more details on the format of the variables *coef* and *a*,
%   please read the function help for these functions.
%
%   `plotnsdgt(coef,a,fs)` does the same assuming a sampling rate of
%   *fs* Hz of the original signal.
%
%   `plotnsdgt(coef,a,fs,dynrange)` additionally limits the dynamic range.
%
%   `C=plotnsdgt(...)` returns the processed image data used in the
%   plotting. Inputting this data directly to `imagesc` or similar
%   functions will create the plot. This is useful for custom
%   post-processing of the image data.
%
%   `plotnsdgt` supports all the optional parameters of |tfplot|. Please
%   see the help of |tfplot| for an exhaustive list. In addition, the
%   following parameters may be specified:
%
%     'xres',xres  Approximate number of pixels along x-axis / time.
%                  The default value is 800
%
%     'yres',yres  Approximate number of pixels along y-axis / frequency
%                  The default value is 600
%
%   See also: tfplot, nsdgt, unsdgt, nsdgtreal

%   AUTHOR : Florent Jaillet and Peter L. Søndergaard
%   TESTING: OK
%   REFERENCE: NA

if nargin<2
  error('%s: Too few input parameters.',upper(mfilename));
end;

definput.import={'ltfattranslate','tfplot'};

definput.keyvals.xres=800;
definput.keyvals.yres=600;

[flags,kv,fs]=ltfatarghelper({'fs','dynrange'},definput,varargin);

timepos=cumsum(a)-a(1);

N=length(a);
cwork=zeros(kv.yres,N);

%% -------- Interpolate in frequency ---------------------

for ii=1:N
  column=coef{ii};
  M=length(column);
  cwork(:,ii)=interp1(linspace(0,1,M),column,linspace(0,1,kv.yres),'nearest');
end;

%% --------  Interpolate in time -------------------------

% Time step in next equidistant spacing on the x-axis (in samples)
aplot=timepos(end)/kv.xres;

% Time positions where we want our pixels plotted (in samples)
xr=(0:kv.xres-1)*aplot;

% Move zero frequency to the center and Nyquist frequency to the top.
if rem(kv.yres,2)==0
  cwork=circshift(cwork,kv.yres/2-1);
else
  cwork=circshift(cwork,(kv.yres-1)/2);
end;

coef=zeros(kv.yres,kv.xres);
for ii=1:kv.yres
  data=interp1(timepos,cwork(ii,:).',xr,'nearest').';
  coef(ii,:)=data;
end;

yr=[-1+2/kv.yres,1];

coef=tfplot(coef,aplot,yr,'argimport',flags,kv);

if nargout<1
    clear coef;
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% Flush out the MATLAB.
close all;
clc;
clear all;

% Read the desired image file.
ImageData = imread('ct_scan.pnm');

% Add noise to oriiginal image.
ImageData = imnoise(ImageData,'Salt & Pepper', 0.04);

% Display the original image.
figure,imshow(ImageData);
title(' Original Image with noise: ');

% Take the input of the angle from user.
Angle = input('Enter the Angle = ');

% Convert the degree to radian.
AngleR = Angle*pi/180;

% Create the conversion matrix.
CMatrix = [+cos(AngleR) +sin(AngleR); -sin(AngleR) +cos(AngleR)];

% Calculate the size of the image.
[X,Y,Z] = size(ImageData);

Temp = round( [1 1; 1 Y; X 1; X Y]*CMatrix );

Temp = bsxfun(@minus, Temp, min(Temp)) + 1;

OutputImage = zeros([max(Temp) Z],class(ImageData));

% Implementation of rotaton function.
for a = 1:size(OutputImage,1)

    for b = 1:size(OutputImage,2)
        Rotation = ([a b]-Temp(1,:))*CMatrix.';

        if all(Rotation >= 1) && all(Rotation <= [X Y])
            CL = ceil(Rotation);

            FL = floor(Rotation);

            A = [...
                ((CL(2)-Rotation(2))*(CL(1)-Rotation(1))),...
                ((Rotation(2)-FL(2))*(Rotation(1)-FL(1)));

                ((CL(2)-Rotation(2))*(Rotation(1)-FL(1))),...
                ((Rotation(2)-FL(2))*(CL(1)-Rotation(1)))];

            Color = bsxfun(@times, A, double(ImageData(FL(1):CL(1),FL(2):CL(2),:)));

            OutputImage(a,b,:) = sum(sum(Color),2);
        end

    end

end

% Display the output image.
figure, imshow(OutputImage);
title(' Final Image: ');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [anchors, im_scales] = proposal_locate_anchors(conf, im_size, target_scale, feature_map_size)
% [anchors, im_scales] = proposal_locate_anchors(conf, im_size, target_scale, feature_map_size)
% --------------------------------------------------------
% Faster R-CNN
% Copyright (c) 2015, Shaoqing Ren
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------
% generate anchors for each scale

    % only for fcn
    if ~exist('feature_map_size', 'var')
        feature_map_size = [];
    end

    func = @proposal_locate_anchors_single_scale;

    if exist('target_scale', 'var')
        [anchors, im_scales] = func(im_size, conf, target_scale, feature_map_size);
    else
        [anchors, im_scales] = arrayfun(@(x) func(im_size, conf, x, feature_map_size), ...
            conf.scales, 'UniformOutput', false);
    end
end

function [anchors, im_scale] = proposal_locate_anchors_single_scale(im_size, conf, target_scale, feature_map_size)
    if isempty(feature_map_size)
        im_scale = prep_im_for_blob_size(im_size, target_scale, conf.max_size);
        img_size = round(im_size * im_scale);
        output_size = cell2mat([conf.output_height_map.values({img_size(1)}), conf.output_width_map.values({img_size(2)})]);
    else
        im_scale = prep_im_for_blob_size(im_size, target_scale, conf.max_size);
        output_size = feature_map_size;
    end

    shift_x = [0:(output_size(2)-1)] * conf.feat_stride;
    shift_y = [0:(output_size(1)-1)] * conf.feat_stride;
    [shift_x, shift_y] = meshgrid(shift_x, shift_y);

    % concat anchors as [channel, height, width], where channel is the fastest dimension.
    anchors = reshape(bsxfun(@plus, permute(conf.anchors, [1, 3, 2]), ...
        permute([shift_x(:), shift_y(:), shift_x(:), shift_y(:)], [3, 1, 2])), [], 4);
%   equals to
%     anchors = arrayfun(@(x, y) single(bsxfun(@plus, conf.anchors, [x, y, x, y])), shift_x, shift_y, 'UniformOutput', false);
%     anchors = reshape(anchors, [], 1);
%     anchors = cat(1, anchors{:});

end

author: jasjeetIM

repo: Mask-RCNN

save_path: github-repos/MATLAB/jasjeetIM-Mask-RCNN


text:
function [SRCC,PLCC] = test_TID2013(net,testDatabase, testMethod, Trainingproportion,seed)

    if ~exist(fullfile('result', testDatabase, testMethod), 'dir')
        mkdir(fullfile('result', testDatabase, testMethod));
    end

    fprintf('Results on TID2013...');
    Path = fullfile('databases', testDatabase);
    type_scores = cell(24,1);
    type_mos = cell(24,1);

    fp = fopen(fullfile(Path,'mos_with_names.txt'),'r');
    for i = 1:3000
        temp = fgetl(fp);
        mos(i) = str2num(temp(1:7));
        classes(i) = str2num(temp(10:11));
        disImagePath{i} = fullfile(Path,'distorted_images',temp(9:end));
        type(i) = str2num(temp(13:14));
        if i <= 2880
            refImagePath{i} = fullfile(Path,'reference_images',['I' temp(10:11) '.BMP']);
        else
            refImagePath{i} = fullfile(Path,'reference_images',['i' temp(10:11) '.bmp']);
        end
    end
    fclose(fp);

    if strcmp(testMethod,'CNN')
        [~, testingSet] = generateTrainingSet(classes,seed,Trainingproportion);
        for i = 1:numel(testingSet)
            scores(i) = processOneImage_CNN(net,disImagePath{testingSet(i)});
        end
        SRCC = corr(scores', mos(testingSet)', 'type', 'Spearman');
        PLCC = RegressionTID2013(scores', mos(testingSet)');
        fprintf('SRCC = %.4f;  PLCC = %.4f  \n', abs(SRCC), abs(PLCC))
    else
        if ~exist(fullfile('result', testDatabase, testMethod, 'scores.mat'), 'file')
            for i = 1:numel(disImagePath)
                fprintf('processing image %d / %d \n',i,numel(disImagePath));
                scores(i) = CallingTheSourceCodeForScorePrediction(disImagePath{i},testMethod);
            end
            save(fullfile('result', testDatabase, testMethod, 'scores.mat'),'scores');
        else
            load(fullfile('result', testDatabase, testMethod, 'scores.mat'),'scores');
        end
        for seed = 1:10
            [~, testingSet] = generateTrainingSet(classes,seed,Trainingproportion);
            SRCC(seed) = corr(scores(testingSet)', mos(testingSet)', 'type', 'Spearman');
            PLCC(seed) = RegressionTID2013(scores(testingSet)', mos(testingSet)');
        end
        fprintf('SRCC = %.4f, std = %.4f;  PLCC = %.4f, std = %.4f \n', abs(median(SRCC)), std(SRCC), abs(median((PLCC))), std((PLCC)));
    end

end

author: HuiZeng


text:
function X = coeffs2vals( X )
%COEFFS2VALS Convert an array of Chebyshev--Fourier--Fourier of
%            coefficients to an array of values.
%
%   COEFFS2VALS( CFS ) computes the values on a tensor-product doubled-up
%   spherical grid of a function that is represented by a
%   Chebyshev-Fourier-Fourier expansion with expansion coefficients CFS.
%
% See also CHEBFUN.COEFFS2VALS, VALS2COEFFS.

% Copyright 2019 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

if isempty( X )
    return
end

[m, n, p] = size( X );

% Old slow approach:
% for k = 1:p
%     X(:,:,k) = chebtech2.coeffs2vals( X(:,:,k) );
%     X(:,:,k) = trigtech.coeffs2vals( X(:,:,k).' ).';
% end
% for j = 1:n
%     vj = reshape( X(:,j,:), m, p );
%     vj = trigtech.coeffs2vals( vj.' ).';
%     X(:,j,:) = reshape( vj, m, 1, p );
% end

% Faster approach, but code less readable:
if ( m > 1 )
    X(2:m-1, :, :) = X(2:m-1, :, :)/2;
    X = fft( vertcat(X, X(m-1:-1:2,:,:)), [], 1 );
    X = X(m:-1:1, :, :);
end

scl_p = (n*p)*even_odd_fix( p );
scl_n = even_odd_fix( n );
Enp = reshape(scl_n.'*scl_p, [1 n p]);

X = ifft(ifft(ifftshift(ifftshift(X.*repmat( Enp, m, 1, 1),2),3), [], 2),[],3);

end

function scl = even_odd_fix( n )

if ( mod(n, 2) )
    scl = (-1).^(-(n-1)/2:(n-1)/2);
else
    scl = (-1).^((-n/2):(n/2-1));
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function G = make_G_matrix(T,g,varargin)

if length(g) == 1 && g < 0
    g=0;
end

G = spdiags(ones(T,1)*[-flipud(g(:))',1],-length(g):0,T,T);
if nargin == 3
    sl = [0;cumsum(varargin{1}(:))];
    for i = 1:length(sl)-1
        G(sl(i)+1,sl(i+1))=0;
    end
end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function l4mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% L4MAT_PRINT_SOME prints some of an L4MAT.
%
%  Discussion:
%
%    An L4MAT is an array of L4 values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, logical A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, a title.
%
  incx = 35;

  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );

  for j2lo = max ( jlo, 1 ) : incx : min ( jhi, n )

    j2hi = j2lo + incx - 1;
    j2hi = min ( j2hi, n );
    j2hi = min ( j2hi, jhi );

    inc = j2hi + 1 - j2lo;

    fprintf ( 1, '\n' );

    if ( 100 <= j2hi )
      fprintf ( 1, '      ' );
      for j = j2lo : j2hi
        fprintf ( 1, ' %1d', floor ( j / 100 ) );
      end
      fprintf ( 1, '\n' );
    end

    if ( 10 <= j2hi )
      fprintf ( 1, '      ' );
      for j = j2lo : j2hi
        fprintf ( 1, ' %1d', mod ( floor ( j / 10 ), 10 ) );
      end
      fprintf ( 1, '\n' );
    end

    fprintf ( 1, '  Col ' );
    for j = j2lo : j2hi
      fprintf ( 1, ' %1d', mod ( j, 10 ) );
    end
    fprintf ( 1, '\n' );

    fprintf ( 1, '  Row\n' );
    fprintf ( 1, '\n' );

    i2lo = max ( ilo, 1 );
    i2hi = min ( ihi, m );

    for i = i2lo : i2hi
      fprintf ( 1, '%5d:', i );
      for j = j2lo : j2hi
        fprintf ( 1, ' %1d', a(i,j) );
      end
      fprintf ( 1, '\n' );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [emp_cov, dist] = empCov2(x_obs, y_obs, data_obs, n_classes)
% SYNTAX:
%   [emp_cov, dist] = empCov2(x_obs, y_obs, data_obs, n_classes)
%
% DESCRIPTION:
%   return the array of the empirical covariance funtion
%
% EXAMPLE:
%   [emp_cov, dist] = empCov2(e_obs(id_td), n_obs(id_td), td_obs);
%   figure; plot(dist, emp_cov, 'b'); hold on;

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       Andrea Gatti
%  Contributors:     Andrea Gatti ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

    if nargin < 4
        n_classes = 50;
    end

    [x_mesh, y_mesh] = meshgrid(x_obs, y_obs);
    d_obs = sqrt(abs(x_mesh - x_mesh').^2 + abs(y_mesh - y_mesh').^2);
    classes = ceil(d_obs / max(d_obs(:)) * (n_classes-1))+1;

    emp_cov = zeros(max(classes(:)), size(data_obs,2));
    for i = 1 : size(data_obs,2)
        corr = (data_obs(:,i)-mean(data_obs(:,i))) * (data_obs(:,i)-mean(data_obs(:,i)))';
        for c = 1 : max(classes)
            emp_cov(c, i) = mean(corr(serialize(triu(classes)==c)));
        end
    end
    emp_cov = mean(emp_cov, 2);
    dist = (0 : (n_classes - 1))' * max(d_obs(:))/n_classes;
end

author: goGPS-Project


text:
%  Figure 10.67      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% fig10_67.m is a script to generate Fig. 10.67 the linear RTP response
% PI controller
% RTP chamber demo Example
clf;
a=[-0.068208813989939 0.014929776357245 0.000000065782442;...
   0.045809672136430 -0.118134773528570 0.021802006306129;...
   0.000000433637498 0.046839194867347 -0.100884399149872];
b3=[0.37873508304055 0.110575403544586 0.022912893467038;...
   0.000000002974222 0.449046982554739 0.073572900808161;...
   0.000000027324116 0.000702342292121 0.417797228783230];
c3=eye(3,3);
d3=0*eye(3,3);

% Combine 3 lamps into a single actuator
b=b3(:,1)+b3(:,2)+b3(:,3);
% Select center temperature
c=c3(2,:)
%
d=[0];
% PI controller
ac=[0];
bc=[1.0];
cc=[0.0527];
dc=[1];
sysG=ss(a,b,c,d);
sysD=ss(ac,bc,cc,dc);
sysL=series(sysD,sysG);
sysH=tf(1,1);
sysCL=feedback(sysL,sysH);

%CL Step Response
t=0:.1:100;
R=[0:.1:25, 25*ones(1,500), 0*ones(1,250)];
R11=[R'];
[yy,t]=lsim(sysCL,R11,t);
plot(t,yy,'--');
grid;
hold on;
plot(t,R11,'-');
xlabel('Time (sec)');
ylabel('Temperature (K)');
title('Fig. 10.67 (a) PI controller: temperature tracking response');
%legend('y')
pause
hold off
% Control effort
sysCLu=feedback(sysD,sysG);
[uuu,t]=lsim(sysCLu,R11,t);
plot(t(1:753,:),uuu(1:753,:));
hold on;
plot(t(752:818,:),0*ones(67,3));
hold on;
plot(t(753:818,:),uuu(753:818,:),'--');
hold on;
plot(t(819:1001,:),uuu(819:1001,:));
xlabel('Time (sec)');
ylabel('Lamp voltage');
title('Fig. 10.67 (b) PI controller :control effort');
legend('u')
grid;
hold off;

author: Sable

repo: mcbench-benchmarks


text:
function out = correl_compare_dep_permtest(y1,y2,varargin)
% Compare dependent correlations between pairs of vectors in y1 and y2
%
% :Usage:
% ::
%
%     out = correl_compare_dep_permtest(y1,y2,['alpha',myalpha],['rank'],['table'])
%
% :PERMUTATION TEST: for correl_compare_dep
%
% Repeats dep. correl. analysis for each pair of columns
% Returns results in correlation matrix form, where number of rows and
% cols. are the number of pairs [y1(:,i) y2(:,i)]
%
% myalpha is 2-tailed alpha value; p-values are 2-tailed
% FDR correction is at .01, 2-tailed
%
% Based on Steiger, 1980, tests for comparing dependent correlations.
%
% :Examples:
% ::
%
%    for i = 1:length(cl), y1(:,i) = cl.CONTRAST.data(:,2); y2(:,i) = cl.CONTRAST.data(:,1); end
%    for i = 1:length(cl), y1(:,i) = cl(i).CONTRAST.data(:,2); y2(:,i) = cl(i).CONTRAST.data(:,1); end
%    y1 is matrix of obs x data vectors for condition 1
%    y2 is matrix of obs x data vectors for condition 2
%    out = correl_compare_dep(y1,y2)
%
%    figure('Color','w');nmdsfig(c.GroupSpace,c.ClusterSolution.classes, ...
%    c.names,out.sig,1,{'Pos' 'Neg'});
%    nmdsfig_legend(c.ClusterSolution.X,c.r)
%
%    % compare correlations on cluster averages
%    c_compare = correl_compare_dep(y1avg,y2avg,'alpha',.05,'rank','table','names',c.APPLY_CLUSTER.names);

    myalpha = .05;
    dorankdata = 0;
    dotable = 0;
    names = [];

    for i = 1:length(varargin)
        if isstr(varargin{i})
            switch varargin{i}
                % functional commands
                case 'alpha', myalpha = varargin{i+1};
                case 'rank', dorankdata = 1;
                case 'table', dotable = 1;
                case 'names', names = varargin{i+1};
                otherwise, warning(['Unknown input string option:' varargin{i}]);
            end
        end
    end

    [N,npairs] = size(y1);

    [rows,cols,ncorr] = corrcoef_indices(npairs);

    if nargin < 4, dorankdata = 0; end
    if dorankdata
        str = sprintf('Ranking data: Nonparametric correlations'); fprintf(1,str);

        for i = 1:npairs
            y1(:,i) = rankdata(y1(:,i));
            y2(:,i) = rankdata(y2(:,i));
        end
    else
        str = sprintf('Assuming continuous data (no ranks).'); fprintf(1,str);
    end

    erase_string(str);

    i = 1; j = 2; k = 3; h = 4;

    % ----------------------------------------
    %%% Correct Permutation %%%
    % ----------------------------------------
    str = sprintf('Computing differences among correlations %04d',0); fprintf(1,str);

    [diffr,Zstar2,rr,pvec] = compare_dep_subfcn(ncorr,y1,y2,rows,cols,i,j,k,h,N);

    erase_string(str);

    % ----------------------------------------
    %%% All output stuff %%%
    % ----------------------------------------
    r1 = reconstruct(rr(:,1),npairs,ncorr,rows,cols);
    r2 = reconstruct(rr(:,2),npairs,ncorr,rows,cols);

    Z = reconstruct(Zstar2,npairs,ncorr,rows,cols);
    p = reconstruct(pvec,npairs,ncorr,rows,cols);

    dat = [rr diffr Zstar2 pvec];
    dat = [rows cols dat];
    dat = dat(pvec <= myalpha,:);

    diffr = reconstruct(diffr,npairs,ncorr,rows,cols);

    sig = (p <= myalpha - eye(size(p))) .* sign(diffr);

    % FDR corrected
    pthr = FDR(pvec,.05);
    if isempty(pthr), pthr = 0; end

    sigfdr = (p <= pthr) .* sign(diffr);

    out = struct('alpha',myalpha,'r1',r1,'r2',r2,'diffr',diffr, ...
        'Z',Z,'p',p,'sig',sig,'pthr',pthr,'sigfdr',sigfdr,'sigstats',dat);

    % output table...
    if dotable
        if isempty(names)
            disp(['No names entered; try ''names'' keyword to add them.']);
            for ii = 1:npairs, names{ii} = ['R' num2str(ii)]; end
        end

        disp(['Uncorrected, p < ' num2str(myalpha)])
        maketable(out,'sig',names);
        disp('')

        disp(['FDR corrected, p < ' num2str(myalpha)])
        maketable(out,'sigfdr',names);
        disp('')
    end

    % ----------------------------------------
    %%% Permutations %%%
    % ----------------------------------------
    niter = 5000;

    [out.nalpha,out.n01,out.nalpha_posvsneg,out.n01_posvsneg] = do_perms;
    out.crit_n_at_alpha = prctile(out.nalpha,95);
    out.crit_n_at_01 = prctile(out.n01,95);
    out.crit_sumsigned_at_alpha = prctile(out.nalpha_posvsneg,97.5);
    out.crit_sumsigned_at_01 = prctile(out.n01_posvsneg,97.5);

    out.numsig_at_alpha = sum(squareform(out.p) <= out.alpha);
    out.numsig_at_01 = sum(squareform(out.p) <= .01);

    out.sumsigned_at_alpha = sum(  (squareform(out.p) <= out.alpha) .* sign(squareform(out.diffr))  );
    out.sumsigned_at_01 = sum(  (squareform(out.p) <= .01) .* sign(squareform(out.diffr))  );

    if dotable
        fprintf(1,'Alpha: %3.4f, Num. sig: %3.0f, Critical num sig: %3.0f\n',out.alpha,out.numsig_at_alpha,out.crit_n_at_alpha);
        fprintf(1,'Alpha: %3.4f, Num. sig: %3.0f, Critical num sig: %3.0f\n',.01,out.numsig_at_01,out.crit_n_at_01);

        fprintf(1,'Alpha: %3.4f, Sum of signed sig.: %3.0f, Critical num 2-tailed: %3.0f\n',out.alpha,out.sumsigned_at_alpha, out.crit_sumsigned_at_alpha);
        fprintf(1,'Alpha: %3.4f, Sum of signed sig.: %3.0f, Critical num 2-tailed: %3.0f\n',.01,out.sumsigned_at_01, out.crit_sumsigned_at_01);

    end

    % nested function
    % get number of significant effects under null hypothesis at alpha and at
    % .01
    function [nalpha,n01,nalpha_posvsneg,n01_posvsneg] = do_perms

        nalpha = zeros(1,niter);
        n01 = nalpha;
        nalpha_posvsneg = nalpha;
        n01_posvsneg = nalpha;

        fprintf(1,'Getting distribution for number of significant correlations.\n  Running %3.0f permutations: %04d',niter,0);

        for ii = 1:niter
            fprintf(1,'\b\b\b\b%04d',ii);

            % permute data
            for jj = 2:npairs
                wh = randperm(N);
                y1(:,jj) = y1(wh,jj);
                y2(:,jj) = y2(wh,jj);
            end

            % test correlations
            [diffr,Zstar2,rr,pvec] = compare_dep_subfcn(ncorr,y1,y2,rows,cols,i,j,k,h,N,0);

            nalpha(ii) = sum(pvec <= myalpha);
            n01(ii) = sum(pvec <= .01);

            nalpha_posvsneg(ii) = sum(  (pvec <= myalpha) .* sign(diffr)  );
            n01_posvsneg(ii) = sum(  (pvec <= .01) .* sign(diffr)  );

        end
        fprintf(1,'\n');

    end

    end   %%% end main function

    function [diffr,Zstar2,rr,pvec] = compare_dep_subfcn(ncorr,y1,y2,rows,cols,i,j,k,h,N,dotext)

    diffr = zeros(ncorr,1);     % correlation difference
    Zstar2 = zeros(ncorr,1);    % from Steiger, the Z-score
    rr = zeros(ncorr,2);        % the correl values for y1 and y2

    if ~exist('dotext','var'), dotext = 1; end

    for cc = 1:ncorr

        if dotext, fprintf(1,'\b\b\b\b%04d',cc); end

        % get full correlation matrix for the 4 variables involved
        dat = [y1(:,[rows(cc) cols(cc)]) y2(:,[rows(cc) cols(cc)])];

        % get differences between correlation z-values (estimate)
        [diffr(cc),r,rr(cc,:),diffrz,z] = correlation_diffs(dat,i,j,k,h);

        s = covcorr(r,i,j,k,h);         % covariance of corr coeffs

        Zstar2(cc) = diffrz * sqrt( (N-3) ./ (2-(2*s)) );

        % bootstrap
        %vals = bootstrp(5000,@correlation_diffs,dat,i,j,k,h);
        %Zboot(cc) = mean(vals) ./ std(vals);

        %if Zstar2(cc) > 2, keyboard, end
    end

    pvec = 2 * (1 - normcdf(abs(Zstar2)));

    end

function [diffr,r,rr,diffrz,z] = correlation_diffs(dat,i,j,k,h)
    % get differences between correlation z-values
    r = corrcoef(dat);

    rr = [r(i,j) r(k,h)];                % correls to be compared
    diffr = -diff(rr);              % negative sign means we get (1) - (2)

    if nargout > 3
        z = .5 .* log( (1+rr) ./ (1-rr) );
        diffrz = -diff(z);               % diff btwn correl z values
    end

    end

function [rows,cols,ncorr] = corrcoef_indices(npairs)
    % upper triangle only
    tmp = triu(ones(npairs));
    tmp = tmp - eye(npairs);
    [rows,cols] = find(tmp);
    ncorr = length(rows);
    end

function s = covcorr(R,i,j,k,h)

    % pool correl. coeff for more stable var est.
    % as they are equal under Ho.  Steiger, 1980, eq. 14 for Z*
    pooledr = mean([R(i,j) R(k,h)]);
    s = pearsonf(R,i,j,k,h) ./  ( corvar(pooledr) );

    end

function cv = corvar(r)
    % variance of a correlation coefficient
    % simplified form (special case) of pearsonf
    cv = (1 - r^2) ^ 2;
    end

    % % function cv = co(R,x,y)
    % % % x and y are 2-vector indices of matrix R to compute covariance for
    % % % cv is covariance
    % % cv = pearsonf(R,x(1),x(2),y(1),y(2));
    % % end

function pf = pearsonf(R,i,j,k,h)

    % Pearson-Filon: covariance (or var) of element i,j with element k,h
    % depending on correlation values

    % for variance of 1 correl, works as well, and reduces to:
    % (1 - rr(1)^2)^2

    pf = (1/2) .* R(i,j) .* R(k,h) .* ...
        ( R(i,k).^2 + R(i,h).^2 + R(j,k).^2 + R(j,h).^2 ) + ...
        R(i,k) .* R(j,h) + R(i,h) .* R(j,k) - ...
        R(i,j) .* ( R(j,k) .* R(j,h) + R(i,k) .* R(i,h) ) - ...
        R(k,h) .* ( R(j,k) .* R(i,k) + R(j,h) .* R(i,h) );

    end

    % function pf = pearsonf2(r,j,k,h,m)
    %
    % pf = .5 * ( (r(j,h) - r(j,k)*r(k,h)) * (r(k,m) - r(k,h)*r(h,m))  + ...
    %     (r(j,m) - r(j,h)*r(h,m)) * (r(k,h) - r(k,j)*r(j,h)) + ...
    %     (r(j,h) - r(j,m)*r(m,h)) * (r(k,m) - r(k,j)*r(j,m)) + ...
    %     (r(j,m) - r(j,k)*r(k,m)) * (r(k,h) - r(k,m)*r(m,h))    );
    %
    % end

function valmat = reconstruct(vals,npairs,ncorr,rows,cols)

    valmat = zeros(npairs);
    for i = 1:ncorr
        valmat(rows(i),cols(i)) = vals(i);
    end
    valmat = valmat + valmat';

    end

function erase_string(str1)
    fprintf(1,repmat('\b',1,length(str1))); % erase string
    end

% function maketable(c_compare,whfield,names)
%     [rows,cols] = find(triu(c_compare.(whfield)));
%     if isempty(rows)
%         disp('No significant results.')
%     else
%         fprintf(1,'Name1\tName2\trow\tcol.\t+ or -\tZ\tp\n');
%         for i = 1:length(rows)
%             fprintf(1,'%s\t%s\t%3.0f\t%3.0f\t%3.0f\t%3.2f\t%3.4f\n',names{rows(i)}, ...
%                 names{cols(i)},rows(i),cols(i),c_compare.(whfield)(rows(i),cols(i)),c_compare.Z(rows(i),cols(i)),c_compare.p(rows(i),cols(i)));
%         end
%     end
%     fprintf(1,'\n');
%     end

function maketable(c_compare,whfield,names)

    str = {'-' '+'};
    [rows,cols] = find(triu(c_compare.(whfield)));
    if isempty(rows)
        disp('No significant results.')
    else
        fprintf(1,'Name1\tName2\trow\tcol.\t+ or -\tr1\tr2\tZ\tp\n');
        for i = 1:length(rows)
            fprintf(1,'%s\t%s\t%3.0f\t%3.0f\t%s\t%3.3f\t%3.3f\t%3.2f\t%3.4f\n', ...
                names{rows(i)}, names{cols(i)},rows(i),cols(i), ...
                str{1.5 + (.5.*c_compare.(whfield)(rows(i),cols(i)))}, ...          % + or - sign
                c_compare.r1(rows(i),cols(i)),c_compare.r2(rows(i),cols(i)), ...   % correlations
                c_compare.Z(rows(i),cols(i)),c_compare.p(rows(i),cols(i)));        % Z and p
        end
    end
    fprintf(1,'\n');
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function times = LoadSpikeTimes(filename,rate)

%LoadSpikeTimes - Load spike times from file.
%
%  USAGE
%
%    times = LoadSpikeTimes(filename,rate)
%
%    filename            spike file name (either .clu or .res)
%    rate                sampling rate
%
%  OUTPUT
%
%    The output is a list of (timestamp,group,cluster) t-uples.
%
%  SEE
%
%    See also GetSpikeTimes, PlotTicks.

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

warning('This function is now deprecated from /buzcode, try using bz_GetSpikes.m')

[path,basename,extension] = fileparts(filename);
if isempty(path), path = '.'; end

electrodeGroup = str2num(extension(2:end));
[unused,basename,unused] = fileparts(basename);

% Load .res file
filename = [path '/' basename '.res.' int2str(electrodeGroup)];
if ~exist(filename),
	error(['File ''' filename ''' not found.']);
end
res = load(filename);

% Load .clu file
filename = [path '/' basename '.clu.' int2str(electrodeGroup)];
if ~exist(filename),
	error(['File ''' filename ''' not found.']);
end
clu = load(filename);

times = [res/rate electrodeGroup*ones(size(res)) clu(2:end)];

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% plotmesh() - plot mesh defined by faces and vertex
%
% Usage:
%     plotmesh(faces, vertex);
%
% Input:
%   faces   - array of N x 3. Each row defines a triangle. The 3 points
%             in each row are row indices in the matrix below.
%   vertex  - array of M x 3 points, (x = first colum; y=second colum
%             z=3rd column). Each row defines a point in 3-D.
%
% Optional input:
%   normal  - normal orientation for each face (for better lighting)
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, 2003

% Copyright (C) May 6, 2003 Arnaud Delorme, SCCN/INC/UCSD,
% arno@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function p1 = plotmesh(faces, vertex, normal, newfig)

    if nargin < 2
        help plotmesh;
        return;
    end;

    FaceColor  = [.8 .55 .35]*1.1; % ~= ruddy Caucasian - pick your complexion!

    if any(any(faces == 0)), faces = faces+1; end;
    %vertex(:,3) = -vertex(:,3);
    %FCmap = [jet(64); FaceColor; FaceColor; FaceColor];
    %colormap(FCmap)
    %W = ones(1,size(vertex,1))*(size(FCmap,1)-1);
    %W = ones(1,size(vertex,1))' * FaceColor;
    %size(W)
    if nargin < 4
        figure;
    end;
    if nargin < 3
        normal = [];
    end;
    if isempty(normal)
        p1 = patch('vertices', vertex, 'faces', faces, ...
                   'facecolor', [1,.75,.65]);
    else
        p1 = patch('vertices', vertex, 'faces', faces, ...
                   'facecolor', [1,.75,.65], 'vertexnormals', normal);
    end;
        %           'FaceVertexCdata',W(:), 'FaceColor','interp', 'vertexnormals', normal);
    set(p1,'EdgeColor','none')

    % Lights
    %Lights = [-125  125  80; ...
    %          125  125  80; ...
    %          125 -125 125; ...
    %          -125 -125 125];    % default lights at four corners
    %for i = 1:size(Lights,1)
    %    hl(i) = light('Position',Lights(i,:),'Color',[1 1 1],...
    %                  'Style','infinite');
    %end
    %camlight left;
    lightangle(45,30);
    lightangle(45+180,30);
    %set(gcf, 'renderer', 'zbuffer'); % cannot use alpha then
    %set(hcap, 'ambientstrength', .6);
    set(p1, 'specularcolorreflectance', 0, 'specularexponent',50);

    set(p1,'DiffuseStrength',.6,'SpecularStrength',0,...
           'AmbientStrength',.4,'SpecularExponent',5);
    axis equal
    view(18,8);
    rotate3d
    axis off;
    lighting phong;

author: PatternRecognition

repo: OpenBMI


text:
function [err_rel,err_abs] = derivcheck(f,x,flag)
%DERIVCHECK Check analytical vs numerical differentiation for a function

if nargin < 3 || isempty(flag); flag = false; end

tic
if flag
    dy_num = fgrad(f,x,'five-points');
else
    dy_num = gradest(f,x);
end
toc
tic
[y,dy_ana] = f(x);
toc

if size(dy_num,1) == size(dy_num,2)
    dy_num = sum(dy_num,1);
end

% Reshape to row vectors
dy_num = dy_num(:)';
dy_ana = dy_ana(:)';

fprintf('Relative errors:\n');
err_rel = (dy_num(:)' - dy_ana(:)')./dy_num(:)'

fprintf('Absolute errors:\n');
err_abs = dy_num(:)' - dy_ana(:)'

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function nn = nnff(nn, x, y)
%NNFF performs a feedforward pass
% nn = nnff(nn, x, y) returns an neural network structure with updated
% layer activations, error and loss (nn.a, nn.e and nn.L)

    n = nn.n;
    m = size(x, 1);

    nn.a{1} = x;

    %feedforward pass
    for i = 2 : n-1
        nn.a{i} = sigm(repmat(nn.b{i - 1}', m, 1) + nn.a{i - 1} * nn.W{i - 1}');
        if(nn.dropoutFraction > 0)
            if(nn.testing)
                nn.a{i} = nn.a{i}.*(1 - nn.dropoutFraction);
            else
                nn.a{i} = nn.a{i}.*(rand(size(nn.a{i}))>nn.dropoutFraction);
            end
        end
        %calculate running exponential activations for use with sparsity
        if(nn.nonSparsityPenalty>0)
            nn.p{i} = 0.99 * nn.p{i} + 0.01 * mean(nn.a{i}, 1);
        end
    end
    switch nn.output
        case 'sigm'
            nn.a{n} = sigm(repmat(nn.b{n - 1}', m, 1) + nn.a{n - 1} * nn.W{n - 1}');
        case 'linear'
            nn.a{n} = repmat(nn.b{n - 1}', m, 1) + nn.a{n - 1} * nn.W{n - 1}';
        case 'softmax'
            nn.a{n} = repmat(nn.b{n - 1}', m, 1) + nn.a{n - 1} * nn.W{n - 1}';
            nn.a{n} = exp(bsxfun(@minus, nn.a{n}, max(nn.a{n},[],2)));
            nn.a{n} = bsxfun(@rdivide, nn.a{n}, sum(nn.a{n}, 2));
    end

    %error and loss
    nn.e = y - nn.a{n};
    switch nn.output
        case {'sigm','linear'}
            nn.L = 1/2 * sum(sum(nn.e .^ 2)) / m;
        case 'softmax'
            nn.L = -sum(sum(y .* log(nn.a{n}))) / m;
    end
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function [irm,ibm] = idealMasks(st,si,q,a)
%IDEALMASKS Calculate ideal time-frequency masks from STFTs
%
%   IRM = IOSR.BSS.IDEALMASKS(ST,SI) calculates the ideal ratio mask (IRM)
%   using target STFT ST and interferer STFT SI. ST and SI must be the same
%   size; IRM is the same size as ST and SI.
%
%   [IRM,IBM] = IOSR.BSS.IDEALMASKS(...) returns the ideal binary mask
%   (IBM).
%
%   IRM = IOSR.BSS.IDEALMASKS(ST,SI,Q) uses the exponent Q to create an
%   ideal sigmoidal mask by raising each of the time-frequency powers to Q.
%   With Q=1 (default) the mask is the IRM. As Q->Inf the IRM will tend
%   towards a binary mask. As Q->0 the mask will tend towards 0.5.
%
%   [IRM,IBM] = IOSR.BSS.IDEALMASKS(ST,SI,Q,A) uses the threshold A to
%   calculate the IBM. The threshold is in terms of the ratio of time
%   frequency powers. With A=1 (default), the IBM is 1 when the target
%   power is greater than the interferer power. Setting A=2, for example,
%   requires the target power to be 6dB greater than the interference
%   power.
%
%   See also IOSR.BSS.APPLYMASK, IOSR.BSS.APPLYIDEALMASKS.

%   Copyright 2016 University of Surrey.

    %% check input

    assert(isequal(size(st),size(si)), 'iosr:idealMasks:invalidInputs', 'ST and SI must be the same size')

    % check sigmoid
    if nargin<3
        q = 1;
    else
        assert(isscalar(a), 'iosr:idealMasks:invalidQ', 'Q must be an scalar')
    end

    % check threshold
    if nargin<4
        a = 1;
    else
        assert(isscalar(a), 'iosr:idealMasks:invalidA', 'A must be an scalar')
    end

    %% calculate masks

    % powers
    St = abs(st).^2;
    Si = abs(si).^2;

    % ideal ratio mask
    irm = (St.^q)./((St.^q)+(Si.^q));
    irm(isnan(irm)) = 0;

    % ideal binary mask
    if nargout>1
        ibm = +(St./Si>a);
        ibm(isnan(ibm)) = 0;
    end

end

author: IoSR-Surrey

repo: MatlabToolbox


text:
% This function performs the binary decision tree training.

% Copyright (C) 2012 Quan Wang <wangq10@rpi.edu>,
% Signal Analysis and Machine Perception Laboratory,
% Department of Electrical, Computer, and Systems Engineering,
% Rensselaer Polytechnic Institute, Troy, NY 12180, USA
%
% You are free to use this software for academic purposes if you cite our paper:
% Q. Wang, Y. Ou, A.A. Julius, K.L. Boyer, M.J. Kim,
% Tracking tetrahymena pyriformis cells using decision trees,
% in: 2012 International Conference on Pattern Recognition, Tsukuba Science City, Japan.
%
% For commercial use, please contact the authors.

function T=create01Tree(X,label,Depth,Splits,MinNode)

% X: N*M training data matrix, each row is one sample
% label: N*1 labels of training data, each entry takes value 0 or 1
% Depth: maximal depth of decision tree
% Splits: number of candidate thresholds at each node
% MinNode: minimal size of a non-leaf node
% T: resulting decision tree

T=struct('feature',zeros(2^(Depth-1)-1,1),...
    'threshold',zeros(2^(Depth-1)-1,1),...
    'leaf',zeros(2^(Depth-1),2),...
    'leaf_prob',zeros(2^(Depth-1),1),...
    'depth',Depth);

T=get01Tree(1,T,X,label,Depth,Splits,MinNode);


text:
function im1 = at_imageresize(im1)

% isz = size(im1(:,:,1));
% if (1920*1440) < prod(isz)
%   if isz(1) > isz(2)
%     im1 = imresize(im1,[1920 NaN]);
%   else
%     im1 = imresize(im1,[NaN 1920]);
%   end
% end

isz = size(im1(:,:,1));
if (1600*1200) < prod(isz)
  scale = 1600/max(isz);
  im1 = imresize(im1, scale);
%   if isz(1) > isz(2)
%     im1 = imresize(im1,[640 NaN]);
%   else
%     im1 = imresize(im1,[NaN 640]);
%   end
end

% isz = size(im1(:,:,1));
% if (1920*1440) < prod(isz)
%   im1 = imresize(im1,sqrt((1920*1440)/prod(isz)));
% end

author: HajimeTaira


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD 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.
%
% TLD 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 TLD.  If not, see <http://www.gnu.org/licenses/>.

function tld = tldUpdateDetectorConservative(tld,I)

bb    = tld.bb(:,I);
img   = tld.img{I};

% Consistency =============================================================

pPatt  = tldGetPattern(img,bb,tld.model.patchsize);
[pConf1,~,pIsin] = tldNN(pPatt,tld);

% Inconsistency: colildes with negative
if pConf1 < 0.5
%     disp('Fast change.');
    tld.valid(I) = 0;
    return;
end
if pIsin(3) == 1 || var(pPatt) < tld.var
%     disp('In negative data.');
    return;
end

% If strong position, identify indexes that will be removed and updated
pRank  = pIsin(2);
nRank  = []; % indexes to patches that will be removed from PEX
nEx    = []; % patches that will be added to NEX
if ~isnan(pRank)

    id_far = bb_overlap(bb,tld.dt{I}.bb) < 0.2;
    id_inP = tld.dt{I}.isin(1,:) == 1;
    nRank  = tld.dt{I}.isin(2,id_far & id_inP);

    % Inconsistency: negative patterns is in PEX with lower rank
    if any(nRank<=pRank),
        disp('inconsistent 2');
        return;
    end

    nEx = tld.dt{I}.patch(:,id_far);
end
tld.nEx{I} = nEx;
% Fern ====================================================================

% Generate positive patterns from current frame
overlap  = bb_overlap(bb,tld.grid);
[pX,pEx] = tldGeneratePositiveData(tld,overlap,img,tld.p_par_update);
pY       = ones(1,size(pX,2));

% Always use negative patterns from the first frame
% nX       = tld.X{1}(:,tld.Y{1}==0);
% idx      = pseudorandom_indexes(size(nX,2),10);
% nX       = nX(:,idx);

% If the current frame is valid, use negativ data from surrounding
% if pConf2 > tld.model.thr_nn_valid
idx = overlap < tld.n_par.overlap & tld.tmp.conf >= 0;

    %idx = pseudorandom_indexes(size(nX,2),10);
    %nX = nX(:,idx);

% Update fern
fern(2,[pX tld.tmp.patt(:,idx)],[pY zeros(1,sum(idx))],tld.model.thr_fern,2);

% Debug
if 0
    figure(10); clf;
    subplot(121);
    imshow(img.input);bb_draw(tld.dt{I}.bb); title('before update');
    subplot(122);
    fern(4,img,tld.control.maxbbox,tld.var,tld.tmp.conf,tld.tmp.patt);
    id_dt  = tld.tmp.conf > tld.model.thr_fern*tld.model.num_trees;

    imshow(img.input); bb_draw(tld.grid(:,id_dt));
end

% Nearest neighbour =======================================================

% Remove patterns from PEX
tld.pex(:,nRank) = [];

% Always use negative patches from the first frame
idx = pseudorandom_indexes(size(tld.nEx{1},2),10);
nEx = [nEx tld.nEx{1}(:,idx) tld.nEx{I-1}];

% Update NN
% TODO sort
tld = tldTrainNN(pEx,nEx,tld);

% sound(tld.gong(1:200));

repo: OpenTLD


text:
function ccff_asymptotic ( n_min, n_inc, n_max, f, f_integral )

%*****************************************************************************80
%
%% CCFF_ASYMPTOTIC computes asymptotic errors for the Legendre integral.
%
%  Discussion:
%
%    The Legendre integral being approximated is:
%      integral ( -1 <= x <= +1 ) f(x) dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N_MIN, N_INC, N_MAX, the minimum, increment, and maximum
%    for the number of points in the rule.
%
%    Input, function pointer F(x), returns the value of the integrand
%    at the N points X.
%
%    Input, real F_INTEGRAL, the exact value of the integral.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N  |Quad error|\n' );
  fprintf ( 1, '\n' );
  for n = n_min : n_inc : n_max
    [ x, w ] = ccff ( n );
    fx = f ( x );
    q = w' * fx;
    e = abs ( f_integral - q );
    fprintf ( 1, '  %4d  %8.2e\n', n, e );
  end

  return
end


text:
function [endpos, maxlevel] = match_bracket(str,startpos,brackets)
%
% [endpos, maxlevel] = match_bracket(str,startpos,brackets)
%
% Looking for the position of a closing bracket token in a string
%
% authors:Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%      str: the full string to be searched
%      startpos: the index in the string as the start position to search; the
%               startpos must be at least 1 greater than the opening bracket position
%      brackets: (optional), a string of length 2, with the first character
%               being the opening token and the 2nd being the closing token.
%               if not given, brackets is set to '[]' to find matching square-brackets;
%               for example, '{}' looks for a matching closing curly-bracket in
%               the string key(pos(startpos,:end))
%
% output:
%      endpos: if a matching bracket is found, return its position in the original
%              string
%      maxlevel: return the depth of the enclosed brackets between the searched pair,
%              includig the searching pair. For example, the matching closing-bracket
%              of the 1st square bracket (startpos=2) in  '[[[]],[]]' returns a
%              position of 9, with a maximum depth of 3; searching for the closing
%              bracket for the 2nd square bracket (startpos=3) returns a position of
%              5 and max-depth of 2.
%
% example:
%      str='[[ [1,2], 1], 10, [5,10] ]';
%      [p1,dep]=match_bracket(str,3)
%      [p2,dep]=match_bracket(str,2)
%      [p3,dep]=match_bracket(str,3)
%
% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

if(nargin<3)
    brackets='[]';
end
count = str(startpos:end);
flag=cumsum(count==brackets(1))-cumsum(count==brackets(2))+1;
endpos = find(flag==0,1);
maxlevel=max(flag(1:endpos));
endpos = endpos + startpos-1;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [timeqrs,indexes,tol,messages]=qrscandidatesnew(timeqrs1,timeqrs2,timeqrs3,heasig,tol,messages)
% [timeqrs,indexes]=qrscandidates(timeqrs1,timeqrs2,timeqrs3,heasig,tol)
%
% Construction QRS positions candidates vector from 3 single detections vectors
% a QRS is admitted as candidate if:
%   - it was found in 3 leads and the middle mark differ less than tol from the others
%   - it was found in 2 leads and the marks differ less than tol
%
% marks that not fullfill the above are replaced by NaN and a QRS positions
% candidates vector is constructed with the detections in each lead in one
% line, one QRS in each column
%
% INPUT:
% timeqrs1,timeqrs2,timeqrs3  -  line vector swith single lead detections
%                                (in samples)
% heasig - header information
% tol - maximum distance to be admited as same complex (in sec)
%       by default half of the default refractary period tol = 0.275/2 sec
%
% OUTPUT:
% timeqrs - QRS positions candidates vector
% indexes - beats corresponding to the QRS positions candidates vector
%
% Rute Almeida  14.APR.2005
% Last update: 26JUL2011
%
% MATLAB Version R13
if nargin<6
    messages.status=1;
end
if ~isfield(messages,'warnings'), messages.warnings=[]; end
spf=heasig.spf_ecg; %24.MAR.09
if nargin==4
    tol=ceil(0.275/2*heasig.freq*spf); %24.MAR.09
elseif nargin<4
    messages.errors=[messages.errors {'Fatal error in qrscandidatesnew: not enough inputs.'}];
    warning(char(messages.errors(end)))
    messages.errors_desc=[messages.errors_desc 'Mandatory inputs not defined.'];
    messages.status=0;
    return
end
indexes=[];
indexes1=1:length(timeqrs1);
indexes2=1:length(timeqrs2);
indexes3=1:length(timeqrs3);
timeqrs=[];
lengthaux=([length(timeqrs1) length(timeqrs2) length(timeqrs3)]);

for g=1:min(lengthaux)
    [aux auxi]=sort([timeqrs1(g) timeqrs2(g) timeqrs3(g)]);
    if aux(1)>(aux(2)-tol-1)
        if (aux(2)+tol+1)>aux(3)
            timeqrs=[timeqrs [timeqrs1(g);timeqrs2(g);timeqrs3(g)]]; %#ok<AGROW>
            indexes=[indexes [indexes1(g);indexes2(g);indexes3(g)]]; %#ok<AGROW>
        else
            if auxi(3)==1
                timeqrs1(g+1:end+1)=timeqrs1(g:end);
                indexes1(g+1:end+1)= indexes1(g:end);
                indexes1(g)=NaN;
                timeqrs1(g)=NaN;
            else if auxi(3)==2
                    timeqrs2(g+1:end+1)=timeqrs2(g:end);
                    indexes2(g+1:end+1)= indexes2(g:end);
                    timeqrs2(g)=NaN;
                    indexes2(g)=NaN;
                else
                    timeqrs3(g+1:end+1)=timeqrs3(g:end);
                    indexes3(g+1:end+1)= indexes3(g:end);
                    timeqrs3(g)=NaN;
                    indexes3(g)=NaN;
                end
            end
            timeqrs=[timeqrs [timeqrs1(g);timeqrs2(g);timeqrs3(g)]]; %#ok<AGROW>
            indexes=[indexes [indexes1(g);indexes2(g);indexes3(g)]]; %#ok<AGROW>
        end
    else

        if auxi(1)==1
            g1=g+1;
            flag=0;
            while  g1<=length(timeqrs1) && flag==0 && timeqrs1(g1)<(aux(2)+tol+1)
                if (timeqrs1(g1)>(aux(2)-tol-1))
                    flag=1;
                    timeqrs1(g:end)=[timeqrs1(g1:end) NaN*ones(1,(g1-g))];
                    indexes1(g:end)=[indexes1(g1:end) NaN*ones(1,(g1-g))];
                end
                g1=g1+1;
            end
            if flag==0
                if g1>length(timeqrs1)
                    timeqrs1(g:end)= NaN*ones(1,length(timeqrs1(g:end)));
                    indexes1(g:end)= NaN*ones(1,length(timeqrs1(g:end)));
                else
                    timeqrs1(g:end)=[timeqrs1((g1-1):end) NaN*ones(1,(g1-1-g))];
                    indexes1(g:end)=[indexes1((g1-1):end) NaN*ones(1,(g1-1-g))];
                    timeqrs1(g)=NaN;
                    indexes1(g)=NaN;
                end
            end
        elseif auxi(1)==2
            g1=g+1;
            flag=0;
            while  g1<=length(timeqrs2) && flag==0 && timeqrs2(g1)<(aux(2)+tol+1)
                if (timeqrs2(g1)>(aux(2)-tol-1))
                    flag=1;
                    timeqrs2(g:end)=[timeqrs2(g1:end) NaN*ones(1,(g1-g))];
                    indexes2(g:end)=[indexes2(g1:end) NaN*ones(1,(g1-g))];
                end
                g1=g1+1;
            end
            if flag==0
                if g1>length(timeqrs1)
                    timeqrs2(g:end)= NaN*ones(1,length(timeqrs2(g:end)));
                    indexes2(g:end)= NaN*ones(1,length(timeqrs2(g:end)));
                else
                    timeqrs2(g:end)=[timeqrs2((g1-1):end) NaN*ones(1,(g1-1-g))];
                    indexes2(g:end)=[indexes2((g1-1):end) NaN*ones(1,(g1-1-g))];
                    timeqrs2(g)=NaN;
                    indexes2(g)=NaN;
                end
            end
        else
            g1=g+1;
            flag=0;
            while  g1<=length(timeqrs3) && flag==0 && timeqrs3(g1)<(aux(2)+tol+1)
                if (timeqrs3(g1)>(aux(2)-tol-1))
                    flag=1;
                    timeqrs3(g:end)=[timeqrs3(g1:end) NaN*ones(1,(g1-g))];
                    indexes3(g:end)=[indexes3(g1:end) NaN*ones(1,(g1-g))];
                end
                g1=g1+1;
            end
            if flag==0
                if g1>length(timeqrs1)
                    timeqrs3(g:end)= NaN*ones(1,length(timeqrs3(g:end)));
                    indexes3(g:end)= NaN*ones(1,length(timeqrs3(g:end)));
                else
                    timeqrs3(g:end)=[timeqrs3((g1-1):end) NaN*ones(1,(g1-1-g))];
                    indexes3(g:end)=[indexes3((g1-1):end) NaN*ones(1,(g1-1-g))];
                    timeqrs3(g)=NaN;
                    indexes3(g)=NaN;
                end
            end
        end

        if aux(3)-aux(2)>(tol+1)
            if auxi(3)==1
                timeqrs1(g+1:end+1)=timeqrs1(g:end);
                indexes1(g+1:end+1)= indexes1(g:end);
                indexes1(g)=NaN;
                timeqrs1(g)=NaN;
            else if auxi(3)==2
                    timeqrs2(g+1:end+1)=timeqrs2(g:end);
                    indexes2(g+1:end+1)= indexes2(g:end);
                    timeqrs2(g)=NaN;
                    indexes2(g)=NaN;
                else
                    timeqrs3(g+1:end+1)=timeqrs3(g:end);
                    indexes3(g+1:end+1)= indexes3(g:end);
                    timeqrs3(g)=NaN;
                    indexes3(g)=NaN;
                end
            end
        else
            flag=1;
        end
        if flag==1
            timeqrs=[timeqrs [timeqrs1(g);timeqrs2(g);timeqrs3(g)]]; %#ok<AGROW>
            indexes=[indexes [indexes1(g);indexes2(g);indexes3(g)]]; %#ok<AGROW>
        else
            timeqrs1(g)=NaN;indexes1(g)=NaN;
            timeqrs2(g)=NaN;indexes2(g)=NaN;
            timeqrs3(g)=NaN;indexes3(g)=NaN;
        end
    end
end
lengthaux=([length(timeqrs1) length(timeqrs2) length(timeqrs3)]);
if isempty(g)
    g=0;
end
aux1=find(lengthaux-g >0);
%remains beats in all 3 leads
while length(aux1)==3
%     intervalaux=(g+1):(min(lengthaux(aux1)));
    for g=(g+1):(min(lengthaux(aux1)))
        [aux auxi]=sort([timeqrs1(g) timeqrs2(g) timeqrs3(g)]);
        if aux(1)>(aux(2)-tol-1)
            if (aux(2)+tol+1)>aux(3)
                timeqrs=[timeqrs [timeqrs1(g);timeqrs2(g);timeqrs3(g)]];  %#ok<AGROW>
                indexes=[indexes [indexes1(g);indexes2(g);indexes3(g)]];  %#ok<AGROW>
            else
                if auxi(3)==1
                    timeqrs1(g+1:end+1)=timeqrs1(g:end);
                    indexes1(g+1:end+1)= indexes1(g:end);
                    indexes1(g)=NaN;
                    timeqrs1(g)=NaN;
                else if auxi(3)==2
                        timeqrs2(g+1:end+1)=timeqrs2(g:end);
                        indexes2(g+1:end+1)= indexes2(g:end);
                        timeqrs2(g)=NaN;
                        indexes2(g)=NaN;
                    else
                        timeqrs3(g+1:end+1)=timeqrs3(g:end);
                        indexes3(g+1:end+1)= indexes3(g:end);
                        timeqrs3(g)=NaN;
                        indexes3(g)=NaN;
                    end
                end
                timeqrs=[timeqrs [timeqrs1(g);timeqrs2(g);timeqrs3(g)]];  %#ok<AGROW>
                indexes=[indexes [indexes1(g);indexes2(g);indexes3(g)]];  %#ok<AGROW>
            end
        else
            if auxi(1)==1
                g1=g+1;
                flag=0;
                while  g1<=length(timeqrs1) && flag==0 && timeqrs1(g1)<(aux(2)+tol+1)
                    if (timeqrs1(g1)>(aux(2)-tol-1))
                        flag=1;
                        timeqrs1(g:end)=[timeqrs1(g1:end) NaN*ones(1,(g1-g))];
                        indexes1(g:end)=[indexes1(g1:end) NaN*ones(1,(g1-g))];
                    end
                    g1=g1+1;
                end
                if flag==0
                    if g1>length(timeqrs1)
                        timeqrs1(g:end)=NaN*ones(1,length(timeqrs1(g:end)));
                        indexes1(g:end)=NaN*ones(1,length(timeqrs1(g:end)));
                    else
                        timeqrs1(g:end)=[timeqrs1((g1-1):end) NaN*ones(1,(g1-1-g))];
                        indexes1(g:end)=[indexes1((g1-1):end) NaN*ones(1,(g1-1-g))];
                        timeqrs1(g)=NaN;
                        indexes1(g)=NaN;
                    end
                end
            elseif auxi(1)==2
                g1=g+1;
                flag=0;
                while  g1<=length(timeqrs2) && flag==0 && timeqrs2(g1)<(aux(2)+tol+1)
                    if (timeqrs2(g1)>(aux(2)-tol-1))
                        flag=1;
                        timeqrs2(g:end)=[timeqrs2(g1:end) NaN*ones(1,(g1-g))];
                        indexes2(g:end)=[indexes2(g1:end) NaN*ones(1,(g1-g))];
                    end
                    g1=g1+1;
                end
                if flag==0
                    if g1>length(timeqrs1)
                        timeqrs2(g:end)=NaN*ones(1,length(timeqrs2(g:end)));
                        indexes2(g:end)=NaN*ones(1,length(timeqrs2(g:end)));
                    else
                        timeqrs2(g:end)=[timeqrs2((g1-1):end) NaN*ones(1,(g1-1-g))];
                        indexes2(g:end)=[indexes2((g1-1):end) NaN*ones(1,(g1-1-g))];
                        timeqrs2(g)=NaN;
                        indexes2(g)=NaN;
                    end
                end
            else
                g1=g+1;
                flag=0;
                while  g1<=length(timeqrs3) && flag==0 && timeqrs3(g1)<(aux(2)+tol+1)
                    if (timeqrs3(g1)>(aux(2)-tol-1))
                        flag=1;
                        timeqrs3(g:end)=[timeqrs3(g1:end) NaN*ones(1,(g1-g))];
                        indexes3(g:end)=[indexes3(g1:end) NaN*ones(1,(g1-g))];
                    end
                    g1=g1+1;
                end
                if flag==0
                    if g1>length(timeqrs1)
                        timeqrs3(g:end)=NaN*ones(1,length(timeqrs3(g:end)));
                        indexes3(g:end)=NaN*ones(1,length(timeqrs3(g:end)));
                    else
                        timeqrs3(g:end)=[timeqrs3((g1-1):end) NaN*ones(1,(g1-1-g))];
                        indexes3(g:end)=[indexes3((g1-1):end) NaN*ones(1,(g1-1-g))];
                        timeqrs3(g)=NaN;
                        indexes3(g)=NaN;
                    end
                end
            end
            if aux(3)-aux(2)>(tol+1)
                if auxi(3)==1
                    timeqrs1(g+1:end+1)=timeqrs1(g:end);
                    indexes1(g+1:end+1)= indexes1(g:end);
                    indexes1(g)=NaN;
                    timeqrs1(g)=NaN;
                else if auxi(3)==2
                        timeqrs2(g+1:end+1)=timeqrs2(g:end);
                        indexes2(g+1:end+1)= indexes2(g:end);
                        timeqrs2(g)=NaN;
                        indexes2(g)=NaN;
                    else
                        timeqrs3(g+1:end+1)=timeqrs3(g:end);
                        indexes3(g+1:end+1)= indexes3(g:end);
                        timeqrs3(g)=NaN;
                        indexes3(g)=NaN;
                    end
                end
            else
                flag=1;
            end
            if flag==1
                timeqrs=[timeqrs [timeqrs1(g);timeqrs2(g);timeqrs3(g)]];  %#ok<AGROW>
                indexes=[indexes [indexes1(g);indexes2(g);indexes3(g)]];  %#ok<AGROW>
            else
                timeqrs1(g)=NaN;indexes1(g)=NaN;
                timeqrs2(g)=NaN;indexes2(g)=NaN;
                timeqrs3(g)=NaN;indexes3(g)=NaN;
            end
        end
    end
    lengthaux=([length(timeqrs1) length(timeqrs2) length(timeqrs3)]);
    aux1=find(lengthaux-g >0);
end
[min1i,min1]=min(find(~isnan(timeqrs1(end:-1:1)))); %#ok<NASGU,MXFND>
[min2i,min3]=min(find(~isnan(timeqrs2(end:-1:1)))); %#ok<NASGU,MXFND>
[min3i,min3]=min(find(~isnan(timeqrs3(end:-1:1)))); %#ok<NASGU,MXFND>
timeqrs1(length(timeqrs1)-min1i+2:end)=[];
timeqrs2(length(timeqrs2)-min2i+2:end)=[];
timeqrs3(length(timeqrs3)-min3i+2:end)=[];
lengthaux=([length(timeqrs1) length(timeqrs2) length(timeqrs3)]);
aux1=find(lengthaux-g >0);
%remaining beats in 2 leads
if length(aux1)>1
    for intervalauxi=(g+1):(min(lengthaux(aux1)));
        if aux1(1)==1
            if aux1(2)==2
                [aux auxi]=sort([timeqrs1(intervalauxi) timeqrs2(intervalauxi)]);
                if aux(1)>(aux(2)-tol-1)
                    timeqrs=[timeqrs [timeqrs1(intervalauxi);timeqrs2(intervalauxi);NaN]];  %#ok<AGROW>
                    indexes=[indexes [indexes1(intervalauxi);indexes2(intervalauxi);NaN]];  %#ok<AGROW>
                else
                    if auxi(2)==2
                        timeqrs2(intervalauxi+1:end+1)=timeqrs2(intervalauxi:end);
                        indexes2(intervalauxi+1:end+1)= indexes2(intervalauxi:end);
                        timeqrs2(intervalauxi)=NaN;
                        indexes2(intervalauxi)=NaN;
                    elseif auxi(2)==1
                        timeqrs1(intervalauxi+1:end+1)=timeqrs1(intervalauxi:end);
                        indexes1(intervalauxi+1:end+1)= indexes1(intervalauxi:end);
                        timeqrs1(intervalauxi)=NaN;
                        indexes1(intervalauxi)=NaN;
                    end
                end
            elseif aux1(2)==3
                [aux auxi]=sort([timeqrs1(intervalauxi) timeqrs3(intervalauxi)]);
                if aux(1)>(aux(2)-tol-1)
                    timeqrs=[timeqrs [timeqrs1(intervalauxi);NaN;timeqrs3(intervalauxi)]];  %#ok<AGROW>
                    indexes=[indexes [indexes1(intervalauxi);NaN;indexes3(intervalauxi)]];  %#ok<AGROW>
                else
                    if auxi(2)==2 %25SET08
                        timeqrs3(intervalauxi+1:end+1)=timeqrs3(intervalauxi:end);
                        indexes3(intervalauxi+1:end+1)= indexes3(intervalauxi:end);
                        timeqrs3(intervalauxi)=NaN;
                        indexes3(intervalauxi)=NaN;
                    elseif auxi(2)==1
                        timeqrs1(intervalauxi+1:end+1)=timeqrs1(intervalauxi:end);
                        indexes1(intervalauxi+1:end+1)= indexes1(intervalauxi:end);
                        timeqrs1(intervalauxi)=NaN;
                        indexes1(intervalauxi)=NaN;
                    end
                end
            end
        elseif aux1(1)==2
            [aux auxi]=sort([timeqrs2(intervalauxi) timeqrs3(intervalauxi)]);
            if aux(1)>(aux(2)-tol-1)
                timeqrs=[timeqrs [NaN;timeqrs2(intervalauxi);timeqrs3(intervalauxi)]];  %#ok<AGROW>
                indexes=[indexes [NaN;indexes2(intervalauxi);indexes3(intervalauxi)]];  %#ok<AGROW>
            else
                if auxi(2)==1%25SET08
                    timeqrs2(intervalauxi+1:end+1)=timeqrs2(intervalauxi:end);
                    indexes2(intervalauxi+1:end+1)= indexes2(intervalauxi:end);
                    timeqrs2(intervalauxi)=NaN;
                    indexes2(intervalauxi)=NaN;
                elseif auxi(2)==2%25SET08
                    timeqrs3(intervalauxi+1:end+1)=timeqrs3(intervalauxi:end);
                    indexes3(intervalauxi+1:end+1)= indexes3(intervalauxi:end);
                    timeqrs3(intervalauxi)=NaN;
                    indexes3(intervalauxi)=NaN;
                end
            end
        end
    end
    lengthaux=([length(timeqrs1) length(timeqrs2) length(timeqrs3)]);
    if isempty(intervalauxi);
        intervalauxi=g;
    end
    aux1=find(lengthaux-intervalauxi >0);
    %remains beats in 2 leads
    while length(aux1)==2
        intervalauxi=intervalauxi+1;
        if aux1(1)==1
            if aux1(2)==2
                [aux auxi]=sort([timeqrs1(intervalauxi) timeqrs2(intervalauxi)]);
                if aux(1)>(aux(2)-tol-1)
                    timeqrs=[timeqrs [timeqrs1(intervalauxi);timeqrs2(intervalauxi);NaN]];  %#ok<AGROW>
                    indexes=[indexes [indexes1(intervalauxi);indexes2(intervalauxi);NaN]];  %#ok<AGROW>
                else
                    if auxi(2)==2
                        timeqrs2(intervalauxi+1:end+1)=timeqrs2(intervalauxi:end);
                        indexes2(intervalauxi+1:end+1)= indexes2(intervalauxi:end);
                        timeqrs2(intervalauxi)=NaN;
                        indexes2(intervalauxi)=NaN;
                    elseif auxi(2)==1
                        timeqrs1(intervalauxi+1:end+1)=timeqrs1(intervalauxi:end);
                        indexes1(intervalauxi+1:end+1)= indexes1(intervalauxi:end);
                        timeqrs1(intervalauxi)=NaN;
                        indexes1(intervalauxi)=NaN;
                    end
                end
            elseif aux1(2)==3
                [aux auxi]=sort([timeqrs1(intervalauxi) timeqrs3(intervalauxi)]);
                if aux(1)>(aux(2)-tol-1)
                    timeqrs=[timeqrs [timeqrs1(intervalauxi);NaN;timeqrs3(intervalauxi)]];  %#ok<AGROW>
                    indexes=[indexes [indexes1(intervalauxi);NaN;indexes3(intervalauxi)]];  %#ok<AGROW>
                else
                    if auxi(2)==3
                        timeqrs3(intervalauxi+1:end+1)=timeqrs3(intervalauxi:end);
                        indexes3(intervalauxi+1:end+1)= indexes3(intervalauxi:end);
                        timeqrs3(intervalauxi)=NaN;
                        indexes3(intervalauxi)=NaN;
                    elseif auxi(2)==1
                        timeqrs1(intervalauxi+1:end+1)=timeqrs1(intervalauxi:end);
                        indexes1(intervalauxi+1:end+1)= indexes1(intervalauxi:end);
                        timeqrs1(intervalauxi)=NaN;
                        indexes1(intervalauxi)=NaN;
                    end
                end
            end
        elseif aux1(1)==2
            [aux auxi]=sort([timeqrs2(intervalauxi) timeqrs3(intervalauxi)]);
            if aux(1)>(aux(2)-tol-1)
                timeqrs=[timeqrs [NaN; timeqrs2(intervalauxi);timeqrs3(intervalauxi)]];  %#ok<AGROW>
                indexes=[indexes [NaN; indexes2(intervalauxi);indexes3(intervalauxi)]];  %#ok<AGROW>
            else
                if auxi(2)==2
                    timeqrs2(intervalauxi+1:end+1)=timeqrs2(intervalauxi:end);
                    indexes2(intervalauxi+1:end+1)= indexes2(intervalauxi:end);
                    timeqrs2(intervalauxi)=NaN;
                    indexes2(intervalauxi)=NaN;
                elseif auxi(2)==3
                    timeqrs3(intervalauxi+1:end+1)=timeqrs3(intervalauxi:end);
                    indexes3(intervalauxi+1:end+1)= indexes3(intervalauxi:end);
                    timeqrs3(intervalauxi)=NaN;
                    indexes3(intervalauxi)=NaN;
                end
            end
        end
        lengthaux=([length(timeqrs1) length(timeqrs2) length(timeqrs3)]);
        aux1=find(lengthaux-intervalauxi >0);
    end
    lengthaux=([length(timeqrs1) length(timeqrs2) length(timeqrs3)]);
    g=intervalauxi;
    aux1=find(lengthaux-g >0);
end
%NOTE THAT ONLY ALIGNED ONES ARE CONSIDERED!!!!!!!!!!!!!!!!!!
intervalaux=(g+1):(min(lengthaux(aux1)));
if length(aux1)>1
    timeqrs1((min(lengthaux(aux1)))+1:end)=[];
    timeqrs2((min(lengthaux(aux1)))+1:end)=[];
    timeqrs3((min(lengthaux(aux1)))+1:end)=[];
    indexes1((min(lengthaux(aux1)))+1:end)=[];
    indexes2((min(lengthaux(aux1)))+1:end)=[];
    indexes3((min(lengthaux(aux1)))+1:end)=[];
    if ~ismember(1,aux1)
        timeqrs1(intervalaux)=NaN*ones(size(intervalaux));
        indexes1(intervalaux)=NaN*ones(size(intervalaux));
        aux=[timeqrs2(intervalaux); timeqrs3(intervalaux) ];
    end
    if ~ismember(2,aux1)
        timeqrs2(intervalaux)=NaN*ones(size(intervalaux));
        indexes2(intervalaux)=NaN*ones(size(intervalaux));
        aux=[timeqrs1(intervalaux); timeqrs3(intervalaux) ];
    end
    if ~ismember(3,aux1)
        timeqrs3(intervalaux)=NaN*ones(size(intervalaux));
        indexes3(intervalaux)=NaN*ones(size(intervalaux));
        aux=[timeqrs1(intervalaux); timeqrs2(intervalaux) ];
    end
    timeqrs=[timeqrs [timeqrs1(intervalaux(abs(diff(aux))<(tol+1)));timeqrs2(intervalaux(abs(diff(aux))<(tol+1)));timeqrs3(intervalaux(abs(diff(aux))<(tol+1)))]];
    indexes=[indexes [indexes1(intervalaux(abs(diff(aux))<(tol+1)));indexes2(intervalaux(abs(diff(aux))<(tol+1)));indexes3(intervalaux(abs(diff(aux))<(tol+1)))]];
%     indexes1(intervalaux(abs(diff(aux))>=(tol+1) | isnan(diff(aux))))=[];
%     indexes2(intervalaux(abs(diff(aux))>=(tol+1) | isnan(diff(aux))))=[];
%     indexes3(intervalaux(abs(diff(aux))>=(tol+1) | isnan(diff(aux))))=[];
%     timeqrs1(intervalaux(abs(diff(aux))>=(tol+1) | isnan(diff(aux))))=[];
%     timeqrs2(intervalaux(abs(diff(aux))>=(tol+1) | isnan(diff(aux))))=[];
%     timeqrs3(intervalaux(abs(diff(aux))>=(tol+1) | isnan(diff(aux))))=[];
else
%     timeqrs1((min(lengthaux(aux1))):end)=[];
%     timeqrs2((min(lengthaux(aux1))):end)=[];
%     timeqrs3((min(lengthaux(aux1))):end)=[];
%     indexes1((min(lengthaux(aux1))):end)=[];
%     indexes2((min(lengthaux(aux1))):end)=[];
%     indexes3((min(lengthaux(aux1))):end)=[];
end
if ~isempty(indexes) %17ABRIL08
indexes=indexes(:,~isnan(indexes(1,:))| ~isnan(indexes(2,:))|~isnan(indexes(3,:)));
end
timeqrs(:,(size(indexes,2)+1):end)=[];

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function output = temporal(curr,prev,alpha,T)
curr=double(curr);
prev=double(prev);
diff = curr-prev;
m = sign(abs(diff)).*(abs(diff)<T);
av = alpha.*curr+(1-alpha).*prev;
output = curr;
m=m&m;
output(m)=av(m);
output=uint8(output);
end

author: zhangqianqianQQ


text:
% Codes for CVPR-15 work `Face Alignment by Coarse-to-Fine Shape Searching'
% Any question please contact Shizhan Zhu: zhshzhutah2@gmail.com
% Released on July 25, 2015

if exist('config','var'), error('Please remove variable config mannally!'); end;

addpath(genpath('./codes_release'));
addpath(genpath('./external'));

config.stageTot = 3;
config.win_size = 250;

% From Pr to sub-region center
regsInfo.iterTot = [3 3 3];
regsInfo.samplingOffset = [NaN,1500,100];
regsInfo.trainSampleTot = [10 10 10];
regsInfo.testSampleTot = [10 10 10];
regsInfo.aug_eyes_id = [37 46];
regsInfo.win_size = config.win_size;
regsInfo.SIFTscale = [12 6 6 6];
regsInfo.mirror = [17:-1:1 27:-1:18 28:31 36:-1:32 ...
    46 45 44 43 48 47 40 39 38 37 42 41 ...
    55:-1:49 60:-1:56 65:-1:61 68:-1:66];
regsInfo.mirror = [regsInfo.mirror regsInfo.mirror+68];
regsInfo.regressorInfo.trainMethod = @getLR_lcScale;
regsInfo.regressorInfo.lambda = [50*25 50*30 50*20 50*20];
regsInfo.regressorInfo.times = 20;

regsInfo.dominantIterTot = [100 100 100];

config.regs = regsInfo;

% From sub-region center to Pr
probsInfo.semantic_id = [18 22 23 27 37 40 43 46 49 55 51 52 53 58 9];
probsInfo.fix_id = [3:7 11:15];
probsInfo.sigmaCutoff = [0.25 0.25];

probsInfo.SVCradius = [20 10];
probsInfo.SVCthre = {[0.3 0.6],[0.3 0.6]};
probsInfo.probSamplings = [5 5];
probsInfo.pyramidScale = [6 12 18];
probsInfo.representativeNum1 = [12 12];
probsInfo.representativeNum2 = [8 8];

probsInfo.acceptThre = 0.9;
probsInfo.gamma_current = 1/50;
probsInfo.gamma_search = 1/50;
probsInfo.gamma_fix = 1/50;
probsInfo.gamma_prior = [1/50,1];

config.probs = probsInfo;

% Prior model
priorsInfo.augTimes = 10;
priorsInfo.maxRoll = 45;
priorsInfo.nose_id = 34;
priorsInfo.noseUpLength = 15;
priorsInfo.rotationCenter = [0.5 0.35]; %(x,y)
priorsInfo.win_size = config.win_size;
priorsInfo.rotatorLength = 15;
priorsInfo.extractWindow = [0.2 0.8 0.1 0.7]; %(xX,xD,yX,yD)
priorsInfo.extractCellNum = 3;
priorsInfo.BaggerNum = 100;

config.priors = priorsInfo;

clear regsInfo probsInfo priorsInfo;

%% Generate testConf
testConf.stageTot = config.stageTot;
testConf.win_size = config.win_size;
testConf.priors.win_size = config.win_size;
testConf.priors.rotationCenter = config.priors.rotationCenter;
testConf.priors.rotatorLength = config.priors.rotatorLength;
testConf.priors.extractWindow = config.priors.extractWindow;
testConf.priors.extractCellNum = config.priors.extractCellNum;
testConf.priors.predictedVoteThre = 10;
testConf.regs.samplingTot = [10 10 10];
testConf.regs.dominantIterTot = config.regs.dominantIterTot;
testConf.regs.SIFTscale = config.regs.SIFTscale;
testConf.regs.iterTot = config.regs.iterTot;
testConf.probs.probSamplings = config.probs.probSamplings;
testConf.probs.semantic_id = config.probs.semantic_id;
testConf.probs.fix_id = config.probs.fix_id;
testConf.probs.pyramidScale = config.probs.pyramidScale;
testConf.probs.acceptThre = config.probs.acceptThre;
testConf.probs.gamma_current = config.probs.gamma_current;
testConf.probs.gamma_search = config.probs.gamma_search;
testConf.probs.gamma_fix = config.probs.gamma_fix;
testConf.probs.gamma_prior = config.probs.gamma_prior;

author: zhusz


text:
classdef MOEADDE_F7 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP for testing MOEA/D-DE

%------------------------------- Reference --------------------------------
% H. Li and Q. Zhang, Multiobjective optimization problems with complicated
% Pareto sets, MOEA/D and NSGA-II, IEEE Transactions on Evolutionary
% Computation, 2009, 13(2): 284-302.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 10; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            J1 = 3 : 2 : obj.D;
            J2 = 2 : 2 : obj.D;
            Y  = X - repmat(X(:,1),1,obj.D).^repmat((1+3*((1:obj.D)-2)/(obj.D-2))/2,size(X,1),1);
            PopObj(:,1) = X(:,1)         + 2*mean(4*Y(:,J1).^2-cos(8*Y(:,J1)*pi)+1,2);
            PopObj(:,2) = 1-sqrt(X(:,1)) + 2*mean(4*Y(:,J2).^2-cos(8*Y(:,J2)*pi)+1,2);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - R(:,1).^0.5;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function res = bf_inverse_deflect(BF, S)
% Used DeFleCT framework to compute spatial filters.
% Copyright (C) 2013 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak, using the code from Matti Stenroos and Olaf Hauk
% http://imaging.mrc-cbu.cam.ac.uk/meg/AnalyzingData/DeFleCT_SpatialFiltering_Tools
%
% Please cite:
% Hauk O, Stenroos M.
% A framework for the design of flexible cross-talk functions for spatial filtering of EEG/MEG data: DeFleCT.
% Human Brain Mapping 2013
% $Id: bf_inverse_deflect.m 7703 2019-11-22 12:06:29Z guillaume $

%--------------------------------------------------------------------------
if nargin == 0

    list = cfg_entry;
    list.tag = 'list';
    list.name = 'List vertex indices';
    list.strtype = 'n';
    list.num = [1 Inf];
    list.help = {'Specify sources of interest by listing vertex indices.'};

    pos = cfg_entry;
    pos.tag = 'pos';
    pos.name = 'MNI coordinates';
    pos.strtype = 'r';
    pos.num = [1 3];
    pos.help = {'Locations for the VOI center in MNI coordinates'};
    pos.val = {};

    radius = cfg_entry;
    radius.tag = 'radius';
    radius.name = 'Radius';
    radius.strtype = 'r';
    radius.num = [1 1];
    radius.val = {0};
    radius.help = {'Radius (in mm) for the VOI (leave 0 for closest point)'};

    voi = cfg_branch;
    voi.tag = 'voi';
    voi.name = 'VOI';
    voi.val = {pos, radius};

    passband = cfg_repeat;
    passband.tag = 'passband';
    passband.name = 'Passband sources';
    passband.num  = [1 Inf];
    passband.values = {voi, list};

    stopband = cfg_repeat;
    stopband.tag = 'stopband';
    stopband.name = 'Stopband sources';
    stopband.num  = [0 Inf];
    stopband.values = {voi, list};

    svdpassband = cfg_entry;
    svdpassband.tag = 'svdpassband';
    svdpassband.name = 'SVD passband';
    svdpassband.strtype = 'w';
    svdpassband.num = [1 1];
    svdpassband.val = {0};
    svdpassband.help = {'Number of components to summarise the passband in.,',...
        'Leave at zero for no SVD.'};

    svdstopband = cfg_entry;
    svdstopband.tag = 'svdstopband';
    svdstopband.name = 'SVD stopband';
    svdstopband.strtype = 'w';
    svdstopband.num = [1 1];
    svdstopband.val = {0};
    svdstopband.help = {'Number of components to summarise the stopband in.,',...
        'Leave at zero for no SVD.'};

    forcepassband         = cfg_menu;
    forcepassband.tag     = 'forcepassband';
    forcepassband.name    = 'Force passband';
    forcepassband.help    = {'Forces the output for all passband components '};
    forcepassband.labels  = {'yes', 'no'};
    forcepassband.values  = {1, 0};
    forcepassband.val = {0};

    label = cfg_entry;
    label.tag = 'label';
    label.name = 'Label';
    label.strtype = 's';
    label.help = {'Label for the output source'};

    usecov         = cfg_menu;
    usecov.tag     = 'usecov';
    usecov.name    = 'Use covariance matrix';
    usecov.help    = {'Use covariance matrix for pre-whitening'};
    usecov.labels  = {'yes', 'no'};
    usecov.values  = {1, 0};
    usecov.val = {1};

    filter = cfg_branch;
    filter.tag = 'filter';
    filter.name = 'Filter';
    filter.val = {label, passband,  svdpassband, forcepassband, stopband, svdstopband, usecov};

    filters = cfg_repeat;
    filters.tag = 'filters';
    filters.name = 'Filters';
    filters.num  = [1 Inf];
    filters.values = {filter};

    snr = cfg_entry;
    snr.tag = 'snr';
    snr.name = 'SNR';
    snr.strtype = 'r';
    snr.num = [1 1];
    snr.val = {5};
    snr.help = {'The assumed ratio of variances of signal and noise,',...
        'used for setting the regularisation parameter.'};

    trunc = cfg_entry;
    trunc.tag = 'trunc';
    trunc.name = 'Truncation parameter';
    trunc.strtype = 'w';
    trunc.num = [1 1];
    trunc.val = {0};
    trunc.help = {'The number of (smallest) singular values of the covariance matrix that are set to ',...
        'zero before making the whitener. For example, if the data has been SSP-projected, it needs to be at least the number of ',...
        'components projected away.'};

    deflect     = cfg_branch;
    deflect.tag  = 'deflect';
    deflect.name = 'DeFleCT';
    deflect.val  = {filters, snr, trunc};
    deflect.help = {'DeFleCT spatial filter design framework by Matti Stenroos and Olaf Hauk',...
        'http://imaging.mrc-cbu.cam.ac.uk/meg/AnalyzingData/DeFleCT_SpatialFiltering_Tools',...
        'Please cite:'...
        'Hauk O, Stenroos M.',...
        'A framework for the design of flexible cross-talk functions for spatial filtering of EEG/MEG data: DeFleCT.',...
        'Human Brain Mapping 2013'};
    res = deflect;

    return
elseif nargin < 2
    error('Two input arguments are required');
end

C     =  BF.features.(S.modality).C;
U     =  BF.features.(S.modality).U;

L = [];
Li = {};
for i = 1:numel(S.L)
    cL    = U'*S.L{i};
    Li{i} = size(L, 2)+(1:size(cL, 2));
    L  =[L cL];
end

mnipos = spm_eeg_inv_transform_points(BF.data.transforms.toMNI,  BF.sources.pos);

nfilt = numel(S.filter);

W     = cell(1, nfilt);
label = {};

spm('Pointer', 'Watch');drawnow;
spm_progress_bar('Init', nfilt, ['Computing ' S.modality ' filters']); drawnow;
if nfilt > 100, Ibar = floor(linspace(1, nfilt,100));
else Ibar = 1:nfilt; end

for i = 1:nfilt
    filter = S.filter(i);
    if ~filter.svdpassband
        filter.svdpassband = [];
    end

     if ~filter.svdstopband
        filter.svdstopband = [];
    end

    passband = get_vertices(filter.passband, mnipos);
    stopband = get_vertices(filter.stopband, mnipos);

    passband = cat(2, Li{passband});

    stopband = cat(2, Li{stopband});

    if filter.usecov
        W{i} = DeFleCT(passband, filter.svdpassband, filter.forcepassband, stopband, filter.svdstopband,...
            L, C, S.snr, S.trunc);
    else
        W{i} = DeFleCT(passband, filter.svdpassband, filter.forcepassband, stopband, filter.svdstopband,...
            L, [], S.snr, eye(size(L, 1)));
    end

    label{i} = filter.label;

    if ismember(i, Ibar)
        spm_progress_bar('Set', i); drawnow;
    end
end

spm_progress_bar('Clear');

disp('Note: it is crucial to visualise the filter-cosstalk functions to verify correct design.');

res.W     = W;
res.label = label;

function vertices = get_vertices(band, mnipos)

vertices = [];
for j = 1:numel(band)
    switch char(fieldnames(band{j}))
        case 'voi'
            pnt = band{j}.voi.pos;
            rad = band{j}.voi.radius;

            dist = sqrt(sum((mnipos - repmat(pnt, size(mnipos, 1), 1)).^2, 2));

            if ~rad
                [mindist, ind] = min(dist);
                if mindist>30
                    error('There are no sources within 3cm of the specified location');
                end
            else
                ind = find(dist<=rad);
            end
        case 'list'
            ind = band{j}.list;
    end

    vertices = [vertices; ind(:)];
end

author: spm


text:
% SP_GET_CELLS: Compute the indices of the cells within the support of a list of tensor-product B-spline function.
%
% [cell_indices, indices_per_function] = sp_get_cells (space, msh, indices)
%
% INPUT:
%    space:   object defining the space of discrete functions (see sp_multipatch)
%    msh:     object defining the domain partition and the quadrature rule (see msh_multipatch)
%    indices: indices of the functions.
%
% OUTPUT:
%    cell_indices: indices of the cells within the support of the basis functions.
%    indices_per_function: indices of the cells within the support of each basis function.
%
% Copyright (C) 2015, 2016 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function [cell_indices, indices_per_function] = sp_get_cells (space, msh, fun_indices)

fun_indices = fun_indices(:).';

indices_per_function = cell (numel (fun_indices), 1);
cell_indices = [];

Nelem = cumsum ([0 msh.nel_per_patch]);
for iptc = 1:space.npatch
  [~, patch_indices, local_funs] = intersect (space.gnum{iptc}, fun_indices);
  if (~isempty (patch_indices))
    patch_indices = patch_indices(:).';
    [aux_cell_indices, ind_per_fun] = sp_get_cells (space.sp_patch{iptc}, msh.msh_patch{iptc}, patch_indices);

    cell_indices = union (cell_indices, Nelem(iptc)+aux_cell_indices);

    if (nargout == 2)
      local_funs = local_funs(:).';
      for ifun = 1:numel(patch_indices)
        indices_per_function{local_funs(ifun)} = union (indices_per_function{local_funs(ifun)}, Nelem(iptc)+ind_per_fun{ifun});
      end
    end
  end
end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function ppma_io_test02 ( )

%*****************************************************************************80
%
%% TEST02 tests PPMA_READ.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 May 2011
%
%  Author:
%
%    John Burkardt
%
  file_name = 'test02.ascii.ppm';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  PPMA_READ reads an ASCII PPMA file.\n' );

  ppma_write_test ( file_name );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PPMA_WRITE_TEST created data,\n' );
  fprintf ( 1, '  and wrote it to the file "%s".\n', file_name );

  [ r, g, b ] = ppma_read ( file_name );
  [ nrow, ncol ] = size ( r );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PPMA_READ has read the file "%s"\n', file_name );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of rows of data =    %d\n', nrow );
  fprintf ( 1, '  Number of columns of data = %d\n', ncol );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample values:\n' );
  fprintf ( 1, '\n' );

  for k = 1 : 10
    i = floor ( ( ( 10 - k ) * 1 + ( k - 1 ) * nrow ) / ( 10 - 1 ) );
    j = floor ( ( ( 10 - k ) * 1 + ( k - 1 ) * ncol ) / ( 10 - 1 ) );
    fprintf ( 1, '  %4d  %4d  %6d  %6d  %6d\n', i, j, r(i,j), g(i,j), b(i,j) );
  end

  return
end

author: johannesgerer


text:
function [ly] = nautmi2ly(nautmi)
% Convert length from nautical miles to light years.
% Chad A. Greene 2012
ly = nautmi*1.957565544614e-13;


text:
function epsilon = QPSK_Alamouti_eps(Mr, k, L, nc, fp, SNR_dB, spa)
% QPSK_ALAMOUTI_EPS
%
% Inputs
%   Mr:     no of rx antennas
%   k:      number of information bits
%   L:      no of coherence blocks
%   nc:     length of coherence block
%   fp:     fraction of block for training
%   SNR_dB: SNR in dB
%   spa:    flag for saddlepoint approximation
%

%% Parameters and initialisation

% Simulation parameters
RELSPREAD_MAX = 0.40;           % max relative spread of outcomes

% Parallelisation
ncores = feature('numcores');   % number of cores on current machine
nworkers = ncores;              % use all cores for the algorithm
poolobj = gcp('nocreate');      % if no pool do not create new one
if isempty(poolobj)
    parpool(nworkers); % initialise the pool
end

% System parameters
n = L * nc;             % blocklength
R = k/n;                % rate (bits)
np = floor(fp*nc);      % number of pilots

nd = nc - np;
if mod(nd,2) == 1
    np = np + 1;
end

%% Compute error probability

rho_dB = SNR_dB;
rho = 10^(rho_dB/10);

eps_trial = eps_RCUs_Alamouti_stable(R, Mr, L, np, nc, rho, spa, RELSPREAD_MAX);
epsilon = median( eps_trial );

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function  [F, G] = F_update(upd, DCMbn, imu)
% F_update: updates F and G matrices before the execution of Kalman filter.
%
% INPUT
%   upd, 1x8 vector with data from the INS.
%   DCMbn, DCM body-to-nav.
%   imu, IMU data structure.
%
% OUTPUT
%   F,  15x15 state transition matrix.
%   G,  15x12 control-input matrix.
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% References:
%
%   Groves, P.D. (2013), Principles of GNSS, Inertial, and
% Multisensor Integrated Navigation Systems (2nd Ed.). Artech House.
% Matrix F from Eq. 14.63.
%
% 	Farrell, J. (2008). Aided Navigation: GPS With High Rate
% Sensors. McGraw-Hill Professional, USA. Matrix G from Eq. 11.108.
%
% Version: 010
% Date:    2022/03/06
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

Vn =  upd(1);
Ve =  upd(2);
Vd =  upd(3);
lat = upd(4);
h  =  upd(5);

fn =  upd(6:8);
% wn =  upd(9:11);

Om = 7.292115e-5;
I = eye(3);
Z = zeros(3);

[RM,RN] = radius(lat);

RO = sqrt(RN*RM) + h;

%% ATTITUDE MATRICES

a11 = 0;
a12 = -( (Om * sin(lat)) + (Ve / RO * tan(lat)) );
a13 = Vn / RO;
a21 = (Om * sin(lat)) + (Ve / RO * tan(lat));
a22 = 0 ;
a23 = (Om * cos(lat)) + (Ve / RO) ;
a31 = -Vn / RO;
a32 = -Om * cos(lat) - (Ve / RO);
a33 = 0;
F11 = [a11 a12 a13; a21 a22 a23; a31 a32 a33;];

% Groves, 14.64
% F11 = skewm(wn);

a11 = 0;
a12 = 1 / RO;
a13 = 0;
a21 = -1 / RO;
a22 = 0;
a23 = 0;
a31 = 0;
a32 = -tan(lat) / RO;
a33 = 0;
F12 = [a11 a12 a13; a21 a22 a23; a31 a32 a33;];

a11 = -Om * sin(lat);
a12 = 0;
a13 = -Ve / (RO^2);
a21 = 0 ;
a22 = 0 ;
a23 = Vn / (RO^2);
a31 =  -Om * cos(lat) - (Ve / ((RO) * (cos(lat))^2));
a32 = 0 ;
a33 = (Ve * tan(lat)) / (RO^2) ;
F13 = [a11 a12 a13; a21 a22 a23; a31 a32 a33;];

%% VELOCITY MATRICES

F21 = skewm(fn);

a11 = Vd / RO;
a12 = -2 * ((Om * sin(lat)) + ((Ve / RO) * tan(lat))) ;
a13 = Vn / RO ;
a21 = (2 * Om * sin(lat)) + ( (Ve / RO) * tan(lat) );
a22 = (1 / RO) * ((Vn * tan(lat)) + Vd) ;
a23 = 2 * Om * cos(lat) + (Ve / RO);
a31 = (-2 * Vn) / RO;
a32 = -2 * (Om * cos(lat) +  (Ve / RO)) ;
a33 = 0;
F22 = [a11 a12 a13; a21 a22 a23; a31 a32 a33;];

e = 0.0818191908425;        % WGS84 eccentricity
res = RN * sqrt( cos(lat)^2 + (1-e^2)^2 * sin(lat)^2);
g = gravity(lat,h);
g0 = g(3);

a11 = -Ve * ((2 * Om * cos(lat)) + (Ve / (RO * (cos(lat))^2)));
a12 = 0 ;
a13 = (1 / RO^2) * ( (Ve^2 * tan(lat)) - (Vn * Vd) );
a21 = 2 * Om * ( (Vn * cos(lat)) - (Vd * sin(lat)) ) + ( (Vn * Ve) / (RO * (cos(lat))^2) ) ;
a22 = 0 ;
a23 = -(Ve / RO^2) * (Vn * tan(lat) + Vd);
a31 = 2 * Om * Ve * sin(lat);
a32 = 0;
% a33 = (1 / RO^2) * (Vn^2 + Ve^2);
a33 = Ve^2 / (RN+h)^2 + Vn^2 / (RM+h)^2 - 2 * g0 / res;
F23 = [a11 a12 a13; a21 a22 a23; a31 a32 a33;];

%% POSITIONING MATRICES

F31 = zeros(3);

a11 = 1 / RO;
a12 = 0;
a13 = 0;
a21 = 0;
a22 = 1 / (RO * cos(lat));
a23 = 0;
a31 = 0;
a32 = 0;
a33 = -1;
F32 = [a11 a12 a13; a21 a22 a23; a31 a32 a33;];

a11 = 0;
a12 = 0;
a13 = -Vn / RO^2;
a21 = (Ve * tan(lat)) / (RO * cos(lat));
a22 = 0;
a23 = -Ve / (RO^2 * cos(lat));
a31 = 0;
a32 = 0;
a33 = 0;
F33 = [a11 a12 a13; a21 a22 a23; a31 a32 a33;];

Fbg = I;
Fba = I;

if (isinf(imu.gb_corr))
    Fgg = Z;
else
    Fgg = -diag( 1./ imu.gb_corr);
    %     Fbg = -diag(sqrt (2 ./ imu.gb_corr .* imu.gb_dyn.^2));
end

if (isinf(imu.ab_corr))
    Faa = Z;
else
    Faa = -diag(1 ./ imu.ab_corr);
    %     Fba = -diag(sqrt (2 ./ imu.ab_corr .* imu.ab_dyn.^2));
end

% Eq. 14.63 from Groves
F = [F11 F12 F13 DCMbn Z  ;
    F21  F22 F23 Z     DCMbn  ;
    F31  F32 F33 Z     Z      ;
    Z    Z   Z   Fgg   Z      ;
    Z    Z   Z   Z     Faa    ;
    ];

% Eq. 11.108 from Farrell
G = [DCMbn Z     Z   Z ;
    Z      DCMbn Z   Z ;
    Z      Z     Z   Z ;
    Z      Z     Fbg Z ;
    Z      Z     Z   Fba ;
    ];

end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function [K] = spm_sptop(sigma,q,c)
% Sparse Toeplitz convolution matrix given convolution kernel
% FORMAT [K] = spm_sptop(sigma,q,c)
%
% sigma - of Gaussian kernel K (or kernel itself)
% q     - order of matrix
% c     - kernel index at t = 0 {default c = length(sigma)/2)
% K     - q x q sparse convolution matrix
%_______________________________________________________________________
%
% Returns a q x q sparse convolution matrix.  If sigma is a scalar then
% a symmetrical Gaussian convolution matrix is returned with kernel width
% = sigma.  If sigma is a vector than sigma constitutes the kernel.  To
% obtain an assymmetrical convolution matrix (i.e. implement a phase shift
% set c = 1.
%
% Boundary handling: The row-wise sum of K is set to unity (kernel truncation)
%
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_sptop.m 1143 2008-02-07 19:33:33Z spm $

% if sigma = 0, return identity matrix; if q = 1, return 1.
%-----------------------------------------------------------------------
if ~any(sigma); K = speye(q); return; end
if q == 1;      K = 1;        return; end

% otherwise get kernel function
%-----------------------------------------------------------------------
if length(sigma) == 1
    E  = ceil(3*sigma);
    x  = [-E:E];
    k  = exp(-x.^2/(2*sigma^2));
else
    if nargin < 3
        c = length(sigma)/2;
    end
    E  = length(sigma);
    x  = [1:E] - fix(c);
    k  = sigma;
end

% and create convolution matrix
%-----------------------------------------------------------------------
K  = k(:)*ones(1,q);
j  = ones(length(k),1)*[1:q];
i  = x(:)*ones(1,q) + j;

% setting the row-wise sum to unity
%-----------------------------------------------------------------------
Q  = find((i >= 1) & (i <= q));
K  = sparse(i(Q),j(Q),K(Q));
Q  = sum(K');
K  = inv(diag(Q + (~Q)))*K;


text:
function area = triangle_area ( t )

%*****************************************************************************80
%
%% TRIANGLE_AREA returns the area of a triangle.
%
%  Location:
%
%    http://people.sc.fsu.edu/~jburkardt/m_src/triangle_integrals/triangle_area.m
%
%  Discussion:
%
%    If the vertices are given in counter clockwise order, the area
%    will be positive.
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    18 April 2015
%
%  Author:
%
%    John Burkardt.
%
%  Parameters:
%
%    Input, real T(2,3), the vertices of the triangle.
%
%    Output, real AREA, the area of the triangle.
%
  area = 0.5 * ...
    ( ...
        ( t(1,2) - t(1,1) ) * ( t(2,3) - t(2,1) ) ...
      - ( t(1,3) - t(1,1) ) * ( t(2,2) - t(2,1) ) ...
    );

  return
end


text:
function Pr = gbm_upcrossing_prob(X0, K, T, mu, sigma)
%gbm_upcrossing_prob - Compute probability of an upcrossing of a geometric Brownian motion.
%
%	Pr = gbm_upcrossing_prob(X0, K, T, mu, sigma);
%
% Inputs:
%	X0 - Almost surely known initial value of geometric Brownian motion process [scalar].
%	K - Level for upcrossing of geometric Brownian motion during a specified period [scalar].
%	T - Duration of the period of observation in a specified unit of time [scalar].
%	mu - Drift of geometric Brownian motion in rate of change per specified unit of time [scalar].
%	sigma - Diffusion of geometric Brownian motion in rate of change per specified unit of time
%		[scalar].
%
% Outputs:
%	Pr - Probability of an upcrossing of a geometric Brownian motion during the period T in the
%		specified units of time [scalar].
%
% Comments:
%	1) This function computes the probability that a geometric Brownian motion crosses at or above a
%		specified level at any time during a specified fixed period of time.
%	2) The term "unit of time" indicates that all times and rates must be in the same units of time.

% Copyright (C) 2012 The MathWorks, Inc.

if K <= X0
	Pr = 1;
else
	d0 = log(X0/K);
	d1 = (d0 + (mu - 0.5*sigma^2)*T)/(sigma*sqrt(2*T));
	d2 = (d0 - (mu - 0.5*sigma^2)*T)/(sigma*sqrt(2*T));
	d3 = -0*(mu - 0.5*sigma^2)/(0.5*sigma^2);

	Pr = 0.5*(1 + erf(d1)) + 0.5*exp(-d3)*(1 + erf(d2));
end

author: Sable


text:
function value = r4_besi0e ( x )

%*****************************************************************************80
%
%% R4_BESI0E evaluates the exponentially scaled Bessel function I0(X).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the exponentially scaled Bessel function I0(X).
%
  persistent ai02cs
  persistent ai0cs
  persistent bi0cs
  persistent ntai0
  persistent ntai02
  persistent nti0
  persistent xsml

  if ( isempty ( nti0 ) )
    ai02cs = [ ...
      0.05449041101410882E+00, ...
      0.00336911647825569E+00, ...
      0.00006889758346918E+00, ...
      0.00000289137052082E+00, ...
      0.00000020489185893E+00, ...
      0.00000002266668991E+00, ...
      0.00000000339623203E+00, ...
      0.00000000049406022E+00, ...
      0.00000000001188914E+00, ...
     -0.00000000003149915E+00, ...
     -0.00000000001321580E+00, ...
     -0.00000000000179419E+00, ...
      0.00000000000071801E+00, ...
      0.00000000000038529E+00, ...
      0.00000000000001539E+00, ...
     -0.00000000000004151E+00, ...
     -0.00000000000000954E+00, ...
      0.00000000000000382E+00, ...
      0.00000000000000176E+00, ...
     -0.00000000000000034E+00, ...
     -0.00000000000000027E+00, ...
      0.00000000000000003E+00 ]';
    ai0cs = [ ...
      0.07575994494023796E+00, ...
      0.00759138081082334E+00, ...
      0.00041531313389237E+00, ...
      0.00001070076463439E+00, ...
     -0.00000790117997921E+00, ...
     -0.00000078261435014E+00, ...
      0.00000027838499429E+00, ...
      0.00000000825247260E+00, ...
     -0.00000001204463945E+00, ...
      0.00000000155964859E+00, ...
      0.00000000022925563E+00, ...
     -0.00000000011916228E+00, ...
      0.00000000001757854E+00, ...
      0.00000000000112822E+00, ...
     -0.00000000000114684E+00, ...
      0.00000000000027155E+00, ...
     -0.00000000000002415E+00, ...
     -0.00000000000000608E+00, ...
      0.00000000000000314E+00, ...
     -0.00000000000000071E+00, ...
      0.00000000000000007E+00 ]';
    bi0cs = [ ...
     -0.07660547252839144951E+00, ...
      1.927337953993808270E+00, ...
      0.2282644586920301339E+00, ...
      0.01304891466707290428E+00, ...
      0.00043442709008164874E+00, ...
      0.00000942265768600193E+00, ...
      0.00000014340062895106E+00, ...
      0.00000000161384906966E+00, ...
      0.00000000001396650044E+00, ...
      0.00000000000009579451E+00, ...
      0.00000000000000053339E+00, ...
      0.00000000000000000245E+00 ]';
    nti0 = r4_inits ( bi0cs, 12, 0.1 * r4_mach ( 3 ) );
    ntai0 = r4_inits ( ai0cs, 21, 0.1 * r4_mach ( 3 ) );
    ntai02 = r4_inits ( ai02cs, 22, 0.1 * r4_mach ( 3 ) );
    xsml = sqrt ( 4.0 * r4_mach ( 3 ) );
  end

  y = abs ( x );

  if ( y <= xsml )
    value = 1.0;
  elseif ( y <= 3.0 )
    value = exp ( - y ) * ( 2.75 + ...
      r4_csevl ( y * y / 4.5 - 1.0, bi0cs, nti0 ) );
  elseif ( y <= 8.0 )
    value = ( 0.375 + r4_csevl ...
      ( ( 48.0 / y - 11.0 ) / 5.0, ai0cs, ntai0 ) ) / sqrt ( y );
  else
    value = ( 0.375 + r4_csevl ...
      ( 16.0 / y - 1.0, ai02cs, ntai02 ) ) / sqrt ( y );
  end

  return
end