text:
function num = cpu_cores()
    num = feature('numcores');
end

repo: R-FCN


text:
close all;
clear all;
clc;
rng('default');
sample_no = 101;
filename = sprintf('sample_%d.mat', sample_no);
load(filename);
Phi = spx.dict.MatrixOperator(PhiMtx);

% Solve the sparse recovery problem using OMP
matching_mode = 4;
options.norm_factor = 2;
options.VERBOSE = true;
result = ar_omp(PhiMtx, K, y, matching_mode, options);
% Solution vector
z = result.z;
stats = spx.commons.sparse.recovery_performance(Phi, K, y, x, z);
spx.commons.sparse.print_recovery_performance(stats);
fprintf('total_matched_atoms_count: %d \n',result.total_matched_atoms_count);
fprintf('Average atom index: %.2f \n',result.atom_index_average);
fprintf('\n\n\n\n');

% Solve the sparse recovery problem using OMP
matching_mode = 2;
options.VERBOSE = true;
result = ar_omp(PhiMtx, K, y, matching_mode, options);
% Solution vector
z = result.z;
stats = spx.commons.sparse.recovery_performance(Phi, K, y, x, z);
spx.commons.sparse.print_recovery_performance(stats);
fprintf('total_matched_atoms_count: %d \n',result.total_matched_atoms_count);
fprintf('Average atom index: %.2f \n',result.atom_index_average);

if 0
mf = spx.graphics.Figures();
mf.new_figure('OMP solution');
subplot(411);
stem(x, '.');
title('Sparse vector');
subplot(412);
stem(z, '.');
title('Recovered sparse vector');
subplot(413);
stem(abs(x - z), '.');
title('Recovery error');
subplot(414);
stem(y, '.');
title('Measurement vector');
end

author: indigits

repo: sparse-plex


text:
function varargout = flipud(varargin)
%FLIPUD   Flip/reverse a CHEBFUN2 in the y-direction.
%   G = FLIPUD(F) returns a CHEBFUN2 G with the same domain as F but reversed;
%   that is, G(x,y) = F(x, c+d-y), where the domain is [a, b, c, d].
%
% See also FLIPLR, FLIPDIM.

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

[varargout{1:nargout}] = flipud@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% Chapter 14 - Poincare Maps and Nonautonomous Systems in the Plane.
% Programs 14a - Solving an initial value problem.
% Copyright Birkhauser 2013. Stephen Lynch.

% Solve a differential equation (Example 1).
r=dsolve('Dr=-r^2','r(0)=1');

% List the first eight returns on the segment {y=0, 0<x<1}.
% There may be small inaccuracies due to the numerical solution.
deq=inline('[-(r(1))^2]','t','r');
options=odeset('RelTol',1e-6,'AbsTol',1e-6);
[t,returns]=ode45(deq,0:2*pi:16*pi,1);
returns

% End of Programs 14a.

author: Sable

repo: mcbench-benchmarks


text:
function [sift, SIFTparam] = LMdenseSift(D, HOMEIMAGES, SIFTparam, HOMESIFT)
%
% Computes dense SIFT features.
% The SIFT grid will be defined by the parameters:
%    SIFTparam.grid_spacing = 1; % distance between grid centers
%    SIFTparam.patch_size = 16;  % size of patch from which to compute SIFT
%    descriptor (it has to be a factor of 4)
%
% Run demoSIFT.m to see an example of how it works.
%
% The SIFT descriptor at each location has 128 dimensions.
%
% This function can be called as:
%
% [sift, param] = LMdenseSift(D(n), HOMEIMAGES, param);
% [sift, param] = LMdenseSift(filename, HOMEIMAGES, param);
% [sift, param] = LMdenseSift(filename, HOMEIMAGES, param, HOMESIFT);
% LMdenseSift(D, HOMEIMAGES, param, HOMESIFT);
%
% 'sift' corresponds to the features of the last image. So, call it passing
% just one image. But you can precompute the SIFT features for a set of
% images: When calling LMdenseSift with a fourth argument it will store the sift descriptors in a
% new folder structure mirroring the folder structure of the images. Then,
% when called again, if the sift files already exist, it will just read
% them without recomputing them.
%
% Antonio Torralba, 2008

if nargin==4
    precomputed = 1;
    % get list of folders and create non-existing ones
    %listoffolders = {D(:).annotation.folder};
else
    precomputed = 0;
    HOMESIFT = '';
end

if nargin<3
    % Default parameters
    SIFTparam.grid_spacing = 1; % distance between grid centers
    SIFTparam.patch_size = 16; % size of patch from which to compute SIFT descriptor (it has to be a factor of 4)
end
SIFTparam.w = SIFTparam.patch_size/2; % boundary
Nfeatures = 128;

if isstruct(D)
    % [gist, param] = LMdenseSift(D, HOMEIMAGES, param);
    Nscenes = length(D);
    typeD = 1;
end

if iscell(D)
    % [gist, param] = LMdenseSift(filename, HOMEIMAGES, param);
    Nscenes = length(D);
    typeD = 2;
end

if isnumeric(D)
    % [gist, param] = LMdenseSift(img, HOMEIMAGES, param);
    Nscenes = size(D,4);
    typeD = 3;
end

if Nscenes >1
    fig = figure;
end

% Loop: Compute SIFT features for all scenes
sift = zeros([Nscenes Nfeatures], 'single');
for n = 1:Nscenes
    g = [];
    todo = 1;

    % if SIFT has already been computed, just read the file
    if precomputed==1
        filesift = fullfile(HOMESIFT, D(n).annotation.folder, [D(n).annotation.filename(1:end-4) '.mat']);
        if exist(filesift, 'file')
            load(filesift, 'sift', 'SIFTparam');
            todo = 0;
        end
    end

    % otherwise compute SIFT
    if todo==1
        disp([n Nscenes])

        % load image
        try
            switch typeD
                case 1
                    img = LMimread(D, n, HOMEIMAGES);
                case 2
                    img = imread(fullfile(HOMEIMAGES, D{n}));
                case 3
                    img = D(:,:,:,n);
            end
        catch
            disp(D(n).annotation.folder)
            disp(D(n).annotation.filename)
            rethrow(lasterror)
        end

        % get SIFT descriptors
        [sift, SIFTparam.grid_x, SIFTparam.grid_y] = dense_sift(img, SIFTparam);

        if isfield(SIFTparam, 'edges')
            % 'dont-compute': default if field not present
            % 'siftrepeat'
            w = SIFTparam.w-1;
            switch lower(SIFTparam.edges)
                case 'siftrepeat'
                    sift = [repmat(sift(1,:,:),[w 1 1]); sift; repmat(sift(end,:,:),[w 1 1])];
                    sift = [repmat(sift(:,1,:),[1 w 1]), sift, repmat(sift(:,end,:),[1 w 1])];
                otherwise
                    error('Unknown edges method')
            end
        end

        % save SIFT if a HOMESIFT file is provided
        if precomputed
            mkdir(fullfile(HOMESIFT, D(n).annotation.folder))
            save (filesift, 'sift', 'SIFTparam')
        end

        if Nscenes >1
            figure(fig);
            subplot(121)
            imshow(uint8(img))
            subplot(122)
            showColorSIFT(sift)
        end
    end

    drawnow
end

function [sift_arr, grid_x, grid_y] = dense_sift(I, SIFTparam)
% Original script by Svetlana Lazebnick
% Antonio Torralba: modified using convolutions to speed up the
% computations.

grid_spacing = SIFTparam.grid_spacing;
patch_size = SIFTparam.patch_size;

I = double(I);
I = mean(I,3);
I = I /max(I(:));

% parameters
num_angles = 8;
num_bins = 4;
num_samples = num_bins * num_bins;
alpha = 9; %% parameter for attenuation of angles (must be odd)

if nargin < 5
    sigma_edge = 1;
end

angle_step = 2 * pi / num_angles;
angles = 0:angle_step:2*pi;
angles(num_angles+1) = []; % bin centers

[hgt wid] = size(I);

[G_X,G_Y]=gen_dgauss(sigma_edge);

% add boundary:
I = [I(2:-1:1,:,:); I; I(end:-1:end-1,:,:)];
I = [I(:,2:-1:1,:) I I(:,end:-1:end-1,:)];

I = I-mean(I(:));
I_X = filter2(G_X, I, 'same'); % vertical edges
I_Y = filter2(G_Y, I, 'same'); % horizontal edges

I_X = I_X(3:end-2,3:end-2,:);
I_Y = I_Y(3:end-2,3:end-2,:);

I_mag = sqrt(I_X.^2 + I_Y.^2); % gradient magnitude
I_theta = atan2(I_Y,I_X);
I_theta(find(isnan(I_theta))) = 0; % necessary????

% grid
grid_x = patch_size/2:grid_spacing:wid-patch_size/2+1;
grid_y = patch_size/2:grid_spacing:hgt-patch_size/2+1;

% make orientation images
I_orientation = zeros([hgt, wid, num_angles], 'single');

% for each histogram angle
cosI = cos(I_theta);
sinI = sin(I_theta);
for a=1:num_angles
    % compute each orientation channel
    tmp = (cosI*cos(angles(a))+sinI*sin(angles(a))).^alpha;
    tmp = tmp .* (tmp > 0);

    % weight by magnitude
    I_orientation(:,:,a) = tmp .* I_mag;
end

% Convolution formulation:
weight_kernel = zeros(patch_size,patch_size);
r = patch_size/2;
cx = r - 0.5;
sample_res = patch_size/num_bins;
weight_x = abs((1:patch_size) - cx)/sample_res;
weight_x = (1 - weight_x) .* (weight_x <= 1);

for a = 1:num_angles
    %I_orientation(:,:,a) = conv2(I_orientation(:,:,a), weight_kernel, 'same');
    I_orientation(:,:,a) = conv2(weight_x, weight_x', I_orientation(:,:,a), 'same');
end

% Sample SIFT bins at valid locations (without boundary artifacts)
% find coordinates of sample points (bin centers)
[sample_x, sample_y] = meshgrid(linspace(1,patch_size+1,num_bins+1));
sample_x = sample_x(1:num_bins,1:num_bins); sample_x = sample_x(:)-patch_size/2;
sample_y = sample_y(1:num_bins,1:num_bins); sample_y = sample_y(:)-patch_size/2;

sift_arr = zeros([length(grid_y) length(grid_x) num_angles*num_bins*num_bins], 'single');
b = 0;
for n = 1:num_bins*num_bins
    sift_arr(:,:,b+1:b+num_angles) = I_orientation(grid_y+sample_y(n), grid_x+sample_x(n), :);
    b = b+num_angles;
end
clear I_orientation

% Outputs:
[grid_x,grid_y] = meshgrid(grid_x, grid_y);
[nrows, ncols, cols] = size(sift_arr);

% normalize SIFT descriptors

%sift_arr = reshape(sift_arr, [nrows*ncols num_angles*num_bins*num_bins]);
%sift_arr = normalize_sift(sift_arr);
%sift_arr = reshape(sift_arr, [nrows ncols num_angles*num_bins*num_bins]);

ct = .1;
sift_arr = sift_arr + ct;
tmp = sqrt(sum(sift_arr.^2, 3));
sift_arr = sift_arr ./ repmat(tmp, [1 1 size(sift_arr,3)]);

function [GX,GY]=gen_dgauss(sigma)

% laplacian of size sigma
%f_wid = 4 * floor(sigma);
%G = normpdf(-f_wid:f_wid,0,sigma);
%G = G' * G;
G = gen_gauss(sigma);
[GX,GY] = gradient(G);

GX = GX * 2 ./ sum(sum(abs(GX)));
GY = GY * 2 ./ sum(sum(abs(GY)));

function G=gen_gauss(sigma)

if all(size(sigma)==[1, 1])
    % isotropic gaussian
	f_wid = 4 * ceil(sigma) + 1;
    G = fspecial('gaussian', f_wid, sigma);
%	G = normpdf(-f_wid:f_wid,0,sigma);
%	G = G' * G;
else
    % anisotropic gaussian
    f_wid_x = 2 * ceil(sigma(1)) + 1;
    f_wid_y = 2 * ceil(sigma(2)) + 1;
    G_x = normpdf(-f_wid_x:f_wid_x,0,sigma(1));
    G_y = normpdf(-f_wid_y:f_wid_y,0,sigma(2));
    G = G_y' * G_x;
end

author: CSAILVision

repo: LabelMeToolbox


text:
%

% @
% Copyright (C) 2017 Jonas Ruesch
%
% 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/>.
% @

% The tornado installation directory.
% Must contain the aircraft geometry file and the status used below.
tornado_root_directory = 'F:\svn\dev\matlab\tornado\T135_export';

% The state:
% Velocity 12.5 m/s , ~45km/h (average velocity as measured in real flight)
% Velocity doesn't matter much for coefficient computation.
% Coefficients are largely velocity independent (tested)).
state.AS=      12.5;      %Airspeed m/s
state.alpha=   0.0;       %Angle of attack, radians
state.betha=   0;         %Angle of sideslip, radians
state.P=       0;         %Rollrate, rad/s
state.Q=       0;         %pitchrate, rad/s
state.R=       0;         %yawrate, rad/s
state.adot=    0;         %Alpha time derivative rad/s
state.bdot=    0;         %Betha time derivative rad/s
state.ALT=     0;         %Altitude, m
state.rho=     1.225;     %Desity, kg/m^3
state.pgcorr=  0;         %Apply prandtl glauert compressibility correction

