text:
% DIPOLEDENSITY - compute and optionally plot a measure of the 3-D spatial
%                   (in)homogeneity of a specified (large) set of 1- or 2-dipole
%                   component models, either as physical dipole density or as
%                   dipole-position entropy across subjects. In either case,
%                   take into account either all the dipoles, or only the nearest
%                   dipole from each of the subjects. If no output arguments,
%                   or if 'plot','on', paints a 3-D density|entropy brain image
%                   on slices of the Montreal Neurological Institute (MNI) mean
%                   MR brain image ('standard_BESA/avg152t1.mat'). Calls
%                   DIPPLOT,
%                   MRI3DPLOT, and Fieldtrip function FT_INSIDE_HEADMODEL.
% Usage:
%               >> [dens3d mri] = dipoledensity( dipoles, 'key',val, ... );
%
% Inputs:
%    dipoles - this may be either the same dipole structure given as input to
%              the DIPPLOT function, a 3 by n array of dipole localization or
%              a cell array containing arguments for the dipplot function. Note that
%              the 'coordformat' option below defines the coordinate space for these
%              dipoles (default is MNI). See help dipplot for more information.
%
% Optional 'key', val input pairs:
%  'mri'        - [string or struct] mri file (matlab format or file format read
%                 by fcdc_read_mri). See dipplot.m help for more information.
%  'method'     - ['alldistance'|'distance'|'entropy'|'relentropy'] method for
%                            computing density:
%                 'alldistance' - {default} take into account the gaussian-weighted
%                            distances from each voxel to all the dipoles. See
%                            'methodparam' (below) to specify a standard deviation
%                            (in mm) for the gaussian weight kernel.
%                 'distance' - take into account only the distances to the nearest
%                              dipole for each subject. See 'methodparam' (below).
%                 'entropy' - taking into account only the nearest dipole to each
%                             voxel for each subject. See 'methodparam' below.
%                 'relentropy' - as in 'entropy,' but take into account all the
%                             dipoles for each subject.
% 'methodparam' - [number] for 'distance'|'alldistance' methods (see above), the
%                 standard deviation (in mm) of the 3-D gaussian smoothing kernel.
%                 For 'entropy'|'relentropy' methods, the number of closest dipoles
%                 to include {defaults: 20 mm | 20 dipoles }
% 'subsample'   - [integer] subsampling of native MNI image {default: 2 -> 2x2x2}
% 'weight'      - [(1,ncomps) array] for 'distance'|'alldistance' methods, the
%                 relative weight of each component dipole {default: ONES}
% 'coordformat' - ['mni'|'spherical'] coordinate format if dipole location or
%                 a structure is given as input. Default is 'mni'.
% 'subjind'     - [(1,ncomps) array] subject index for each dipole model. If two
%                 dipoles are in one component model, give only one subject index.
% 'nsessions'   - [integer] for 'alldistance' method, the number of sessions to
%                 divide the output values by, so that the returned measure is
%                 dipole density per session {default: 1}
% 'plot'        - ['on'|'off'] force plotting dipole density|entropy
%                 {default: 'on' if no output arguments, else 'off'}
% 'dipplot'     - ['on'|'off'] plot the dipplot image (used for converting
%                 coordinates (default is 'off')
% 'plotargs'    - {cell array} plotting arguments for MRI3DPLOT function.
% 'volmesh_fname' - [string] precomputed mesh volume file name. If not
%                 given as input the function will recompute it (it can take from
%                 five to 20 minutes). By default this function save the volume file
%                 mesh into a file named volmesh_local.mat in the current
%                 folder.
% 'norm2JointProb' - ['on'|'off'] Use joint probability (i.e. sum of all
%                    voxel values == 1) instead of number of dipoles/cm^3.
%                    Should be used for group comparison. (default 'off')
%
% Outputs:
%  dens3d       - [3-D num array] density in dipoles per cubic centimeter. If output
%                 is returned, no plot is produced unless 'plot','on' is specified.
%  mri          - {MRI structure} used in MRI3DPLOT.
%
% Example:
%         >> fakedipoles = (rand(3,10)-0.5)*80;
%         >> [dens3d mri] = dipoledensity( fakedipoles, 'coordformat', 'mni');
%         >> mri3dplot(dens3d,mri); % replot if no output is given above
%                                   % function is called automatically
%
% ------------------------------------
% NOTES:   to do multiple subject causal-weighted density map,
% (1) concatenate dipplot coord matrices for all subject
% (2) make g.subjind vector [ones(1,ncompsS1) 2*ones(1,ncompsS2) ... N*ones(1,ncompssN)]
% (3) concatenate normalized outflows for all subjects to form weight vector
% (4) call dipoledensity function with method = 'entropy' or 'relentropy'
% ------------------------------------
%
% See also:
%           EEGLAB: DIPPLOT, MRI3DPLOT, Fieldtrip: FT_INSIDE_HEADMODEL
%
% Authors: Arnaud Delorme & Scott Makeig SCCN, INC, UCSD
% Modified by: Makoto Miyakoshi
%              Ramon Martinez-Cancino
% Copyright (C) Arnaud Delorme & Scott Makeig, SCCN/INC/UCSD, 2003-
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [prob3d, mri] = dipoledensity(dipplotargs, varargin)

    % TO DO: return in dipplot() the real 3-D location of dipoles (in posxyz)
    %        FIX the dimension order here

prob3d = []; mri = [];
if nargin < 1
   help dipoledensity
   return
end

g = finputcheck(varargin, { 'subjind'     'integer'  []               [];
                            'method' 'string' { 'relentropy','entropy','distance','alldistance' } 'alldistance';
                            'methodparam' 'real'     []               20;
                            'weight'      { 'real','cell' }  []               [];
                            'smooth'      'real'     []               0;
                            'nsessions'   'integer'  []               1;
                            'subsample'   'integer'  []               2;
                            'plotargs'    'cell'     []               {};
                            'plot'        'string'  { 'on','off' }    fastif(nargout == 0, 'on', 'off');
                            'dipplot'     'string'  { 'on','off' }   'off';
                            'coordformat' 'string'  { 'mni','spherical' }   'mni';
                            'normalization' 'string'  { 'on','off' } 'on';
                            'volmesh_fname' 'string'  []  'volmesh_local.mat';
                            'mri'         { 'struct','string' } [] '';
                            'norm2JointProb' 'string'  { 'on','off' } 'off'});
if ischar(g), error(g); end
if ~strcmpi(g.method, 'alldistance') && isempty(g.subjind)
    error('Subject indices are required for this method');
end
if ~iscell(g.weight), g.weight = { g.weight }; end

% plotting dipplot
% ----------------
if ~iscell(dipplotargs) % convert input
    if ~isstruct(dipplotargs)
        if size(dipplotargs,1) == 3, dipplotargs = dipplotargs';
        elseif size(dipplotargs,2) ~= 3
            error('If an array of dipoles is given as entry, there must be 3 columns or 3 rows for x y z');
        end
        model = [];
        for idip = 1:length(dipplotargs)
            model(idip).posxyz = dipplotargs(idip,:);
            model(idip).momxyz = [1 0 0];
            model(idip).rv = 0.5;
        end
        dipplotargs = model;
    end
    dipplotargs = { dipplotargs 'coordformat' g.coordformat };
else
    dipplotargs = { dipplotargs{:} 'coordformat' g.coordformat };
end
struct = dipplot(dipplotargs{:}, 'plot', g.dipplot, 'density', 'off');
if nargout == 0
    drawnow;
end

% retrieve coordinates in MNI space
% ---------------------------------
if 0 % deprecated
     % find dipoles
     % ------------
    hmesh = findobj(gcf, 'tag', 'mesh');
    if isempty(hmesh), error('Current figure must contain dipoles'); end
    hh = [];
    disp('Finding dipoles...');
    dips = zeros(1,200);
    for index = 1:1000
        hh = [ hh(:); findobj(gcf, 'tag', ['dipole' int2str(index) ]) ];
        dips(index) = length(findobj(gcf, 'tag', ['dipole' int2str(index) ]));
    end

    disp('Retrieving dipole positions ...');
    count = 1;
    for index = 1:length(hh)
        tmp = get(hh(index), 'userdata');
        if length(tmp) == 1
            allx(count) = tmp.eleccoord(1,1);
            ally(count) = tmp.eleccoord(1,2);
            allz(count) = tmp.eleccoord(1,3);
            alli(count) = index;
            count = count + 1;
        end
    end
end

% check weights
% -------------
if ~isempty(g.weight{1})
    if ~iscell(g.weight)
        if length(g.weight) ~= length(struct)
            error('There must be as many elements in the weight matrix as there are dipoles')
        end
    else
        if length(g.weight{1}) ~= length(struct) || length(g.weight{1}) ~= length(g.weight{end})
            error('There must be as many elements in the weight matrix as there are dipoles')
        end
    end
else
    g.weight = { ones( 1, length(struct)) };
end
if ~isempty(g.subjind)
    if length(g.subjind) ~= length(struct)
        error('There must be as many element in the subject matrix as there are dipoles')
    end
else
    g.subjind = ones( 1, length(struct));
end

% decoding dipole locations
% -------------------------
disp('Retrieving dipole positions ...');
count = 1;
for index = 1:length(struct)
    dips = size(struct(index).eleccoord,1);
    for dip = 1:dips
        allx(count) = struct(index).eleccoord(dip,1);
        ally(count) = struct(index).eleccoord(dip,2);
        allz(count) = struct(index).eleccoord(dip,3);
        alli(count) = index;
        allw1(count) = g.weight{1}(  index)/dips;
        allw2(count) = g.weight{end}(index)/dips;
        alls(count) = g.subjind(index);
        count = count + 1;
    end
end
g.weight{1}    = allw1;
g.weight{end}  = allw2;
g.subjind = alls;

% read MRI file
% -------------
if isempty(g.mri) % default MRI file
    dipfitdefs;
    load('-mat', template_models(1).mrifile); % load mri variable
    g.mri = mri;
end
if ischar(g.mri)
    try
        mri = load('-mat', g.mri);
        mri = mri.mri;
    catch
        disp('Failed to read Matlab file. Attempt to read MRI file using function read_fcdc_mri');
        try
            warning off;
            mri = read_fcdc_mri(g.mri);
            mri.anatomy = round(gammacorrection( mri.anatomy, 0.8));
            mri.anatomy = uint8(round(mri.anatomy/max(reshape(mri.anatomy, prod(mri.dim),1))*255));
            % WARNING: if using double instead of int8, the scaling is different
            % [-128 to 128 and 0 is not good]
            % WARNING: the transform matrix is not 1, 1, 1 on the diagonal, some slices may be
            % misplaced
            warning on;
        catch
            error('Cannot load file using read_fcdc_mri');
        end
    end
    g.mri = mri; % output the anatomic mri image
end

% reserve array for density
% -------------------------
prob3d = {zeros(ceil(g.mri.dim/g.subsample)) };
for i = 2:length(g.weight), prob3d{i} = prob3d{1}; end

% compute voxel size
% ------------------
point1 = g.mri.transform * [ 1 1 1 1 ]';
point2 = g.mri.transform * [ 2 2 2 1 ]';
voxvol = sum((point1(1:3)-point2(1:3)).^2)*g.subsample^3; % in mm

% compute global subject entropy if necessary
% -------------------------------------------
vals   = unique_bc(g.subjind); % the unique subject indices
if strcmpi(g.method, 'relentropy') || strcmpi(g.method, 'entropy') %%%%% entropy %%%%%%%
    newind = zeros(size(g.subjind));
    for index = 1:length(vals) % foreach subject in the cluster
        tmpind = find(g.subjind == vals(index)); % dipoles for the subject
        totcount(index) = length(tmpind); % store the number of subject dipoles
        newind(tmpind) = index; % put subject index into newind
    end
    g.subjind = newind;
    gp = totcount/sum(totcount);
    globent = -sum(gp.*log(gp));
end

% compute volume inside head mesh
% -------------------------------
dipfitdefs; % get the location of standard BEM volume file
tmp = load('-mat',DIPOLEDENSITY_STDBEM); % load MNI mesh

if ~exist(g.volmesh_fname)
    % Checking for Fieldtrip
    if exist('ft_electroderealign', 'file')~=2,error('dipoledensity: Fieldtrip toolbox is required'); end

    disp('Computing volume within head mesh...');
    [X Y Z]           = meshgrid(g.mri.xgrid(1:g.subsample:end)+g.subsample/2, ...
                                 g.mri.ygrid(1:g.subsample:end)+g.subsample/2, ...
                                 g.mri.zgrid(1:g.subsample:end)+g.subsample/2);
    [indX indY indZ ] = meshgrid(1:length(g.mri.xgrid(1:g.subsample:end)), ...
                                 1:length(g.mri.ygrid(1:g.subsample:end)), ...
                                 1:length(g.mri.zgrid(1:g.subsample:end)));
    allpoints = [ X(:)'    ; Y(:)'   ; Z(:)' ];
    allinds   = [ indX(:)' ; indY(:)'; indZ(:)' ];
    allpoints = g.mri.transform * [ allpoints ; ones(1, size(allpoints,2)) ];
    allpoints(4,:) = [];

    olddir = pwd;
    tmppath = which('ft_electroderealign');
    tmppath = fullfile(fileparts(tmppath), 'private');
    cd(tmppath);
    inside = ft_inside_headmodel(allpoints', tmp.vol);
    Inside = find(inside); Outside = find(~inside);
    cd(olddir);
    disp('Done.');

    if 0 % old code using Delaunay %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        P = tmp.vol.bnd(1).pnt;
        T = delaunayn(P); % recompute triangularization (the original one is not compatible
                          % with tsearchn) get coordinates of all points in the volume
        % search for points inside or outside the volume (takes about 14 minutes!)
        IO = tsearchn(P, T, allpoints');
        Inside        = find(isnan(IO));
        Outside       = find(~isnan(IO));
        disp('Done.');
    end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    try
        save('-mat', g.volmesh_fname, 'allpoints', 'allinds', 'Inside', 'Outside');
        disp('Saving file containing inside/outide voxel indices...');
    catch, end
else
    disp('Loading file containing inside/outide voxel indices...');
    load('-mat',g.volmesh_fname);
end
InsidePoints  = allpoints(:, Inside);
InsideIndices = allinds(:, Inside);

% scan grid and compute entropy at each voxel
% -------------------------------------------
edges = [0.5:1:length(vals)+0.5];

if ~strcmpi(g.method, 'alldistance')
    fprintf('Computing (of %d):', size(InsideIndices,2));
    % entropy calculation: have to scan voxels
    % ----------------------------------------
    for i = 1:size(InsideIndices,2)

        alldists = (InsidePoints(1,i) - allx).^2 ...
                 + (InsidePoints(2,i) - ally).^2 ...
                 + (InsidePoints(3,i) - allz).^2;
        [tmpsort indsort] = sort(alldists); % sort dipoles by distance
        tmpweights{1}   = g.weight{1}(  indsort);
        tmpweights{end} = g.weight{end}(indsort);

        if strcmpi(g.method, 'relentropy') || strcmpi(g.method, 'entropy') %%%%% entropy %%%%%%%

            subjs  = g.subjind(indsort(1:g.methodparam)); % get subject indices of closest dipoles
            p      = histc(subjs, edges);
            if strcmpi(g.method, 'relentropy')
                p      = p(1:end-1)./totcount;
                % this should be uniform if p conforms to global count for all subjects
            end
            p      = p/sum(p);
            p(find(p == 0)) = [];
            for tmpi = 1:length(g.weight)
                prob3d{1}(InsideIndices(1,i), InsideIndices(2,i), InsideIndices(3,i)) = -sum(p.*log(p));
            end
        else
            % distance to each subject
            ordsubjs  = g.subjind(indsort);
            for index = 1:length(vals) % for each subject
                tmpind = find(ordsubjs == vals(index));
                if strcmpi(g.method,'distance')
                    use_dipoles(index) = tmpind(1); % find their nearest dipole
                end
            end
            for tmpi = 1:length(g.weight)
                prob3d{tmpi}(InsideIndices(1,i), InsideIndices(2,i), InsideIndices(3,i)) = ...
                   sum(tmpweights{tmpi}(use_dipoles).*exp(-tmpsort(use_dipoles)/ ...
                           (2*g.methodparam^2))); % 3-D gaussian smooth
            end
        end
        if mod(i,100) == 0, fprintf('%d ', i); end
    end
else % 'alldistance'
    % distance calculation: can scan dipoles instead of voxels (since linear)
    % --------------------------------------------------------
    %alldists = allx.^2 + ally.^2 + allz.^2;
    %figure; hist(alldists); return; % look at distribution of distances

    fprintf('Computing (of %d):', size(allx,2));
    for tmpi=1:length(g.weight)
        tmpprob{tmpi} = zeros(1, size(InsidePoints,2));
    end
    if length(g.weight) > 1, tmpprob2 = tmpprob; end
    for i = 1:size(allx,2)
        alldists = (InsidePoints(1,:) - allx(i)).^2 + ...
                   (InsidePoints(2,:) - ally(i)).^2 + ...
                   (InsidePoints(3,:) - allz(i)).^2;
%         alldists = 1;  % TM
        for tmpi=1:length(g.weight)
            tmpprob{tmpi} = tmpprob{tmpi} + g.weight{tmpi}(i)*exp(-alldists/(2*g.methodparam^2)); % 3-D gaussian smooth
            if any(isinf(tmpprob{tmpi})), error('Infinite value in probability calculation'); end
        end
        if mod(i,50) == 0, fprintf('%d ', i); end
    end
    % copy values to 3-D mesh
    % -----------------------
    for i = 1:length(Inside)
        pnts = allinds(:,Inside(i));
        for tmpi = 1:length(g.weight)
            prob3d{tmpi}(pnts(1), pnts(2), pnts(3)) = tmpprob{tmpi}(i);
        end
    end

end
fprintf('\n');

% normalize for points inside and outside the volume
% norm2JointProb is applied before plotting
% --------------------------------------------------
if strcmpi(g.method, 'alldistance') && strcmpi(g.normalization,'on')
    for i =1:length(g.weight)
        disp('Normalizing to dipole/mm^3');
        if any(prob3d{i}(:)<0)
            fprintf('WARNING: Some probabilities are negative, this will likely cause problems when normalizing probabilities.\n');
            fprintf('It is highly recommended to turn normaliziation off by using ''normalization'' key to ''off''.\n');
        end
        totval = sum(prob3d{i}(:));                                     % total values in the head
        totdip = size(allx,2);                                          % number of dipoles
        prob3d{i} = (prob3d{i}/totval*totdip/voxvol*1000)/g.nsessions;  % time 1000 to get cubic centimeters
    end
end

% resample matrix
% ----------------
if g.subsample ~= 1
    for i =1:length(g.weight)
        prob3d{i} = prob3d{i}/g.subsample;
        newprob3d = zeros(g.mri.dim);
        X = ceil(g.mri.xgrid/g.subsample);
        Y = ceil(g.mri.ygrid/g.subsample);
        Z = ceil(g.mri.zgrid/g.subsample);
        for index = 1:size(newprob3d,3)
            newprob3d(:,:,index) = prob3d{i}(X,Y,Z(index));
        end
        prob3d{i} = newprob3d;
    end
end

% 3-D smoothing
% -------------
if g.smooth ~= 0
    disp('Smoothing...');
    for i =1:length(g.weight)
        prob3d{i} = smooth3d(prob3d{i}, g.smooth);
    end
end

% Perform normalization so that the total sum of joint prob == 1
if strcmpi(g.norm2JointProb, 'on')
    prob3d{i} = prob3d{i}/sum(prob3d{i}(:));
end

% plotting
% --------
if strcmpi(g.plot, 'off')
    close gcf;
else
    mri3dplot( prob3d, g.mri, g.plotargs{:}); % plot the density using mri3dplot()
end
return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function mne_write_cov(fid,cov)
%
%
%   mne_write_cov(fid,cov)
%
%   Write a covariance matrix to an open file
%
%   fid     - an open file id
%   cov     - the covariance matrix to write
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.4  2007/12/10 01:02:55  msh
%   Fixed writing of a diagonal covariance matrix
%
%   Revision 1.3  2006/05/05 03:50:40  msh
%   Added routines to compute L2-norm inverse solutions.
%   Added mne_write_inverse_sol_stc to write them in stc files
%   Several bug fixes in other files
%
%   Revision 1.2  2006/05/03 18:53:06  msh
%   Approaching Matlab 6.5 backward compatibility
%
%   Revision 1.1  2006/04/29 12:44:10  msh
%   Added covariance matrix writing routines.
%
%

me='MNE:mne_write_cov';

global FIFF;
if isempty(FIFF)
   FIFF = fiff_define_constants();
end

fiff_start_block(fid,FIFF.FIFFB_MNE_COV);
%
%   Dimensions etc.
%
fiff_write_int(fid,FIFF.FIFF_MNE_COV_KIND,cov.kind);
fiff_write_int(fid,FIFF.FIFF_MNE_COV_DIM,cov.dim);
if cov.nfree > 0
    fiff_write_int(fid,FIFF.FIFF_MNE_COV_NFREE,cov.nfree);
end
%
%   Channel names
%
if ~isempty(cov.names)
    fiff_write_name_list(fid,FIFF.FIFF_MNE_ROW_NAMES,cov.names);
end
%
%   Data
%
if cov.diag
   fiff_write_double(fid,FIFF.FIFF_MNE_COV_DIAG,cov.data);
else
   if issparse(cov.data)
      fiff_write_float_sparse_rcs(fid,FIFF.FIFF_MNE_COV,cov.data);
   else
      q = 1;
      vals = zeros(cov.dim*(cov.dim+1)/2,1);
      for j = 1:cov.dim
          for k = 1:j
              vals(q) = cov.data(j,k);
              q = q + 1;
          end
      end
      fiff_write_double(fid,FIFF.FIFF_MNE_COV,vals);
   end
end
%
%   Eigenvalues and vectors if present
%
if ~isempty(cov.eig) && ~isempty(cov.eigvec)
    fiff_write_float_matrix(fid,FIFF.FIFF_MNE_COV_EIGENVECTORS,cov.eigvec);
    fiff_write_double(fid,FIFF.FIFF_MNE_COV_EIGENVALUES,cov.eig);
end
%
%   Projection operator
%
fiff_write_proj(fid,cov.projs);
%
%   Bad channels
%
if ~isempty(cov.bads)
    fiff_start_block(fid,FIFF.FIFFB_MNE_BAD_CHANNELS);
    fiff_write_name_list(fid,FIFF.FIFF_MNE_CH_NAME_LIST,cov.bads);
    fiff_end_block(fid,FIFF.FIFFB_MNE_BAD_CHANNELS);
end
%
%   Done!
%
fiff_end_block(fid,FIFF.FIFFB_MNE_COV);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function w = dwt3D(x, J, af)

% 3-D Discrete Wavelet Transform
%
% USAGE:
%   w = dwt3D(x, stages, af)
% INPUT:
%   x - N1 by N2 by N3 matrix
%       1) Ni all even
%       2) min(Ni) >= 2^(J-1)*length(af)
%   J - number of stages
%   af  - analysis filters
% OUTPUT:
%   w - cell array of wavelet coefficients
% EXAMPLE:
%   [af, sf] = farras;
%   x = rand(128,64,64);
%   w = dwt3D(x,3,af);
%   y = idwt3D(w,3,sf);
%   err = x-y;
%   max(max(max(abs(err))))
%
% WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
% http://taco.poly.edu/WaveletSoftware/

for k = 1:J
    [x w{k}] = afb3D(x, af, af, af);
end
w{J+1} = x;

author: thomaskuestner


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcImPlay(Temp,Event,handles)
handles = guidata(handles.MU_matrix_display);

try
    TMatrix=handles.TMatrix;
    TMatrix=TMatrix-min(TMatrix(:));
    TMatrix=uint8((double(TMatrix)/double(max(TMatrix(:))))*255);
    implay(TMatrix);
    colormap(handles.V.Color_map);
catch me
    error_msg{1,1}='ERROR!!! implay is not working in this Matlab version.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
% StackExchange Signal Processing Q86094
% https://dsp.stackexchange.com/questions/86094
% Analyzing 2 2D Kernels Which Approximates a Gaussian Kernel
% References:
%   1.
% Remarks:
%   1.  B
% TODO:
% 	1.  C
% Release Notes Royi Avital RoyiAvital@yahoo.com
% - 1.0.000     10/01/2023
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Constants

KERNEL_A = 1; %<! The 1st kernel
KERNEL_B = 2; %<! The 2nd Kernel

%% Parameters

mK1 = (1 / 16) * [1, 2, 1; 2, 4, 2; 1, 2, 1];
mK2 = (1 / 48) * [0, 1, 2, 1, 0; 1, 2, 4, 2, 1; 2, 4, 8, 4, 2; 1, 2, 4, 2, 1; 0, 1, 2, 1, 0];

%% Generate / Load Data

%% Analysis

% Padding the 1st array to have the same size as the 2nd
mK1 = padarray(mK1, [1, 1], 0, 'both');

% SVD Decomposition of the kernels
[mU1, mS1, mV1] = svd(mK1);
[mU2, mS2, mV2] = svd(mK2);

% Separable (Approximation)
mD1 = mS1(1, 1) * mU1(:, 1) * mV1(:, 1).';
mD2 = mS2(1, 1) * mU2(:, 1) * mV2(:, 1).';

vU1 = sqrt(mS1(1, 1)) * mU1(:, 1);
vU2 = sqrt(mS2(1, 1)) * mU2(:, 1);

vU1U2 = vU1 - vU2;
vU2U1 = vU2 - vU1;

%% Display Results

figureIdx = figureIdx + 1;

maxVal = max(max(mK1(:)), max(mK2(:)));

hF = figure('Position', [100, 100, 900, 400]);
hA   = subplot(1, 2, 1);
hImgObj = imagesc(mK1);
set(hA, 'CLim', [0, maxVal]);
set(hA, 'DataAspectRatio', [1, 1, 1]);
set(get(hA, 'Title'), 'String', {['Kernel A']}, ...
    'FontSize', fontSizeTitle);

hA   = subplot(1, 2, 2);
hImgObj = imagesc(mK2);
set(hA, 'CLim', [0, maxVal]);
set(hA, 'DataAspectRatio', [1, 1, 1]);
set(get(hA, 'Title'), 'String', {['Kernel B']}, ...
    'FontSize', fontSizeTitle);

if(generateFigures == ON)
    % saveas(hF, ['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
    print(hF, ['Figure', num2str(figureIdx, figureCounterSpec), '.png'], '-dpng', '-r0'); %<! Saves as Screen Resolution
end

figureIdx = figureIdx + 1;

maxVal = max(max(mD1(:)), max(mD2(:)));

hF = figure('Position', [100, 100, 900, 400]);
hA   = subplot(1, 2, 1);
hImgObj = imagesc(mD1);
set(hA, 'CLim', [0, maxVal]);
set(hA, 'DataAspectRatio', [1, 1, 1]);
set(get(hA, 'Title'), 'String', {['Separable Approximation of Kernel A']}, ...
    'FontSize', fontSizeTitle);

hA   = subplot(1, 2, 2);
hImgObj = imagesc(mD2);
set(hA, 'CLim', [0, maxVal]);
set(hA, 'DataAspectRatio', [1, 1, 1]);
set(get(hA, 'Title'), 'String', {['Separable Approximation of Kernel B']}, ...
    'FontSize', fontSizeTitle);

if(generateFigures == ON)
    % saveas(hF, ['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
    print(hF, ['Figure', num2str(figureIdx, figureCounterSpec), '.png'], '-dpng', '-r0'); %<! Saves as Screen Resolution
end

figureIdx = figureIdx + 1;

hF = figure('Position', figPosLarge);
hA   = axes(hF);
set(hA, 'NextPlot', 'add');
hLineObj = plot(-mU1(:, 1), 'DisplayName', 'Kernel A');
set(hLineObj, 'LineWidth', lineWidthNormal);
hLineObj = plot(-mU2(:, 1), 'DisplayName', 'Kernel B');
set(hLineObj, 'LineWidth', lineWidthNormal);
set(get(hA, 'Title'), 'String', {['Separable Filters of the Kernels']}, ...
    'FontSize', fontSizeTitle);
hLegend = ClickableLegend();

if(generateFigures == ON)
    % saveas(hF, ['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
    print(hF, ['Figure', num2str(figureIdx, figureCounterSpec), '.png'], '-dpng', '-r0'); %<! Saves as Screen Resolution
end

figureIdx = figureIdx + 1;

hF = figure('Position', [100, 100, 1200, 800]);
hA   = subplot(1, 2, 1);
[~, ~, hLineObj] = PlotDft(vU1U2, 1, 'plotTitle', 'DFT of Kernal A - Kernel B in 1D Seprable Approximation', 'numFreqBins', 100);
set(hLineObj, 'LineWidth', lineWidthNormal);

hA   = subplot(1, 2, 2);
[~, ~, hLineObj] = PlotDft(vU2U1, 1, 'plotTitle', 'DFT of Kernal B - Kernel A in 1D Seprable Approximation', 'numFreqBins', 100);
set(hLineObj, 'LineWidth', lineWidthNormal);

if(generateFigures == ON)
    % saveas(hF, ['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
    print(hF, ['Figure', num2str(figureIdx, figureCounterSpec), '.png'], '-dpng', '-r0'); %<! Saves as Screen Resolution
end

% figureIdx = figureIdx + 1;
%
% hFigure = figure('Position', figPosLarge);
% hAxes   = axes(hFigure);
% set(hAxes, 'NextPlot', 'add');
% hLineObj = plot(vTheta, 10 * log10(abs(vM)));
% set(hLineObj, 'LineWidth', lineWidthNormal);
%
% set(get(hAxes, 'Title'), 'String', {['MUSIC Pseudo Spectrum']}, ...
%     'FontSize', fontSizeTitle);
% set(get(hAxes, 'XLabel'), 'String', {['Spatial Angle [Deg]']}, ...
%     'FontSize', fontSizeAxis);
% set(get(hAxes, 'YLabel'), 'String', {['Value [dB]']}, ...
%     'FontSize', fontSizeAxis);
%
% hLineObj = plot(vTheta(vPeakIdx(1:numSig)), 10 * log10(abs(vPeakVal(1:numSig))));
% set(hLineObj, 'LineStyle', 'none', 'LineWidth', lineWidthNormal, 'Marker', 'x', 'MarkerSize', markerSizeLarge, 'Color', 'r');
% for ii = 1:numSig
%     hLineObj = xline(vTheta(vPeakIdx(ii)), '-', {['Angle: ', num2str(vTheta(vPeakIdx(ii))), ' [Deg]']});
%     set(hLineObj, 'LineStyle', ':', 'LineWidth', lineWidthThin, 'Color', 'r');
% end
%
% % hLegend = ClickableLegend({['Ground Truth'], ['Input Noisy Samples'], ['TV Estimation']});
%
% if(generateFigures == ON)
%     % saveas(hFigure, ['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
%     print(hFigure, ['Figure', num2str(figureIdx, figureCounterSpec), '.png'], '-dpng', '-r0'); %<! Saves as Screen Resolution
% end

%% Auxilizary Functions

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function pointLocations = markImagePoints(imgax,varname,clr)
% Pass in a handle to an image-containing axes
% (optional; default is returned by imgca);
%
% Brett Shoelson, PhD
% brett.shoelson@mathworks.com
% 3/25/2011
%
% Copyright 2011 MathWorks, Inc.

if nargin == 0 || isempty(imgax)
    imgax = imgca;
end
axes(imgax);%#ok
imghandle = imhandles(imgax);
if isempty(imghandle)
    error('MarkImagePoints: specified axes does not contain an image!')
end
imghandle = imghandle(1); % Just in case there are multiple images in the axes
parentfig = ancestor(imgax,'figure');
pointLocationsRequested = nargin > 1;
pointLocations = [];
if pointLocationsRequested
    requestedvar = varname;
    if isempty(requestedvar)
        requestedvar = 'MarkedPoints';
    end
end
if nargin < 3
    clr = [1 0 0]; %default = red
end
oldTitleProps = get(get(imgax,'title'));
tmp       = title(imgax,'Click to define object(s). Press <ENTER> to finish selection.');
set(tmp,'color',clr,'fontsize',12, 'fontweight','b','tag','markImagePoints', 'Visible', 'on');
currBDF   = get(imgax,'ButtonDownFcn');
currKPF   = get(parentfig,'KeyPressFcn');
set(imgax,'ButtonDownFcn','');
set(imghandle,'ButtonDownFcn',@placePoint);
set(parentfig,'KeyPressFcn',  @noMorePoints);

    function roi = impoint2(varargin)
        % impoint2: improved impoint object
        roi      = impoint(varargin{:});
        % Add a context menu for adding points
        l        = findobj(roi,'type','hggroup');
        uic      = unique( get(l,'UIContextMenu') );
        for u = 1:numel(uic)
            uimenu( uic(u), 'Label', 'Delete', 'Callback', @deleteROI )
        end

        function deleteROI(src,evt) %#ok
            delete(roi);
        end

    end

    function pointLocations = noMorePoints(~,evt)
        finished  = strcmpi(evt.Key,'return');
        if finished
            % Delete title, reset original functionality
            %delete(findall(parentfig,'tag','markImagePoints'));
            set(findall(parentfig,'tag','markImagePoints'),...
                'string',oldTitleProps.String,...
                'color',oldTitleProps.Color,...
                'tag',oldTitleProps.Tag,...
                'fontsize',oldTitleProps.FontSize,...
                'fontweight',oldTitleProps.FontWeight,...
                'visible',oldTitleProps.Visible);
            set(imghandle,'ButtonDownFcn','');
            set(parentfig,'KeyPressFcn','');
            set(imgax,'ButtonDownFcn',currBDF);
            set(parentfig,'KeyPressFcn',currKPF);
            if pointLocationsRequested
                roi = findall(imgax,'type','hggroup');
                pointLocations = NaN(numel(roi),2);
                for ii = 1:numel(roi)
                    tmp = get(roi(ii),'children');
                    pointLocations(ii,1) = get(tmp(1),'xdata');
                    pointLocations(ii,2) = get(tmp(1),'ydata');
                end
                assignin('base',requestedvar,pointLocations);
            end
        end

    end

    function placePoint(varargin)
        point_loc = get(imgax,'CurrentPoint');
        point_loc = point_loc(1,1:2);
        impoint2(imgax,point_loc);
    end

end

author: Sable

repo: mcbench-benchmarks


text:
function nMC=grMinVerCover(E,d)
% Function nMC=grMinVerCover(E,d) solve the minimal vertex cover problem.
% Input parameters:
%   E(m,2) - the edges of graph;
%     1st and 2nd elements of each row is numbers of vertexes;
%     m - number of edges.
%   d(n) (optional) - the weights of vertexes,
%     n - number of vertexes.
%     If we have only 1st parameter E, then all d=1.
% Output parameter:
%   nMC - the list of the numbers of vertexes included
%     in the minimal (weighted) vertex cover.
% Uses the reduction to integer LP-problem.
% Required the Optimization Toolbox v.3.0.1 or over.
% Author: Sergiy Iglin
% e-mail: siglin@yandex.ru
% personal page: http://iglin.exponenta.ru

% ============= Input data validation ==================
if nargin<1,
  error('There are no input data!')
end
[m,n,E] = grValidation(E); % E data validation
if nargin<2, % we may only 1st parameter
  d=ones(n,1); % all weights =1
else
  d=d(:); % reshape to vector-column
  if length(d)<n, % the poor length
    error('The length of the vector d is poor!')
  else
    d=d(1:n); % First n numbers
  end
end

% ============= Parameters of integer LP problem ==========
A=zeros(n,m); % for incidence matrix
A(E(:,1:2)+repmat(([1:m]'-1)*n,1,2))=1; % we fill the incidence matrix
options=optimset('bintprog'); % the default options
options.Display='off'; % we change the output

% ============= We solve the MILP problem ==========
xmin=bintprog(d,-A',-ones(m,1),[],[],[],options);
nMC=find(round(xmin)); % the answer - numbers of vertexes
return

author: HuangCongQing


text:
%% test modulate for all oasis functions.
col = {[0 114 178],[0 158 115], [213 94 0],[230 159 0],...
    [86 180 233], [204 121 167], [64 224 208], [240 228 66]}; % colors
plot_cvx = false;

%% threshold foopsi, convolution kernel
g = [1.7, -0.712];         % AR coefficient
noise = 1;
T = 3000;
framerate = 30;
firerate = 0.5;
b = 0;              % baseline
N = 20;              % number of trials
seed = 3;          % seed for genrating random variables
[Y, trueC, trueS] = gen_data(g, noise, T, framerate, firerate, b, N, seed);
y = Y(1,:);
true_c = trueC(1,:);  %#ok<*NASGU>
true_s = trueS(1,:);
temp = roots([1, -g(1), -g(2)]);
d = max(temp);
r = min(temp);
w = 200;
ht = (exp(log(d)*(1:w)) - exp(log(r)*(1:w))) / (d-r); % convolution kernel

% case 1: all parameters are known, the kernel is the sum of two
% exponential functions
smin = 0.5;
pars = [d, r];
[c_oasis, s_oasis] = deconvolveCa(y, 'exp2', pars, 'thresholded', 'smin', smin);  %#ok<*ASGLU>
figure('name', 'threshold, exp2, known: taur, taud, smin', 'papersize', [15, 4]);
show_results;

% case 1: all parameters are known
smin = 0.5;
[c_oasis, s_oasis] = deconvolveCa(y, 'kernel', ht, 'thresholded', 'smin', smin);  %#ok<*ASGLU>
figure('name', 'threshold, kernel, known: kernel, smin', 'papersize', [15, 4]);
show_results;
%%%%%%%%%%%%%%  END %%%%%%%%%%%%%%%%%%

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function im_scale = prep_im_for_blob_size(im_size, target_size, max_size)

    im_size_min = min(im_size(1:2));
    im_size_max = max(im_size(1:2));
    im_scale = double(target_size) / im_size_min;

    % Prevent the biggest axis from being more than MAX_SIZE
    if round(im_scale * im_size_max) > max_size
        im_scale = double(max_size) / double(im_size_max);
    end
end

author: jasjeetIM

repo: Mask-RCNN

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


text:
function outsig=rampup(L,wintype)
%RAMPUP  Rising ramp function
%   Usage: outsig=rampup(L);
%
%   `rampup(L)` will return a rising ramp function of length *L*. The
%   ramp is a sinusoide starting from zero and ending at one. The ramp
%   is centered such that the first element is always 0 and the last
%   element is not quite 1, such that the ramp fits with following ones.
%
%   `rampup(L,wintype)` will use another window for ramping. This may be any
%   of the window types from |firwin|. Please see the help on |firwin| for
%   more information. The default is to use a piece of the Hann window.
%
%   See also: rampdown, rampsignal, firwin

complainif_argnonotinrange(nargin,1,2,mfilename);

if nargin==1
  wintype='hann';
end;

win=firwin(wintype,2*L,'inf');
outsig=win(L+1:2*L);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [ n_data, x, fx ] = cinh_values ( n_data )

%*****************************************************************************80
%
%% CINH_VALUES returns some values of the alternate hyperbolic cosine integral function.
%
%  Discussion:
%
%    The alternate hyperbolic cosine integral is defined by
%
%      CINH(X) =integral ( 0 <= T < X ) ( cosh ( T ) - 1 ) / T  dT
%
%    In Mathematica, the function can be evaluated by:
%
%      Integrate [ ( Cosh[t] - 1 ) / t, { t, 0, x } ]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Cambridge University Press, 1999,
%    ISBN: 0-521-64314-7,
%    LC: QA76.95.W65.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 17;

  fx_vec = [ ...
     0.00000000000000000, ...
     0.06315467070191883, ...
     0.09136085223843649, ...
     0.1250284547325902, ...
     0.1643278712460683, ...
     0.2094587379417273, ...
     0.2606512760786754, ...
     0.3823047024751071, ...
     0.5318061742668980, ...
     0.7122865135136963, ...
     0.9275748842583805, ...
     1.182304077185436, ...
     2.030919091578478, ...
     3.284564141195967, ...
     5.129213294250493, ...
     7.850037532801762, ...
    11.88451858691463 ];

  x_vec = [ ...
     0.0, ...
     0.5, ...
     0.6, ...
     0.7, ...
     0.8, ...
     0.9, ...
     1.0, ...
     1.2, ...
     1.4, ...
     1.6, ...
     1.8, ...
     2.0, ...
     2.5, ...
     3.0, ...
     3.5, ...
     4.0, ...
     4.5 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y=H_obj(x)

y=[8-2*x(2) -2*x(1);-2*x(1) 2];

repo: Optimization-Theory-and-Algorithm


text:
function [ events ] = getEventsFromPlugData(dataset, household, appliance_id, evaluation_days, granularity, filter)

    % get ground truth from plug-level data
    events = [];

    p_edgeThreshold = filter.p_edgeThreshold;
    p_filtering = filter.p_filtering;
    p_filtLength = filter.p_filtLength;

    % detect events in plug data
    appliance_consumption = read_plug_data(dataset, household, appliance_id, evaluation_days, granularity);
    function_handle = str2func(p_filtering);
    appliance_consumption_filtered = function_handle(appliance_consumption, p_filtLength);
    diff_consumption = diff(appliance_consumption_filtered);
    edges = abs(diff_consumption) > p_edgeThreshold;
    events_start_time = find(diff(edges) == 1) + 1;
    events_end_time = find(diff(edges) == -1) + 1;

%     % skip appliance if the number of events is too low
%     if length(events_start_time) < 2 || length(events_end_time) < 2
%         continue;
%     end

    % make sure that the events' end and start time are correct
    if events_end_time(1) < events_start_time(1)
        events_end_time = events_end_time(2:end);
    end
    if events_end_time(end) < events_start_time(end)
        events_start_time = events_start_time(1:end-1);
    end

    % only slect significant events
    idx_siginificant_events = find(events_start_time(2:end) - events_end_time(1:end-1) > 2);
    events_start_time = events_start_time([1, 1+idx_siginificant_events]);
    events_end_time = events_end_time([idx_siginificant_events, end]);

    % compute change in power caused by events
    power_change = appliance_consumption_filtered(1, events_end_time) - appliance_consumption_filtered(1, events_start_time);

    % ignore events that last too long ( > 60 seconds)
    idx_valid_events = events_end_time - events_start_time < 60;

    % store the ground truth (events in plug data)
    events = [events; events_start_time(idx_valid_events)', events_end_time(idx_valid_events)', power_change(idx_valid_events)'];
end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function [c,ceq] = traj_cnstr(opt_vars, traj_par)
% -------------------------------------------------------------------------
% The function computes constraints on trajectory for trajectoty
% optimization needed for dynamic parameter identification
%
% Constraints include joint limits, maximum velocities and accelerations
% ------------------------------------------------------------------------
% Trajectory parameters
N = traj_par.N;
wf = traj_par.wf;
T = traj_par.T;
t = traj_par.t;

% As paramters of the trajectory are in a signle vector we reshape them as
% to feed the function that computes the trajectory
ab = reshape(opt_vars,[12,N]);
a = ab(1:6,:); % sin coeffs
b = ab(7:12,:); % cos coeffs

% To guarantee that positions, velocities and accelerations are zero in the
% beginning and at time T, we add fifth order polynomial to fourier
% series. The parameters of the polynomial depends on the parameters of
% fourier series. Here we compute them.
c_pol = getPolCoeffs(T, a, b, wf, N, traj_par.q0);

% Compute trajectory (Fouruer series + fifth order polynomail)
[q,qd,q2d] = mixed_traj(t, c_pol, a, b, wf, N);

% Inequality constraints
c(1:6) = traj_par.q_min - min(q,[],2); % upper joint limit constraint
c(7:12) = max(q,[],2) - traj_par.q_max; % lower joint limit constraint
c(13:18) = max(abs(qd),[],2) - traj_par.qd_max; % max joint velocity const
c(19:24) = max(abs(q2d),[],2) - traj_par.q2d_max; % max joint acceleration constr

% Equality contrsints
ceq = [];

author: shamilmamedov


text:
function [assigns, totvalue] = auction(valuematrix)

% [assigns, totvalue] = auction(valuematrix)
%
% First column is null assignment - any number of objects can be assigned
% this
%
% Author: Paul Horridge

bidincrement = 1e-6; % choose better later
[nobj, nhyps] = size(valuematrix); % nhyps includes null hypothesis
prices = zeros(1, nhyps);
assigns = ones(nobj, 1);
hyps2obj = zeros(1, nhyps);
happy = false(nobj, 1);

% Repeat while someone unhappy
while any(~happy)
    % Find unhappy person
    i = find(~happy, 1, 'first');
    % Get their net values
    netvals = valuematrix(i,:) - prices;
    [bestval, bestj] = max(netvals);
    if netvals(assigns(i)) >= bestval - bidincrement
        % Check if happy
        happy(i) = true;
    else
        % Get second best hypothesis
        nextval = max(netvals([1:bestj-1 bestj+1:end]));
        % Increase price if required
        if bestj > 1
            prices(bestj) = prices(bestj) + bestval - nextval + bidincrement;
            % De-assign original holder if required
            oldi = hyps2obj(bestj);
            if oldi > 0
                assigns(oldi) = 1;
                happy(oldi) = false;
            end
        end
        % Assign best hypothesis to i
        assigns(i) = bestj;
        hyps2obj(bestj) = i;
        happy(i) = true;
    end
end
totvalue = sum(valuematrix((assigns'-1)*nobj + (1:nobj)));

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function obj = addContact(obj, contact, fric_coef, geometry, load_path)
    % Adds contact constraints for the robot manipulator
    %
    % Parameters:
    %  contact: the contact frame object @type ContactFrame
    %  fric_coef: the friction coefficient @type struct
    %  geometry: the geometry constant of the contact @type struct
    %
    % Optional fields of geometry:
    %  la: the distance from the origin to the rolling edge along
    %  the positive x-axis  @type double
    %  lb: the distance from the origin to the rolling edge along
    %  the negative x-axis  @type double
    %  La: the distance from the origin to the rolling edge along
    %  the negative y-axis  @type double
    %  Lb: the distance from the origin to the rolling edge along
    %  the positive y-axis  @type double
    %  RefFrame: The frame describing the ground contact zero configuration.
    %
    % Optional fields of fric_coef:
    %  mu: the (static) coefficient of friction. @type double
    %  gamma: the coefficient of torsional friction @type double

    assert(isa(contact, 'ContactFrame'),...
        'The contact must be given as an object of ContactFrame class.');

    if nargin < 3
        fric_coef = [];
    end
    % Equivalent to grasp map in Murray Ch 5. Maps the wrench base into the
    % nominal contact reference frame.
    if nargin < 4
        geometry = [];
        if ((~isfield(geometry, 'RefFrame')) && ~isempty(geometry))
            geometry.RefFrame = eye(3);
        end
    end

    if nargin < 5
       load_path = [];
    end

    if isempty(load_path)
        if isempty(geometry) || ~isfield(geometry, 'RefFrame')
            ref = eye(3);
        else
            ref = geometry.RefFrame;
        end
        G = [eye(3), zeros(3,3); zeros(3,3), ref] * contact.WrenchBase;

        % compute the spatial position (cartesian position + Euler angles)
        pos = getCartesianPosition(obj, contact);
        rpy = getRelativeEulerAngles(obj, contact, ref);
        %         rpy = getEulerAngles(obj, contact);

        h = transpose([pos, rpy]); %effectively as transpose
        % extract the contrained elements
        constr =  G' * h;
        % compute the body jacobian
        jac = getBodyJacobian(obj, contact);
        %         if strcmp(contact.Type,'PointContactWithFriction') || strcmp(contact.Type,'PointContactWithoutFriction')
        %             jac_pos = jacobian(pos,obj.States.x); % directly use partial derivatives for position
        %             jac = [jac_pos;jac_rot(4:6,:)];
        %         else
        %             jac = jac_rot;
        %         end
        % extract the contrained elements
        idx = sum(contact.WrenchBase,2);
        constr_jac = jac(find(idx),:); %#ok<FNDSB>

        % label for the holonomic constraint
        label_full = cellfun(@(x)[contact.Name,x],...
            {'PosX','PosY','PosZ','Roll','Pitch','Yaw'},'UniformOutput',false);
        for i=size(contact.WrenchBase,2):-1:1
            label{i} = label_full{find(contact.WrenchBase(:,i))};         %#ok<FNDSB>
        end

        % create a holonomic constraint object
        contact_constr = HolonomicConstraint(obj,...
            constr, contact.Name,...
            'Jacobian',constr_jac,...
            'ConstrLabel',{label},...
            'DerivativeOrder',2);
        % add as a set of holonomic constraints
        obj = addHolonomicConstraint(obj, contact_constr);

    else
        % label for the holonomic constraint
        label_full = cellfun(@(x)[contact.Name,x],...
            {'PosX','PosY','PosZ','Roll','Pitch','Yaw'},'UniformOutput',false);
        for i=size(contact.WrenchBase,2):-1:1
            label{i} = label_full{find(contact.WrenchBase(:,i))};         %#ok<FNDSB>
        end
        % create an empty holonomic constraint object first with correct
        % name
        contact_constr = HolonomicConstraint(obj, ...
            [], contact.Name,...
            'LoadPath',load_path,...
            'ConstrLabel',{label},...
            'DerivativeOrder',2);
        % add as a set of holonomic constraints
        obj = addHolonomicConstraint(obj, contact_constr, load_path);
    end

    % the contact wrench input vector
    f_name = contact_constr.InputName;
    f = obj.Inputs.ConstraintWrench.(f_name);

    % if the friction coefficients are given, enforce friction cone
    % constraints
    if ~isempty(fric_coef)
        % get the friction cone constraint
        [friction_cone, fc_label, auxdata] = getFrictionCone(contact, f, fric_coef);

        % create an unilateral constraint object
        fc_cstr = UnilateralConstraint(obj, friction_cone,...
            ['fc' contact.Name], f_name, ...
            'ConstrLabel',{fc_label(:)'},...
            'AuxData',auxdata);
        % add as a set of unilateral constraitns
        obj = addUnilateralConstraint(obj, fc_cstr);
    end

    % if the contact geometry is given, enforce zero moment point
    % constraints
    if ~isempty(geometry)
        % get the friction cone constraint
        [zmp, zmp_label, auxdata] = getZMPConstraint(contact, f, geometry);

        % create an unilateral constraint object
        if ~isempty(zmp)
            zmp_cstr = UnilateralConstraint(obj, zmp,...
                ['zmp' contact.Name], f_name, ...
                'ConstrLabel',{zmp_label(:)'},...
                'AuxData',auxdata);
            % add as a set of unilateral constraitns
            obj = addUnilateralConstraint(obj, zmp_cstr);
        end
    end
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function model = SRModel
%SRMODEL Super-resolution model parameter structure.

    % Magnification factor used for super-resolution (default: 2)
    model.magFactor = 2.0;

    % Width of isotropic Gaussian PSF (default: 0.4)
    model.psfWidth = 0.4;

    % Motion parameters used for super-resolution. This can be either a
    % 2-D homography to model parametric motion or displacement vector
    % fields in case of local motion estimation (optical flow).
    %   - 2-D homography:
    %       Motion is given by a cell array of 2-D homographies where the
    %       k-th element models the motion for the k-th low-resolution
    %       frame:
    %           motionParams{1} = [h11 h12 h13; h21 h22 h23; h31 h32 h33]
    %           motionParams{2} = ...
    %   - Displacement fields:
    %       Motion is given by cell array of displacement fields where the
    %       k-th element models the motion for the k-th low-resolution.
    %       A displacement field is given by a structure with matrices to
    %       model displacements in x- and y-direction:
    %           v.vx = [...]    Displacements in x-direction (pixel-wise)
    %           v.vy = [...]    Displacements in y-direction (pixel-wise)
    %           motionParams{1} = v;
    %           motionParams{2} = ...
    model.motionParams = [];

    % Photometric parameters used for super-resolution (optional).
    % see also SRPhotometricParams
    model.photometricParams = [];

    % Error (noise) model assumed for the low-resolution observations
    % (optional). Can be either 'l2NormErrorModel' (least square, default)
    % in case  of Gaussian noise or 'l1NormErrorModel' (least absolute
    % deviation) in case of Laplacian noise. Additionally, user-defined
    % error models can be specified, e.g. a Lorentzian M-estimator:
    %   sigma = 0.05;
    %   lorentzian = @(r) deal(log(1 + 0.5*(r/sigma).^2), ... Function
    %                          (2*r) ./ (r.^2 + 2*sigma^2));  Gradient
    %   model.errorModel = lorentzian;
    model.errorModel = 'l2NormErrorModel';

    % Image prior used for regularization of super-resolved images.
    % see also SRPrior
    model.imagePrior = [];

    % Confidence map for the low-resolution observations (optional). A 3-D
    % array of size M x N x K can be used to weight K low-resolution frames
    % of size M x N pixel-wise.
    model.confidence = [];

    % Initial guess for the super-resolved image (optional). This can be
    % used as starting point for numerical optimization algorithms for
    % image reconstruction.
    model.SR = [];

author: thomas-koehler

repo: SupER


text:
function [c, s, kernel, iter] = deconvCa(y, kernel, smin, fit_gt, debug_on, sn, maxIter, theta, lambda)
%% deconvolve calcium traces to infer spike counts
%% inputs:
%   y:  1*T vector, observed calcium traces
%   kernel: struct variable with two fields {'fhandle', 'pars', nMax}. kernel
%   deterines the convolution kernel
%      fhandle: function handle, the function form of the response function
%       pars: p*1 vector, parameters for fhandle
%       nMax: scalar, maximum number of frames required by calcium indicators
%           to return resting states
%   smin: scalar, minimum number of nonzero spike count within each bin
%   fit_gt: true or false. iteratively fit gt or not
%   debug_on: play and save video of the whole procedure.
%   sn:     noise power
%   maxIter: scalar, maximum iterations for running OASIS, default (3)
%   theta: 1*T vecotr,  weight vector at each time point
%   lambda: scalar, tuning parameter

%% outputs:
%   c: 1*T vector, inferred calcium trace
%   s: 1*T vector, inferred spike train
%   kernel: convolution kernel
%   iters: number of iterations for updating irf

%% Author: Pengcheng Zhou, Carnegie Mellon University, 2016
% This work is based on one NIPS paper by Johannes Friedrich & Liam
% Paninski

%% input arugments
y = reshape(y, 1, []);  % convert data into row vector
if ~exist('sn', 'var') || isempty(sn)
    sn = get_noise_fft(y);
end
T = length(y);          % number of frames

% get convolution kernel
if ~exist('kernel', 'var') || isempty(kernel)
    kernel = create_kernel('exp2', ar2exp(estimate_time_constant(y,2)));
end
fhandle = kernel.fhandle;
nMax = kernel.nMax;

y = [y, zeros(1, nMax)];  % add few more elements for computation convenience

% threshold of the spike
if ~exist('smin', 'var') || isempty(smin)
    smin = 3*sn;               % min spike count
else
    smin = smin * sn;
end

% fit response function or not
if ~exist('fit_gt', 'var') || isempty(fit_gt)
    fit_gt = false;
else
    kernel.pars = ar2exp(estimate_time_constant(y,2));
end
thresh = 1e-3;

% debug mode
if ~exist('debug_on', 'var') || isempty(debug_on)
    debug_on = false;
end

% maximum iterations for running OASIS
if ~exist('maxIter', 'var') || isempty(maxIter)
    maxIter = 5;
end

% tuning parameter for enforcing sparsity
if ~exist('lambda', 'var') || isempty(lambda)
    lambda = 0;   % tuning parameter
end

% weight for each frame
if ~exist('theta', 'var')|| isempty(theta) || (length(theta)==1)     % weight vector
    theta = ones(1, T);
elseif length(theta)==T
    theta = reshape(theta, 1, T);
else
    disp('elements in theta should be equal to elements in y');
    return;
end

% tuning parameter for enforcing sparsity
if ~exist('lambda', 'var') || isempty(lambda)
    lambda = 0;   % tuning parameter
end

%% running OASIS
if debug_on
    figure('position', [1, 1 1500, 250]);
    plot(y); hold on;
    a = plot(y, 'r');
    avi_file = VideoWriter('example.avi');
    avi_file.open();
    xlim([1, T]);
end

f0 = inf;       % objective function
y0 = y;         % backup the raw trace
iter = 1;       % iteratiosn
for iter = 1:maxIter
    if debug_on
        title(sprintf('Iter %d, norm(residual) = %.4f', iter, f0));
    end
    % normalize gt to make its maximum value to be 1
    gt = fhandle(kernel.pars, 1:nMax); % 1*nMax, response function of calcium transients
    ind0 = 1; %correct the response function of the first event
    gt = [reshape(gt, 1, nMax), zeros(1, T)]; % add few more elements to gt for computation convenience

    gt_max = max(gt);
    gt = gt/gt_max;
    gt1 = [gt(ind0:end), zeros(1, ind0-1)]; % kernel for the first frame

    %% initialize values for the results
    s = zeros(size(y));     % spike count within each bin
    v = zeros(size(y));     % sum(theta_t*y_t) - lambda
    w = zeros(size(y));     % sum(theta_t^2 * g(t-ti+1)^2))
    pool_ti = zeros(size(y)); % first frame of each pool

    % initialize the first pool
    v(1) = theta(1)*y(1)*gt1(1) - lambda;
    w(1) = (theta(1)*gt1(1))^2;
    s(1) = max(0, v(1) / w(1));
    %     if s(1)<smin
    %         s(1) = 0;
    %     end
    pool_ti(1) = 1;
    tpre = 1;   % time of the last event

    %% start to initialize the next pool
    ii = 2;     % start to initialize the next pool
    t = 2;      % frame to begin
    frame_valid = true;
    while t<=T
        % plot results
        if debug_on
            c = zeros(1,T);
            c(1:T) = s(1) * gt1(1:T);
            for m=2:(ii-1)
                t0 = pool_ti(m);
                c(t0:T) = c(t0:T)+s(t0) * gt(1:(T-t0+1));
            end
            delete(a);
            a = plot(c(1:t), 'r', 'linewidth', 2);
            drawnow;
            temp = getframe();
            temp.cdata = imresize(temp.cdata, [200, 1500]);
            avi_file.writeVideo(temp);
        end
        %
        % find the spike count that minimizes the objective function
        if (t-tpre)>=nMax  % no influences from the previous events
            vi = theta(t)*y(t)*gt(1) - lambda;
        elseif tpre>1        % estimate vi and subtract the effect of the previous event
            vi = theta(t) * (y(t)-s(tpre)*gt(t-tpre+1)) * gt(1) - lambda;
        else   % special treatment with the first event
            vi = theta(t) * (y(t)-s(tpre)*gt1(t-tpre+1)) * gt(1) - lambda;
        end
        wi = (theta(t) * gt(1))^2;
        si = vi/wi;

        % check the violatoin of si
        if si>smin && frame_valid % no violation, create a new pool to save the result and move to the next new pool
            % peel off the previous event
            if (t-tpre)<nMax
                if tpre>1
                    y(tpre+(1:nMax)-1) = y(tpre+(1:nMax)-1) - s(tpre)*gt(1:nMax);
                else
                    y(tpre+(1:nMax)-1) = y(tpre+(1:nMax)-1) - s(tpre)*gt1(1:nMax);
                end
            end
            pool_ti(ii) = t;   % create a new pool
            tpre = t;
            v(t) = vi;
            w(t) = wi;
            s(t) = si;
            ii = ii+1;  % move to the next pool
            t = t+1;
        else  % with violation
            frame_valid = true; % allows the next frame to be valid
            if (t-tpre)>=nMax
                % ignore this frame directly and move to the next frame
                t = t + 1;
                continue;
            elseif tpre>1         % merge it to the current pool
                v(tpre) = v(tpre) + theta(t)*y(t)*gt(t-tpre+1);
                w(tpre) = w(tpre) + (theta(t)*gt(t-tpre+1))^2;
            else
                v(tpre) = v(tpre) + theta(t)*y(t)*gt1(t-tpre+1);
                w(tpre) = w(tpre) + (theta(t)*gt1(t-tpre+1));
            end
            s(tpre) = v(tpre)/w(tpre);  % update the current event

            % check the violation of the current pool
            if s(tpre)>smin      % the previous event is still avaiable
                t = t+1;
                continue;
            elseif ii==2  %the previous pool is not available anymore, but it's the first pool
                s(tpre) = max(0, s(tpre));
                t = t+1;
            else   % not available, then delete the current pool and force its first frame to be
                % event-free.
                t = tpre;  % go back to the first frame of the current pull

                frame_valid = false; % force the frame t to be invalid
                ii = ii-1;
                tpre = pool_ti(ii-1);

                % add back the signal of the previous pool
                if (t-tpre)>nMax
                    t = t+1;
                    frame_valid = true;
                    continue;
                elseif tpre>1
                    y(tpre+(1:nMax)-1) = y(tpre+(1:nMax)-1) + s(tpre)*gt(1:nMax);
                else
                    y(tpre+(1:nMax)-1) = y(tpre+(1:nMax)-1) + s(tpre)*gt1(1:nMax);
                end
            end
        end
    end
    %% collect the results
    pool_ti(ii:end) = [];
    temp = s(pool_ti);
    s = zeros(1, T);
    s(pool_ti) = temp;
    temp = s;
    temp(1) = 0;
    c = conv(temp, gt(1:nMax));
    c = c(1:T) + s(1)*gt1(1:T);
    f1 = norm(y(1:T)-c, 2);
    s = s*sum(gt);
    if debug_on
        delete(a);
        a = plot(c, 'r');
        drawnow();
        temp = getframe();
        temp.cdata = imresize(temp.cdata, [200, 1500]);
        avi_file.writeVideo(temp);
    end

    %% break the while loop
    if ~fit_gt  % don't iterate gt
        break;
    elseif (f0-f1)/f0 <= thresh % improvement is small, stop
        break;
    else  % move to the next iteration
        y = y0;
        if strcmpi(kernel.type, 'exp2')
            [kernel, ~] = update_kernel_exp2(y(1:T), s, kernel);
        else
            break;
        end
        f0 = f1;
        iter = iter + 1;
        %             return;
    end
end

if debug_on
    avi_file.close();
    close(gcf);
end
end

%% estimate the noise power
function [sn,psdx,ff] = get_noise_fft(Y,options)
% Written by:
% Eftychios A. Pnevmatikakis, Simons Foundation, 2015
% with minor adaption by Pengcheng Zhou, Carnegie Mellon University, 2015
options.noise_range = [.25, .5];
range_ff = options.noise_range;
options.noise_method = 'logmexp';
method = options.noise_method;
options.block_size = [64, 64];
block_size = options.block_size;
options.split_data = false;
split_data = options.split_data;
options.max_timesteps = 3000;

dims = ndims(Y);
sizY = size(Y);
N = min(sizY(end),options.max_timesteps);
if N < sizY(end)
    %Y = reshape(Y,prod(sizY(1:end-1)),[]);
    Y(prod(sizY(1:end-1))*N+1:end) = [];
    Y = reshape(Y,[sizY(1:end-1),N]);
end

Fs = 1;
ff = 0:Fs/N:Fs/2;
indf=ff>range_ff(1);
indf(ff>range_ff(2))=0;
if dims > 1
    d = prod(sizY(1:dims-1));
    Y = reshape(Y,d,N);
    Nb = prod(block_size);
    SN = cell(ceil(d/Nb),1);
    PSDX = cell(ceil(d/Nb),1);
    if ~split_data
        for ind = 1:ceil(d/Nb);
            xdft = fft(Y((ind-1)*Nb+1:min(ind*Nb,d),:),[],2);
            xdft = xdft(:,1: floor(N/2)+1); % FN: floor added.
            psdx = (1/(Fs*N)) * abs(xdft).^2;
            psdx(:,2:end-1) = 2*psdx(:,2:end-1);
            %SN{ind} = mean_psd(psdx(:,indf),method);
            switch method
                case 'mean'
                    SN{ind}=sqrt(mean(psdx(:,indf)/2,2));
                case 'median'
                    SN{ind}=sqrt(median(psdx(:,indf)/2),2);
                case 'logmexp'
                    SN{ind} = sqrt(exp(mean(log(psdx(:,indf)/2),2)));
            end
            PSDX{ind} = psdx;
        end
    else
        nc = ceil(d/Nb);
        Yc = mat2cell(Y,[Nb*ones(nc-1,1);d-(nc-1)*Nb],N);
        parfor ind = 1:ceil(d/Nb);
            xdft = fft(Yc{ind},[],2);
            xdft = xdft(:,1:floor(N/2)+1);
            psdx = (1/(Fs*N)) * abs(xdft).^2;
            psdx(:,2:end-1) = 2*psdx(:,2:end-1);
            Yc{ind} = [];
            switch method
                case 'mean'
                    SN{ind}=sqrt(mean(psdx(:,indf)/2,2));
                case 'median'
                    SN{ind}=sqrt(median(psdx(:,indf)/2),2);
                case 'logmexp'
                    SN{ind} = sqrt(exp(mean(log(psdx(:,indf)/2),2)));
            end

        end
    end
    sn = cell2mat(SN);
else
    xdft = fft(Y);
    xdft = xdft(:,1:floor(N/2)+1);
    psdx = (1/(Fs*N)) * abs(xdft).^2;
    psdx(:,2:end-1) = 2*psdx(:,2:end-1);
    switch method
        case 'mean'
            sn = sqrt(mean(psdx(:,indf)/2,2));
        case 'median'
            sn = sqrt(median(psdx(:,indf)/2),2);
        case 'logmexp'
            sn = sqrt(exp(mean(log(psdx(:,indf)/2),2)));
    end
end
psdx = cell2mat(PSDX);
if dims > 2
    sn = reshape(sn,sizY(1:dims-1));
end
end

author: zhoupc


text:
function [ xstar, seed ] = rk2_ti_step ( x, t, h, q, fi, gi, seed )

%*****************************************************************************80
%
%% RK2_TI_STEP takes one step of a stochastic Runge Kutta scheme.
%
%  Discussion:
%
%    The Runge-Kutta scheme is second-order, and suitable for time-invariant
%    systems.
%
%    d/dx X(t,xsi) = F ( X(t,xsi) ) + G ( X(t,xsi) ) * w(t,xsi)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Jeremy Kasdin,
%    Runge-Kutta algorithm for the numerical integration of
%    stochastic differential equations,
%    Journal of Guidance, Control, and Dynamics,
%    Volume 18, Number 1, January-February 1995, pages 114-120.
%
%    Jeremy Kasdin,
%    Discrete Simulation of Colored Noise and Stochastic Processes
%    and 1/f^a Power Law Noise Generation,
%    Proceedings of the IEEE,
%    Volume 83, Number 5, 1995, pages 802-827.
%
%  Parameters:
%
%    Input, real X, the value at the current time.
%
%    Input, real T, the current time.
%
%    Input, real H, the time step.
%
%    Input, real Q, the spectral density of the input white noise.
%
%    Input, external real FI, the name of the deterministic
%    right hand side function.
%
%    Input, external real GI, the name of the stochastic
%    right hand side function.
%
%    Input/output, integer SEED, a seed for the random
%    number generator.
%
%    Output, real XSTAR, the value at time T+H.
%
  a21 = 1.0;
  a31 = 0.5;
  a32 = 0.5;

  q1 = 2.0;
  q2 = 2.0;

  t1 = t;
  x1 = x;
  [ n1, seed ] = r8_normal_01 ( seed );
  w1 = n1 * sqrt ( q1 * q / h );
  k1 = h * fi ( x1 ) + h * gi ( x1 ) * w1;

  t2 = t1 + a21 * h;
  x2 = x1 + a21 * k1;
  [ n2, seed ] = r8_normal_01 ( seed );
  w2 = n2 * sqrt ( q2 * q / h );
  k2 = h * fi ( x2 ) + h * gi ( x2 ) * w2;

  xstar = x1 + a31 * k1 + a32 * k2;

  return
end


text:
function obj = setRotationMatrix(obj,r)
    % set the rotation matrix (R) of the frame w.r.t to the
    % reference frame
    %
    % Parameters:
    % r: the rotation matrix or Euler angles (in radian) @type
    % rowvec|matrix

    if isvector(r) && length(r) == 3
        %         rpy = rad2deg(r);
        obj.R = rotz(r(3)) * roty(r(2)) * rotx(r(1));
    elseif all(size(r)==[3,3])
        if isnumeric(r)
            assert(abs(det(r))-1 <= 1e-6,...
                'The determinant of the rotation matrix must equal 1.');
        end
        obj.R = r;
    end
    % remove small numbers generated from the rotation matrix
    if isnumeric(r)
        obj.R = roundn(obj.R,-6);
    end
    % update the homogeneous transformation matrix
    obj.computeHomogeneousTransform();
end

function R = rotx(alpha)
%rotz  rotate around X by ALPHA
%
%	R = rotx(alpha)
%
% See also: roty, rotz
% Author: Jake Reher: jreher@caltech.edu

R = [1 0 0; ...
     0 cos(alpha) -sin(alpha); ...
     0 sin(alpha)  cos(alpha)];

end

function R = roty(alpha)
%roty  rotate around Y by ALPHA
%
%	R = roty(alpha)
%
% See also: rotx, rotz
% Author: Jake Reher: jreher@caltech.edu

R = [cos(alpha)  0 sin(alpha); ...
             0   1          0; ...
    -sin(alpha)  0 cos(alpha)];

end

function R = rotz(alpha)
%rotz  rotate around Z by ALPHA
%
%	R = rotz(alpha)
%
% See also: ROTX, ROTY, ROT, POS.
% Author: Jake Reher: jreher@caltech.edus

R = [cos(alpha) -sin(alpha) 0; ...
     sin(alpha)  cos(alpha) 0; ...
              0           0 1];

end


text:
%compute dtailinflmag
function [data,units]=compute_dtailinflmag(trx,n)

%
larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
dtailinflmag=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    dtailinflmag{1,i}=(trx(larva).tailinflmag(2:end)-trx(larva).tailinflmag(1:end-1))./trx(larva).dt;
end

units=parseunits('mm/s');
data=dtailinflmag;

author: kristinbranson

repo: JAABA


text:
function y=expSlope2(x)
%%EXPSLOPE2  Evaluate the function 2*(exp(x)-1-x)/x^2 avoiding problems at
%            the point x=0, where the limit equals 1. This function arises
%            in range enclosure algorithms for interval arithmetic.
%
%INPUTS: x A vector of matrix of real values at which the function should
%          be evaluated.
%
%OUTPUTS: y The value(s) of the function 2*(exp(x)-1-x)/x^2 at the point(s)
%           in x.
%
%For values of x having magnitude over 1e-2 or for complex values, the
%function is evaluated as written. For smaller values, a Taylor series
%expansion is used to avoid the singularity near the origin. The Taylor
%series expansion is simple, because the value of the nth derivative at 0
%is just 1/(n+1).
%
%This is one of the slope function in Table 10.6 of Section 10.6.2 of [1].
%This file implements the function for non-intervals. See the Interval
%class for an implementation for Interval arithmetic.
%
%REFERENCES:
%[1] IEEE Standard for Interval Arithmetic," in IEEE Std 1788-2015,
%    pp.1-97, June 30 2015.
%
%November 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Allocate space
y=zeros(size(x));

%For large enough values, use the standard function.
sel=abs(x)>1e-2|~isreal(x);
x1=x(sel);
y(sel)=2*(expm1(x1)-x1)./x1.^2;

%For smaller magnitude values, use a sixth-order Taylor series expansion.
%The relative error at 10^-2 is on the order of 1e-20, which is less than
%the precision of a double floating point number at that value. The series
%is given in Horner form.
x1=x(~sel);
y(~sel)=1+x1.*(1/3+x1.*(1/12+x1.*(1/60+x1.*(1/360+(1/2520+x1./20160).*x1))));

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [coeff,RMSout]=FitLaw_GA(X,y,fh,Int,LB,UB)
%
% function [coeff,RMSout]=FitLaw_GA(X,y,fh,Int,LB,UB)
%
% The script fits a given general analytical law using Genetic Algorithms.
% Such law can depend on several input variables and various parameters,
% but it must have only one output, to be matched to the input variable y.
% The optimum coefficients are derived using an optimization based on the
% Genetic Algorithm.
%
% INPUT
%   X = matrix containing the values of the indipendent variables of the
%       function. Each column refers to a variable
%   y = target values to be matched by the function (row vector)
%   fh = handle of the function to be fitted to y and whose parameters have
%        to be optimized
%   Int = initial interval of the parameters to be optmized. Column i of
%         the matrix refer to parameter i
%   LB,UB = vectors defining the lower and upper bounds never to be
%           exceeded by the parameters. Position i of the vectors refer to
%           parameter i.
%
% OUTPUT
%   coeff = parameters providing the best fitting
%   RMSout = lowest RMS of the fitting error
%
% EXAMPLE OF USAGE
%   Suppose one wants to fit to y the function p=x.*(1-exp(-a.*w+b.*z)),
%   where (x,w,z) are the indipendent variables and (a,b) are the
%   parameters to be determined. It is first necessary to define the handle
%   to the function p before calling the optimization function. This is
%   done as follows:
%
%      fh=@(X,par) X(:,1)'.*(1-exp(-par(1).*X(:,2)'+par(2).*X(:,3)'))
%
%   Afterwards, the callback of the function would be:
%
%   [coeff,RMSout]=FitLaw_GA([x' w' z'],y,fh,[-5 -5; 5 5],[-10 10],[10 10])
%
%   Note that since fh has been defined considering that each variable is
%   assigned to each column of the matrix X, then X in the callback has to
%   be in such form (as in the example). y is a row vector containing the
%   target values of the function to be fitted. [-5 -5; 5 5] defines the
%   initial intervals in which the parameters are sampled when the
%   algorithm starts. Specifically, column i of the matrix Int contains the
%   lower and upper limits for parameter i. However, as the algorithm
%   proceeds, the parameters may exceed such intervals. That is why lower
%   and upper bounds can be introduced to limit the acceptable values of
%   the parameters (in this example, both paramaters are such that
%   -10 < a,b < 10).
%
%   To improve the fitting power, increase the variables PopSize and Iter
%   defined below in the cell GENETIC ALGORITHM OPTIONS. Specifically,
%   a high value of PopSize provides a higher probability to avoid local
%   minima in the optimization procedure, while a high value of Iter
%   permits a proper convergence of the algorithm.
%   Obviously, high values of those variables imply longer calculation
%   times.
%
% By: L.Luini

%% GENETIC ALGORITHM OPTIONS
PopSize=200;   % population size (default==200)
Iter=100;   % number of iterations of the algorithm (default==100)
MigrInt=Iter/10;   % migration interval
Fig=1;   % if Fig==1, plot optimization results

%% RUN THE GENETIC ALGORITHM
numPar=size(Int,2);
options=gaoptimset('PlotFcns',{@gaplotbestf,@gaplotbestindiv}, ...
    'PopInitRange',Int, ...
    'PopulationSize',PopSize,'MigrationInterval',MigrInt, ...
    'StallGenLimit',Inf,'StallTimeLimit',Inf,'Generations',Iter);
coeff=ga(@FitFcn,numPar,[],[],[],[],LB,UB,[],options);
est=fh(X,coeff);
error=100.*(est-y)./y;
RMSout=sqrt(mean(error).^2+std(error).^2);

% plot results if the function depends only on 1 variable
if Fig==1&&numPar==1
    figure
    plot(X,y,'b','LineWidth',1.5)
    hold on;grid on;
    plot(X,est,'r','LineWidth',1.5)
    legend('Input data','Fitting curve')
    grid
end

    function fitness=FitFcn(param)

        % Fitness function definition: RMS of the error
        est=fh(X,param);
        error=100.*(est-y)./y;   % definition of the error
        fitness=sqrt(mean(error).^2+std(error).^2);
    end

end


text:
function setdt (delt)

% this function allows for the specification of the delta-t value
% (delta-t = tt - ut1) to be used in the calculation of sidereal
% time and the terrestrial-to-celestial transformation.  it allows
% these calculations to be performed, correctly, using ut1 as the
% time argument for the earth rotation angle and tdb as the time
% argument for the precession and nutation components.  this
% function, if used, should be called before any function
% related to earth rotation (e.g., sidtim or tercel) for a given
% date.  the value of delta-t specified here will be used until
% explicitly changed.

%  delt = value of delta-t (tt-ut1) in seconds (in)

% note 1:  the computed value of sidereal time, and the equivalent
% earth orientation angles, are relatively insensitive to the value
% of delta-t: up to only ~3 microarcseconds per second of delta-t.
% therefore, for many applications, this function either need not
% be called at all, or can be called just once for a wide range of
% dates (e.g., a year).  if this call is not used, a default
% value of delta-t of 64 seconds is used, which is appropriate to
% 2000.0.

% note 2:  the input time arguments to sidtim and tercel (tjdh and
% tjdl) are expressed in ut1 regardless of whether this call is
% used.

% ported from NOVAS 3.1

%%%%%%%%%%%%%%%%%%%%%%%

global dt

dt = delt / 86400.0d0;

end

author: goGPS-Project


text:
function f = real(f)
%REAL   Real part of a CHEBTECH.
%   REAL(F) is the real part of F.
%
% See also ISREAL, IMAG, CONJ.

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

% Compute the real part of the coefficients:
f.coeffs = real(f.coeffs);

if ( ~any(f.coeffs(:)) )
    % Input was imaginary, so output a zero CHEBTECH:
    f = f.make(zeros(1, size(f.coeffs, 2)));
    f.ishappy = 1;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [x,zo,xs]=v_estnoisem(yf,tz,pp)
%V_ESTNOISEM - estimate noise spectrum using minimum statistics
%
% Usage:    ninc=round(0.016*fs);   % frame increment [fs=sample frequency]
%           ovf=2;                  % overlap factor
%           f=v_rfft(v_enframe(s,v_windows(2,ovf*ninc,'l'),ninc),ovf*ninc,2);
%           f=f.*conj(f);           % convert to power spectrum
%           x=v_estnoisem(f,ninc/fs); % estimate the noise power spectrum
%
% Inputs:
%   yf      input power spectra (one row per frame)
%   tz      frame increment in seconds
%           Alternatively, the input state from a previous call (see below)
%   pp      algorithm parameters [optional]
%
% Outputs:
%   x       estimated noise power spectra (one row per frame)
%   zo      output state
%   xs      estimated std error of x (one row per frame)
%           xs seems often to be an underestimate by a factor of 2 or 3
%
% The algorithm parameters are defined in reference [1] from which equation
% numbers are given in parentheses. They are as follows:
%
%        pp.taca      % (11): smoothing time constant for alpha_c [0.0449 seconds]
%        pp.tamax     % (3): max smoothing time constant [0.392 seconds]
%        pp.taminh    % (3): min smoothing time constant (upper limit) [0.0133 seconds]
%        pp.tpfall    % (12): time constant for P to fall [0.064 seconds]
%        pp.tbmax     % (20): max smoothing time constant [0.0717 seconds]
%        pp.qeqmin    % (23): minimum value of Qeq [2]
%        pp.qeqmax    % max value of Qeq per frame [14]
%        pp.av        % (23)+13 lines: fudge factor for bc calculation  [2.12]
%        pp.td        % time to take minimum over [1.536 seconds]
%        pp.nu        % number of subwindows to use [3]
%        pp.qith      % Q-inverse thresholds to select maximum noise slope [0.03 0.05 0.06 Inf ]
%        pp.nsmdb     % corresponding noise slope thresholds in dB/second   [47 31.4 15.7 4.1]
%
% Example use:      y=v_enframe(s,w,ni);                  % divide speech signal s(n) into
%                                                       % overlapping frames using window w(n)
%                   yf=v_rfft(y,nf,2);                    % take fourier transform
%                   dp=v_estnoisem(yf.*conj(yf),tinc);    % estimate the noise
%
% If convenient, you can call v_estnoisem in chunks of arbitrary size. Thus the following are equivalent:
%
%                   (a) dp=v_estnoisem(yp(1:300),tinc);
%
%                   (b) [dp(1:100),z]=v_estnoisem(yp(1:100),tinc);
%                       [dp(101:200),z]=v_estnoisem(yp(101:200),z);
%                       [dp(201:300),z]=v_estnoisem(yp(201:300),z);

% This is intended to be a precise implementation of [1] with Table III
% replaced by the updated table 5 from [2]. The only deliberate algorithm
% change is the introduction of a minimum value for 1/Qeq in equation (23).
% This change only affects the first few frames and improves the
% convergence of the algorithm. A minor improveemnt was reported in [3] but
% this has not yet been included.
%
% Refs:
%    [1] Rainer Martin.
%        Noise power spectral density estimation based on optimal smoothing and minimum statistics.
%        IEEE Trans. Speech and Audio Processing, 9(5):504-512, July 2001.
%    [2] Rainer Martin.
%        Bias compensation methods for minimum statistics noise power spectral density estimation
%        Signal Processing, 2006, 86, 1215-1229
%    [3] Dirk Mauler and Rainer Martin
%        Noise power spectral density estimation on highly correlated data
%        Proc IWAENC, 2006

%	   Copyright (C) Mike Brookes 2008
%      Version: $Id: v_estnoisem.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nr,nrf]=size(yf);          % number of frames and freq bins
x=zeros(nr,nrf);            % initialize output arrays
xs=zeros(nr,nrf);           % will hold std error in the future
if isempty(yf) && isstruct(tz)             % no real data
    zo=tz;              % just keep the same state
else
    if isstruct(tz)       % take parameters from a previous call
        nrcum=tz.nrcum;
        p=tz.p;          % smoothed power spectrum
        ac=tz.ac;               % correction factor (9)
        sn2=tz.sn2;              % estimated noise power
        pb=tz.pb;               % smoothed noisy speech power (20)
        pb2=tz.pb2;
        pminu=tz.pminu;
        actmin=tz.actmin;   % Running minimum estimate
        actminsub=tz.actminsub;           % sub-window minimum estimate
        subwc=tz.subwc;                   % force a buffer switch on first loop
        actbuf=tz.actbuf;  % buffer to store subwindow minima
        ibuf=tz.ibuf;
        lminflag=tz.lminflag;      % flag to remember local minimum
        tinc=tz.tinc;     % frame increment
        qq=tz.qq;         % parameter structure
    else
        tinc = tz;          % second argument is frame increment
        nrcum=0;            % no frames so far
        % default algorithm constants

        qq.taca=0.0449;    % smoothing time constant for alpha_c = -tinc/log(0.7) in equ (11)
        qq.tamax=0.392;    % max smoothing time constant in (3) = -tinc/log(0.96)
        qq.taminh=0.0133;    % min smoothing time constant (upper limit) in (3) = -tinc/log(0.3)
        qq.tpfall=0.064;   % time constant for P to fall (12)
        qq.tbmax=0.0717;   % max smoothing time constant in (20) = -tinc/log(0.8)
        qq.qeqmin=2;       % minimum value of Qeq (23)
        qq.qeqmax=14;      % max value of Qeq per frame
        qq.av=2.12;             % fudge factor for bc calculation (23 + 13 lines)
        qq.td=1.536;       % time to take minimum over
        qq.nu=8;           % number of subwindows
        qq.qith=[0.03 0.05 0.06 Inf]; % noise slope thresholds in dB/s
        qq.nsmdb=[47 31.4 15.7 4.1];

        if nargin>=3 && ~isempty(pp)
            qqn=fieldnames(qq);
            for i=1:length(qqn)
                if isfield(pp,qqn{i})
                    qq.(qqn{i})=pp.(qqn{i});
                end
            end
        end
    end

    % unpack parameter structure

    taca=qq.taca;    % smoothing time constant for alpha_c = -tinc/log(0.7) in equ (11)
    tamax=qq.tamax;    % max smoothing time constant in (3) = -tinc/log(0.96)
    taminh=qq.taminh;    % min smoothing time constant (upper limit) in (3) = -tinc/log(0.3)
    tpfall=qq.tpfall;   % time constant for P to fall (12)
    tbmax=qq.tbmax;   % max smoothing time constant in (20) = -tinc/log(0.8)
    qeqmin=qq.qeqmin;       % minimum value of Qeq (23)
    qeqmax=qq.qeqmax;      % max value of Qeq per frame
    av=qq.av;             % fudge factor for bc calculation (23 + 13 lines)
    td=qq.td;       % time to take minimum over
    nu=qq.nu;           % number of subwindows
    qith=qq.qith; % noise slope thresholds in dB/s
    nsmdb=qq.nsmdb;   % maximum permitted +ve noise slope in dB/s

    % derived algorithm constants

    aca=exp(-tinc/taca); % smoothing constant for alpha_c in equ (11) = 0.7
    acmax=aca;          % min value of alpha_c = 0.7 in equ (11) also = 0.7
    amax=exp(-tinc/tamax); % max smoothing constant in (3) = 0.96
    aminh=exp(-tinc/taminh); % min smoothing constant (upper limit) in (3) = 0.3
    bmax=exp(-tinc/tbmax); % max smoothing constant in (20) = 0.8
    snrexp = -tinc/tpfall;
    nv=round(td/(tinc*nu));    % length of each subwindow in frames
    if nv<4            % algorithm doesn't work for miniscule frames
        nv=4;
        nu=max(round(td/(tinc*nv)),1);
    end
    nd=nu*nv;           % length of total window in frames
    [md,hd]=mhvals(nd); % calculate the constants M(D) and H(D) from Table III
    [mv,hv]=mhvals(nv); % calculate the constants M(D) and H(D) from Table III
    nsms=10.^(nsmdb*nv*tinc/10);  % [8 4 2 1.2] in paper
    qeqimax=1/qeqmin;  % maximum value of Qeq inverse (23)
    qeqimin=1/qeqmax; % minumum value of Qeq per frame inverse

    if isempty(yf)      % provide dummy initialization
        ac=1;               % correction factor (9)
        subwc=nv;                   % force a buffer switch on first loop
        ibuf=0;
        p=x;          % smoothed power spectrum
        sn2=p;              % estimated noise power
        pb=p;               % smoothed noisy speech power (20)
        pb2=pb.^2;
        pminu=p;
        actmin=repmat(Inf,1,nrf);   % Running minimum estimate
        actminsub=actmin;           % sub-window minimum estimate
        actbuf=repmat(Inf,nu,nrf);  % buffer to store subwindow minima
        lminflag=zeros(1,nrf);      % flag to remember local minimum
    else

        if ~nrcum       % initialize values for first frame
            p=yf(1,:);          % smoothed power spectrum
            ac=1;               % correction factor (9)
            sn2=p;              % estimated noise power
            pb=p;               % smoothed noisy speech power (20)
            pb2=pb.^2;
            pminu=p;
            actmin=repmat(Inf,1,nrf);   % Running minimum estimate
            actminsub=actmin;           % sub-window minimum estimate
            subwc=nv;                   % force a buffer switch on first loop
            actbuf=repmat(Inf,nu,nrf);  % buffer to store subwindow minima
            ibuf=0;
            lminflag=zeros(1,nrf);      % flag to remember local minimum
        end

        % loop for each frame

        for t=1:nr              % we use t instead of lambda in the paper
            yft=yf(t,:);        % noise speech power spectrum
            acb=(1+(sum(p)./sum(yft)-1).^2).^(-1);  % alpha_c-bar(t)  (9)
            ac=aca*ac+(1-aca)*max(acb,acmax);       % alpha_c(t)  (10)
            ah=amax*ac.*(1+(p./sn2-1).^2).^(-1);    % alpha_hat: smoothing factor per frequency (11)
            snr=sum(p)/sum(sn2);
            ah=max(ah,min(aminh,snr^snrexp));       % lower limit for alpha_hat (12)

            p=ah.*p+(1-ah).*yft;            % smoothed noisy speech power (3)
            b=min(ah.^2,bmax);              % smoothing constant for estimating periodogram variance (22 + 2 lines)
            pb=b.*pb + (1-b).*p;            % smoothed periodogram (20)
            pb2=b.*pb2 + (1-b).*p.^2;     	% smoothed periodogram squared (21)

            qeqi=max(min((pb2-pb.^2)./(2*sn2.^2),qeqimax),qeqimin/(t+nrcum));   % Qeq inverse (23)
            qiav=sum(qeqi)/nrf;             % Average over all frequencies (23+12 lines) (ignore non-duplication of DC and nyquist terms)
            bc=1+av*sqrt(qiav);             % bias correction factor (23+11 lines)
            bmind=1+2*(nd-1)*(1-md)./(qeqi.^(-1)-2*md);      % we use the simplified form (17) instead of (15)
            bminv=1+2*(nv-1)*(1-mv)./(qeqi.^(-1)-2*mv);      % same expression but for sub windows
            kmod=bc*p.*bmind<actmin;        % Frequency mask for new minimum
            if any(kmod)
                actmin(kmod)=bc*p(kmod).*bmind(kmod);
                actminsub(kmod)=bc*p(kmod).*bminv(kmod);
            end
            if subwc>1 && subwc<nv              % middle of buffer - allow a local minimum
                lminflag=lminflag | kmod;    	% potential local minimum frequency bins
                pminu=min(actminsub,pminu);
                sn2=pminu;
            else
                if subwc>=nv                    % end of buffer - do a buffer switch
                    ibuf=1+rem(ibuf,nu);     	% increment actbuf storage pointer
                    actbuf(ibuf,:)=actmin;    	% save sub-window minimum
                    pminu=min(actbuf,[],1);
                    i=find(qiav<qith);
                    nsm=nsms(i(1));          	% noise slope max
                    lmin=lminflag & ~kmod & actminsub<nsm*pminu & actminsub>pminu;
                    if any(lmin)
                        pminu(lmin)=actminsub(lmin);
                        actbuf(:,lmin)=repmat(pminu(lmin),nu,1);
                    end
                    lminflag(:)=0;
                    actmin(:)=Inf;
                    subwc=0;
                end
            end
            subwc=subwc+1;
            x(t,:)=sn2;
            qisq=sqrt(qeqi);
            % empirical formula for standard error based on Fig 15 of [2]
            xs(t,:)=sn2.*sqrt(0.266*(nd+100*qisq).*qisq/(1+0.005*nd+6/nd)./(0.5*qeqi.^(-1)+nd-1));
        end
    end
    if nargout>1    % we need to store the state for next time
        zo.nrcum=nrcum+nr;      % number of frames so far
        zo.p=p;          % smoothed power spectrum
        zo.ac=ac;               % correction factor (9)
        zo.sn2=sn2;              % estimated noise power
        zo.pb=pb;               % smoothed noisy speech power (20)
        zo.pb2=pb2;
        zo.pminu=pminu;
        zo.actmin=actmin;   % Running minimum estimate
        zo.actminsub=actminsub;           % sub-window minimum estimate
        zo.subwc=subwc;                   % force a buffer switch on first loop
        zo.actbuf=actbuf;  % buffer to store subwindow minima
        zo.ibuf=ibuf;
        zo.lminflag=lminflag;      % flag to remember local minimum
        zo.tinc=tinc;     % must be the last one
        zo.qq=qq;
    end
    if ~nargout
        clf;
        subplot(212);
        plot((1:nr)*tinc,10*log10([sum(yf,2) sum(x,2)]))
        ylabel('Frame Energy (dB)');
        xlabel(sprintf('Time (s)   [%d ms frame incr]',round(tinc*1000)));
        v_axisenlarge([-1 -1.05]);
        legend('input','noise','Location','Best');
        subplot(211);
        plot(1:nrf,10*log10([sum(yf,1)'/nr sum(x,1)'/nr]))
        ylabel('Power (dB)');
        xlabel('Frequency bin');
        v_axisenlarge([-1 -1.05]);
        legend('input','noise','Location','Best');
    end
end

function [m,h,d]=mhvals(d)
% Values are taken from Table 5 in [2]
%[2] R. Martin,"Bias compensation methods for minimum statistics noise power
%               spectral density estimation", Signal Processing Vol 86, pp1215-1229, 2006.

% approx: plot(d.^(-0.5),[m 1-d.^(-0.5)],'x-'), plot(d.^0.5,h,'x-')
persistent dmh
if isempty(dmh)
    dmh=[
        1   0       0;
        2   0.26    0.15;
        5   0.48    0.48;
        8   0.58    0.78;
        10  0.61    0.98;
        15  0.668   1.55;
        20  0.705   2;
        30  0.762   2.3;
        40  0.8     2.52;
        60  0.841   3.1;
        80  0.865   3.38;
        120 0.89    4.15;
        140 0.9     4.35;
        160 0.91    4.25;
        180 0.92    3.9;
        220 0.93    4.1;
        260 0.935   4.7;
        300 0.94    5];
end

if nargin>=1
    i=find(d<=dmh(:,1));
    if isempty(i)
        i=size(dmh,1);
        j=i;
    else
        i=i(1);
        j=i-1;
    end
    if d==dmh(i,1)
        m=dmh(i,2);
        h=dmh(i,3);
    else
        qj=sqrt(dmh(i-1,1));    % interpolate using sqrt(d)
        qi=sqrt(dmh(i,1));
        q=sqrt(d);
        h=dmh(i,3)+(q-qi)*(dmh(j,3)-dmh(i,3))/(qj-qi);
        m=dmh(i,2)+(qi*qj/q-qj)*(dmh(j,2)-dmh(i,2))/(qi-qj);
    end
else
    d=dmh(:,1);
    m=dmh(:,2);
    h=dmh(:,3);
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function rgb = black

% returns a predefined color as [red green blue] values
%
% skin_surface         = [255 213 119]/255;
% outer_skull_surface  = [140  85  85]/255;
% inner_skull_surface  = [202 100 100]/255;
% cortex = [255 213 119]/255;
% black  = [0   0   0  ]/255;
% white  = [255 255 255]/255;
% red    = [255 0   0  ]/255;
% green  = [0   192 0  ]/255;
% blue   = [0   0   255]/255;
% yellow = [255 255 0  ]/255;
% cortex_light = [199 194 169]/255;
% cortex_dark  = [100 97 85]/255;

rgb = [0   0   0  ]/255;

author: spm