% The ExperimentalCarrier geometry.
aircraft_geometry_name = 'ExperimentalCarrier';

% Center of gravity
centerOfGravity = 0.092; % Tornado x-axis extends aft (!).

% The sweep ranges and sampling resolution for Simulink modelling
alphaStart = -10;
alphaEnd = 20;
numAlphas = 31;
betaStart = -10;
betaEnd = 10;
numBetas = 21;

% % Quick Test
% alphaStart = -10;
% alphaEnd = 20;
% numAlphas = 2;
% betaStart = -10;
% betaEnd = 10;
% numBetas = 2;

% Coefficient computation
[alpha, beta, CX, CY, CZ, Cl, Cm, Cn, CX_d, CY_d, CZ_d, Cl_d, Cm_d, Cn_d, ...
     CX_P, CY_P, CZ_P, Cl_P, Cm_P, Cn_P, CX_Q, CY_Q, CZ_Q, Cl_Q, Cm_Q, Cn_Q, CX_R, CY_R, CZ_R, Cl_R, Cm_R, Cn_R]=computeAerodynamicCoefficients(...
     aircraft_geometry_name, state, tornado_root_directory, alphaStart, ...
     alphaEnd, numAlphas, betaStart, betaEnd, numBetas, centerOfGravity);

plotCoefficients(alpha, beta, CX, CY, CZ, Cl, Cm, Cn, CX_d, CY_d, CZ_d, Cl_d, Cm_d, Cn_d, ...
    CX_P, CY_P, CZ_P, Cl_P, Cm_P, Cn_P, CX_Q, CY_Q, CZ_Q, Cl_Q, Cm_Q, Cn_Q, CX_R, CY_R, CZ_R, Cl_R, Cm_R, Cn_R);

plotAircraftGeometry(aircraft_geometry_name, state, tornado_root_directory);

save;

author: jrgenerative

repo: fixed-wing-sim


text:
function Wmat = bst_shepards(destLoc, srcLoc, nbNeighbors, excludeParam, expDistance)
% BST_SHEPARDS: 3D nearest-neighbor interpolation using Shepard's weighting.
%
% USAGE:  Wmat = bst_shepards(destLoc, srcLoc, nbNeighbors=8, excludeParam=0, expDistance=2)
%
% INPUT:
%    - srcLoc       : Nx3 array of original locations, or tesselation structure (Faces,Vertices,VertConn)
%    - destLoc      : NNx3 array of locations onto original data will be interpolated, or tesselation structure (Faces,Vertices,VertConn)
%    - nbNeighbors  : Number of nearest neighbors to be considered in the interpolation (default is 8)
%    - excludeParam : If > 0, the source points that are two far away from the destination surface are ignored.
%                     Excluded points #i that have: (minDist(i) > mean(minDist) + excludeParam * std(minDist))
%                     where minDist represents the minimal distance between each source point and the destination surface
%                     If < 0, exclude the vertices that are further from the absolute distance excludeParam  (in millimeters)
%    - expDistance  : Distance exponent (if higher, influence of a value decreases faster)
%
% OUTPUT:
%    - Wmat : Interpolation matrix

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2010-2017

%% ===== PARSE INPUTS =====
% Check number of arguments
if (nargin < 2)
    error('Usage: Wmat = bst_shepards(destLoc, srcLoc, nbNeighbors, excludeParam, expDistance)');
end
% Check matrices orientation
if ((size(destLoc, 2) ~= 3) || (size(srcLoc, 2) ~= 3)) && ((size(destLoc, 2) ~= 2) || (size(srcLoc, 2) ~= 2))
    error('destLoc and srcLoc must have 2 or 3 columns.');
end
% Argument: Number of neighbors for interpolation
if (nargin < 3) || isempty(nbNeighbors)
    nbNeighbors = 8;
end
% Argument: excludeParam
if (nargin < 4) || isempty(excludeParam)
    excludeParam = 0;
end
% Argument: expDistance
if (nargin < 5) || isempty(expDistance)
    expDistance = 2;
end

%% ===== SHEPARDS INTERPOLATION =====
% Allocate interpolation matrix
nDest = size(destLoc,1);
nSrc  = size(srcLoc,1);
% Maximum number of neighbors = number of electrodes
if (nbNeighbors > nSrc)
    nbNeighbors = nSrc;
end

% Find nearest neighbors
[I,dist] = bst_nearest(srcLoc, destLoc, nbNeighbors, 1);
% Square the distance matrix
dist = dist .^ 2;
% Eliminate zeros in distance matrix for stability
dist(dist == 0) = eps;

% One neighbor
if (nbNeighbors == 1)
    Wmat = sparse(1:nDest, I(:)', ones(1,nDest), nDest, nSrc);

% More complicated cases
elseif (nbNeighbors > 1)
    % Interpolation weights from Shepards method
    W = (bst_bsxfun(@minus, dist(:,nbNeighbors), dist) ./ bst_bsxfun(@times, dist(:,nbNeighbors), dist)) .^ expDistance;
    sumW = sum(W(:,1:nbNeighbors-1),2);
    % Correct zero values: points overlap exactly => take only the first point
    iZeroW = find(sumW == 0);
    if ~isempty(iZeroW)
        sumW(iZeroW) = 1;
        W(iZeroW, 1:nbNeighbors-1) = ones(length(iZeroW),1) * [1,zeros(1,nbNeighbors-2)];
    end
    W = W(:,1:nbNeighbors-1) ./ (sumW * ones(1,nbNeighbors-1));
    % Create sparse matrix with those weights
    i = repmat((1:nDest)', nbNeighbors-1, 1);
    j = reshape(I(:, 1:nbNeighbors-1), [], 1);
    Wmat = sparse(i, j, W(:), nDest, nSrc);
end

%% ===== IGNORE VERTICES TOO FAR AWAY =====
% Set to zero the weights of the vertices that are too far away from the sources
% EEG: Distance relative to the mean distance between sensors
if (excludeParam > 0)
    % Find vertices that are too far from their nearest neighbors
    iTooFarVertices = (dist(:,1) > mean(dist(:,1)) + excludeParam * std(dist(:,1)));
    % Remove them from the interpolation matrix
    Wmat(iTooFarVertices, :) = 0;
% SEEG/ECOG: Absolute distance
elseif (excludeParam < 0)
    % Find vertices that are too far from their nearest neighbors (in millimeters)
    iTooFarVertices = (sqrt(dist(:,1)) >  abs(excludeParam));
    % Remove them from the interpolation matrix
    Wmat(iTooFarVertices, :) = 0;
end

author: brainstorm-tools


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Q=SELECT_CONFIGURATION(ROBOT, Qinv, CONF)
% Returns the joint coordinates Q that comply with the axes configuration
% vector CONF, given a set of solutions of the inverse kinematic problem Qinv.
% For 6DOF or less manipulators, the variable CONF={CF1, CF4, CF6, CFX}
% specifies univoquely only one of the solutions.
%
% See also:
%   COMPUTE_CONFIGURATION, GET_CONF_DATA
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function q=select_configuration(robot, qinv, conf)

q=qinv(:,1); %zeros(robot.DOF,1);

for i=1:size(qinv,2),
    confi=compute_configuration(robot, qinv(:,i))
    if isequal(conf(1:4), confi(1:4))
       %if the same configuration is found, store the joint values and
       %return
       q = qinv(:,i);
       return;
    end
end

disp('ERROR: RAPID/select_configuration: No solutions complies with the exact specified configuration ');

disp('WARNING: Selecting now the closest configuration');

q=select_closest_configuration(robot, qinv, conf);

repo: ARTE


text:
function Z = mldivide(X,Y)
%MLDIVIDE Slash left division for tensors.
%
%   MLDIVIDE(A,B) is called for the syntax 'A \ B' when A is a scalar and B
%   is a tensor.
%
%   Example
%   X = tenrand([4 3 2],5);
%   3 \ X
%
%   See also TENSOR, TENSOR/LDIVIDE.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if isscalar(X)
    Z = tenfun(@ldivide,X,Y);
    return;
end

error('MLDIVIDE only supports the scalar case for tensors');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function c8_tanh_test ( )

%*****************************************************************************80
%
%% C8_TANH_TEST tests C8_TANH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2015
%
%  Author:
%
%    John Burkardt
%
  seed = 123456678;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8_TANH_TEST\n' );
  fprintf ( 1, '  C8_TANH computes the hyperbolic sine of a C8.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '       C1=C8_UNIFORM_01          C2 = C8_TANH(C1)           C3 = C8_ATANH(C1)\n' );
  fprintf ( 1, '     ---------------------     ---------------------     ---------------------\n' );
  fprintf ( 1, '\n' );

  for test = 1 : 10

    [ c1, seed ] = c8_uniform_01 ( seed );
    c2 = c8_tanh ( c1 );
    c3 = c8_atanh ( c2);

    fprintf ( 1, '  (%12f  %12f)  (%12f  %12f)  (%12f  %12f)\n', ...
      real ( c1 ), imag ( c1 ), real ( c2 ), imag ( c2 ), real ( c3 ), imag ( c3 ) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function plotSection(ori,varargin)
% plot orientations to ODF sections
%
% Input
%  ori - @orientation
%
% Options
%  sections - number of sections
%  points   - number of orientations to be plotted
%  all      - plot all orientations
%
% Flags
%  phi2      - phi2 sections (default)
%  phi1      - phi1 sections
%  sigma     - sigma = phi1 ~ phi2 sections
%  axisAngle - rotational angle sections
%
% See also
% vector3d/scatter saveFigure

if check_option(varargin,{'contour','contourf','smooth'})
  warning('Using the options contour, contourf or smooth in orientation sections plots is not recommented. Computing an ODF and plotting it is usually much better');
end

if ori.antipodal, varargin = [varargin,'antipodal']; end
oS = newODFSectionPlot(ori.CS,ori.SS,varargin{:});

oS.plot(ori,varargin{:});

author: mtex-toolbox

repo: mtex

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


text:
function [inplane, isEmpty] = vol2ipCurROI(volume,inplane)
%
% [inplane, isEmpty] = vol2ipCurROI(volume,inplane)
%
% Calls vol2ipROI with the currently selected ROI.  Called from
% callback in xformInplaneMenu.
%
% Also returns isEmpty, a logical flag indicating whether the xformed
% ROI has coordinates within the inplanes or not.
%
% djh, 8/98
% djh, 2/2001, replaced globals with local variables
if ~volume.selectedROI
    myErrorDlg('Must have a selected ROI in the Volume window before it can be transformed to the Inplane window.');
end

volROI = volume.ROIs(volume.selectedROI);
ipROI = vol2ipROI(volROI,volume,inplane);

% if no coords exist in the inplanes, offer user the option of not xforming
% this ROI
if isempty(ipROI.coords)
    isEmpty = true;
    q = 'This ROI lies outside the inplane range. Xform empty ROI?';
    resp = questdlg(q, ipROI.name);
    if ~isequal(resp, 'Yes'); % cancel or close dialog
        fprintf('%s: user aborted. \n', mfilename);
        return
    end
else
    isEmpty = false;
end

inplane = addROI(inplane, ipROI, 1);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef TestInitUndistortRectifyMap
    %TestInitUndistortRectifyMap

    methods (Static)
        function test_1
            camMat = eye(3);
            distCoeffs = zeros(1,4);
            sz = [50 100];  % [w,h]
            [map1,map2] = cv.initUndistortRectifyMap(...
                camMat, distCoeffs, sz, 'M1Type','int16', ...
                'NewCameraMatrix',camMat, 'R',eye(3));
            validateattributes(map1, {'int16'}, {'size',[sz(2) sz(1) 2]});
            validateattributes(map2, {'uint16'}, {'size',[sz(2) sz(1)]});
        end

        function test_2
            camMat = eye(3);
            distCoeffs = zeros(1,4);
            sz = [50 100];  % [w,h]
            [map1,map2] = cv.initUndistortRectifyMap(...
                camMat, distCoeffs, sz, 'M1Type','single1');
            validateattributes(map1, {'single'}, {'size',[sz(2) sz(1)]});
            validateattributes(map2, {'single'}, {'size',[sz(2) sz(1)]});
        end

        function test_3
            camMat = eye(3);
            distCoeffs = zeros(1,4);
            sz = [50 100];  % [w,h]
            [map1,map2] = cv.initUndistortRectifyMap(...
                camMat, distCoeffs, sz, 'M1Type','single2');
            validateattributes(map1, {'single'}, {'size',[sz(2) sz(1) 2]});
            assert(isempty(map2));
        end

        function test_error_argnum
            try
                cv.initUndistortRectifyMap();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function b = wraprad(a)
%WRAPRAD Map angles measured in radians to the interval [-pi,pi).
%
%   B = WRAPRAD(A) maps the angles in A to their equivalent in the interval
%   [-pi,pi) by adding or subtracting the appropriate multiple of 2*pi.
%
%   See also WRAPDEG, WRAPGRAD, UNWRAPDEG, UNWRAP.

%   Author:      Peter J. Acklam
%   Time-stamp:  2003-10-13 14:27:25 +0200
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   % check number of input arguments
   error(nargchk(1, 1, nargin));

   PI = pi;
   TWOPI = 2*PI;

   b = a - TWOPI * floor((a + PI) / TWOPI);

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function area = trimeshSurfaceArea(v, e, f)
%TRIMESHSURFACEAREA Surface area of a triangular mesh
%
%   S = trimeshSurfaceArea(V, F)
%   S = trimeshSurfaceArea(V, E, F)
%   Computes the surface area of the mesh specified by vertex array V and
%   face array F. Vertex array is a NV-by-3 array of coordinates.
%   Face array is a NF-by-3, containing vertex indices of each face.
%
%   Example
%     % Compute area of an octahedron (equal to 2*sqrt(3)*a*a, with
%     % a = sqrt(2) in this case)
%     [v f] = createOctahedron;
%     trimeshSurfaceArea(v, f)
%     ans =
%         6.9282
%
%     % triangulate a compute area of a unit cube
%     [v f] = createCube;
%     f2 = triangulateFaces(f);
%     trimeshSurfaceArea(v, f2)
%     ans =
%         6
%
%   See also
%   meshes3d, meshSurfaceArea, trimeshMeanBreadth, triangulateFaces

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2011-08-26,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% check input number
if nargin == 2
    f = e;
end

% compute two direction vectors, using first vertex of each face as origin
v1 = v(f(:, 2), :) - v(f(:, 1), :);
v2 = v(f(:, 3), :) - v(f(:, 1), :);

% area of each triangle is half the cross product norm
vn = vectorNorm3d(crossProduct3d(v1, v2));

% sum up and normalize
area = sum(vn) / 2;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function vw = makeROIfromIndices(vw,ind,name,select,color,comments)
% makeROIfromIndices - make ROI from indices
%
% vw = makeROIfromIndices(vw,ind)
%
% 2010 SD & WZ: wrote it.

if notDefined('vw'),        error('Need view struct');   end
if notDefined('ind'),       error('Need indices');       end
if notDefined('name'),      name    = [];           end
if notDefined('select'),    select  = true;         end
if notDefined('color'),     color   = [];           end
if notDefined('comments'),  comments= [];           end

% define coordinates from indices
coords = viewGet(vw,'coords');
coords = coords(:,logical(ind));

% make new ROI
vw = newROI(vw,name,select,color,coords, comments);

return


text:
function out = vscale(F)
%VSCALE   Vertical scale of a DISKFUNV.
%   VSCL = VSCALE(F) returns the maximal vertical scale of the components of a
%   DISKFUNV F as determined by evaluating F on a coarse tensor-product grid.

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

out = max([ vscale(F.components{1}), vscale(F.components{2}) ]);

end


text:
function img = discreteTrefoil(varargin)
%DISCRETETREFOIL  Discretize a trefoil curve
%
%   IMG = discreteTrefoil(LX, LY, TREFOIL);
%   Computes the discretisation of a trefoil curve. Inputs are LX and LY,
%   two row vectors specifying pixel positions along each direction, and
%   TREFOIL = [XC YC ROUT RIN THETA], where (XC YC) is the center of the
%   trefoil, ROUT and RIN are the outer and the inner radius, and THETA is
%   the orientation of the first lobe, in degrees counted counter
%   clockwise.
%
%   Example
%     trefoil = [50.12 50.23 45 15 10];
%     img = discreteTrefoil(1:100, 1:100, trefoil);
%     imshow(img)
%
%   See also
%   imShapes, discreteEllipse, discreteStarfish
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-07-13,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% compute coordinate of image pixels
[lx, ly, varargin] = parseGridArgs(varargin{:});
[x, y]   = meshgrid(lx, ly);

% process input parameters
var = varargin{1};
center  = var(:, 1:2);
rOut    = var(:, 3);
rIn     = var(:, 4);
theta   = var(:, 5);

% compute middle radius, and radius extent
rc = (rOut + rIn) / 2;
dr = (rOut - rIn) / 2;

% transforms pixels according to shape position
tra     = createTranslation(-center);
rot     = createRotation(-deg2rad(theta));

[x, y]   = transformPoint(x, y, rot*tra);

% convert to polar coordinates
[th, rho] = cart2pol(x(:), y(:));

% compute theoretical polar distance
rhoTh = rc + dr * cos(3 * th);

% create image
img = false(size(x));
img(rho < rhoTh) = 1;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
 function g = gibbsKernDiagGradient(kern, x, covDiag)

% GIBBSKERNDIAGGRADIENT Compute the gradient of the GIBBS kernel's diagonal wrt parameters.
% FORMAT
% DESC computes the gradient of functions of the diagonal of the
% Mark Gibbs's non-stationary kernel matrix with respect to the parameters of the kernel. The
% parameters' gradients are returned in the order given by the
% gibbsKernExtractParam command.
% ARG kern : the kernel structure for which the gradients are
% computed.
% ARG x : the input data for which the gradient is being computed.
% ARG factors : partial derivatives of the function of interest with
% respect to the diagonal elements of the kernel.
% RETURN g : gradients of the relevant function with respect to each
% of the parameters. Ordering should match the ordering given in
% gibbsKernExtractParam.
%
% SEEALSO : gibbsKernParamInit, kernDiagGradient, gibbsKernExtractParam, gibbsKernGradient
%
% COPYRIGHT : Neil D. Lawrence, 2006

% KERN

g = zeros(1, kern.nParams);
g(end) = sum(covDiag);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text: %PLUS Dataset overload

author: marianux

repo: ecg-kit

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


text:
function plotQualityMeasuresForMotionTypes(srResults_globalMotion, srResults_mixedMotion, srResults_localMotion, qualityMeasure, binningFactors, srMethods, normalized, outputDir)

    figure_global = plotQualityMeasuresAsBarplot(srResults_globalMotion, qualityMeasure, binningFactors, srMethods, normalized);
    figure_mixed = plotQualityMeasuresAsBarplot(srResults_mixedMotion, qualityMeasure, binningFactors, srMethods, normalized);
    figure_local = plotQualityMeasuresAsBarplot(srResults_localMotion, qualityMeasure, binningFactors, srMethods, normalized);
    LinkFigures([figure_global figure_mixed figure_local]);

    if normalized
        qualityMeasure = [qualityMeasure, '_norm'];
    end

    % Save to Tikz.
    figure(figure_global);
    matlab2tikz([outputDir, '/', 'baselineDatasets_globalMotion_', qualityMeasure, '.tikz'], ...
        'height', '\figureheight', 'width', '\figurewidth', 'showInfo', false, ...
        'extraaxisoptions', ['xlabel near ticks,', ...
        'ylabel near ticks,', ...
        'scaled y ticks=false,', ...
        'yticklabel style={/pgf/number format/fixed, /pgf/number format/precision=2},']);

    figure(figure_mixed);
    matlab2tikz([outputDir, '/', 'baselineDatasets_mixedMotion_', qualityMeasure, '.tikz'], ...
        'height', '\figureheight', 'width', '\figurewidth', 'showInfo', false, ...
        'extraaxisoptions', ['xlabel near ticks,', ...
        'ylabel near ticks,', ...
        'scaled y ticks=false,', ...
        'yticklabel style={/pgf/number format/fixed, /pgf/number format/precision=2},']);

    figure(figure_local);
    matlab2tikz([outputDir, '/', 'baselineDatasets_localMotion_', qualityMeasure, '.tikz'], ...
        'height', '\figureheight', 'width', '\figurewidth', 'showInfo', false, ...
        'extraaxisoptions', ['xlabel near ticks,', ...
        'ylabel near ticks,', ...
        'scaled y ticks=false,', ...
        'yticklabel style={/pgf/number format/fixed, /pgf/number format/precision=2},']);

function figureNumber = plotQualityMeasuresAsBarplot(srResults, qualityMeasure, binningFactors, srMethods, normalized)

    % Compute mean of the quality measure for different SR methods and
    % binning factors.
    barMeans = zeros(length(binningFactors), length(srMethods));
    for binningFactorIdx = 1:length(binningFactors)
        for srMethodIdx = 1:length(srMethods)
            % Get quality measures for given SR method and binning factor.
            qualityMeasureVals = [srResults([srResults(:).binningFactor] == binningFactors(binningFactorIdx) &  strcmp(srMethods{srMethodIdx}, {srResults(:).srMethod}) > 0).(qualityMeasure)];
            if normalized
                % Normalize the quality measure with the input
                % low-resolution image (after NN interpolation).
                qualityMeasureVals_nn = [srResults([srResults(:).binningFactor] == binningFactors(binningFactorIdx) &  strcmp('nn', {srResults(:).srMethod}) > 0).(qualityMeasure)];
                qualityMeasureVals = (qualityMeasureVals - qualityMeasureVals_nn) ./ abs(qualityMeasureVals_nn);
            end

            % Compute the mean of the quality measure.
            barMeans(binningFactorIdx, srMethodIdx) = mean(qualityMeasureVals);
        end
    end

    % Plot mean quality measures as bar plot.
    figure;
    figureNumber = get(gcf, 'Number');
    bar_handle = bar(barMeans);
    for k = 1:length(srMethods)
        set(bar_handle(k), 'FaceColor', getFaceColorForSRMethod(srMethods{k}));
    end
    hold off;
    pbaspect([4 1.25 1]);
    set(gca, 'XTickLabel', {'2x', '3x', '4x'}, 'XTick', 1:(3*length(srMethods)));
    ymax = max(1.05*max(barMeans(:)));
    ymin = max(0.99*min(barMeans(:)), -0.25*ymax);
    if ymax < ymin
        tmp = ymin;
        ymin = ymax;
        ymax = tmp;
    end
    ylim([ymin ymax]);
    grid on;

author: thomas-koehler

repo: SupER


text:
%============================================================================
% Copyright (C) 2015, Heikki Hyyti
%
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the "Software"),
% to deal in the Software without restriction, including without limitation
% the rights to use, copy, modify, merge, publish, distribute, sublicense,
% and/or sell copies of the Software, and to permit persons to whom the
% Software is furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
% FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
% DEALINGS IN THE SOFTWARE.
%============================================================================

classdef DCM_IMU < handle
% DCM_IMU Implementation of Hyyti's IMU algorithm
%
%   If you use the algorithm in any scientific context, please cite:
%   Heikki Hyyti and Arto Visala, "A DCM Based Attitude Estimation Algorithm for Low-Cost MEMS IMUs,"
%   International Journal of Navigation and Observation, vol. 2015, Article ID 503814, 18 pages, 2015.
%   http://dx.doi.org/10.1155/2015/503814
%
%   Date          Author          Notes
%   1/12/2015     Heikki Hyyti    Initial release

    %% Public properties
    properties (Access = public)
        g0 = 9.8189;                % gravitation around Helsinki, Finland (change according to your area)
        state = [0 0 1 0 0 0]';     % States are lowest row of rotation matrix and gyroscope x y and z biases
                                    % (C_31, C_32, C_33, w_b1, w_b2, w_b3)
        q_dcm2 = 0.1^2;             % estimated variance of dcm states (gyro variance per second)
        q_gyro_bias2 = 0.0001^2;    % very small number to make bias change slowly
        r_acc2 = 0.5^2;             % variance of calibrated accelerometer (g-component)
        r_a2 = 10^2;                % large variance for some unknown acceleration (acc = a + g)
        q_dcm2_init = 1^2;          % initial variance of dcm states (for attitude estimation)
        q_gyro_bias2_init = 0.1^2;  % initial variance of bias states (for bias estimator)
        a = zeros(3,1);             % estimated non-gravitational accelerations
        yaw = 0;                    % Yaw angle around z axis (in ZYX convention)
        pitch = 0;                  % Pitch angle around y axis
        roll = 0;                   % Roll angle around x axis
        P = [];                     % estimate covariance (these are initialized in constructor below)
        H = [];                     % observation model (static)
        Q = [];                     % proces noise covariance (static part)
        first_row = [1 0 0]';       % first row of of the rotation matrix (for yaw angle estimate)
    end

    %% Public methods
    methods (Access = public)
        function obj = DCM_IMU(varargin)
            updateP = true;
            for i = 1:2:nargin
                if  strcmp(varargin{i}, 'Gravity'), obj.g0 = varargin{i+1};
                elseif  strcmp(varargin{i}, 'State'), obj.state = varargin{i+1};
                elseif  strcmp(varargin{i}, 'Covariance'), obj.P = varargin{i+1}; updateP = false;
                elseif  strcmp(varargin{i}, 'DCMVariance'), obj.q_dcm2 = varargin{i+1};
                elseif  strcmp(varargin{i}, 'BiasVariance'), obj.q_gyro_bias2 = varargin{i+1};
                elseif  strcmp(varargin{i}, 'InitialDCMVariance'), obj.q_dcm2_init = varargin{i+1};
                elseif  strcmp(varargin{i}, 'InitialBiasVariance'), obj.q_gyro_bias2_init = varargin{i+1};
                elseif  strcmp(varargin{i}, 'MeasurementVariance'), obj.r_acc2 = varargin{i+1};
                elseif  strcmp(varargin{i}, 'MeasurementVarianceVariableGain'), obj.r_a2 = varargin{i+1};
                else error('Invalid argument');
                end
            end;

            if (updateP), obj.P = [obj.q_dcm2_init*eye(3), zeros(3,3); zeros(3,3), obj.q_gyro_bias2_init*eye(3)]; end;
            obj.H = [eye(3)*obj.g0, zeros(3,3)];
            obj.Q = [obj.q_dcm2*eye(3), zeros(3,3); zeros(3,3) obj.q_gyro_bias2*eye(3)];
        end
        function obj = UpdateIMU(obj, Gyroscope, Accelerometer, SamplePeriod)
            x = obj.state;
            x_last = x;
            Q_ = SamplePeriod^2 * obj.Q; %Process noise covariance with time dependent noise

            % control input (angular velocities from gyroscopes)
            if (size(Gyroscope,1) == 3), u = Gyroscope;
            else u = Gyroscope';
            end

            % "rotation operators"
            C3X = [0 -x(3) x(2); x(3) 0 -x(1); -x(2) x(1) 0];
            UX = [0 -(u(3)-x(6)) u(2)-x(5);
                u(3)-x(6) 0 -(u(1)-x(4));
                -(u(2)-x(5)) u(1)-x(4) 0];

            % Model generation
            A = [zeros(3,3) -SamplePeriod*C3X; zeros(3,6)];
            B = [SamplePeriod*C3X; zeros(3,3)];
            F = eye(6) + [-SamplePeriod*UX, -SamplePeriod*C3X; zeros(3,6)];

            % Kalman a priori prediction
            x_predict = x + A*x + B*u;
            P_predict = F * obj.P * F' + Q_;

            % measurements/observations (acceleromeres)
            if (size(Accelerometer,1) == 3), z = Accelerometer;
            else z = Accelerometer';
            end

            % recompute R using the error between acceleration and the model of g
            % (estimate of the magnitude of a0 in a = a0 + g)
            a_predict = z - x_predict(1:3)*obj.g0;
            a_len = sqrt(a_predict'*a_predict);
            R = (a_len*obj.r_a2 + obj.r_acc2)*eye(3);

            % Kalman innovation
            y = z - obj.H*x_predict;
            S = obj.H * P_predict * obj.H' + R;

            % Kalman gain
            K = P_predict * obj.H' / S;

            % update a posteriori
            x = x_predict + K * y;

            % update a posteriori covariance
            IKH = eye(6) - K*obj.H;
            obj.P = IKH * P_predict * IKH' + K * R * K'; % for using any K

            % normalization of x & P (divide by DCM vector length)
            dcm_vector_length = sqrt(x(1)^2 + x(2)^2 + x(3)^2);
            J_33 = [x(2)^2 + x(3)^2,    -x(1)*x(2),         -x(1)*x(3); ...
                    -x(1)*x(2),         x(1)^2 + x(3)^2,    -x(2)*x(3); ...
                    -x(1)*x(3),         -x(2)*x(3),         x(1)^2 + x(2)^2];
            J = [ J_33 / (dcm_vector_length^3), zeros(3,3); zeros(3,3), eye(3)];

            % Laplace approximation of normalization function for x to P, J = Jacobian(f,x)
            % P_new = E[J*(x-x0)*(x-x0)'*J'] = J*E[(x-x0)*(x-x0)']*J' = J*P*J'
            obj.P = J*obj.P*J';
            x(1:3) = x(1:3) ./ dcm_vector_length;
            obj.state = x;

            % compute Euler angles (not exactly a part of the extended Kalman filter)
            % yaw integration through full rotation matrix
            u_nb = u - x(4:6);
            if (true)
                % Fill rotation matrix from angular values

                % cy = cos(obj.yaw); %old angles (last state before integration)
                % sy = sin(obj.yaw);
                % cp = cos(obj.pitch);
                % sp = sin(obj.pitch);
                % cr = cos(obj.roll);
                % sr = sin(obj.roll);

                % % compute needed parts of rotation matrix R
                % R11 = cy*cp;
                % R12 = cy*sp*sr-sy*cr;
                % R13 = cy*sp*cr+sy*sr;
                % R21 = sy*cp;
                % R22 = sy*sp*sr+cy*cr;
                % R23 = sy*sp*cr-cy*sr;

                % compute needed parts of rotation matrix R using state x and yaw
                cy = cos(obj.yaw); %old yaw angle (last state before integration)
                sy = sin(obj.yaw);
                d = sqrt(x_last(2)^2 + x_last(3)^2);
                d_inv = 1 / d;

                % compute needed parts of rotation matrix R (state and angle based version, equivalent with the commented version above)
                R11 = cy * d;
                R12 = -(x_last(3)*sy + x_last(1)*x_last(2)*cy) * d_inv;
                R13 = (x_last(2)*sy - x_last(1)*x_last(3)*cy) * d_inv;
                R21 = sy * d;
                R22 = (x_last(3)*cy - x_last(1)*x_last(2)*sy) * d_inv;
                R23 = -(x_last(2)*cy + x_last(1)*x_last(3)*sy) * d_inv;

                % update needed parts of R for yaw computation
                R11_new = R11 + SamplePeriod*(u_nb(3)*R12 - u_nb(2)*R13);
                R21_new = R21 + SamplePeriod*(u_nb(3)*R22 - u_nb(2)*R23);

                obj.yaw = atan2(R21_new,R11_new);
            else
                % alternative method estimating the whole rotation matrix
                % integrate full rotation matrix (using first row estimate in memory)
                x1 = obj.first_row + SamplePeriod*UX'*obj.first_row; %rotate x1 by x1 x u_nb
                x2 = C3X * x1; %second row x2 = (state x x1)
                x2 = x2 ./ sqrt(x2(1)^2 + x2(2)^2 + x2(3)^2); % normalize length of the second row
                x1 = C3X' * x2; %recalculate first row x1 = (x2 * state) (ensure perpendicularity)
                obj.first_row = x1 ./ sqrt(x1(1)^2 + x1(2)^2 + x1(3)^2); % normalize length
                obj.yaw = atan2(x2(1),obj.first_row(1));
            end

            %compute new pitch and roll angles from a posteriori states
            obj.pitch = asin(-x(1));
            obj.roll = atan2(x(2),x(3));

            % save the estimated non-gravitational acceleration
            obj.a = z - x(1:3)*obj.g0; % acceleration estimate (g reduced)
        end
    end
end

author: hhyyti

repo: dcm-imu

save_path: github-repos/MATLAB/hhyyti-dcm-imu


text:
function newFunction(fxnname)
%CODETEMPLATE creates a new m-file with proper help layout
%
% SYNOPSIS: newFunction(fxnname)
%
% INPUT fxnname: name of the new function. CodeTemplate will check that the
%                function name has not been used yet, and it will ask for
%                additional input, such as H1-line via inputdlg
%
% OUTPUT None. The function creates a new m-file and opens it in the editor
%
% REMARKS
%
% created with MATLAB ver.: 7.10.0.499 (R2010a) on Microsoft Windows 7 Version 6.1 (Build 7600)
%
% created by: Jonas Dorn
% DATE: 31-Mar-2010
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%==============
%% DEFAULTS
%==============

% add common superclass for value and handle classes, respectively
%superclass = {'myClass','myHandle'};
superclass = {[],'handle'};

%===============
%% CHECK INPUT
%===============

% if no input argument, assign empty to functionname
if nargin < 1 || isempty(fxnname) || ~ischar(fxnname)
    fxnname = [];
end

% check whether function already exists - only if given as input
if ~isempty(fxnname)
    fxnCheck = which(fxnname);
    if ~isempty(fxnCheck)
        error('%s already exists: %s',fxnname,fxnCheck)
    end
end

%==============

%================
%% GATHER DATA
%================

% data-gathering stage begins below
% getenv() and version are used to derive environmental variables

% find username. This will return empty on Linux, but we have a line for
% the username in the input dialogue, anyway.
username = getenv('username');

try %#ok<TRYNC>
    % try to resolve user name - I use a function I called username2name to
    % change e.g.  'jdorn' into 'Jonas Dorn'
    username = username2name(username);
end

% get version, OS, date
vers    = version;
datetoday = date;

% getenv('OS') does not work on Mac. Query the OS like ver.m
% find platform OS
if ispc
    platform = [system_dependent('getos'),' ',system_dependent('getwinsys')];
elseif ismac
    [fail, input] = unix('sw_vers');
    if ~fail
        platform = strrep(input, 'ProductName:', '');
        platform = strrep(platform, sprintf('\t'), '');
        platform = strrep(platform, sprintf('\n'), ' ');
        platform = strrep(platform, 'ProductVersion:', ' Version: ');
        platform = strrep(platform, 'BuildVersion:', 'Build: ');
    else
        platform = system_dependent('getos');
    end
else
    platform = system_dependent('getos');
end
os = platform;

% ask for the rest of the input with inputdlg

% set up inputdlg
inPrompt = {'Username',...
    'Function name',...
    'Description: ''FUNCTIONNAME does ...'' (captitalized function name)',...
    'Synopsis: ''[output1, output2] = functionname(input1, input2)',...
    'Description of input arguments (use \n for additional line breaks)',...
    'Description of output arguments (use \n for additional line breaks)',...
    'Function: 1, Value Class: 2, Handle Class: 3, Subclass: superclass-name'};
inTitle = 'Please describe your function';
numLines = repmat([1,100],7,1);
numLines([5,6],1) = 5;
% assign defaultAnswer
[defaultAnswer{1:7}] = deal('');
% assign username in default answers
defaultAnswer{1} = username;
defaultAnswer{7} = '1';
% if we have a function name already, all will be simpler
if ~isempty(fxnname)
    defaultAnswer{2} = fxnname;
    defaultAnswer{3} = sprintf('%s ...', upper(fxnname));
end

% loop till description is ok
ok = false;

while ~ok

    % get input
    description = inputdlg(inPrompt, inTitle, numLines, defaultAnswer);

    % check for user abort
    if isempty(description)
        error('description cancelled by user');
    else
        ok = true; % hope for the best
    end

    % read description. Functionname, description and synopsis are required
    fxnname = description{2};
    if isempty(fxnname) || isempty(description{2})...
            || isempty(description{3}) || isempty(description{4}) ||...
            any(findstr(description{3},'...')) || isempty(regexpi(description{4},fxnname))
        h = errordlg('Username, function name, description and synopsis are required inputs!');
        uiwait(h);
        ok = false;
    end

    % check whether function already exists (again)
    fxnCheck = which(fxnname);
    if ~isempty(fxnCheck)
        h = errordlg('%s already exists: %s',fxnname,fxnCheck);
        uiwait(h);
        ok = false;
    end

    % check for ok and update defaultAnswer with description if necessary
    if ~ok
        defaultAnswer = description;
    end

end % while

% read other input
username = description{1};
desc = description{3};
synopsis = description{4};
% if there are line breaks in inputtext and outputtext: make sure that
% these lines will still be commented!
inputtext = description{5};
% add line breaks if necessary, turn into single line of text
if size(inputtext,1) > 1
    inputtext = [inputtext,repmat('\n',size(inputtext,1),1)];
    inputtext(end,end-1:end) = ' ';
    inputtext = inputtext';
    inputtext = inputtext(:)';
    % kill some white space
    inputtext = regexprep(inputtext,'(\s*)\\n','\\n');
end
if strfind(inputtext,'\n')
    inputtext = regexprep(inputtext,'\\n','\\n%%\\t\\t');
end
outputtext = description{6};
if size(outputtext,1) > 1
    outputtext = [outputtext,repmat('\n',size(outputtext,1),1)];
    outputtext(end,end-1:end) = ' ';
    outputtext = outputtext';
    outputtext = outputtext(:)';
    outputtext = regexprep(outputtext,'(\s*)\\n','\\n');
end
if strfind(outputtext,'\n')
    outputtext = regexprep(outputtext,'\\n','\\n%%\\t\\t\\t');
end

% ask for directory to save the function
dirName = uigetdir(pwd,'select save directory');
if dirName == 0
    error('directory selection cancelled by user')
end
% end of data-gathering stage

%=================

%=================
%% WRITE FUNCTION
%=================

% create the filename based on the function name
fsuffix = '.m';
filename = fullfile(dirName,[fxnname fsuffix]);
% end filename creation

% check for function vs. class and add default superclasses if necessary
d7 = str2double(description{7});
if isnan(d7)
    superclass = description{7};
    isClass = true;
elseif d7 == 2 || d7 == 3
    isClass = true;
    superclass = superclass{d7-1};
else
    isClass = false;
end

if ~isClass % it's a function
    % beginning of file-printing stage
    fid = fopen(filename,'wt');
    fprintf(fid,'function %s\n',synopsis);
    fprintf(fid,'%%%s\n',desc);
    fprintf(fid,'%%\n');
    fprintf(fid,'%% SYNOPSIS: %s\n',synopsis);
    fprintf(fid,'%%\n');
    fprintf(fid,['%% INPUT ',inputtext,'\n']);
    fprintf(fid,'%%\n');
    fprintf(fid,['%% OUTPUT ',outputtext,'\n']);
    fprintf(fid,'%%\n');
    fprintf(fid,'%% REMARKS\n');
    fprintf(fid,'%%\n');
    fprintf(fid,'%% created with MATLAB ver.: %s on %s\n',vers,os);
    fprintf(fid,'%%\n');
    fprintf(fid,'%% created by: %s\n',username);
    fprintf(fid,'%% DATE: %s\n',datetoday);
    fprintf(fid,'%%\n');
    fprintf(fid,'%s\n',repmat('%',1,75));
    fclose(fid);
    % end of file-printing stage
else % it's a class
    fid = fopen(filename,'wt');
    fprintf(fid,'%%%s\n',desc);
    fprintf(fid,'%%\n');
    fprintf(fid,'%% CONSTRUCTOR: %s\n',synopsis);
    fprintf(fid,['%%   IN : ',inputtext,'\n']);
    fprintf(fid,['%%   OUT: ',outputtext,'\n']);
    fprintf(fid,'%%\n');
    fprintf(fid,'%% PROPERTIES\n');
    fprintf(fid,'%%   #property1:\n');
    fprintf(fid,'%%\n');
    fprintf(fid,'%% METHODS\n');
    fprintf(fid,'%%   #method1: short description\n');
    fprintf(fid,'%%        out = method1(obj,in)\n');
    fprintf(fid,'%%        Description of input and output\n');
    fprintf(fid,'%%\n');
    fprintf(fid,'%% REMARKS\n');
    fprintf(fid,'%%\n');
    fprintf(fid,'%% created with MATLAB ver.: %s on %s\n',vers,os);
    fprintf(fid,'%%\n');
    fprintf(fid,'%% created by: %s\n',username);
    fprintf(fid,'%% DATE: %s\n',datetoday);
    fprintf(fid,'%%\n');
    fprintf(fid,'%s\n',repmat('%',1,75));
    % subclass myClass or myHandle by default
    if ~isempty(superclass)
        fprintf(fid,'\nclassdef %s < %s\n\nend',fxnname,superclass); %#ok<PFCEL>
    else
        fprintf(fid,'\nclassdef %s\n\nend',fxnname);
    end
    fclose(fid);
end

% pop up the newly generated file
edit(filename);


text:
function [qrs_amp_raw,qrs_i_raw,delay]=pan_tompkin(ecg,fs,gr)

%% function [qrs_amp_raw,qrs_i_raw,delay]=pan_tompkin(ecg,fs)
% Complete implementation of Pan-Tompkins algorithm

%% Inputs
% ecg : raw ecg vector signal 1d signal
% fs : sampling frequency e.g. 200Hz, 400Hz and etc
% gr : flag to plot or not plot (set it 1 to have a plot or set it zero not
% to see any plots
%% Outputs
% qrs_amp_raw : amplitude of R waves amplitudes
% qrs_i_raw : index of R waves
% delay : number of samples which the signal is delayed due to the
% filtering
%% Method :

%% PreProcessing
% 1) Signal is preprocessed , if the sampling frequency is higher then it is downsampled
% and if it is lower upsampled to make the sampling frequency 200 Hz
% with the same filtering setups introduced in Pan
% tompkins paper (a combination of low pass and high pass filter 5-15 Hz)
% to get rid of the baseline wander and muscle noise.

% 2) The filtered signal
% is derivated using a derivating filter to high light the QRS complex.

% 3) Signal is squared.4)Signal is averaged with a moving window to get rid
% of noise (0.150 seconds length).

% 5) depending on the sampling frequency of your signal the filtering
% options are changed to best match the characteristics of your ecg signal

% 6) Unlike the other implementations in this implementation the desicion
% rule of the Pan tompkins is implemented completely.

%% Decision Rule
% At this point in the algorithm, the preceding stages have produced a roughly pulse-shaped
% waveform at the output of the MWI . The determination as to whether this pulse
% corresponds to a QRS complex (as opposed to a high-sloped T-wave or a noise artefact) is
% performed with an adaptive thresholding operation and other decision
% rules outlined below;

% a) FIDUCIAL MARK - The waveform is first processed to produce a set of weighted unit
% samples at the location of the MWI maxima. This is done in order to localize the QRS
% complex to a single instant of time. The w[k] weighting is the maxima value.

% b) THRESHOLDING - When analyzing the amplitude of the MWI output, the algorithm uses
% two threshold values (THR_SIG and THR_NOISE, appropriately initialized during a brief
% 2 second training phase) that continuously adapt to changing ECG signal quality. The
% first pass through y[n] uses these thresholds to classify the each non-zero sample
% (CURRENTPEAK) as either signal or noise:
% If CURRENTPEAK > THR_SIG, that location is identified as a QRS complex
% candidate?and the signal level (SIG_LEV) is updated:
% SIG _ LEV = 0.125 CURRENTPEAK + 0.875?SIG _ LEV

% If THR_NOISE < CURRENTPEAK < THR_SIG, then that location is identified as a
% Noise peak?and the noise level (NOISE_LEV) is updated:
% NOISE _ LEV = 0.125CURRENTPEAK + 0.875?NOISE _ LEV
% Based on new estimates of the signal and noise levels (SIG_LEV and NOISE_LEV,
% respectively) at that point in the ECG, the thresholds are adjusted as follows:
% THR _ SIG = NOISE _ LEV + 0.25 ?(SIG _ LEV-NOISE _ LEV )
% THR _ NOISE = 0.5?(THR _ SIG)
% These adjustments lower the threshold gradually in signal segments that are deemed to
% be of poorer quality.

% c) SEARCHBACK FOR MISSED QRS COMPLEXES - In the thresholding step above, if
% CURRENTPEAK < THR_SIG, the peak is deemed not to have resulted from a QRS
% complex. If however, an unreasonably long period has expired without an abovethreshold
% peak, the algorithm will assume a QRS has been missed and perform a
% searchback. This limits the number of false negatives. The minimum time used to trigger
% a searchback is 1.66 times the current R peak to R peak time period (called the RR
% interval). This value has a physiological origin - the time value between adjacent
% heartbeats cannot change more quickly than this. The missed QRS complex is assumed
% to occur at the location of the highest peak in the interval that lies between THR_SIG and
% THR_NOISE. In this algorithm, two average RR intervals are stored,the first RR interval is
% calculated as an average of the last eight QRS locations in order to adapt to changing heart
% rate and the second RR interval mean is the mean
% of the most regular RR intervals . The threshold is lowered if the heart rate is not regular
% to improve detection.

% d) ELIMINATION OF MULTIPLE DETECTIONS WITHIN REFRACTORY PERIOD - It is
% impossible for a legitimate QRS complex to occur if it lies within 200ms after a previously
% detected one. This constraint is a physiological one ?due to the refractory period during
% which ventricular depolarization cannot occur despite a stimulus[1]. As QRS complex
% candidates are generated, the algorithm eliminates such physically impossible events,
% thereby reducing false positives.

% e) T WAVE DISCRIMINATION - Finally, if a QRS candidate occurs after the 200ms
% refractory period but within 360ms of the previous QRS, the algorithm determines
% whether this is a genuine QRS complex of the next heartbeat or an abnormally prominent
% T wave. This decision is based on the mean slope of the waveform at that position. A slope of
% less than one half that of the previous QRS complex is consistent with the slower
% changing behaviour of a T wave ?otherwise, it becomes a QRS detection.
% Extra concept : beside the points mentioned in the paper, this code also
% checks if the occured peak which is less than 360 msec latency has also a
% latency less than 0,5*mean_RR if yes this is counted as noise

% f) In the final stage , the output of R waves detected in smoothed signal is analyzed and double
% checked with the help of the output of the bandpass signal to improve
% detection and find the original index of the real R waves on the raw ecg
% signal

%% References :

%[1]PAN.J, TOMPKINS. W.J,"A Real-Time QRS Detection Algorithm" IEEE
%TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. BME-32, NO. 3, MARCH 1985.

%% Author : Hooman Sedghamiz
% Linkoping university
% email : hoose792@student.liu.se
% hooman.sedghamiz@medel.com

% Any direct or indirect use of this code should be referenced
% Copyright march 2014
%%
if ~isvector(ecg)
  error('ecg must be a row or column vector');
end

if nargin < 3
    gr = 1;   % on default the function always plots
end
ecg = ecg(:); % vectorize

%% Initialize
qrs_c =[]; %amplitude of R
qrs_i =[]; %index
SIG_LEV = 0;
nois_c =[];
nois_i =[];
delay = 0;
skip = 0; % becomes one when a T wave is detected
not_nois = 0; % it is not noise when not_nois = 1
selected_RR =[]; % Selected RR intervals
m_selected_RR = 0;
mean_RR = 0;
qrs_i_raw =[];
qrs_amp_raw=[];
ser_back = 0;
test_m = 0;
SIGL_buf = [];
NOISL_buf = [];
THRS_buf = [];
SIGL_buf1 = [];
NOISL_buf1 = [];
THRS_buf1 = [];

%% Plot differently based on filtering settings
if gr
 if fs == 200
  figure,  ax(1)=subplot(321);plot(ecg);axis tight;title('Raw ECG Signal');
 else
  figure,  ax(1)=subplot(3,2,[1 2]);plot(ecg);axis tight;title('Raw ECG Signal');
 end
end
%% Noise cancelation(Filtering) % Filters (Filter in between 5-15 Hz)
if fs == 200
%% Low Pass Filter  H(z) = ((1 - z^(-6))^2)/(1 - z^(-1))^2
b = [1 0 0 0 0 0 -2 0 0 0 0 0 1];
a = [1 -2 1];
h_l = filter(b,a,[1 zeros(1,12)]);
ecg_l = conv (ecg ,h_l);
ecg_l = ecg_l/ max( abs(ecg_l));
delay = 6; %based on the paper
if gr
ax(2)=subplot(322);plot(ecg_l);axis tight;title('Low pass filtered');
end
%% High Pass filter H(z) = (-1+32z^(-16)+z^(-32))/(1+z^(-1))
b = [-1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 32 -32 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1];
a = [1 -1];
h_h = filter(b,a,[1 zeros(1,32)]);
ecg_h = conv (ecg_l ,h_h);
ecg_h = ecg_h/ max( abs(ecg_h));
delay = delay + 16; % 16 samples for highpass filtering
if gr
ax(3)=subplot(323);plot(ecg_h);axis tight;title('High Pass Filtered');
end
else
%% bandpass filter for Noise cancelation of other sampling frequencies(Filtering)
f1=5; %cuttoff low frequency to get rid of baseline wander
f2=15; %cuttoff frequency to discard high frequency noise
Wn=[f1 f2]*2/fs; % cutt off based on fs
N = 3; % order of 3 less processing
[a,b] = butter(N,Wn); %bandpass filtering
ecg_h = filtfilt(a,b,ecg);
ecg_h = ecg_h/ max( abs(ecg_h));
if gr
ax(3)=subplot(323);plot(ecg_h);axis tight;title('Band Pass Filtered');
end
end
%% derivative filter H(z) = (1/8T)(-z^(-2) - 2z^(-1) + 2z + z^(2))
h_d = [-1 -2 0 2 1]*(1/8);%1/8*fs
ecg_d = conv (ecg_h ,h_d);
ecg_d = ecg_d/max(ecg_d);
delay = delay + 2; % delay of derivative filter 2 samples
if gr
ax(4)=subplot(324);plot(ecg_d);axis tight;title('Filtered with the derivative filter');
end
%% Squaring nonlinearly enhance the dominant peaks
ecg_s = ecg_d.^2;
if gr
ax(5)=subplot(325);plot(ecg_s);axis tight;title('Squared');
end

%% Moving average Y(nt) = (1/N)[x(nT-(N - 1)T)+ x(nT - (N - 2)T)+...+x(nT)]
ecg_m = conv(ecg_s ,ones(1 ,round(0.150*fs))/round(0.150*fs));
delay = delay + 15;

if gr
ax(6)=subplot(326);plot(ecg_m);axis tight;title('Averaged with 30 samples length,Black noise,Green Adaptive Threshold,RED Sig Level,Red circles QRS adaptive threshold');
axis tight;
end

%% Fiducial Mark
% Note : a minimum distance of 40 samples is considered between each R wave
% since in physiological point of view no RR wave can occur in less than
% 200 msec distance
[pks,locs] = findpeaks(ecg_m,'MINPEAKDISTANCE',round(0.2*fs));

%% initialize the training phase (2 seconds of the signal) to determine the THR_SIG and THR_NOISE
THR_SIG = max(ecg_m(1:2*fs))*1/3; % 0.25 of the max amplitude
THR_NOISE = mean(ecg_m(1:2*fs))*1/2; % 0.5 of the mean signal is considered to be noise
SIG_LEV= THR_SIG;
NOISE_LEV = THR_NOISE;

%% Initialize bandpath filter threshold(2 seconds of the bandpass signal)
THR_SIG1 = max(ecg_h(1:2*fs))*1/3; % 0.25 of the max amplitude
THR_NOISE1 = mean(ecg_h(1:2*fs))*1/2; %
SIG_LEV1 = THR_SIG1; % Signal level in Bandpassed filter
NOISE_LEV1 = THR_NOISE1; % Noise level in Bandpassed filter
%% Thresholding and online desicion rule

for i = 1 : length(pks)

   %% locate the corresponding peak in the filtered signal
    if locs(i)-round(0.150*fs)>= 1 && locs(i)<= length(ecg_h)
          [y_i x_i] = max(ecg_h(locs(i)-round(0.150*fs):locs(i)));
       else
          if i == 1
            [y_i x_i] = max(ecg_h(1:locs(i)));
            ser_back = 1;
          elseif locs(i)>= length(ecg_h)
            [y_i x_i] = max(ecg_h(locs(i)-round(0.150*fs):end));
          end

     end

  %% update the heart_rate (Two heart rate means one the moste recent and the other selected)
    if length(qrs_c) >= 9

        diffRR = diff(qrs_i(end-8:end)); %calculate RR interval
        mean_RR = mean(diffRR); % calculate the mean of 8 previous R waves interval
        comp =qrs_i(end)-qrs_i(end-1); %latest RR
        if comp <= 0.92*mean_RR || comp >= 1.16*mean_RR
            % lower down thresholds to detect better in MVI
                THR_SIG = 0.5*(THR_SIG);
                %THR_NOISE = 0.5*(THR_SIG);
               % lower down thresholds to detect better in Bandpass filtered
                THR_SIG1 = 0.5*(THR_SIG1);
                %THR_NOISE1 = 0.5*(THR_SIG1);

        else
            m_selected_RR = mean_RR; %the latest regular beats mean
        end

    end

      %% calculate the mean of the last 8 R waves to make sure that QRS is not
       % missing(If no R detected , trigger a search back) 1.66*mean

       if m_selected_RR
           test_m = m_selected_RR; %if the regular RR availabe use it
       elseif mean_RR && m_selected_RR == 0
           test_m = mean_RR;
       else
           test_m = 0;
       end

    if test_m
          if (locs(i) - qrs_i(end)) >= round(1.66*test_m)% it shows a QRS is missed
              [pks_temp,locs_temp] = max(ecg_m(qrs_i(end)+ round(0.200*fs):locs(i)-round(0.200*fs))); % search back and locate the max in this interval
              locs_temp = qrs_i(end)+ round(0.200*fs) + locs_temp -1; %location

              if pks_temp > THR_NOISE
               qrs_c = [qrs_c pks_temp];
               qrs_i = [qrs_i locs_temp];

               % find the location in filtered sig
               if locs_temp <= length(ecg_h)
                [y_i_t x_i_t] = max(ecg_h(locs_temp-round(0.150*fs):locs_temp));
               else
                [y_i_t x_i_t] = max(ecg_h(locs_temp-round(0.150*fs):end));
               end
               % take care of bandpass signal threshold
               if y_i_t > THR_NOISE1

                      qrs_i_raw = [qrs_i_raw locs_temp-round(0.150*fs)+ (x_i_t - 1)];% save index of bandpass
                      qrs_amp_raw =[qrs_amp_raw y_i_t]; %save amplitude of bandpass
                      SIG_LEV1 = 0.25*y_i_t + 0.75*SIG_LEV1; %when found with the second thres
               end

               not_nois = 1;
               SIG_LEV = 0.25*pks_temp + 0.75*SIG_LEV ;  %when found with the second threshold
             end

          else
              not_nois = 0;

          end
    end

    %%  find noise and QRS peaks
    if pks(i) >= THR_SIG

                 % if a QRS candidate occurs within 360ms of the previous QRS
                 % ,the algorithm determines if its T wave or QRS
                 if length(qrs_c) >= 3
                      if (locs(i)-qrs_i(end)) <= round(0.3600*fs)
                        Slope1 = mean(diff(ecg_m(locs(i)-round(0.075*fs):locs(i)))); %mean slope of the waveform at that position
                        Slope2 = mean(diff(ecg_m(qrs_i(end)-round(0.075*fs):qrs_i(end)))); %mean slope of previous R wave
                             if abs(Slope1) <= abs(0.5*(Slope2))  % slope less then 0.5 of previous R
                                 nois_c = [nois_c pks(i)];
                                 nois_i = [nois_i locs(i)];
                                 skip = 1; % T wave identification
                                 % adjust noise level in both filtered and
                                 % MVI
                                 NOISE_LEV1 = 0.125*y_i + 0.875*NOISE_LEV1;
                                 NOISE_LEV = 0.125*pks(i) + 0.875*NOISE_LEV;
                             else
                                 skip = 0;
                             end

                      end
                 end

        if skip == 0  % skip is 1 when a T wave is detected
        qrs_c = [qrs_c pks(i)];
        qrs_i = [qrs_i locs(i)];

        % bandpass filter check threshold
         if y_i >= THR_SIG1
                        if ser_back
                           qrs_i_raw = [qrs_i_raw x_i];  % save index of bandpass
                        else
                           qrs_i_raw = [qrs_i_raw locs(i)-round(0.150*fs)+ (x_i - 1)];% save index of bandpass
                        end
                           qrs_amp_raw =[qrs_amp_raw y_i];% save amplitude of bandpass
          SIG_LEV1 = 0.125*y_i + 0.875*SIG_LEV1;% adjust threshold for bandpass filtered sig
         end

        % adjust Signal level
        SIG_LEV = 0.125*pks(i) + 0.875*SIG_LEV ;
        end

    elseif THR_NOISE <= pks(i) && pks(i)<THR_SIG

         %adjust Noise level in filtered sig
         NOISE_LEV1 = 0.125*y_i + 0.875*NOISE_LEV1;
         %adjust Noise level in MVI
         NOISE_LEV = 0.125*pks(i) + 0.875*NOISE_LEV;

    elseif pks(i) < THR_NOISE
        nois_c = [nois_c pks(i)];
        nois_i = [nois_i locs(i)];

        % noise level in filtered signal
        NOISE_LEV1 = 0.125*y_i + 0.875*NOISE_LEV1;
        %end

         %adjust Noise level in MVI
        NOISE_LEV = 0.125*pks(i) + 0.875*NOISE_LEV;

    end

    %% adjust the threshold with SNR
    if NOISE_LEV ~= 0 || SIG_LEV ~= 0
        THR_SIG = NOISE_LEV + 0.25*(abs(SIG_LEV - NOISE_LEV));
        THR_NOISE = 0.5*(THR_SIG);
    end

    % adjust the threshold with SNR for bandpassed signal
    if NOISE_LEV1 ~= 0 || SIG_LEV1 ~= 0
        THR_SIG1 = NOISE_LEV1 + 0.25*(abs(SIG_LEV1 - NOISE_LEV1));
        THR_NOISE1 = 0.5*(THR_SIG1);
    end

% take a track of thresholds of smoothed signal
SIGL_buf = [SIGL_buf SIG_LEV];
NOISL_buf = [NOISL_buf NOISE_LEV];
THRS_buf = [THRS_buf THR_SIG];

% take a track of thresholds of filtered signal
SIGL_buf1 = [SIGL_buf1 SIG_LEV1];
NOISL_buf1 = [NOISL_buf1 NOISE_LEV1];
THRS_buf1 = [THRS_buf1 THR_SIG1];

 skip = 0; %reset parameters
 not_nois = 0; %reset parameters
 ser_back = 0;  %reset bandpass param
end

if gr
hold on,scatter(qrs_i,qrs_c,'m');
hold on,plot(locs,NOISL_buf,'--k','LineWidth',2);
hold on,plot(locs,SIGL_buf,'--r','LineWidth',2);
hold on,plot(locs,THRS_buf,'--g','LineWidth',2);
if ax(:)
linkaxes(ax,'x');
zoom on;
end
end

%% overlay on the signals
if gr
figure,az(1)=subplot(311);plot(ecg_h);title('QRS on Filtered Signal');axis tight;
hold on,scatter(qrs_i_raw,qrs_amp_raw,'m');
hold on,plot(locs,NOISL_buf1,'LineWidth',2,'Linestyle','--','color','k');
hold on,plot(locs,SIGL_buf1,'LineWidth',2,'Linestyle','-.','color','r');
hold on,plot(locs,THRS_buf1,'LineWidth',2,'Linestyle','-.','color','g');
az(2)=subplot(312);plot(ecg_m);title('QRS on MVI signal and Noise level(black),Signal Level (red) and Adaptive Threshold(green)');axis tight;
hold on,scatter(qrs_i,qrs_c,'m');
hold on,plot(locs,NOISL_buf,'LineWidth',2,'Linestyle','--','color','k');
hold on,plot(locs,SIGL_buf,'LineWidth',2,'Linestyle','-.','color','r');
hold on,plot(locs,THRS_buf,'LineWidth',2,'Linestyle','-.','color','g');
az(3)=subplot(313);plot(ecg-mean(ecg));title('Pulse train of the found QRS on ECG signal');axis tight;
line(repmat(qrs_i_raw,[2 1]),repmat([min(ecg-mean(ecg))/2; max(ecg-mean(ecg))/2],size(qrs_i_raw)),'LineWidth',2.5,'LineStyle','-.','Color','r');
linkaxes(az,'x');
zoom on;
end
end

author: Aiwiscal

repo: ECG-ML-DL-Algorithm-Matlab


text:
% Fig. 5.28  Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
% script for right side of Figure 5.28

clf
n=[1 2];
d=conv([1 1 0],[1 13]);
nc=[1 .05];
dc=[1 .01];
nol=conv([0 0 n],nc);
dol=conv(d,dc);
sysOL=tf(nol,dol);
K=0:.001:10;
rlocus(sysOL,K);
 hold on
 title('Fig.5.28b Root locus for lead plus lag')
  axis([-.15 .05 -.1 .1])
 z=0:.1:.9;
 wn=2:2:19;
 sgrid(z, wn)
 dcl=91*nol+dol;
 r=roots(dcl) % shows that extra  root from Lag comp is almost right on the zero.
 plot(r,'*')
 hold off

author: Sable

repo: mcbench-benchmarks


text:
function varargout = PlotPaths(varargin)

% Begin initialization code
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @GUI_OpeningFcn, ...
                   'gui_OutputFcn',  @GUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code

% --- Executes just before GUI is made visible.
function GUI_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to GUI (see VARARGIN)

% Choose default command line output for GUI
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes GUI wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = GUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
kappa=2;
theta=1;
sigma=.4;
X0=1;
n=1000;
m=7;
W=randn(n,m);
X=zeros(n,m);
X(1,:)=X0;

t=[0:1/n:1];

for j=1:m
    for i=1:n
        X(i+1,j)=X(i,j)+kappa*(theta-X(i,j))/n+sigma*W(i,j)/sqrt(n);
    end
end
plot(t,X)
xlabel('$t$','interpreter','latex','FontSize',14)
ylabel('$r_t$','interpreter','latex','FontSize',14)
title('$Some~Vasicek~Paths~:~~dr_{t}=\kappa(\theta-r_{t})dt+\sigma dW_{t}$','interpreter','latex','FontSize',14)
axis([0 1 0 2])

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
kappa=2;
theta=1;
sigma=.4;
X0=1;
n=1000;
m=7;
W=randn(n,m);
X=zeros(n,m);
X(1,:)=X0;
t=[0:1/n:1];
for j=1:m
    for i=1:n
        X(i+1,j)=X(i,j)+kappa*(theta-X(i,j))/n+sigma*sqrt(X(i,j))*W(i,j)/sqrt(n);
    end
end
plot(t,X)
xlabel('$t$','interpreter','latex','FontSize',14)
ylabel('$r_t$','interpreter','latex','FontSize',14)
title('$Some~Cox.Ingersoll.Ross~Paths~:~~dr_t=\kappa(\theta-r_t)dt+\zeta \sqrt{r_t} dW_t$','interpreter','latex','FontSize',14)
axis([0 1 0 2])

% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

n=1000;
m=50;
X0=ones(1,m);
r=0.35;
sigma=.25;
X=[X0;1+r/n+sigma*sqrt(1/n)*randn(n,m)];
t=[0:1/n:1];
Xt=cumprod(X);
plot(t,Xt)
xlabel('$t$','interpreter','latex','FontSize',13)
ylabel('$X_t$','interpreter','latex','FontSize',13)
title('$Geometric ~Brownian ~Motion~:~~dX_t=r X_t dt+ \sigma X_t dW_t$','interpreter','latex','FontSize',13)
axis([0 1 0 3])

% --- Executes on button press in pushbutton4.
function pushbutton4_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

n=1000;
m=50;
X0=zeros(1,m);
r=0.05;
sigma=.2;
X=[X0;r/n+sigma*sqrt(1/n)*randn(n,m)];
Xt=cumsum(X);
t=[0:1/n:1];
plot(t,Xt)
xlabel('$t$','interpreter','latex','FontSize',13)
ylabel('$X_t$','interpreter','latex','FontSize',13)
title('$Arithmetic ~Brownian ~Motion~:~~dX_t=r dt+ \sigma dW_t$','interpreter','latex','FontSize',13)
axis([0 1 -1 1])

% --- Executes on button press in pushbutton5.
function pushbutton5_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

n=1000;
m=20;
X0=zeros(1,m);
sigma=.35;
X=[X0;sigma*sqrt(1/n)*randn(n,m)];
Xt=cumsum(X);
Y=[X0;repmat(Xt(n+1,:),n,1)/n];
Yt=cumsum(Y);
Z=Xt-Yt;
t=[0:1/n:1];
plot(t,Z)
xlabel('$t$','interpreter','latex','FontSize',13)
ylabel('$B_t$','interpreter','latex','FontSize',13)
title('$~Brownian ~Bridge~:~~B(t)=W(t)- \frac{t}{T} W(T)$','interpreter','latex','FontSize',13)
axis([0 1 -1 1])

% --- Executes on button press in pushbutton6.
function pushbutton6_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

n=1000;
m=5;
X0=zeros(1,m);
sigma=.3;
X=[X0;sigma*sqrt(1/n)*randn(n,m)];
Y=[X0;sigma*sqrt(1/n)*randn(n,m)];
Xt=cumsum(X);
Yt=cumsum(Y);
t=[0:1/n:1];
plot(Xt,Yt)
xlabel('$W_x(t)$','interpreter','latex','FontSize',13)
ylabel('$W_y(t)$','interpreter','latex','FontSize',13)
title('$2~Dimension ~Brownian ~Motion~:~~\vec{W}(t)=(W_x(t)~~W_y(t))$','interpreter','latex','FontSize',13)
axis([-1 1 -1 1])

% --- Executes on button press in pushbutton7.
function pushbutton7_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton7 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

n=1000;
m=5;
X0=zeros(1,m);
sigma=.4;
X=[X0;sigma*sqrt(1/n)*randn(n,m)];
Y=[X0;sigma*sqrt(1/n)*randn(n,m)];
Z=[X0;sigma*sqrt(1/n)*randn(n,m)];
Xt=cumsum(X);
Yt=cumsum(Y);
Zt=cumsum(Z);
t=[0:1/n:1];
plot3(Xt,Yt,Zt)
grid on
xlabel('$W_x(t)$','interpreter','latex','FontSize',13)
ylabel('$W_y(t)$','interpreter','latex','FontSize',13)
zlabel('$W_z(t)$','interpreter','latex','FontSize',13)
title('$3~Dimension ~Brownian ~Motion~:~~\vec{W}(t)=(W_x(t)~~W_y(t)~~W_z(t))$','interpreter','latex','FontSize',13)
axis([-1 1 -1 1 -1 1])

% Rodolphe Sitter


text:
% Bezier parameters
% Copyright by Nguyen Quoc Duan - EMMC11

% Purpose : to compute the parameter terms in Bezier curves
% n : the polynomial degree of Bezier curve

function B_para=B_para(i,n,t)
B_para=(1-t)^(n-i)*t^i;


text:
function sphere_delaunay_test ( )

%*****************************************************************************80
%
%% SPHERE_DELAUNAY_TEST tests the SPHERE_DELAUNAY library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 May 2010
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_DELAUNAY_TEST\n' );
  fprintf ( 1, '  MATLAB version:\n' );
  fprintf ( 1, '  Test the SPHERE_DELAUNAY library.\n' );

  sphere_delaunay_test01 ( );
  sphere_delaunay_test02 ( );
  sphere_delaunay_test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_DELAUNAY_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function speedup_fg(sp, host, port)
% SPEEDUP_FG(SP, HOST, PORT)
%
% Speeds up Flightgear's simulation speed by a factor SP over socket
% connection (HOST, PORT).
% Run Flightgear on HOST (e.g. '127.0.0.1' = localhost) with
% parameter "--props = PORT".
%
% SP = 1: Realtime
% SP = 2: Realtime speeded up by factor 2
% ...

import java.net.Socket
import java.io.*

fg_socket = Socket(host, port);
out = PrintWriter(fg_socket.getOutputStream(), true);

out.println(sprintf('set /sim/speed-up %2.1f\r\n', sp));

out.close();
fg_socket.close();

end


text:
function	[xdata,ydata,Wout,Weff,Wrdn,Win,Wirr] = ...
			make_train_data(parm,T,Wout,Wrdn,Win,Wirr)
%  Make training and test data set
%   [xdata,ydata,Wout,Weff,Wrdn,Win,Wirr] = ...
%			make_train_data(parm,T,Wout,Wrdn,Win,Wirr)
% - Output
% xdata : Xdim x Tsample x Ntrial
% ydata : Ydim x Tsample x Ntrial
% Wout  : Output weight for time embedded input (Ydim,Xdim*Dtau)
% Weff  : nonzero output weight for (Ydim,Meff*Dtau)
%         effective input dimension : [Meff x Dtau]
% Win   : Effective input mixing matrix
% Wrdn  : Redundant input mixing matrix
% Wirr  : Irrelevant input mixing matrix
% - Input
% parm.Ntrial =  Number of trials
% parm.Ydim   =  Output dim
% parm.Xdim   =  Input dim
% parm.Meff   =  Effective Input dim
% parm.Mrdn   =  Redundant Input dim [= fix(Meff/2)]
% parm.sy     =  Output noise variance
% parm.Tau    =  Lag time steps
% parm.Dtau   =  Number of embedding dimension
%
% *** sinusoidal input (this is part of effective input)
% parm.Sdim = 	number of sinusoidal input
% parm.Tmin = 	minimum period
% parm.Tmax =   maximum period
% *** Correlation between effective input and irrelevant input
% data_id : Correlation structure of input
% = 0 : No correlation between input
% = 1 : Correlation within effective input
%     + Correlation within irrelevant input
% = 2 : Correlation within effective input
%     + Correlation within irrelevant input
%     + Redundant input which have correlation with effctive input
% = 3 : Correlation between effective input and irrelevant input
%
% *** Output weight
% parm.WXmax = maximum weight of X-dim
% parm.WYmax = bias weight for Y-dim
% parm.Wtau  = effctive time delay length of weight

if ~exist('T','var'), T = 2000; end;

% Input dim
if isfield(parm,'Xdim')
	Xdim = parm.Xdim ;
else
	Xdim = 100;
end
% Output dim
if isfield(parm,'Ydim')
	Ydim = parm.Ydim ;
else
	Ydim = 1;
end
% Effective Input dim
if isfield(parm,'Meff')
	Meff = parm.Meff ;
else
	Meff = 10; ;
end
% Sinusoidal input dimension
if isfield(parm,'Sdim')
	Sdim = parm.Sdim ;
else
	Sdim = 0 ;
end
% Redundant input dimension
if isfield(parm,'Mrdn')
	Mrdn = parm.Mrdn ;
else
	Mrdn = fix(Meff/2);
end

if Meff > Xdim, Meff = Xdim; end;
if Meff < Sdim, Sdim = Meff; end;
if (Meff + Mrdn) > Xdim, Mrdn = Xdim - Meff; end;

% data_id : Correlation structure of input
% = 0 : No correlation between input
% = 1 : Correlation within effective input
%     + Correlation within irrelevant input
% = 2 : Correlation within effective input
%     + Correlation within irrelevant input
%     + Redundant input which have correlation with effctive input
% = 3 : Correlation between effective input and irrelevant input

if isfield(parm,'data_id')
	data_id = parm.data_id;
else
	data_id = 0;
end
% Input correlation parameter
% # of principal dimension in mixing matrix
if isfield(parm,'Rcor')
	Rcor = parm.Rcor ;
else
	Rcor = 10 ;
end

if isfield(parm,'Ntrial')
	Ntrial = parm.Ntrial;
else
	Ntrial = 1;
end

if isfield(parm,'sx')
	% Sinusoidal input noise variance
	sx = parm.sx   ;
else
	sx = 0.0;
end
if isfield(parm,'sy')
	% output noise variance
	sy = parm.sy   ;
else
	sy = 0.1;
end

if length(sy) == 1 && Ydim > 1, sy = repmat(sy, [Ydim,1]); end;

if isfield(parm,'Tpred')
	Tpred = parm.Tpred ;
else
	Tpred = 0;
end
if isfield(parm,'Dtau')
	Tau   = parm.Tau   ;
	Dtau  = parm.Dtau  ;
else
	Tau   = 1;
	Dtau  = 1;
end

TPRED = Tpred + Tau*(Dtau-1);
Ndata = T + TPRED;

xdata = zeros(Meff,Ndata,Ntrial);

% Sinusoidal signal
if Sdim > 0
	xsin   = make_sin_out(Ndata, Sdim , parm);
	xdata(1:Sdim,:,:) = repmat(xsin, [1 1 Ntrial]) ...
	                  + sx*randn(Sdim,Ndata,Ntrial);
end

% Gaussian random input for effctive dimension
if Meff > Sdim
	% Effctive input
	XSdim = (Meff-Sdim);

	x = randn(XSdim,Ndata*Ntrial) ;

	% Input variable (make correlation among input)
	switch	data_id
	case	0
		% No correlation
		xdata((Sdim+1):Meff,:,:) = reshape(x ,[XSdim,Ndata,Ntrial]);
		Win = [];
	case	{1,2}
		% Mixing matrix
		if ~exist('Win','var')
			Win = make_nn_diag(XSdim , Rcor);
		end
		xdata((Sdim+1):Meff,:,:) = reshape(Win * x ,[XSdim,Ndata,Ntrial]);
	case	{3}
		Xrest = (Xdim-Sdim);
		Nrest = (Xdim-Meff);
		% Mixing matrix
		if ~exist('Win','var')
			Win = make_nn_diag(Xrest , Rcor);
		end
		% Correlation between effctive input and irrelevant input
		x = Win * randn(Xrest,Ndata*Ntrial) ;

		xdata((Sdim+1):Meff,:,:) = ...
			reshape( x(1:XSdim,:), [XSdim,Ndata,Ntrial]);
		xres  = ...
			reshape( x((XSdim+1):Xrest,:), [Nrest,Ndata,Ntrial]);
	end
end

xdata  = normalize_data(xdata);
xembed = embed_data(xdata,[],parm);

% Output weight
if exist('Wout','var')
	Wout = reshape(Wout, [Ydim, Xdim, Dtau]);
	Weff = Wout(:,1:Meff,:);
	Weff = reshape(Weff, [Ydim, Meff * Dtau]);
else
	[Wout,Weff] = make_test_weight(parm);
end

y = Weff * reshape(xembed ,[Meff*Dtau, T*Ntrial]);

ymean = mean(y(:));
yvar  = sqrt( mean((y(:) - ymean).^2) );
noise = repmultiply( randn(Ydim,T*Ntrial), sy*yvar);

ydata = zeros(Ydim,Ndata,Ntrial);
ydata(:,(TPRED+1):Ndata,:) = reshape(y + noise ,[Ydim, T, Ntrial]);

% Irrelevant input
XRdim = (Xdim-Meff);

if Xdim > Meff
	switch	data_id
	case	{0}
		% No correlation with effctive input
		% Gaussian random input
		xres = randn(XRdim,Ndata,Ntrial) ;
		Wrdn = [];
		Wirr = [];
	case	{1}
		% No correlation with effctive input
		% Correlation within irrelevant input (correlation with neighbor)
		if ~exist('Wirr','var')
			% mixing matrix for irrelevant component
			Wirr = make_nn_diag(XRdim,Rcor);
		end
		xres = randn(XRdim,Ndata*Ntrial) ;
		xres = reshape(Wirr * xres ,[XRdim,Ndata,Ntrial]);
		Wrdn = [];
	case	{2}
		% Redundant input: Linear correlation with effctive input
		if ~exist('Wrdn','var')
			% projection matrix for redundant component
			Wrdn = make_rand_orth(Mrdn,Meff);
		end

		% Irrelevant inputs
		Nrest = XRdim - Mrdn;
		if Nrest > 0
			% Correlation within irrelevant input
			if ~exist('Wirr','var')
				% mixing matrix for irrelevant component
				Wirr = make_nn_diag(Nrest,Rcor);
			end
			xres2 = randn(Nrest,Ndata*Ntrial) ;
			xres2 = reshape(Wirr * xres2 ,[Nrest,Ndata,Ntrial]);
			xres = [Wrdn * xdata; xres2];
		else
			% Redundant input: Linear correlation with effctive input
			xres  = Wrdn * xdata;
		end
	case	{3}
		Wrdn = [];
		Wirr = [];
	end

	xres  = normalize_data(xres);
else
	xres = [];
	Wrdn = [];
	Wirr = [];
end

xdata = [xdata; xres];

if Ntrial==1,
	xdata = xdata(:,:,1);
end

%save(parm.datafile,'ydata','xdata','Wout')

author: KamitaniLab

repo: GenericObjectDecoding


text:
function [results] = computeQualityMetrics(image)

    % Anisotrophy Test
    [gray,rgb] = biqaa.blindimagequality(image,8,6,0,'degree');
    results.biqaa_gray = gray;
    results.biqaa_rgb = rgb;

    % BIQI Test
    results.biqi = biqi.biqi(image);

    % Bliinds2 Test
    results.bliinds = bliinds2.bliinds2_score(image);

    % IQVG Test
    addpath('/Users/dsoellinger/Documents/git/uni/Matlab Toolbox/cjlin1_svnlib/libsvm/matlab');
    results.iqvg = iqvg.IQVG(image);
    rmpath('/Users/dsoellinger/Documents/git/uni/Matlab Toolbox/cjlin1_svnlib/libsvm/matlab');

    % NIQE Test
    results.niqe = niqe.niqe(image);

    % DIVINE Test
    addpath('/Users/dsoellinger/Documents/git/uni/Matlab Toolbox/gregfreeman_libsvm/libsvm/matlab');
    results.divine = divine.divine(image);
    rmpath('/Users/dsoellinger/Documents/git/uni/Matlab Toolbox/gregfreeman_libsvm/libsvm/matlab');

    % BRISQUE Test
    % Requires libsvm 3.22. Make sure that it is installed.
    results.brisque = brisque.brisquescore(image);

end

author: dsoellinger


text:
%% Run highPassAndPrepCleanICA on VEP PREPPED data

%% Read in the file and set the necessary parameters
basename = 'vep';
indir = 'O:\ARL_Data\VEP\VEP_Robust_1Hz'; % Input data directory used for this demo
outdir = 'O:\ARL_Data\VEP\VEP_PrepClean_Infomax';

%% Run the pipeline
for k = 1:18
    thisName = sprintf('%s_%02d', basename, k);
    fname = [indir filesep thisName '.set'];
    EEG = pop_loadset(fname);
    EEG = highPassAndPrepCleanICA(EEG, 'detrendCutoff', 1.0, ...
        'icatype', 'runica', 'extended', 0, 'fractionBad', 0.1);
    EEG.setname = [EEG.setname  'HP and PREP Clean then ICA Infomax'];
    fname = [outdir filesep thisName '.set'];
    pop_saveset(EEG, 'filename', [thisName '.set'], 'filepath', outdir, ...
                'version', '7.3');
end

author: VisLab

repo: EEG-Clean-Tools


text:
function data = data2struct(data,T,options)

if length(options.embeddedlags)>1
    L = -min(options.embeddedlags) + max(options.embeddedlags);
else
    L = options.maxorder;
end

if ~isstruct(data), data = struct('X',data); end
if ~isfield(data,'C')
    if options.K>1, data.C = NaN(size(data.X,1)-L*length(T),options.K);
    else data.C = ones(size(data.X,1)-L*length(T),1);
    end
elseif size(data.C,1)==sum(T) && L>0 % we need to trim C
    ind = [];
    for j=1:length(T)
        t0 = sum(T(1:j-1)) - L*(j-1);
        if length(options.embeddedlags)>1
            ind = [ind (t0+1-min(options.embeddedlags)):(t0+T(j)-max(options.embeddedlags)) ];
        else
            ind = [ind (t0+options.maxorder+1:t0+T(j))];
        end
    end
    data.C = data.C(ind,:);
    %warning('C has more rows than it should; the first rows of each trial will be discarded\n')
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function [ x, error_norm, iter, flag ] = gmres ( A, x, b, ...
  M, restart, max_it, tol, n )

%*****************************************************************************80
%
%% GMRES solves a linear system Ax=b using the Generalized Minimal residual.
%
%  Discussion:
%
%    The routine uses restarts; it does NOT include preconditioning.
%
%  Modified:
%
%    19 July 2004
%
%  Reference:
%
%    Richard Barrett, Michael Berry, Tony Chan, James Demmel,
%      June Donato, Jack Dongarra, Victor Eijkhout, Roidan Pozo,
%      Charles Romine, Henk van der Vorst
%    Templates for the Solution of Linear Systems: Building Blocks for
%      Iterative Methods,
%    SIAM Publications, 1993.
%
%  Parameters:
%
%    Input, real A(N,N), the nonsymmetric positive definite matrix.
%
%    Input, real X(N), the initial guess vector.
%
%    Input, real B(N), the right hand side vector.
%
%    Input, real M, not used.
%
%    Input, integer RESTART, the number of iterations between restarts.
%
%    Input, integer MAX_IT, the maximum number of iterations.
%
%    Input, real TOL, an error tolerance.
%
%    Input, integer N, the order of the matrix A.
%
%    Output, real X(N), the solution.
%
%    Output, real ERROR_NORM, the norm of the error.
%
%    Output, integer ITER, the number of iterations performed.
%
%    Output, integer FLAG, the return flag.
%    0 = the solution was found to within the specified tolerance.
%    1 = a satisfactory solution was not found.  The iteration limit
%        was exceeded.
%
  iter = 0;
  flag = 0;
  bnrm2 = norm ( b );

  if ( bnrm2 == 0.0 )
    bnrm2 = 1.0
  end

  r = ( b - A * x );
  error_norm = norm ( r ) / bnrm2;
  errorhist = [ ];
  errorhist(1) = error_norm;

  if ( error_norm < tol )
    return
  end

  m = restart;
  V = sparse(n,m+1);
  H = sparse(m+1,m);
  cs(1:m) = zeros(m,1);
  sn(1:m) = zeros(m,1);
  e1 = zeros(n,1);
  e1(1) = 1.0;

  iter2 = 0;

  for iter = 1: ceil(max_it / m) ;    % changed by L. Foster, this change
                                       % will correspond to one iteration
                                       % for every matrix-vector mult.
%%    r = M \ ( b-A*x );          % changed by L. Foster, done before
                                  %    each loop repetition
    V(:,1) = r / norm( r );
    s = norm ( r ) * e1;
%
%  Construct an orthonormal basis using Gram-Schmidt orthogonalization.
%
    for i = 1 : m

      w = A * V(:,i);

      for k = 1 : i
        H(k,i)= w' * V(:,k);
        w = w - H(k,i) * V(:,k);
      end

      H(i+1,i) = norm ( w );
      V(:,i+1) = w / H(i+1,i);
%
%  Apply the Givens rotation.
%
      for k = 1 : i-1
        temp     =  cs(k) * H(k,i) + sn(k) * H(k+1,i);
        H(k+1,i) = -sn(k) * H(k,i) + cs(k) * H(k+1,i);
        H(k,i)   = temp;
      end
%
%  Form the I-th rotation matrix.
%
      aa = H(i,i);
      bb = H(i+1,i);

      if ( bb == 0.0 )
        cs(i) = 1.0;
        sn(i) = 0.0;
      elseif ( abs ( aa ) < abs ( bb ) )
        temp = aa / bb;
        sn(i) = 1.0 / sqrt( 1.0 + temp^2 );
        cs(i) = temp * sn(i);
      else
        temp = bb / aa;
        cs(i) = 1.0 / sqrt( 1.0 + temp^2 );
        sn(i) = temp * cs(i);
      end
%
%  Approximate the residual norm.
%
      temp   = cs(i) * s(i);
      s(i+1) = -sn(i) * s(i);
      s(i)   = temp;
      H(i,i) = cs(i) * H(i,i) + sn(i) * H(i+1,i);
      H(i+1,i) = 0.0;
      error_norm  = abs ( s(i+1) ) / bnrm2;
      iter2 = iter2 + 1;
      errorhist(iter2+1) = error_norm;
%
%  Update the approximate solution.
%
      if ( error_norm <= tol | max_it <= iter2 )
        y = H(1:i,1:i) \ s(1:i);
        x = x + V(:,1:i) * y;
        break;
      end

    end

    if ( error_norm <= tol | max_it <= iter2 )
      break;
    end

    y = H(1:m,1:m) \ s(1:m);
%
%  Update the approximation.
%
    x = x + V(:,1:m) * y;
%
%  Compute the residual.
%
    r =  ( b-A*x );

    s(i+1) = norm ( r );
    error_norm = s(i+1) / bnrm2;
    iter2 = iter2 + 1;
    errorhist(iter2+1) = error_norm;

    if ( error_norm <= tol | max_it <= iter2 )
      break;
    end

  end

  if ( tol < error_norm )
    flag = 1;
  end;

  error_norm = errorhist;
  iter = iter2;

  return
end

author: johannesgerer


text:
function [connections]=findConnectionMatrix(mesh)
% Finds the sparse connection matrix for a mesh
%
%  connections=findConnectionMatrix(mesh)
%
% The mesh slots must be:
%     .uniqueVertices (nPoints,3)     : 3D coords of all mesh vertices
%     .uniqueFaceIndexList (nFaces,3) : triplets of indices specifies one corner of a face triangle
%
% The returned connections is a sparse matrix (nPoints*nPoints)
%
% AUTHOR: Wade
% DATE : Last modified 020701

nVerts=length(mesh.uniqueVertices);

% faceIndexList defines vertices that are connected to each other:
% If two vertices are in the same row, they're in the same triangle and
% therefore are connected

% The values here are going to be used to create scalars
% (connection_locations) that index into the connections matrix.
% The multiplication by nVerts makes it possible to add together
% a row and column value (see below).
%
r1=(mesh.uniqueFaceIndexList(:,1)-1)*nVerts;
r2=(mesh.uniqueFaceIndexList(:,2)-1);

r3=(mesh.uniqueFaceIndexList(:,3)-1);
r4=(mesh.uniqueFaceIndexList(:,2)-1)*nVerts;

% Calculate the connection_locations
% These are connections between (1,2), (1,3), and (2,3)
%
connect_locations=[ (r1+r2+1) ;(r1+r3+1); (r4+r3+1) ];

[cly,clx]=ind2sub([nVerts,nVerts],connect_locations);
connections = sparse(cly,clx,1,nVerts,nVerts);

%Now make it symmetric because if a is connected to b then b is connected to a.
connections=connections+(connections');
connections=(connections~=0);

return;


text:
%% Misorientation Distribution Function
%
%% TODO: Please help to extend this section
% Let us consider the uncorrelated missorientation ODF corresponding to our
% model ODF.

mtexdata titanium

odf = calcDensity(ebsd.orientations)

%%

% the uncorrelated
mdf = calcMDF(odf)

%%

plotSection(mdf,'axisAngle')

%% Axis / Angle Distribution
% Then we can plot the distribution of the rotation axes of this
% missorientation ODF

plotAxisDistribution(mdf)

%%
% and the distribution of the missorientation angles and compare them to a
% uniform ODF

close all
plotAngleDistribution(mdf)
hold all
plotAngleDistribution(ebsd.CS,ebsd.CS)
hold off
legend('model ODF','uniform ODF')

author: mtex-toolbox


text:
function centroids = computeCentroids(X, idx, K)
%COMPUTECENTROIDS returns the new centroids by computing the means of the
%data points assigned to each centroid.
%   centroids = COMPUTECENTROIDS(X, idx, K) returns the new centroids by
%   computing the means of the data points assigned to each centroid. It is
%   given a dataset X where each row is a single data point, a vector
%   idx of centroid assignments (i.e. each entry in range [1..K]) for each
%   example, and K, the number of centroids. You should return a matrix
%   centroids, where each row of centroids is the mean of the data points
%   assigned to it.
%

% Useful variables
[m n] = size(X);

% You need to return the following variables correctly.
centroids = zeros(K, n);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every centroid and compute mean of all points that
%               belong to it. Concretely, the row vector centroids(i, :)
%               should contain the mean of the data points assigned to
%               centroid i.
%
% Note: You can use a for-loop over the centroids to compute this.
%

for k = 1:K
    xPos = find(idx == k);
    if(size(xPos,1)>0)
      avg = mean(X(xPos,:));
      centroids(k, :) = avg;
    end
end

% =============================================================

end

repo: machine-learning