text:
function W_nullpwd = sphNullformer_pwd(order, beam_dirs)
%SPHNULLFORMER_PWD Beamweights for a PWD beamformer at a specified direction,
%with nulls at others
%
%   For a set of K directions, computes the beamforming weights that creates
%   a PWD beamformer at each of the direction in the set, while placing
%   nulls at the rest.
%
%   Inputs:
%       order:  order of SH signals
%       beam_dirs:  Kx2 [azi elev] directions
%
%   Outputs:
%       W_nullpwd:  (order+1)^2xK matrix of beamweights, one set per
%           column for each beamforming direction
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SPHNULLFORMER_PWD.M - 5/10/2016
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    beam_dirs2 = [beam_dirs(:,1) pi/2-beam_dirs(:,2)]; % convert from azi-elev to azi-incl

    % steering vectors in the SHD
    Y_nullpwd = getSH(order, beam_dirs2, 'real');

    % beamform at each direction and nullform at the rest
    W_nullpwd = pinv(Y_nullpwd);
end

author: polarch

repo: Spherical-Array-Processing


text:
function signal = warmStartWithBadChRemoved(signal, pipe)

channel_mask = true(signal.nbchan, 1);
if any(strcmp(pipe.fltorder, 'flt_selchans'))
    mask = ismember({signal.chanlocs.labels}, pipe.pselchans.channels);
    if ~pipe.pselchans.remove_selection
        mask = ~mask; end
    channel_mask(mask) = 0;
end

%{
%% Bad channel detection
% define parameters
min_corr = 0.5;
ignored_quantile = 0.1;
window_len = 2;
max_broken_time = 0.4;
linenoise_aware = true;
rereferenced = false;
protect_channels = [];

% flag channels
if ~exist('removed_channel_mask','var')
    if max_broken_time > 0 && max_broken_time < 1  %#ok<*NODEF>
        max_broken_time = size(signal.data,2)*max_broken_time;
    else
        max_broken_time = signal.srate*max_broken_time;
    end

    [nChs,nPts] = size(signal.data);
    window_len = window_len*signal.srate;
    wnd = 0:window_len-1;
    offsets = round(1:window_len:nPts-window_len);
    W = length(offsets);
    retained = 1:(nChs-ceil(nChs*ignored_quantile));

    % optionally ignore both 50 and 60 Hz spectral components...
    if linenoise_aware && signal.srate > 110
        if signal.srate <= 130
            B = design_fir(500,[2*[0 45 50 55]/signal.srate 1],[1 1 0 1 1]);
        else
            B = design_fir(500,[2*[0 45 50 55 60 65]/signal.srate 1],[1 1 0 1 0 1 1]);
        end
        for c=signal.nbchan:-1:1
            X(:,c) = filtfilt_fast(B,1,signal.data(c,:)'); end
    else
        X = signal.data';
    end

    % optionally subtract common reference from data
    if rereferenced
        X = bsxfun(@minus,X,mean(X,2)); end

    % for each window, flag channels with too low correlation to any other channel (outside the
    % ignored quantile)
    flagged = zeros(nChs,W);
    for o=1:W
        sortcc = sort(abs(corrcoef(X(offsets(o)+wnd,:))));
        flagged(:,o) = all(sortcc(retained,:) < min_corr);
    end

    % mark all channels for removal which have more flagged samples than the maximum number of
    % ignored samples
    removed_channel_mask = sum(flagged,2)*window_len > max_broken_time;
    fprintf('Removing %i channels...',nnz(removed_channel_mask));
    disp({signal.chanlocs(removed_channel_mask==1).labels});

    % remove the channels in the protect list
    if ~isempty(protect_channels)
        removed_channel_mask(set_chanid(signal,protect_channels)) = true; end
end

% annotate the data with what was removed (for visualization)
if ~isfield(signal.etc,'clean_channel_mask')
    signal.etc.clean_channel_mask = true(1,signal.nbchan); end
signal.etc.clean_channel_mask(signal.etc.clean_channel_mask) = ~removed_channel_mask;
signal.etc.badChIndex = find(removed_channel_mask==1);
signal.etc.badChLabels = {signal.chanlocs(signal.etc.badChIndex).labels};
%}

%% Whitening !!! is this being used?
% rowmeans = mean(signal.data(channel_mask, :), 2);
% data = bsxfun(@minus,signal.data(channel_mask, :), rowmeans);
%
% [E,D] = eig ( data * data' / size(data, 2) );
% signal.icasphere = (sqrtm(D)\eye(sum(channel_mask))) * E';

%%

repo: REST


text:
function [ x, y ] = p09_dat ( data_num )

%*****************************************************************************80
%
%% P09_DAT returns the data vector for problem 9.
%
%  Discussion:
%
%    This example is due to Max Waldmeier.
%
%    This data represents a measure of sunspot activity over the
%    years 1700 to 1960.  The X value is the year, and the Y value
%    is a measure of the sunspot activity, which is usually, but
%    not always, an integer.
%
%    The X data is equally spaced.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Max Waldmeier,
%    The Sunspot-Activity in the Years 1610-1960,
%    Shulthess, Zurich, 1961.
%
%  Parameters:
%
%    Input, integer DATA_NUM, the number of data points.
%
%    Output, real X(DATA_NUM,1), the abscissa data.
%
%    Output, real Y(DATA_NUM,1), the ordinate data.
%
  x(1:261,1) = ( 1700 : 1960 )';

  y(1:261,1) = [ ...
    5.0,  11.0,  16.0,  23.0,  36.0, ...
   58.0,  29.0,  20.0,  10.0,   8.0, ...
    3.0,   0.0,   0.0,   2.0,  11.0, ...
   27.0,  47.0,  63.0,  60.0,  39.0, ...
   28.0,  26.0,  22.0,  11.0,  21.0, ...
   40.0,  78.0, 122.0, 103.0,  73.0, ...
   47.0,  35.0,  11.0,   5.0,  16.0, ...
   34.0,  70.0,  81.0, 111.0, 101.0, ...
   73.0,  40.0,  20.0,  16.0,   5.0, ...
   11.0,  22.0,  40.0,  60.0,  80.9, ...
   83.4,  47.7,  47.8,  30.7,  12.2, ...
    9.6,  10.2,  32.4,  47.6,  54.0, ...
   62.9,  85.9,  61.2,  45.1,  36.4, ...
   20.9,  11.4,  37.8,  69.8, 106.1, ...
  100.8,  81.6,  66.5,  34.8,  30.6, ...
    7.0,  19.8,  92.5, 154.4, 125.9, ...
   84.8,  68.1,  38.5,  22.8,  10.2, ...
   24.1,  82.9, 132.0, 130.9, 118.1, ...
   89.9,  66.6,  60.0,  46.9,  41.0, ...
   21.3,  16.0,   6.4,   4.1,   6.8, ...
   14.5,  34.0,  45.0,  43.1,  47.5, ...
   42.2,  28.1,  10.1,   8.1,   2.5, ...
    0.0,   1.4,   5.0,  12.2,  13.9, ...
   35.4,  45.8,  41.1,  30.1,  23.9, ...
   15.6,   6.6,   4.0,   1.8,   8.5, ...
   16.6,  36.3,  49.6,  64.2,  67.0, ...
   70.9,  47.8,  27.5,   8.5,  13.2, ...
   56.9, 121.5, 138.3, 103.2,  85.7, ...
   64.6,  36.7,  24.2,  10.7,  15.0, ...
   40.1,  61.5,  98.5, 124.7,  96.3, ...
   66.6,  64.5,  54.1,  39.0,  20.6, ...
    6.7,   4.3,  22.7,  54.8,  93.8, ...
   95.8,  77.2,  59.1,  44.0,  47.0, ...
   30.5,  16.3,   7.3,  37.6,  74.0, ...
  139.0, 111.2, 101.6,  66.2,  44.7, ...
   17.0,  11.3,  12.4,   3.4,   6.0, ...
   32.3,  54.3,  59.7,  63.7,  63.5, ...
   52.2,  25.4,  13.1,   6.8,   6.3, ...
    7.1,  35.6,  73.0,  85.1,  78.0, ...
   64.0,  41.8,  26.2,  26.7,  12.1, ...
    9.5,   2.7,   5.0,  24.4,  42.0, ...
   63.5,  53.8,  62.0,  48.5,  43.9, ...
   18.6,   5.7,   3.6,   1.4,   9.6, ...
   47.4,  57.1, 103.9,  80.6,  63.6, ...
   37.6,  26.1,  14.2,   5.8,  16.7, ...
   44.3,  63.9,  69.0,  77.8,  64.9, ...
   35.7,  21.2,  11.1,   5.7,   8.7, ...
   36.1,  79.7, 114.4, 109.6,  88.8, ...
   67.8,  47.5,  30.6,  16.3,   9.6, ...
   33.2,  92.6, 151.6, 136.3, 134.7, ...
   83.9,  69.4,  31.5,  13.9,   4.4, ...
   38.0, 141.7, 190.2, 184.8, 159.0, ...
  112.3 ]';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function dispProgress(name, value, maxvalue)
%DISPPROGRESS display calculation progress
% initialize and update progress bar
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

persistent prop

if(nargin < 3)
    maxvalue = 0;
end

if(~exist('prop','var') || isempty(prop) || strcmp(name,'InitDispProgress'))
    if(isstruct(value))
        prop = value;
        return;
    else
%         error('dispProgress(): Unknown data structure for property');
        return; % shortcut for GA simulation
    end
end

if(prop.flagParallel)
    if(verLessThan('matlab','8.4'))
        if(matlabpool('size') > 0)
            return; % shortcut for parallel computing
        end
    else
        if(~isempty(gcp('nocreate')))
            return % shortcut for parallel computing
        end
    end
end

% if(matlabpool('size') > 0)
%     openPool = true;
%     currpath = fileparts(fileparts(mfilename('fullpath')));
%     a = load([currpath, filesep, 'prop.mat']);
% else
%     openPool = false;
% end

if(prop.flagDisp)
    % GUI

    if(~prop.openPool)
        if(prop.disp.openBar(strcmp(prop.disp.names,name)))
            % progress bar already exists
            if(ischar(value) && strcmp(value,'Close'))
                multiWaitbar(name,value);
                prop.disp.openBar(strcmp(prop.disp.names,name)) = false;
            else
                cval = round(value*100);
                if(cval == 0), cval = max(ceil(value*100),1); end;
                if(cval > 100), cval = 100; end;
                abort = multiWaitbar(name, value, 'Color', prop.disp.colors(cval,:));
                if(abort)
                    multiWaitbar('CLOSEALL');
                    error('User Termination of CS reconstruction!');
                end
            end
        else
            % progress bar needs initialization
            if(maxvalue > 1)
                multiWaitbar(name,value,'CancelFcn', @(a,b) disp( ['Cancel ',a] ));
                prop.disp.openBar(strcmp(prop.disp.names,name)) = true;
            end
        end

    else
        % parallel computation
        if(prop.disp.openBar(strcmp(prop.disp.names,name)))
            % progress bar already exists
            if(ischar(value) && strcmp(value,'Close'))
                try
                    delete(prop.disp.ppm{strcmp(prop.disp.names,name)});
                catch err
                end
                prop.disp.openBar(strcmp(prop.disp.names,name)) = false;
            else
                prop.disp.ppm{strcmp(prop.disp.names,name)}.increment(value);
            end
        else
            % progress bar needs initialization
            if(maxvalue > 1)
                helper = ParforProgressStarter2(name, maxvalue, 1);
                prop.disp.ppm{strcmp(prop.disp.names,name)} = helper;
                prop.disp.openBar(strcmp(prop.disp.names,name)) = true;
            end
        end
%         save([currpath, filesep, 'prop.mat'], '-struct', 'a');
    end

else
    % console

    if(ischar(value) && strcmp(value,'Close'))
        if(any(strcmp(name,{'CG','SENSE - Time','SENSE - Kernel','SENSE - Concatenate','Extracting K-Space','GRAPPA calibration','Partial Fourier','ESPReSSo','Trafo','Log barrier','Newton','Proximal Average', 'ADMM'})))
            if(prop.openPool)
                fprintf([repmat('\b',1,length(name)+1)]);
            else
                fprintf([repmat('\b',1,length(prop.lastOut))]);
            end
        end

    elseif(value == 0)
        % init step
        prop.maxvals(strcmp(prop.disp.names,name)) = maxvalue;
        if(strcmp(name,'FOCUSS'))
%             prop.lastOut = sprintf('%s %3d%%',name,value);
            prop.lastVal = 0;
            fprintf('%s [%s]\n%s [', name, repmat('- ',1,25),repmat(' ',1,length(name)));
        elseif(strcmp(name,'CG'))
%             prop.lastOut = sprintf('%s\t\t%s %3d%%\n',prop.lastOut,name,value);
%             if(~prop.openPool)
%                 fprintf('%s',prop.lastOut);
%             else
%                 prop.lastOut = '';
%             end
            prop.lastVal = 0;
            fprintf('%s [%s]\n%s [', name, repmat('- ',1,25),repmat(' ',1,length(name)));
        elseif(any(strcmp(name,{'Extracting K-Space','GRAPPA calibration','Partial Fourier','ESPReSSo','Trafo','Log barrier','Newton','Proximal Average', 'ADMM'})))
%             if(prop.openPool)
%                 prop.lastOut = sprintf('%s\n',name);
%             else
%                 prop.lastOut = sprintf('%s %3d%%\n',name,value);
%             end
%             fprintf('%s',prop.lastOut);
            prop.lastVal = 0;
            fprintf('%s [%s]\n%s [', name, repmat('- ',1,25),repmat(' ',1,length(name)));
        end

    else
        if(~prop.openPool)
            switch name
                case {'Slices','Channels','Repetitions','POCS','Time','Frequency','Averages'}
                    outstring = sprintf('%s %d/%d\n', name, value*prop.maxvals(strcmp(prop.disp.names,name)), prop.maxvals(strcmp(prop.disp.names,name)));
                    delString = 0;

                case {'Extracting K-Space','GRAPPA calibration','Partial Fourier','ESPReSSo','Trafo','Log barrier','Newton','Line Search','Proximal Average', 'ADMM'}
%                     outstring = sprintf('%s %3d%%\n', name, value*100);
%                     delString = length(prop.lastOut);

                    value = round(value*100);
                    if(mod(value,2) ~= 0)
                        value = value-1;
                    end
                    value = value/2;

                    if(value >= 50)
                        outstring = sprintf('*]\n');
                    else
                        if(value == prop.lastVal)
                            outstring = '';
                        else
                            outstring = '*';
                        end
                    end
                    prop.lastVal = value;
                    delString = 0;
                case {'FOCUSS', 'CG'}

%                     if(strcmp(name,'FOCUSS'))
%                         outstring = sprintf('%s %3.0f%%%s',name,value*100,prop.lastOut(12:end));
%                     else
%                         outstring = sprintf('%s\t\t%s %3.0f%%\n', prop.lastOut(1:11),name,value*100);
%                     end
%
%                     if(~isempty(prop.lastOut))
%                         delString = length(prop.lastOut);
%                     else
%                         delString = 0;
%                     end
                    value = round(value*100);
                    if(mod(value,2) ~= 0)
                        value = value-1;
                    end
                    value = value/2;

                    if(value >= 50)
                        outstring = sprintf('*]\n');
                    else
                        if(value == prop.lastVal)
                            outstring = '';
                        else
                            outstring = '*';
                        end
                    end
                    prop.lastVal = value;
                    delString = 0;

                case 'SENSE - Slice'
                    outstring = sprintf('%s %d/%d\n', name, value*prop.maxvals(strcmp(prop.disp.names,name)), prop.maxvals(strcmp(prop.disp.names,name)));
                    delString = 0;

                case 'SENSE - Time'
                    outstring = sprintf('%s %d/%d\n', name, value*prop.maxvals(strcmp(prop.disp.names,name)), prop.maxvals(strcmp(prop.disp.names,name)));

                    if(~isempty(prop.lastOut))
                        delString = length(prop.lastOut);
                    else
                        delString = 0;
                    end

                case {'SENSE - Kernel', 'SENSE - Concatenate'}
                    outstring = sprintf('%s %3.0f%%\n',name,value);
                    if(~isempty(prop.lastOut))
                        delString = length(prop.lastOut);
                    else
                        delString = 0;
                    end

                otherwise
                    error('dispProgress(): Undetermined display name');
            end

            prop.lastOut = outstring;

            fprintf([repmat('\b',1,delString),'%s'],outstring);
        end
    end
end

end

author: thomaskuestner


text:
function p_ = ViewRealizedVol(X,p)

[J,K]=size(X);

% constrain probabilities to sum to one...
Aeq = ones(1,J);
beq=1;

% ...constrain the median...
V = abs(X(:,1));

[V_Sort I_Sort]=sort(V);
F=cumsum(p(I_Sort));

I_Reference=max(find(F<=3/5));
V_Reference=V_Sort(I_Reference);

I_Select=find(V<=V_Reference);

a=zeros(1,J);
a(I_Select)=1;

A = a;
b = .5;

% ...compute posterior probabilities
p_ = EntropyProg(p,A,b,Aeq ,beq);

author: Sable

repo: mcbench-benchmarks


text:
function opti_Install_Help()
%OPTI_INSTALL_HELP  Build and save the OPTI help database

%You Must be in the OPTI Directory!
builddocsearchdb ([cd '\Help\OPTI_HTML'])

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function ui = mrViewAttachTSeries(ui, tSeries);
% Attach tSeries paths / files to a mrViewer UI, checking
% that the files exist, but not loading them.
%
% ui = mrViewAttachTSeries(ui, tSeries);
%
% ras, 12/05.
if notDefined('ui'), ui = mrViewGet; end
if ishandle(ui), ui = get(ui, 'UserData'); end

% allow loaded tSeries files to be passed -- it may be memory-hogging,
% but it may make it easier in some situations to have the data pre-loaded
if isstruct(tSeries)
    ui.tSeries{end+1} = tSeries;
else
    % either a path to one file, or several paths, were passed: check
    % that the files exist
    if ischar(tSeries)
        % path string
%         if ~exist(tSeries, 'file')
%             error(sprintf('tSeries file %s not found.', tSeries));
%         end
        ui.tSeries{end+1} = tSeries;
    elseif iscell(tSeries)
        for i = 1:length(tSeries)
            ui = mrViewAttachTSeries(ui, tSeries{i});
        end
    end
end

% check if we should make the button for plotting time courses
% in the ROI panel visible:
if isfield(ui.controls,'roiTimeCourse') & ishandle(ui.controls.roiTimeCourse)
    if ~isempty(ui.tSeries) & ~isempty(ui.stim)
        set(ui.controls.roiTimeCourse, 'Visible', 'on');
    else
        set(ui.controls.roiTimeCourse, 'Visible', 'off');
    end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [K] = R2K(R)
% Convert temperature from Rankine to Kelvin.
K = R*5/9;


text:
function events = in_events_xltek(sFile, EventFile)
% IN_EVENTS_XLTEK: Open an XLTEK exported events file (.txt)

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

% Open and read file
fid = fopen(EventFile,'r');
% Skip 6 header lines
tmp = fgetl(fid);
tmp = fgetl(fid);
tmp = fgetl(fid);
tmp = fgetl(fid);
tmp = fgetl(fid);
tmp = fgetl(fid);

% Read start time
strLine = fgetl(fid);
[tmp,tmp2] = strtok(strLine);
starttime = strtok(tmp2);
starttime = sum(sscanf(starttime, '%d:%d:%d') .* [3600;60;1]);

% Loop to read all the events
evtLabel = {};
evtTime  = [];
while 1
    % Get line
    strLine = fgetl(fid);
    if (strLine == -1)
        break;
    end
    % Get label and timing
    [tmp1,tmp2] = strtok(strLine);
    [tmp1,tmp2] = strtok(tmp2);
    evtTime(end+1) = sum(sscanf(tmp1, '%d:%d:%d') .* [3600;60;1]) - starttime;
    evtLabel{end+1} = strtrim(tmp2);
end
% Close file
fclose(fid);

% ===== CONVERT TO BRAINSTORM FORMAT =====
% List of events (keep original order)
[uniqueEvt, I] = unique(evtLabel);
uniqueEvt = evtLabel(sort(I));
% Initialize returned structure
events = repmat(db_template('event'), [1, length(uniqueEvt)]);
% Create events list
for iEvt = 1:length(uniqueEvt)
    % Find all the occurrences of event #iEvt
    iMrk = find(strcmpi(evtLabel, uniqueEvt{iEvt}));
    % Add event structure
    events(iEvt).label      = uniqueEvt{iEvt};
    events(iEvt).times      = unique(round(evtTime(iMrk) .* sFile.prop.sfreq)) ./ sFile.prop.sfreq;
    events(iEvt).epochs     = ones(1, length(events(iEvt).times));
    events(iEvt).reactTimes = [];
    events(iEvt).select     = 1;
    events(iEvt).channels   = [];
    events(iEvt).notes      = [];
end

author: brainstorm-tools


text:
%COMSetSubCodes COMSetSubCodes enumeration for the brick
%
% Notes::
% - COMSetSubCodes can be found in the EV3 documentation and source code
% (bytecodes.h)

classdef COMSetSubCodes < uint8
    enumeration
         SetOnOff (1)
         SetVisible (2)
         SetSearch (3)
         SetPin (5)
         SetPasskey (6)
         SetConnection (7)
         SetBrickName (8)
         SetMoveUp (9)
         SetMoveDown (10)
         SetEncrypt (11)
         SetSSID (12)
         SetMode2 (13)
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
% Copyright (c) Roman Garnett, 2012--2014.

function pdfs = evaluate_pdfs(mus, Ks, x, weights)

  persistent individual_pdfs;

  % reset if no arguments given
  if (nargin == 0)
    individual_pdfs = [];
    return;
  end

  % precompute N(x; mu, K) for every x
  if (isempty(individual_pdfs))
    num_points = size(x, 1);
    num_nodes  = size(mus, 1);

    individual_pdfs = zeros(num_nodes, num_points);
    for i = 1:num_nodes
        individual_pdfs(i, :) = mvnpdf(x, mus(i, :), Ks(:, :, i))';
    end
  end

  pdfs = weights * individual_pdfs;

end

author: muhanzhang

repo: DGCNN

save_path: github-repos/MATLAB/muhanzhang-DGCNN


text:
function Fo = blockframeaccel(F, Lb, varargin)
%BLOCKFRAMEACCEL Precompute structures for block processing
%   Usage: F = blockframeaccel(F,Lb);
%
%   `F=blockframeaccel(F,Lb)` has to be called for each frame object prior to
%   entering the main loop where |blockana| and |blocksyn| are called.
%   The function works entirely like |frameaccel| but in addition, it prepares
%   structures for the processing of a consecutive stream of blocks.
%
%      `'sliwin',sliwin`   : Slicing window. `sliwin` have to be a window
%                            of length *2Lb* or a string accepted
%                            by the |firwin| function. It is used only in
%                            the slicing window approach. The default is
%                            `'hann'`.
%
%      `'zpad',zpad`   : Number of zero samples the block will be padded
%                        after it is windowed by a slicing window. This
%                        does not affect the synthesis windowing.

complainif_notenoughargs(nargin,2,'BLOCKFRAMEACCEL');
complainif_notvalidframeobj(F,'BLOCKFRAMEACCEL');

definput.flags.blockalg = {'naive','sliced','segola'};
definput.keyvals.sliwin = [];
definput.keyvals.zpad = 0;
[flags,kv]=ltfatarghelper({},definput,varargin);

isSliProp = ~isempty(kv.sliwin) || kv.zpad~=0;
assert(~(~flags.do_sliced && isSliProp),...
   sprintf(['%s: Definig slicing window properties without setting the',...
            ' ''sliced'' flag.'], mfilename));

if flags.do_sliced
   if isempty(kv.sliwin)
      kv.sliwin = 'hann';
   end

   if ~isnumeric(kv.sliwin)
      kv.sliwin = fftshift(sqrt(firwin(kv.sliwin,2*Lb)));
   end

   Fo = frameaccel(F,2*Lb+2*kv.zpad);
   Fo.sliwin = kv.sliwin;
   Fo.zpad = kv.zpad;
elseif flags.do_segola
   % Determine window length without calling frameaccel
   % Fo = frameaccel(F,Lb);
   winLen = framefirlen(F);

   if winLen==-1
      error(['%s: Segment overlap cannot be used with this frame.,'...
             ' It does not have FIR windows.'],upper(mfilename));
   end

   switch(F.type)
      case {'fwt'}
         Fo = frameaccel(F,Lb);
         Fo.a = F.g.a(:);
      case {'dgt','dgtreal'}
         Fo = frameaccel(F,Lb+winLen-1+F.a);
%       case {'filterbank','filterbankreal','ufilterbank','ufilterbankreal'}
%          lcma =  filterbanklength(1,F.a(:,1));
%          Fo = frameaccel(F,Lb+winLen-1+lcma);
%          assert(all(Fo.a(:,2)==1), '%s: Fractional subsampling is not supported',upper(mfilename) );
%          Fo.lcma =  lcma;
      case {'dwilt'}
         Fo = frameaccel(F,Lb+winLen-1+2*F.M);
         Fo.a = 2*Fo.M;
      case {'wmdct'}
         Fo = frameaccel(F,Lb+winLen-1+F.M);
         Fo.a = Fo.M;
      otherwise
	 error('%s: Unsupported frame for segola.',upper(mfilename));
   end

   % This is important otherwise we would get 0 coefficients for some
   % blocks.
   assert(max(Fo.a(:,1)) <= Lb ,sprintf(['%s: Time step %i is bigger than the',...
      ' block length %i.'],upper(mfilename),max(Fo.a(:,1)),Lb));

   Fo.winLen = winLen;

elseif flags.do_naive
   Fo = frameaccel(F,Lb);
end

Fo.blockalg = flags.blockalg;

function winLen = framefirlen(F)
%FRAMEFIRLEN Frame window/filter length
%
%   Function returns length of the longest FIR window/filter. The function
%   returns -1 if the frame does not have FIR windows.

winLen = -1;
info = [];
switch(F.type)
      case {'dgt','dgtreal'}
        [~, info] =  gabwin(F.g,F.a,F.M,[],F.kv.lt);
      case {'dwilt','wmdct'}
        [~, info] = wilwin(F.g,F.M,[],upper(mfilename));
      case {'filterbank','ufilterbank'}
        [~, ~,info]  = filterbankwin(F.g,F.a);
      case {'filterbankreal','ufilterbankreal'}
        [~, ~,info]  = filterbankwin(F.g,F.a,'real');
      case 'fwt'
        winLen = (F.g.a(1)^F.J-1)/(F.g.a(1)-1)*(numel(F.g.g{1}.h)-1)+1;
end;

if ~isempty(info) && isfield(info,'isfir') && info.isfir
   if isfield(info,'longestfilter')
      winLen = info.longestfilter;
   else
      winLen = max(info.gl);
   end
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% Replace XXXX by XXXXX for better documentation in MATLAB
% allFiles = dir('*.m')
% for iFile = 1:length(allFiles)
%    replace_in_all_files(allFiles(iFile).name);
%    disp(allFiles(iFile).name);
% end

function replace_in_all_files(fileName)

fid = fopen(fileName, 'r');
if fid == -1
    error('File not found');
end

allStrs = {};
count   = 1;
expression1 = '[_A-Za-z0-9]+\(\)'; % + means at least 1
expression2 = '[A-Za-z0-9]\(\)';
while ~feof(fid)
    strTmp = fgetl(fid);

    if length(strTmp) > 0 && strTmp(1) == '%'
        pos1 = regexp(strTmp, expression1);
        pos2 = regexp(strTmp, expression2)+1;

        if length(pos1) ~= length(pos2)
            error('Issue with length of regular expressions')
        end
        for iPos = length(pos1):-1:1
            strTmp(pos2(iPos):pos2(iPos)+1) = [];
            strTmp(pos1(iPos):pos2(iPos)-1) = upper(strTmp(pos1(iPos):pos2(iPos)-1));
        end
    end

    allStrs{count} = strTmp;
    count = count+1;
end
fclose(fid);

% write new file
fid = fopen(fileName, 'w');
for iRow = 1:length(allStrs)
    fprintf(fid, '%s\n', allStrs{iRow});
end
fclose(fid);

function [twoLines, content] = replace_txt(twoLines, pattern1, pattern2)

indB1 = strfind(twoLines, pattern1{1}); indPat1 = 1;
indB2 = strfind(twoLines, pattern2{1}); indPat2 = 1;
content = '';
if isempty(indB1), indB1 = strfind(twoLines, pattern1{2}); indPat1 = 2; end
if isempty(indB2), indB2 = strfind(twoLines, pattern2{2}); indPat2 = 2; end
if ~isempty(indB1) && ~isempty(indB2)
    content = twoLines(indB1+length(pattern1{indPat1}):indB2-1);
    fprintf('Found: %s\n', content);
    twoLines = [ twoLines(1:indB1-1) '[' content '](http://sccn.ucsd.edu/eeglab/locatefile.php?file=' content ')' twoLines(indB2+length(pattern2{indPat2}):end) ];
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [xx, yy] = chebpts2(nx, ny, D, kind)
%CHEBPTS2   Chebyshev tensor product grid.
%   [XX YY] = CHEBPTS2(N) constructs an N by N grid of Chebyshev tensor points
%   on [-1,1]^2.
%
%   [XX YY] = CHEBPTS2(NX,NY) constructs an NX by NY grid of Chebyshev tensor
%   points on [-1,1]^2.
%
%   [XX YY] = CHEBPTS2(NX,NY,D) constructs an NX by NY grid of Chebyshev tensor
%   points on the rectangle [a,b] x [c,d], where D = [a b c d].
%
%   [XX YY] = CHEBPTS2(NX,NY,D,KIND) constructs a Chebyshev tensor grid of
%   the kind KIND. KIND = 2 is the default.
%
% See also CHEBPTS.

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

if ( nargin > 2 )
   % Third argument should be a domain.
   D = D(:).';  % make a row vector.
   if ( ~all( size( D ) == [1 4] ) )
        error('CHEBFUN:CHEBFUN2:chebpts2:domain', 'Unrecognised domain.');
   end
else  % Default to the canoncial domain.
    D = [-1, 1, -1, 1];
end

if ( nargin == 1 )
    % Make it a square Chebyshev grid if only one input.
    ny = nx;
end
if ( nargin < 4 )
    % default to Chebyshev point of the 2nd kind.
    kind = 2;
end

% Get points:
x = chebpts(nx, D(1:2), kind);
y = chebpts(ny, D(3:4), kind);

% Tensor product.
[xx, yy] = meshgrid(x, y);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef ComplianceAndAdjointPrinter < CompositeResultsPrinter

    methods (Access = public)

        function obj = ComplianceAndAdjointPrinter(d)
            obj.init(d);
        end

    end

    methods (Access = protected)

        function createPrinters(obj,d)
            obj.printers{1} = obj.createCompliancePrinter(d);
            obj.printers{2} = obj.createAdjointPrinter(d);
          %  obj.printers{3} = obj.createRegularizedDensityPrinter(d);
        end

        function storeFieldsToPrint(obj,d)
            obj.storeComplianceFields(d);
            obj.storeAdjointFields(d);
       %     obj.storeRegularizedDensity(d);
        end

        function createHeadPrinter(obj,d,dh)
            phyPr = d.cost.shapeFunctions{1}.getPhysicalProblems();
            d.quad = phyPr{1}.element.quadrature;
            obj.printers{1}.createHeadPrinter(d,dh);
            h = obj.printers{1}.getHeadPrinter();
            obj.headPrinter = h;
        end

    end

    methods (Access = private, Static)

        function p = createCompliancePrinter(d)
            p =  ResultsPrinter.create('Elasticity',d);
            p.setStrVariablesNames('Stress','Strain','Disp');
        end

        function p = createAdjointPrinter(d)
            p =  ResultsPrinter.create('Elasticity',d);
            p.setStrVariablesNames('StressAdj','StrainAdj','DispAdj');
        end

     %   function p = createRegularizedDensityPrinter(d)
     %       p = ResultsPrinter.create('DensityGauss',d);
     %   end

        function storeFieldsToPrintFromPhyPr(printer,phyPr)
            d.fields = phyPr.variables;
            printer.storeFieldsToPrint(d);
        end

    end

    methods (Access = private)

        function storeComplianceFields(obj,d)
            d.fields = d.phyProblems{1}.variables;
            obj.printers{1}.storeFieldsToPrint(d);
        end

        function storeAdjointFields(obj,d)
            d.fields = d.phyProblems{2}.variables;
            obj.printers{2}.storeFieldsToPrint(d);
        end

        function storeRegularizedDensity(obj,d)
     %       d.fields = d.regDensity;
     %       obj.printers{3}.storeFieldsToPrint(d);
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function fmtt
%FMTT 	Unit test for the function FMT.

%	O. Lemoine - May 1996.

N=128;

% Perfect reconstruction with ifmt
fmin=0.1; fmax=0.5;
sig=amgauss(N).*fmconst(N,.3);
[MELLIN,BETA]=fmt(sig,fmin,fmax,N);
X=ifmt(MELLIN,BETA,N);
err=abs(X-sig);
if any(err>1e-7),
 error('fmt test 1 failed');
end

% Energy conservation
x=fmconst(N);
fmin=0.1; fmax=0.4;
FMT=fmt(x,fmin,fmax,N);
SP=fft(x);
indmin = 1+round(fmin*(N-2));
indmax = 1+round(fmax*(N-2));
SPana=SP(indmin:indmax);
nu=(indmin:indmax)'/N;
SPp=SPana./nu;
Es=SPp'*SPana;
Efmt=norm(FMT)^2;
if abs(Es-Efmt)>sqrt(eps),
 error('fmt test 2 failed');
end;

% Unitarity of the MT
x1=amgauss(N).*fmlin(N,.15,.35);
x2=amgauss(N).*fmconst(N);
fmin=0.01; fmax=0.49;
FMT1=fmt(x1,fmin,fmax,2*N);
FMT2=fmt(x2,fmin,fmax,2*N);
indmin = 1+round(fmin*(2*N-2));
indmax = 1+round(fmax*(2*N-2));
SP1=fft(x1); SP2=fft(x2);
nu=(indmin:indmax)'/N/2;
SP1p=SP1(indmin:indmax)./nu;
cor1=SP1p'*SP2(indmin:indmax);
cor2=conj(FMT1*FMT2');
if abs(cor1-cor2)>N*1e-2,
 error('fmt test 3 failed');
end;

% Covariance by dilation
% Property of the Mellin transform used in scale.
% So as scale works, this property is verified.

% MT of a product = convolution of the MT
x1=amgauss(N).*fmlin(N,.15,.35);
x2=amgauss(N).*fmsin(N,.15,.35);
FMT1=fmt(x1,fmin,fmax,2*N);
FMT2=fmt(x2,fmin,fmax,2*N);
FMT=conv(FMT1,FMT2);
FMT=FMT/max(real(FMT));
X1=fft(x1); X2=fft(x2);
X=X1.*X2;
x=fftshift(ifft(X));
FMTp=fmt(x,fmin,fmax,2*N);
FMTp=FMTp/max(real(FMTp));
Diff=FMTp-FMT(N+1:3*N);
if any(abs(Diff)>1e-4),
 error('fmt test 4 failed');
end

N=121;

% Perfect reconstruction with ifmt
fmin=0.1; fmax=0.5;
sig=amgauss(N).*fmconst(N,.3);
[MELLIN,BETA]=fmt(sig,fmin,fmax,N+1);
X=ifmt(MELLIN,BETA,N);
err=abs(X-sig);
if any(err>1e-2),
 error('fmt test 5 failed');
end

% Energy conservation
x=fmconst(N);
fmin=0.1; fmax=0.4;
FMT=fmt(x,fmin,fmax,N+1);
SP=fft(hilbert(real(x)));
indmin = 1+round(fmin*(N-2));
indmax = 1+round(fmax*(N-2));
SPana=SP(indmin:indmax);
nu=(indmin:indmax)'/(N+1);
SPp=SPana./nu;
Es=SPp'*SPana;
Efmt=norm(FMT)^2;
if abs(Es-Efmt)>sqrt(eps),
 error('fmt test 6 failed');
end;

% Unitarity of the MT
x1=amgauss(N).*fmlin(N,.15,.35);
x2=amgauss(N).*fmconst(N);
fmin=0.01; fmax=0.49;
FMT1=fmt(x1,fmin,fmax,2*N);
FMT2=fmt(x2,fmin,fmax,2*N);
indmin = 1+round(fmin*(2*N-2));
indmax = 1+round(fmax*(2*N-2));
SP1=fft(x1); SP2=fft(x2);
nu=(indmin:indmax)'/N/2;
SP1p=SP1(indmin:indmax)./nu;
cor1=SP1p'*SP2(indmin:indmax);
cor2=conj(FMT1*FMT2');
if abs(cor1-cor2)>N*1e-2,
 error('fmt test 7 failed');
end;

% MT of a product = convolution of the MT
x1=amgauss(N).*fmlin(N,.15,.35);
x2=amgauss(N).*fmsin(N,.15,.35);
FMT1=fmt(x1,fmin,fmax,2*N);
FMT2=fmt(x2,fmin,fmax,2*N);
FMT=conv(FMT1,FMT2);
FMT=FMT/max(real(FMT));
X1=fft(x1); X2=fft(x2);
X=X1.*X2;
x=fftshift(ifft(X));
FMTp=fmt(x,fmin,fmax,2*N);
FMTp=FMTp/max(real(FMTp));
Diff=FMTp-FMT(N+1:3*N);
if any(abs(Diff)>1e-4),
 error('fmt test 8 failed');
end

author: HeLiangHIT


text:
 function tk = ir_chebyshev_poly(n)
%function tk = ir_chebyshev_poly(n)
%|
%| based on ChebyshevPoly.m by David Terr, Raytheon, 5-10-04
%|
%| Given nonnegative integer n, compute the
%| Chebyshev polynomial T_n. Return the result as a column vector whose mth
%| element is the coefficient of x^(n+1-m).
%| polyval(ir_chebyshev_poly(n), x) evaluates T_n(x).

if nargin < 1, ir_usage, end
if streq(n, 'test'), ir_chebyshev_poly_test, return, end

if n==0
	tk = 1;
elseif n==1
	tk = [1 0]';
else
	tkm2 = zeros(n+1,1);
	tkm2(n+1) = 1;
	tkm1 = zeros(n+1,1);
	tkm1(n) = 1;

	for k=2:n
		tk = zeros(n+1,1);

		for e=n-k+1:2:n
			tk(e) = 2*tkm1(e+1) - tkm2(e);
		end

		if mod(k,2)==0
			tk(n+1) = (-1)^(k/2);
		end

		if k<n
			tkm2 = tkm1;
			tkm1 = tk;
		end
	end
end

% ir_chebyshev_poly_test()
function ir_chebyshev_poly_test
nlist = [0:4];
x = 1 * linspace(-1,1,101)';
y = zeros(numel(x), numel(nlist));
for in=1:numel(nlist);
	n = nlist(in);
	p = ir_chebyshev_poly(n);
	y(:,in) = polyval(p,x);
end
if im
	clf
	plot(x, y)
	axis([-1 1 -1.1 1.1])
	xtick([-1 0 1]), ytick([-1 0 1])
	titlef('Chebychev Polynomials')
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [O] = forwardcomp(X,W,L)
%FORWARD COMPUTATION.
%
%	Syntax
%     [O] = forwardcomp(X,W,L)
%
%	Description
%
%	  FORWARDCOMP computes the outputs for a multiple perceptron Neural Network
%	  from its net input.
%
%	  FORWARDCOMP(X,W,L) takes:
%	    X - matrix of input patterns including one row of bias = 1.
%	    W - matrix of weights of the Network (cell vectors).
%	    L - vector of hidden layers and neurons per layer [m1 m2...mL].
%	  and returns
%       O - The output of the network
%
%
% Paul Acquatella, 08-02-2009
% Revised 08-02-2009
% Copyright Paul Acquatella
% $Revision: 1.1.6.2 $  $Date: 2009/08/02 18:21:07 $

%Calcular Numero de Salidas y de Capas Ocultas
Q = size(W{end},1);         %Numero de Salidas
H = size(L,2);              %Numero de Capas Ocultas "Depth"

%Inicializar los campos locales inducidos 'v'
v = cell(1,H+1);
v{end} = zeros(Q,1);
for i=1:H;
    v{i} = zeros(L(:,i),1);
end

%Inicializar las salidas 'y' de cada capa
y = cell(1,H+1);
y{end} = zeros(Q,1);
for i=1:H;
    y{i} = [zeros(L(:,i),1); 1];
end

for j = 1:(H+1)
    if j == 1
        v{j} = W{j}*X;                  %Compute the induced local field
        y{j}(1:end-1) = logsig(v{j});   %Compute the layer output
    elseif j == (H+1)
        v{j} = W{j}*y{j-1};             %Compute the induced local field
        y{j} = logsig(v{j});            %Compute the layer output
        O = y{j};
    else
        v{j} = W{j}*y{j-1};             %Compute the induced local field
        y{j}(1:end-1) = logsig(v{j});   %Compute the layer output
    end
end


text:
% --------------------------------------------------------
% MDP Tracking
% Copyright (c) 2015 CVGL Stanford
% Licensed under The MIT License [see LICENSE for details]
% Written by Yu Xiang
% --------------------------------------------------------
function MOT_make_videos

is_save = 1;

opt = globals();
N = numel(opt.mot2d_test_seqs);

for seq_idx = 1:N
    close all;
    hf = figure(1);
    seq_name = opt.mot2d_test_seqs{seq_idx};
    seq_num = opt.mot2d_test_nums(seq_idx);

    % build the dres structure for images
    filename = sprintf('%s/%s_dres_image.mat', opt.results, seq_name);
    if exist(filename, 'file') ~= 0
        object = load(filename);
        dres_image = object.dres_image;
        fprintf('load images from file %s done\n', filename);
    else
        dres_image = read_dres_image(opt, seq_set, seq_name, seq_num);
        fprintf('read images done\n');
        save(filename, 'dres_image', '-v7.3');
    end

    % read tracking results
    filename = sprintf('results_MOT/results_MOT_1/%s.txt', seq_name);
    dres_track = read_mot2dres(filename);
    fprintf('read tracking results from %s\n', filename);
    ids = unique(dres_track.id);
    cmap = colormap(hsv(numel(ids)));
    cmap = cmap(randperm(numel(ids)),:);

    if is_save
        file_video = sprintf('results_MOT/results_MOT_1/%s.avi', seq_name);
        aviobj = VideoWriter(file_video);
        aviobj.FrameRate = 9;
        open(aviobj);
        fprintf('save video to %s\n', file_video);
    end

    for fr = 1:seq_num
        show_dres(fr, dres_image.I{fr}, '', dres_track, 2, cmap);
        if is_save
            writeVideo(aviobj, getframe(hf));
        else
            pause;
        end
    end

    if is_save
        close(aviobj);
    end
end

author: yuxng


text:
function [sLDA WLDA M WPCA]=mylda(data,class,n)
% [sLDA WLDA M WPCA]=mylda(data,class,n)
% this function written by muhammet balcilar
% yildiz technical university computer engineering department
% istanbul turkiye 2011

% this function convert data from its original space to LDA space
% if number of data samples is less than number of diamension, PCA is
% implemented for reducing number of diamension to #samples-1.
% after PCA, LDA is implemented for reducing diamention to n.

% data is consist of M rows(sample size), N cols(dimensions)
% class is consist of M rows(sample size), 1 cols , each element of class
% is shows class number of each data sample
% (class number must be integer 1 to classsize)
% n is the number of outputs data diamensions.(optionally)
% sLDA is consist of M rows(sample size) n cols(new dimensions)
% WPCA is translate matrix which convert to original space to PCA space
% M is the mean vector of training set
% WLDA is the translate matrix which convert to original space to LDA space
% exaple: there are 4 samples which have 5 diamensions.first two samples
% are member of class 1 others are member of class 2.
% Train= [5.6,5.7,5.5,5.7 5.6;
%     5.7,5.3,5.1,5.0 5.2;
%     10.6,9.9,10.4,10.7 10.2;
%     10.7,9.8,9.9,10 10];
% Class=[1;1;2;2];
% [sLDA WLDA M WPCA]=mylda(Train,Class)
% Test= [4.9 5.5 4.8 5.7 5];
% LDATEST = (Test-M)*WPCA*WLDA

usinif=unique(class);
if nargin==2
    n=length(usinif)-1;
end

if size(data,2)>=size(data,1)
    % PCA start
    O=data';
    m=(mean(O'))';
    for i=1:size(O,2)
        mO(:,i)=O(:,i)-m;
    end
    CV=mO*mO';
    [v u]=eig(CV);
    D=v(:,end-size(data,1)+2:end);
    yO=(mO'*D)';
    M=m';
    WPCA=D;
    % PCA finished
else
    yO=data';
    M=zeros(1,size(data,2));
    WPCA=1;
end

% LDA start
mU=(mean(yO'))';
mK=[];
for i=1:length(usinif)
    I=find(class==i);
    ort=(mean(yO(:,I)'))';
    mK=[mK ort];
    for j=1:length(I)
        UU(:,I(j))=yO(:,I(j))-ort;
    end
end
for i=1:length(usinif)
    I=find(class==i);
    S{i}= UU(:,I)*UU(:,I)';
end
SW=S{1};
for i=2:length(usinif)
    SW=SW+S{i};
end

for i=1:length(usinif)
    mmK(:,i)=mK(:,i)-mU;
end
SB=2*mmK*mmK';
[w u]=eig(SB,SW);
u=abs(diag(u));
u=[u [1:length(u)]'];
u=sortrows(u,1);
WLDA=w(:,u(end-n+1:end,2));
sLDA=(yO'*WLDA)';


text:
function plot(obj)
% Plot an fmri_model object
%
% :Usage:
% ::
%
%     plot(obj)

% ..
%    Setup
% ..

TR = obj.xY.RT;

bf = obj.xBF.bf;

% Define sessions and number of conditions
nsess = length(obj.Sess);

% ----------------------------------------------
% Image
% ----------------------------------------------
create_figure('design matrix image');
imagesc(obj.xX.X);
title('Design matrix (X)');
set(gca, 'YDir', 'reverse');
axis tight;
xlabel('Regressors');
ylabel('Images (TRs)');
colormap gray
drawnow

% ----------------------------------------------
% Basis sets
% ----------------------------------------------
nconds = length(obj.xBF);
create_figure('basis sets', 1, nconds);
for i = 1:nconds

    subplot(1, nconds, i);
    xx = linspace(0, obj.xBF(i).length - 1, size(obj.xBF(i).bf, 1));

    plot(xx, obj.xBF(i).bf, 'LineWidth', 2);
    title(sprintf('Basis set, Cond. %3.0f, %s', i, obj.xX.name{i}));
    axis tight
    xlabel('Time (sec)');

end

% ----------------------------------------------
% Plot
% ----------------------------------------------

switch obj.build_method
% ----------------------------------------------
    case 'Separate sessions'
% ----------------------------------------------
        nr = 2; nc = ceil(nsess) ./ 2;

        create_figure('design matrix', nr, nc);

        for i = 1:nsess
            subplot(nr, nc, i)
            set(gca, 'FontSize', 10);

            [Xs, dummy, C, dummy, names] = get_session_X(obj, i);

            plot_subfcn(obj, Xs, C, names);
            title(sprintf('Session %3.0f', i));

        end

% ----------------------------------------------
    case 'Concatenated sessions'
% ----------------------------------------------
        create_figure('design matrix');

        plot_subfcn(obj, obj.xX.X(:, obj.xX.iH), obj.xX.X(:, obj.xX.iC), obj.xX.name);

    otherwise
        error('Unknown build method for fmri_model object. See help for valid strings.');

end

end % main function

% ----------------------------------------------
% ----------------------------------------------

% SUB-functions

% ----------------------------------------------
% ----------------------------------------------

function plot_subfcn(obj, Xs, C, names)

% ----------------------------------------------
% SETUP
% ----------------------------------------------

% conditions for unique color codes
nconds = size(obj.xX.cond_assignments, 2);

% bf = obj.xBF.bf;
% nbf = size(bf, 2);

colors = {'r' 'g' 'b' 'm' [1 .5 0], [0 1 .5], [0 .5 1], [1 0 .5]};
while length(colors) < nconds, colors = [colors colors]; end

disp('Plotting design matrix cols');

% ----------------------------------------------
% Plot
% ----------------------------------------------

han = plot_matrix_cols([Xs C], 'horiz');
axis tight
drawnow

% ***need to fix for Separate Session design

for i = 1:nconds
%     whstart = nbf * (i - 1) + 1;
%     whend = nbf * i;
    wh = find(obj.xX.cond_assignments(:, i));
    set(han(wh), 'Color', colors{i});
end

if ~isempty(C)
    names = [names {'Covariates'}];
end

yvals = linspace(1, size([Xs C], 2), nconds + 2);
yvals = yvals(2:end-1);

set(gca, 'YTick', yvals, 'YTickLabel', names)

xlabel('Time (TRs)')

title(obj.build_method);

drawnow

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%% Function List
%
%% ARMAX-GARCH Toolbox
%
% <matlab:web(fullfile('garch.m')) garch.m>, estimates the ARMAX-GARCH family of models
%
% Inputs: a data vector, the GARCH model to be estimated, the distribution, the order of AR, MA, ARCH and GARCH effects as well as a vector of factors for the mean and volatility process.
%
% Outputs: a vector of estimated parameters, standard errors estimated by the inverse Hessian, the Log-Likelihood value, a vector of conditional variances, a vector of residuals, and a summary of results which includes: model statistics, t-statistics, robust standard errors, scores among others.
%
% <matlab:web(fullfile('garchfind.m')) garchfind>, finds the combination of models and distributions that better fits the data based on a set of criteria (i.e. largest log likelihood value and the smallest AIC and BIC criteria).
%
% Inputs: a data vector, a model and distribution vectors, the order of AR, MA, ARCH and GARCH,
%
% Outputs: the best fit model
%
% % <matlab:web(fullfile('garchplot.m')) garchplot>, Plots the data,volatility and residuals
%
% <matlab:web(fullfile('garchfor.m')) garchfor> & <matlab:web(fullfile('garchfor2.m')) garchfor2>, estimates mean, volatility and kurtosis forecasts
%
% Inputs: a data vector, a vector of residuals, a vector of conditional variances, a set of parameters, the  model, the order of AR, MA, ARCH and GARCH and the number of forecasts
%
% Outputs: a vector of mean, and volatility forecasts as well as a vector of cumulative mean and volatility forecasts.
%
% <matlab:web(fullfile('garchsim.m')) garchsim>, simulates GARCH responses
%
% Inputs: estimated parameters, GARCH model, distribution, order of ARCH and GARCH, number of samples and number of paths. Additionally, a vector of time series of positive pre-sample conditional standard deviations may be provided, which the variance model will initialize.
%
% Outputs: a vector of simulated series and conditional standard deviations
%
% <matlab:web(fullfile('garchvar.m')) garchvar> & <matlab:web(fullfile('garchvar2.m')) garchvar2>, estimates  Value-at-Risk for a given confidence level and horizon period for both long and short positions
%
% Inputs: a data vector, a vector of residuals, a vector of conditional% variances, a set of parameters, the  model, the distribution, the order of AR, MA, ARCH and GARCH, the number of forecasts and confidence level
%
% Outputs: a vector of VaR forecasts
%
% <matlab:web(fullfile('garchvolfor.m')) garchvolfor>, an application in Volatility Forecasting & Value-at-Risk. It allows the comparison of volatility and Value-at-Risk estimates for a data vector and for a variety of GARCH models and distributions as specified in the model & distribution variables and at different forecast periods as well as sort the results according to only a sub-set of forecast periods.
%
% Inputs: a data vector, a vector of models and distributions, order of AR, MA, ARCH and GARCH, max forecasts, the forecasts of interest and a vector of a% VaR losses
%
% Outputs:  volatility and VaR back-testing results, vectors of forecasted estimated of volatility, VaR and returns
%
%% ARMAX-GARCH-K Toolbox
%
% <matlab:web(fullfile('garchk.m')) garchk>, estimates the ARMAX-GARCH-K family of models
%
% Inputs: a data vector, the GARCH model to be estimated, the order of AR, MA, ARCH and GARCH effects as well as a vector of factors for the mean and volatility process.
%
% Outputs: a vector of estimated parameters, standard errors estimated by the inverse Hessian, the Log-Likelihood value, a vector of conditional variances, a vector of conditional kurtosis, a vector of residuals, and a summary of results which includes: model statistics, t-statistics, robust standard errors, scores among others.
%
% <matlab:web(fullfile('garchkplot.m')) garchkplot>, Plots the data, conditional volatility, conditional kurtosis and residuals
%
% <matlab:web(fullfile('garchkfor.m')) garchkfor> & <matlab:web(fullfile('garchkfor2.m')) garchkfor2>, estimates mean, volatility and kurtosis forecasts
%
% Inputs: a data vector, a vector of residuals, a vector of conditional variances, a vector of conditional kurtosis, a set of parameters, the  model, the order of AR, MA, ARCH and GARCH and the number of forecasts
%
% Outputs: a vector of mean, volatility and kurtosis forecasts as well as a vector of cumulative mean and volatility forecasts.
%
% <matlab:web(fullfile('garchksim.m')) garchksim>, simulates volatility and kurtosis responses
%
% Inputs: estimated parameters, GARCH model, distribution, order of ARCH and GARCH, number of samples and number of paths. Additionally, a vector of time series of positive pre-sample conditional standard deviations may be provided, which the variance model will initialize.
%
% Outputs: a vector of simulated series, conditional standard deviations and conditional kurtosis
%
% <matlab:web(fullfile('garchkvar.m')) garchkvar> & <matlab:web(fullfile('garchkvar2.m')) garchkvar2>, estimates  Value-at-Risk for a given confidence level and horizon period for both long and short positions
%
% Inputs: a data vector, a vector of residuals, a vector of conditional variances, a set of parameters, the  model, the order of AR, MA, ARCH and GARCH, the number of forecasts and confidence level
%
% Outputs: a vector of VaR forecasts
%
% <..\readme\readme.html Return to Main>


text:
function s = findspan(n,p,u,U)
% FINDSPAN  Find the span of a B-Spline knot vector at a parametric point
% -------------------------------------------------------------------------
% ADAPTATION of FINDSPAN from C
% -------------------------------------------------------------------------
%
% Calling Sequence:
%
%   s = findspan(n,p,u,U)
%
%  INPUT:
%
%    n - number of control points - 1
%    p - spline degree
%    u - parametric point
%    U - knot sequence
%
%  RETURN:
%
%    s - knot span
%
%  Algorithm A2.1 from 'The NURBS BOOK' pg68

                                                % int findspan(int n, int p, double u, double *U) {

                                                %   int low, high, mid;
                                                %   // special case
if (u==U(n+2)), s=n; return,  end               %   if (u == U[n+1]) return(n);
                                                %
                                                %   // do binary search
low = p;                                        %   low = p;
high = n + 1;                                   %   high = n + 1;
mid = floor((low + high) / 2);                  %   mid = (low + high) / 2;
while (u < U(mid+1) || u >= U(mid+2))           %   while (u < U[mid] || u >= U[mid+1])  {
    if (u < U(mid+1))                           %     if (u < U[mid])
        high = mid;                             %       high = mid;
    else                                        %     else
        low = mid;                              %       low = mid;
    end
    mid = floor((low + high) / 2);              %     mid = (low + high) / 2;
end                                             %   }
                                                %
s = mid;                                        %   return(mid);
                                                %   }


text:
% UTILS Utility package for MTEST unit testing framework
%
% Array Comparison
%   compareFloats            - Compare floating-point arrays using tolerance
%
% Test Case Discovery Functions
%   isTestCaseSubclass       - True for name of TestCase subclass
%
% String Functions
%   containsRegexp           - True if string contains regular expression
%   isSetUpString            - True for string that looks like a setup function
%   isTearDownString         - True for string that looks like teardown function
%   isTestString             - True for string that looks like a test function
%
% Miscellaneous Functions
%   generateDoc              - Publish test scripts in mtest/doc
%   parseFloatAssertInputs   - Common input-parsing logic for several functions

% Undocumented Functions
%   isAlmostEqual        - Floating-point equality test using relative tolerance

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

author: ricedsp


text:
% Figure 6.2: Penalty function approximation
% Section 6.1.2
% Boyd & Vandenberghe "Convex Optimization"
% Original by Lieven Vandenberghe
% Adapted for CVX Argyris Zymnis - 10/2005
%
% Comparison of the ell1, ell2, deadzone-linear and log-barrier
% penalty functions for the approximation problem:
%       minimize phi(A*x-b),
%
% where phi(x) is the penalty function
% Log-barrier will be implemented in the future version of CVX

% Generate input data
randn('state',0);
m=100; n=30;
A = randn(m,n);
b = randn(m,1);

% ell_1 approximation
% minimize   ||Ax+b||_1
disp('ell-one approximation');
cvx_begin
    variable x1(n)
    minimize(norm(A*x1+b,1))
cvx_end

% ell_2 approximation
% minimize ||Ax+b||_2
disp('ell-2');
x2=-A\b;

% deadzone penalty approximation
% minimize sum(deadzone(Ax+b,0.5))
% deadzone(y,z) = max(abs(y)-z,0)
dz = 0.5;
disp('deadzone penalty');
cvx_begin
    variable xdz(n)
    minimize(sum(max(abs(A*xdz+b)-dz,0)))
cvx_end

% log-barrier penalty approximation
%
% minimize -sum log(1-(ai'*x+bi)^2)

disp('log-barrier')

% parameters for Newton Method & line search
alpha=.01; beta=.5;

% minimize linfty norm to get starting point
cvx_begin
    variable xlb(n)
    minimize norm(A*xlb+b,Inf)
cvx_end
linf = cvx_optval;
A = A/(1.1*linf);
b = b/(1.1*linf);

for iters = 1:50

   yp = 1 - (A*xlb+b);  ym = (A*xlb+b) + 1;
   f = -sum(log(yp)) - sum(log(ym));
   g = A'*(1./yp) - A'*(1./ym);
   H = A'*diag(1./(yp.^2) + 1./(ym.^2))*A;
   v = -H\g;
   fprime = g'*v;
   ntdecr = sqrt(-fprime);
   if (ntdecr < 1e-5), break; end;

   t = 1;
   newx = xlb + t*v;
   while ((min(1-(A*newx +b)) < 0) || (min((A*newx +b)+1) < 0))
       t = beta*t;
       newx = xlb + t*v;
   end;
   newf = -sum(log(1 - (A*newx+b))) - sum(log(1+(A*newx+b)));
   while (newf > f + alpha*t*fprime)
       t = beta*t;
       newx = xlb + t*v;
       newf = -sum(log(1-(A*newx+b))) - sum(log(1+(A*newx+b)));
   end;
   xlb = xlb+t*v;
end

% Plot histogram of residuals

ss = max(abs([A*x1+b; A*x2+b; A*xdz+b;  A*xlb+b]));
tt = -ceil(ss):0.05:ceil(ss);  % sets center for each bin
[N1,hist1] = hist(A*x1+b,tt);
[N2,hist2] = hist(A*x2+b,tt);
[N3,hist3] = hist(A*xdz+b,tt);
[N4,hist4] = hist(A*xlb+b,tt);

range_max=2.0;  rr=-range_max:1e-2:range_max;

figure(1), clf, hold off
subplot(4,1,1),
bar(hist1,N1);
hold on
plot(rr, abs(rr)*40/3, '-');
ylabel('p=1')
axis([-range_max range_max 0 40]);
hold off

subplot(4,1,2),
bar(hist2,N2);
hold on;
plot(rr,2*rr.^2),
ylabel('p=2')
axis([-range_max range_max 0 11]);
hold off

subplot(4,1,3),
bar(hist3,N3);
hold on
plot(rr,30/3*max(0,abs(rr)-dz))
ylabel('Deadzone')
axis([-range_max range_max 0 25]);
hold off

subplot(4,1,4),
bar(hist4,N4);
rr_lb=linspace(-1+(1e-6),1-(1e-6),600);
hold on
plot(rr_lb, -3*log(1-rr_lb.^2),rr,2*rr.^2,'--')
axis([-range_max range_max 0 11]);
ylabel('Log barrier'),
xlabel('r')
hold off

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function L = log_prob_node(CPD, self_ev, pev)
% LOG_PROB_NODE Compute prod_m log P(x(i,m)| x(pi_i,m), theta_i) for node i (root)
% L = log_prob_node(CPD, self_ev, pev)
%
% self_ev{m} is the evidence on this node in case m
% pev{i,m} is the evidence on the i'th parent in case m (ignored)
% We always return L = 0.

L = 0;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function centroid = voronoi_centroids ( n, first, list_num, list, xyz, v_num, v_xyz )

%*****************************************************************************80
%
%% VORONOI_CENTROIDS computes the centroids of each polygon in a Voronoi diagram.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 May 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of data points and Voronoi polygons.
%
%    Input, integer FIRST(N+1), for each polygon, points to the location
%    in LIST of the index.
%
%    Input, integer LIST_NUM, the number of items in LIST.
%
%    Input, integer LIST(LIST_NUM), the indices of vertices that form
%    each polygon.
%
%    Input, real XYZ(3,N), the coodinates of the data points.
%
%    Input, integer V_NUM, the number of Voronoi vertices.
%
%    Input, real V_XYZ(3,V_NUM), the coordinates of the Voronoi vertices.
%
%    Output, real CENTROID(3,N), the centroid of each Voronoi polygon.
%
  centroid = zeros ( 3, n );

  r = 1.0;

  for poly = 1 : n
%
%  Compute the area of each triangle formed by one side of the polygon
%  and the Delaunay point.
%
    v3 = list(first(poly+1)-1);

    for l = first(poly) : first(poly+1) - 1

      v2 = v3;
      v3 = list(l);

      stri_area = stri_vertices_to_area ( r, ...
        xyz(1:3,poly), v_xyz(1:3,v2), v_xyz(1:3,v3) );

      stri_centroid = stri_vertices_to_centroid ( r, ...
        xyz(1:3,poly), v_xyz(1:3,v2), v_xyz(1:3,v3) );

      centroid(1:3,poly) = centroid(1:3,poly) + stri_area * stri_centroid(1:3,1);

    end

    centroid(1:3,poly) = centroid(1:3,poly) / norm ( centroid(1:3,poly) );

  end

  return
end

author: johannesgerer


text:
function [doAlign, doResamp] = dtiInitAlign(dwParams,dwDir,doResamp)
%
%  doAlign = dtiInitAlign(dwParams,dwDir)
%
% Based on user selected params decide if we align the raw dwi data to a
% reference image (t1).
%
% WEB resources:
%   https://github.com/vistalab/vistasoft/wiki/DWI-Initialization
%   vistaBrowseGit('dtiInitAlign');
%
% (C) Stanford VISTA, 8/2011 [lmp]
%

%%
doAlign  = false;

% If we are overwriting or an alignment file does not exist then we signal
% to align and resample.
if dwParams.clobber == 1 || ~exist(dwDir.acpcFile,'file')
    doAlign  = true;
    doResamp = true;
else
    % Prompt the user to overwrite if clobber = 'ask'
    if dwParams.clobber == 0
        resp = questdlg([dwDir.acpcFile ' exists- would you like to overwrite it?'],...
            'Clobber AcPc','Overwrite','Use Existing File','Abort','Use Existing File');
        if(strcmpi(resp,'Abort')), error('User aborted.'); end
        if(strcmpi(resp,'Overwrite'))
            doAlign  = true;
            doResamp = true;
        end
    end
end

end


text:
function varargout=glotros(varargin)
%
% For calling details please see v_glotros.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_glotros.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: glotros.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_glotros(varargin{:});
else
    v_glotros(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function y = prod( x, dim )

%   Disciplined geometric programming information for PROD:
%      PROD(X) and PROD(X,DIM) are vectorized versions of multiplication
%      so in most cases it would be incompatible with DCPs. Therefore it
%      has not been implemented to support DCPs. DGPs however support
%      products more liberally. When PROD is used in a DCP, elements in
%      each subvector must satisfy the corresponding combination rule
%      for multiplication (see TIMES). For example, suppose that X looks
%      like this:
%         X = [ log-convex log-concave log-affine  ;
%               log-affine log-concave log-concave ]
%      Then PROD(X,1) would be permittted, but PROD(X,2) would not,
%      because the top row contains the product of log-convex and
%      log-concave terms, in violation of the DGP ruleset.

narginchk(1,2);

%
% Size check
%

try
    ox = x;
    if nargin < 2, dim = []; end
    [ x, sx, sy, zx, zy, nx, nv, perm ] = cvx_reduce_size( x, dim ); %#ok
catch exc
    error( exc.message );
end

%
% Quick exit for easy cases
%

if isempty( x ),
    if ~any( zy ), zy = 1; end
    y = cvx( ones( zy ) );
    return
elseif nx == 1,
    y = ox;
    return
end

%
% Type check
%

persistent remap_1 remap_2 remap_3 remap_0
if isempty( remap_3 ),
    remap_0 = cvx_remap( 'zero' );
    remap_1 = cvx_remap( 'constant' );
    remap_2 = cvx_remap( 'log-convex' );
    remap_3 = cvx_remap( 'log-concave' );
end
vx = cvx_reshape( cvx_classify( x ), sx );
t0 = any( reshape( remap_0( vx ), sx ) );
t1 = all( reshape( remap_1( vx ), sx ) );
t2 = all( reshape( remap_2( vx ), sx ) ) | ...
     all( reshape( remap_3( vx ), sx ) );
t3 = t0 & t1;
ta = t3 + 2 * ( t2 & ~t3 );
nu = sort( ta(:) );
nu = nu([true;diff(nu)~=0]);
nk = length( nu );

%
% Perform the computations
%

if nk > 1,
    y = cvx( [ 1, nv ], [] );
end
for k = 1 : nk,

    if nk == 1,
        xt = x;
    else
        tt = ta == nu( k );
        xt = cvx_subsref( x, ':', tt );
    end

    switch nu( k ),
        case 0,
            error( 'Disciplined convex programming error:\n   Invalid computation: prod( {%s} )', cvx_class( xt, true, true ) );
        case 1,
            yt = cvx( prod( cvx_constant( xt ) ) );
        case 2,
            yt = exp( sum( log( xt ) ) );
        otherwise,
            error( 'Shouldn''t be here.' );
    end

    if nk == 1,
        y = yt;
    else
        y = cvx_subsasgn( y, tt, yt );
    end

end

%
% Reverse the reshaping and permutation steps
%

y = reshape( y, sy );
if ~isempty( perm ),
    y = ipermute( y, perm );
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function j_progress(varargin)
% function fn_progress(prompt,max[,size])
% function fn_progress(prompt,max,'%' or 'p')
% function fn_progress(prompt)
% function fn_progress(i)
% function fn_progress('end')
% function fn_progress('cont')
% from Thomas Deneux, modified by Julien Cohen-Adad 2006-11-24
%
% example: j_progress('Pouf pouf. Please wait...');for i=1:50, j_progress(i/50), pause(.01), end;j_progress('elapsed')
%---
% progress indicator

persistent promptsize
persistent curi
persistent after
persistent size
persistent pflag
persistent max

if ischar(varargin{1}) % INITIALIZATION or special
    prompt = varargin{1};

    if ischar(prompt) % special cases
        switch prompt
            case 'end'
                fprintf(repmat('\b',1,promptsize+1+size+length(after)+1))
                return
            case 'cont'
                disp(sprintf(repmat(' ',1,promptsize+1+size+length(after))))
                return
            case 'elapsed'
                disp([sprintf([repmat('\b',1,size+length(after)+1) 'OK (elapsed time ' num2str(toc) ' seconds)'])])

%                 disp([sprintf(['elapsed ' num2str(toc) ' s']) after])

%                 disp(['elapsed ' num2str(toc) ' s'])
                return
            case 'elapsedmin'
                disp(['elapsed ' num2str(toc/60,'%.1f') 's.'])
                return
        end
    end

    % Input
    promptsize = length(prompt);
    i = 0; curi = i;
    if nargin<2
        pflag = true;
        max = 1;
        size = 3;
    else
        max = varargin{2};
        if ischar(max)
            max = str2num(max);
        end
        pflag = false;
        if nargin>2
            size = varargin{3};
            if strcmp(size,'%') || strcmp(size,'p')
                pflag = true;
                size = 3;
            elseif ischar(size)
                size = str2num(size);
            end
        else
            size = floor(log10(max))+1;
        end
    end
    csize = num2str(size);
    if pflag
        after = sprintf('%%');
    else
        after = sprintf(['/%' csize 'i'],max);
    end
    disp([sprintf([prompt ' %' csize 'i'],i) after])
    tic
else % STATE PROGRESS
    i = varargin{1};
    if pflag, i = floor(i/max*100); end
    if (i == curi), return, end
    curi = i;
    csize = num2str(size);
    disp([sprintf([repmat('\b',1,size+length(after)+1) '%' csize 'i'],i) after])
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [classMeans,clusterIndex,maxIterReached] = prtUtilKmeans(data,nClusters,varargin)
%[classMeans,clusterIndex] = prtUtilKmeans(data,nClusters)
%   Perform k-means clustering on the nObservations x nFeatures matrix data
%   using nClusters.  classMeans is a matrix of size nClusters (or less) x
%   nFeatures representing the classMeans, and clusterIndex is a
%   nObservations x 1 vector indicating the closest mean to each
%   observation in data.
%
%[classMeans,clusterIndex] = prtUtilKmeans(data,nClusters,param1,value1,...)
%   Enables inputs of parameter/value pairs as described below:
%
%   Parameters:
%       initialMeans - string or double matrix of size nClusters x
%       nFeatures.  If initialMeans is a string, it can be 'random' or
%       'plusplus'. random  which uses random samples from the data to
%       initialize the classMeans. plusplus uses the kMeans++ algorithm
%           Arthur, D. and Vassilvitskii, S. (2007).
%               k-means++: The advantages of careful seeding.
%               Proceedings of the eighteenth annual ACM-SIAM
%               symposium on Discrete algorithms. 1027?1035.
%       If initialMeans is a double matrix, the rows of the
%       matrix represetnt the initial class means, and nClusters is
%       ignored.  Default is 'plusplus'.
%
%       distanceMetricFn - prtDistance* function specifying the distance
%       metric to use.  distanceMetricFn(x,x) must return 0.  Default value
%       is @(data,centers)prtDistanceEuclidean(data,centers));
%
%       handleEmptyClusters - string specifying how to handle empty
%       clusters.  Allowed values are 'remove' and 'random'.  Default value
%       is 'remove'.
%
%       plotVisualization - bool value specifying whether to display a plot
%       of the current class centers and the corresponding data on every
%       iteration.  plotVisualization can also be a counting number, in
%       which case it specifies how often to update the plot.  Default
%       value is 'false'.
%
%   Example usage:
%
%       ds = prtDataGenBimodal(100);
%       data = ds.getObservations;
%       close all;
%       [classMeans,clusterIndex] = prtUtilKmeans(data,4,'plotVisualization',4);
%
%       close all;
%       [classMeans,clusterIndex] = prtUtilKmeans(data,4,'plotVisualization',4,'distanceMetricFn',@prtDistanceCityBlock);
%

p = inputParser;

%p.addParamValue('initialMeans',randn(50,2)/10);
p.addParamValue('initialMeans','random');
p.addParamValue('distanceMetricFn',@(data,centers)prtDistanceEuclidean(data,centers));
p.addParamValue('handleEmptyClusters','remove');
p.addParamValue('maxIterations',1000);
p.addParamValue('plotVisualization',false);
p.addParamValue('logicalMeans',false);

p.parse(varargin{:});
inputStructure = p.Results;

[nSamples,nDimensions] = size(data); %#ok<NASGU>
clusterIndexOld = nan(nSamples,1);
clusterIndex = [];

for iter = 1:inputStructure.maxIterations

    if iter == 1; %initialize
        if isa(inputStructure.initialMeans,'char')
            switch lower(inputStructure.initialMeans)
                case 'random'
                    randInds = max(1,ceil(rand(1,nClusters)*nSamples));
                    classMeans = data(randInds,:);
                case 'plusplus'

                    classMeans = nan(nClusters, size(data,2));
                    classMeans(1,:) = data(max(1,ceil(rand*nSamples)),:);
                    initDistanceMat = nan(nSamples,nClusters);
                    for iCluster = 2:nClusters
                        initDistanceMat(:,iCluster-1) = sum(bsxfun(@minus,data,classMeans(iCluster-1,:)).^2,2); % KMeans++ uses a squared euclidean distance.

                        minDistances = min(initDistanceMat,[],2);
                        drawProbabilities = minDistances./sum(minDistances);

                        newClusterMeanIndex = prtRvUtilRandomSample(drawProbabilities, 1);

                        classMeans(iCluster,:) = data(newClusterMeanIndex,:);
                    end

                otherwise
                    error('invalid');
            end
        elseif isnumeric(inputStructure.initialMeans)
            classMeans = inputStructure.initialMeans;
            nClusters = size(classMeans,1);
        else
            error('invalid');
        end
    else
        for clusterInd = 1:nClusters
            classMeans(clusterInd,:) = mean(data(clusterIndex == clusterInd,:),1);
        end
        if inputStructure.logicalMeans
            classMeans = classMeans>0.5;
        end
    end

    if ~mod(iter,inputStructure.plotVisualization);
        prtUtilKmeansPlotVisualization(data,classMeans,clusterIndex,inputStructure,iter);
    end

    distanceMat = inputStructure.distanceMetricFn(data,classMeans);
    [twiddle,clusterIndex] = min(distanceMat,[],2); %#ok<ASGLU>

    %Handle empty clusters:
    nMaxFixSteps = 10;
    for iFix = 1:nMaxFixSteps
        if length(unique(clusterIndex)) ~= nClusters
            invalidClusters = setdiff(1:nClusters,clusterIndex);
            validClusters = intersect(1:nClusters,clusterIndex);
            switch lower(inputStructure.handleEmptyClusters)
                case 'remove'
                    classMeans = classMeans(validClusters,:);
                    nClusters = nClusters - length(invalidClusters);
                    distanceMat = distanceMat(:,validClusters);
                    [twiddle,clusterIndex] = min(distanceMat,[],2); %#ok<ASGLU>
                case 'random'
                    randInds = max(1,ceil(rand(1,length(invalidClusters))*nSamples));
                    classMeans(invalidClusters,:) = data(randInds,:);
                    distanceMat =  inputStructure.distanceMetricFn(data,classMeans);
                    [twiddle,clusterIndex] = min(distanceMat,[],2); %#ok<ASGLU>
                otherwise
                    error('invalid');
            end
        else
            break
        end
    end
    if iFix == nMaxFixSteps
        maxIterReached = true;
        return
    end
    %Check convergence:
    if all(clusterIndexOld == clusterIndex)
        maxIterReached = false;
        return;
    else
        clusterIndexOld = clusterIndex;
    end

end
maxIterReached = true;

function prtUtilKmeansPlotVisualization(data,classMeans,clusterIndex,inputStructure,iter)

[nSamples,nDimensions] = size(data); %#ok<ASGLU>
if iter == 1
    distanceMat =  inputStructure.distanceMetricFn(data,classMeans);
    [twiddle,clusterIndex] = min(distanceMat,[],2); %#ok<ASGLU>
end
if nDimensions > 3
    warning('prt:prtUtilKmeans','plotVisualization is true, but dimensionality of data (%d) is > 3',nDimensions);
    return;
end
ds = prtDataSetClass(data,clusterIndex);

plot(ds);
hold on;
switch nDimensions
    case 1
        plot(classMeans,'b.');
    case 2
        plot(classMeans(:,1),classMeans(:,2),'b.');
    case 3
        plot3(classMeans(:,1),classMeans(:,2),classMeans(:,3),'b.');
end
hold off;
title(sprintf('Iteration %d',iter));
drawnow; %pause;

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
classdef featureMatrices < handle

% featureMatrices object for ECGwrapper
% -------------------------------------
%
% Description:
%
%
% Arguments: (specified as ECGwrapper('arg_name1', arg_val1, ... , 'arg_nameN', arg_valN) )
%
%     a. ECG specification:
%
%         a.1 Recording filename where ECG is stored in a valid format.
%
%           + recording_name: ECG recording to be classified.
%           + recording_format : Valid ECG format. (MIT, ISHNE, AHA, HES, MAT)
%
%     b. Operating modes
%
%
%     c. Modifiers
%
%
% Output:
%
% Examples:
%
%
% Limits and Known bugs:
%   Probably a lot :( ... but dont panic! send me feedback if you need
%   help.
%
% Author: Mariano Llamedo Soria (llamedom at frba.utn.edu.ar)
% Version: 0.1 beta
% Birthdate  : 20/3/2012
% Last update: 20/3/2012

    properties(GetAccess = private, Constant)
        sampling_rate = 360; %Hz
        fnPayload = {'featMat_clust' 'featMat_ldc' 'ECG' 'QRS_locations' 'RR_intervals'};
        fnFMfile = {'featMat_clust' 'featMat_ldc'};
        fnECGfile = {'ECG' 'QRS_locations' 'RR_intervals'};
        SmallValue = 0.0001;
    end

    properties ( Access = private )
        delayLP
        autovec
        q_filters
        LP
        LPsize
        scale_idx
        scales
        global_struct
    end

    properties (SetAccess = private, GetAccess = public)
        %read-only
        name
        class_labeling
        true_labels
        %auxiliary indexes
        file_count = 1;
        File_list = [];
        File_idx = [];
        QRS_idx = [];

    end

    properties
        %read-write
        recording_name
        recording_format
        progress_hdl
        tmp_path
        bCache
        cant_pids
        this_pid
        cant_iteration
        this_iteration
    end

    methods

        function obj = featureMatrices(class_labeling)

            %% Constants and definitions

            obj.name = 'feat_matrix_construction';
            obj.class_labeling = class_labeling;

        end

        function Start(obj)

            %% Calculate accesories signals

            %load low-pass preprocessing filter.
            load( ['LP Fs ' num2str(obj.sampling_rate) 'Hz - Fc 35Hz - 80db stop.mat' ]);
            obj.LPsize = length(LP.numerator);
            obj.delayLP = round((obj.LPsize-1)/2);
            obj.LP = LP;
            clear LP

            % Wavelet transform filter bank design
            obj.scales = 1:6;
            CantScales = length(obj.scales);
            MaxScales = max(obj.scales);
            obj.scale_idx= nan(MaxScales,1);

            filters_cache_filename = ['wt_filters_' num2str(MaxScales) ' scales_' num2str(obj.sampling_rate) ' Hz.mat' ];
            if( exist(filters_cache_filename, 'file') )
                aux = load( filters_cache_filename );
                obj.q_filters = aux.q_filters;
                clear aux
            else
                obj.q_filters = qs_filter_design(MaxScales, obj.sampling_rate);
            end

            for ii = 1:CantScales
                %indice para saber que escala corresponde a cada columna de MyWavelets
                obj.scale_idx(obj.scales(ii)) = ii;
            end

            % Update point
            obj.progress_hdl.checkpoint('PCA projection');

%             obj.autovec = PCA_proj_basis(obj.recording_name, obj.recording_format, obj.q_filters);

            obj.autovec = [-0.696609733483207 -0.403719836632163 -0.593081084444745;-0.487887339191286 -0.339521516938205 0.804171053814316;-0.526023595928265 0.849550135686930 0.0395441965529324;];

            % data structure to keep data through all iterations.
            obj.global_struct = [];

        end

        function payload = Process(obj, ECG, this_header, this_iter_ECG_relative_start_end_idx, this_iter_annotations)

            obj.progress_hdl.checkpoint('Label parsing');

            [QRS_locations, obj.true_labels ] = Annotation_process(this_iter_annotations, obj.recording_format, obj.class_labeling);

            % RR interval sequence
            RR_intervals = diff(QRS_locations, 1);
            RR_intervals = colvec([ RR_intervals(1); RR_intervals ]);

            dRR_intervals = diff(RR_intervals, 1);
            dRR_intervals = colvec([ dRR_intervals(2); dRR_intervals(2); dRR_intervals(2:end) ]);

            % interval resampling.
            RR_intervals = RR_intervals * obj.sampling_rate / this_header.freq;
            dRR_intervals = dRR_intervals * obj.sampling_rate / this_header.freq;

            %% Preprocessing

            % Update point
            obj.progress_hdl.checkpoint('Resampling');

            cant_QRS_locations = length(QRS_locations);
            this_iter_QRS_seq_idx = 1:cant_QRS_locations;

            %resample to obj.sampling_rate Hz
            ECG = resample(ECG, obj.sampling_rate, this_header.freq);
            QRS_locations = colvec(round(QRS_locations * obj.sampling_rate / this_header.freq));
        %         this_iter_true_labels = true_labels(this_iter_QRS_seq_idx);
            this_iter_ECG_resampled_size = size(ECG,1);

            if( this_header.nsig > 2 )
                %multilead approach.
                % project ECG and wtECG to obtain PCA components
                ECG = ECG * obj.autovec(:,1:2);
        %             wtECG = cellfun(@(a)( a * obj.autovec(:,1:2) ), mat2cell(wtECG, this_iter_ECG_resampled_size, obj.ECG_header.nsig, ones(1,CantScales) ), 'UniformOutput', false);
        %             wtECG = cell2mat(wtECG);
            end

            % Update point
            obj.progress_hdl.checkpoint('Filtering');

            %Low pass filtering @ 35Hz => ECG recovery
            ECG = filter(obj.LP, ECG);
            %delay compensation.
            ECG = [ zeros(obj.delayLP, size(ECG,2) ) ;ECG(obj.LPsize:end,:) ; zeros(obj.delayLP, size(ECG,2))];

            %Quito la linea de base.
        %     ECG = BaselineWanderRemovalMedian( ECG, obj.sampling_rate);
            ECG = BaselineWanderRemovalSplines( ECG, QRS_locations, obj.sampling_rate);

            % Update point
            obj.progress_hdl.checkpoint('Wavelet transform calculation');

            wtECG = qs_wt(ECG, obj.scales, obj.sampling_rate, obj.q_filters);

            %% Features Calculation

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            this_iter_seq_idx = 1:length(QRS_locations);

            % log(RRactual)
            %%%%%%%%%%%%%%%
            featMat_ldc = log(colvec(RR_intervals(this_iter_QRS_seq_idx))/obj.sampling_rate);

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % log(RRpost)
            %%%%%%%%%%%%%%%

            featMat_ldc = [ featMat_ldc log(colvec(RR_intervals(min(cant_QRS_locations, this_iter_QRS_seq_idx+1)))/obj.sampling_rate)];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % log(RRmean1)
            %%%%%%%%%%%%%%%

            aux_idx = arrayfun(@(a)( find(  QRS_locations >= (QRS_locations(a) - 1*60*obj.sampling_rate) & ...
                                            QRS_locations <= QRS_locations(a) )), ...
                               this_iter_seq_idx, 'UniformOutput', false);

            aux_featval = cellfun(@(a)(mean(RR_intervals(a))/obj.sampling_rate), aux_idx);
            featMat_ldc = [ featMat_ldc colvec(log(aux_featval)) ];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % log(RRmean20)
            %%%%%%%%%%%%%%%

            aux_idx = arrayfun(@(a)( find(  QRS_locations >= (QRS_locations(a) - 20*60*obj.sampling_rate) & ...
                                            QRS_locations <= QRS_locations(a) )), ...
                               this_iter_QRS_seq_idx, 'UniformOutput', false);

            aux_featval = cellfun(@(a)(mean(RR_intervals(a))/obj.sampling_rate), aux_idx);
            featMat_ldc = [ featMat_ldc colvec(log(aux_featval)) ];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % AutoCorr1PlanoWTZeroCross
            %%%%%%%%%%%%%%%
            % AutoCorr1PlanoWTModMaxPos
            %%%%%%%%%%%%%%%

            aux_idx = arrayfun(@(a)( max(1, QRS_locations(a) - round(0.08*obj.sampling_rate)): ...
                                     min(this_iter_ECG_resampled_size, QRS_locations(a) + round(0.08*obj.sampling_rate))) , ...
                               this_iter_seq_idx, 'UniformOutput', false);

            aux_featval = cellfun(@(a)(squeeze(wtECG(a,:,obj.scale_idx(4)))), aux_idx, 'UniformOutput', false);

            %calculate PCA matrix in this slices for feature
            %AutoCorr1PlanoWTModMaxPos
            autovec_slices = cellfun(@(a)(autovec_calculation(a)), aux_featval, 'UniformOutput', false );

            aux_idx = arrayfun(@(a)( max(1, QRS_locations(a) - round(0.13*obj.sampling_rate)): ...
                                     min(this_iter_ECG_resampled_size, QRS_locations(a) + round(0.2*obj.sampling_rate))) , ...
                               this_iter_seq_idx, 'UniformOutput', false);

            aux_featval = cellfun(@(a,b)(squeeze(wtECG(a,:,obj.scale_idx(4))) * b), aux_idx, autovec_slices, 'UniformOutput', false);
            [ aux_mp aux_zc ] = cellfun(@(a)(CalcModMaxPos(a(:,1))), aux_featval );

            featMat_ldc = [ featMat_ldc colvec(aux_zc) colvec(aux_mp) ];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % AutoCorr2PlanoWTZeroCross
            %%%%%%%%%%%%%%%
            % AutoCorr2PlanoWTModMaxPos
            %%%%%%%%%%%%%%%

            [ aux_mp aux_zc ] = cellfun(@(a)(CalcModMaxPos(a(:,2))), aux_featval );

            featMat_ldc = [ featMat_ldc colvec(aux_zc) colvec(aux_mp) ];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            %clear some auxiliar variables.
            clear aux*

            %calculate clustering features.

            % log(RRactual)
            %%%%%%%%%%%%%%%
            featMat_clust = featMat_ldc(:,1);

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % log(RRant)
            %%%%%%%%%%%%%%%
            featMat_clust = [ featMat_clust log(colvec(RR_intervals(max(1, this_iter_QRS_seq_idx-1)))/obj.sampling_rate)];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % log(Prematuridad_loc)
            %%%%%%%%%%%%%%%

            aux_idx = arrayfun(@(a)(max(1,a-1):min(cant_QRS_locations,a+1)), this_iter_QRS_seq_idx, 'UniformOutput', false);
            %tengo que contemplar los casos extremos
            if( length(aux_idx{1}) < 3 )
                aux_idx{1} = [1 aux_idx{1}];
            end
            if( length(aux_idx{end}) < 3 )
                aux_idx{end} = [aux_idx{end} cant_QRS_locations];
            end

            aux_featval = cellfun(@(a)( RR_intervals(a(2))/sum(RR_intervals(a)) ), aux_idx);
            featMat_clust = [ featMat_clust colvec(log(aux_featval)) ];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % log(AbsRRlocalVar)
            %%%%%%%%%%%%%%%

            aux_idx = arrayfun(@(a)(max(1,a-1):min(cant_QRS_locations,a+1)), this_iter_QRS_seq_idx, 'UniformOutput', false);
            if( length(aux_idx{1}) < 3 )
                aux_idx{1} = [1 aux_idx{1}];
            end
            if( length(aux_idx{end}) < 3 )
                aux_idx{end} = [aux_idx{end} cant_QRS_locations];
            end

            aux_featval = cellfun(@(a)(obj.SmallValue+sum(abs(dRR_intervals(a)))/obj.sampling_rate), aux_idx);

            featMat_clust = [ featMat_clust colvec(log(aux_featval)) ];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % log(RRmean20)
            %%%%%%%%%%%%%%%

            featMat_clust = [ featMat_clust featMat_ldc(:,4) ];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % log(QRS_max_scale_proj12)
            %%%%%%%%%%%%%%%

            aux_idx = arrayfun(@(a)( max(1, QRS_locations(a) - round(0.08*obj.sampling_rate)): ...
                                     min(this_iter_ECG_resampled_size, QRS_locations(a) + round(0.08*obj.sampling_rate))) , ...
                               this_iter_seq_idx, 'UniformOutput', false);

            aux_featval = cellfun(@(a)(wtECG(a,:,:)), aux_idx, 'UniformOutput', false);

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            aux_featval = cellfun(@(a)(squeeze(sum(abs(a(:,1,:))))'), aux_featval, 'UniformOutput', false );
            aux_featval = cellfun(@(a)(obj.scales * colvec(a) ./ sum(a)), aux_featval );

        %     iAreaAbs = squeeze(sum(abs(wtAux)))';
        %     MaxProjArea = obj.scales * iAreaAbs ./ sum(iAreaAbs);

            featMat_clust = [ featMat_clust colvec(log(aux_featval)) ];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % AutoCorr1PlanoWTModMaxPos
            %%%%%%%%%%%%%%%

            aux_idx = arrayfun(@(a)( max(1, QRS_locations(a) - round(0.13*obj.sampling_rate)): ...
                                     min(this_iter_ECG_resampled_size, QRS_locations(a) + round(0.2*obj.sampling_rate))) , ...
                               this_iter_seq_idx, 'UniformOutput', false);

            aux_featval = cellfun(@(a,b)(squeeze(wtECG(a,:,obj.scale_idx(4))) * b), aux_idx, autovec_slices, 'UniformOutput', false);
            aux_featval = cellfun(@(a)(CalcModMaxPos(a(:,1))), aux_featval);

            featMat_clust = [ featMat_clust colvec(aux_featval) ];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % CrossCorr_QRST_ModMaxVal13
            %%%%%%%%%%%%%%%

            aux_idx = arrayfun(@(a)( max(1, QRS_locations(a) - round(0.13*obj.sampling_rate)): ...
                                     min( [ this_iter_ECG_resampled_size, ...
                                            QRS_locations(a) + round(0.4*obj.sampling_rate), ...
                                            QRS_locations(a) + RR_intervals(a) - round(0.13*obj.sampling_rate) ] ) ), ...
                               this_iter_seq_idx, 'UniformOutput', false);

            if( this_header.nsig > 2 )
                % wtECG already projected in obj.autovec.
                aux_featval = cellfun(@(a)( CalcModMaxVal( squeeze(wtECG(a,:,obj.scale_idx(3))) ) ), aux_idx);
            else
                aux_featval = cellfun(@(a)( CalcModMaxVal( squeeze(wtECG(a,:,obj.scale_idx(3))) * obj.autovec ) ), aux_idx);
            end

            featMat_clust = [ featMat_clust colvec(aux_featval) ];

            %save ECG for user interface and cluster analysis
            for ii = 1:length(obj.fnPayload)
                payload.(obj.fnPayload{ii}) = eval(obj.fnPayload{ii});
            end

        end

        function allPayloads = Concatenate(obj, allPayloads, payload )

            if( isempty(allPayloads) )
                for ii = 1:length(obj.fnFMfile)
                     allPayloads.(obj.fnFMfile{ii}) = payload.(obj.fnFMfile{ii});
                end
            else
                for ii = 1:length(obj.fnFMfile)
                     allPayloads.(obj.fnFMfile{ii}) = [allPayloads.(obj.fnFMfile{ii}); payload.(obj.fnFMfile{ii})];
                end
            end

            [~, rec_filename] = fileparts(obj.recording_name);

            for ii = 1:length(obj.fnECGfile)
                 aux_payload.(obj.fnECGfile{ii}) = payload.(obj.fnECGfile{ii});
            end

            save([obj.tmp_path 'tmpfile_' rec_filename '_ECG_cantpids_' num2str(obj.cant_pids) '_' num2str(obj.this_pid) '_' num2str(obj.this_iteration) '_' num2str(obj.cant_iteration)  '.mat'], '-struct', 'aux_payload' );

            aux_ECG_filename = [obj.tmp_path 'tmpfile_' rec_filename '_ECG_cantpids_' num2str(obj.cant_pids) '_thispid_' num2str(obj.this_pid) '_iteration_' num2str(obj.this_iteration) '_of_' num2str(obj.cant_iteration)  '.mat'];
            %do this for cluster supercomputers reasons
            movefile(   [obj.tmp_path 'tmpfile_' rec_filename '_ECG_cantpids_' num2str(obj.cant_pids) '_' num2str(obj.this_pid) '_' num2str(obj.this_iteration) '_' num2str(obj.cant_iteration)  '.mat'], ...
                aux_ECG_filename, 'f' );

            cant_QRS = length(payload.QRS_locations);
            obj.QRS_idx = [obj.QRS_idx; uint16(colvec(1:cant_QRS))];
            obj.File_idx = [obj.File_idx; repmat(uint16(obj.file_count), cant_QRS, 1)];
            obj.File_list = [obj.File_list; cellstr(aux_ECG_filename)];
            obj.file_count = obj.file_count + 1;

        end

    end

end

function autovec = autovec_calculation(wtECGslice)

    mean_wtECGslice = mean(wtECGslice);
    wtECGslice_cov = cov( bsxfun( @minus, wtECGslice, mean_wtECGslice ));
    [autovec autoval] = eig(wtECGslice_cov);
    [~, autoval_idx] = sort(diag(autoval), 'descend');
    autovec = autovec(:,autoval_idx);

end

function [ ModMaxPos ZeroCross ] = CalcModMaxPos( wtSignal )

    ModMaxPos = nan;
    ZeroCross = nan;

    lwtSignal = size(wtSignal, 1);

    if( all( wtSignal == 0 ) )
        return
    end

    ac1 = conv(wtSignal(:,1), flipud(wtSignal(:,1)));

    ac1 = ac1(lwtSignal:end);
    ac1 = 1/ac1(1)*ac1;

    interp_idx = (1:0.1:lwtSignal)';

    ac1_interp = spline( 1:lwtSignal, ac1, interp_idx );
    iZeroCrossPos_ac1 = myzerocros(ac1);

    if( isempty(iZeroCrossPos_ac1))
        return
    else
        if( iZeroCrossPos_ac1(1) < lwtSignal )
            if( sign(ac1(iZeroCrossPos_ac1(1))) ~= sign(ac1(iZeroCrossPos_ac1(1)+1)) )
                ZeroCross = iZeroCrossPos_ac1(1) - ac1(iZeroCrossPos_ac1(1)) / ( ac1(iZeroCrossPos_ac1(1)+1) - ac1(iZeroCrossPos_ac1(1)) );
            else
                ZeroCross = (iZeroCrossPos_ac1(1)-1) - ac1(iZeroCrossPos_ac1(1)-1) / ( ac1(iZeroCrossPos_ac1(1)) - ac1(iZeroCrossPos_ac1(1)-1) );
            end
        else
            return
        end
    end

    ModMaxPos_ac1 = modmax( ac1_interp, 2, 0, 0 );

    if( ~isempty(ModMaxPos_ac1) )
        valid_ac1_max_idx = find(interp_idx(ModMaxPos_ac1) > ZeroCross);
        if( ~isempty(valid_ac1_max_idx) )
            ModMaxPos = interp_idx(ModMaxPos_ac1(valid_ac1_max_idx(1)));
        else
            return
        end
    else
        return
    end

end

function ModMaxVal = CalcModMaxVal( wtSignal )

    % std_wtSignal = std(wtSignal);
    % wtSignal = bsxfun( @times, wtSignal, 1./std_wtSignal);

    cc = conv(wtSignal(:,1), flipud(wtSignal(:,2)));

    [~, max_pos] = max(abs(cc));

    ModMaxVal = cc(max_pos);

end

author: marianux

repo: ecg-kit

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


text:
%book : Signals and Systems Laboratory with MATLAB
%by Alex Palamides & Anastasia Veloni

% Fourier Series properties

%time shifting

t0=0;
T=10;
w=2*pi/T;
syms t
x=t*exp(-5*t)
k=-5:5;
a=(1/T)*int(x*exp(-j*k*w*t),t,t0,t0+T);
a1=eval(a);
subplot(211);
stem(k,abs(a1));
title(' Coefficients of x(t)=te^-^5^t');
legend('Magnitude');
subplot(212);
stem(k,angle(a1));
legend('Angle');

figure
t1=3;
right= exp(-j*k*w*t1).*a;
right =eval(right);
subplot(211);
stem(k,abs(right));
legend('Magnitude');
title('Right part');
subplot(212);
stem(k,angle(right));
legend('Angle');

figure
x=(t-t1).*exp(-5*(t-t1));
a=(1/T)*int(x*exp(-j*k*w*t),t,t0+t1,t0+T+t1);
coe=eval(a);
subplot(211);
stem(k,abs(coe));
legend('Magnitude');
title(' Coefficient of (t-3)exp(-5(t-3)) ');
subplot(212);
stem(k,angle(coe));
legend('Angle');

author: Sable

repo: mcbench-benchmarks


text:
function C = modwt_brick_wall(X, wavelet, N)
%%
%% Purpose:  Sets the first N_j coefficients to NaN, for j = 1,...,J.
%% -------------------------------------------------------------------------
%% Reference: Lindsay et al. (1996).  The Discrete Wavelet Transform
%%            and the Scale Anlaysis of the Surface Properties of Sea
%%            Ice.  IEEE Trans. on Geo. and Rem. Sen., 34(3), pp.
%%            771-787.
%%
%% Input: X        Matrix containing wavelet coefficients with appropriate
%%                 boundary condition
%%        wavelet  Character string; 'haar', 'd4', 'la8', 'la16'
%%        N        Length of original vector of observations
%%
%% Output: C  Matrix containing wavelet coefficients where ones affected
%%            by boundary conditions are replaced with NaNs
%%
[h, g, l] = myfilter(wavelet);
[I, J] = size(X);

C = X;
for j = 1:(J-1)
  n = (2.^j - 1) * (l - 1);
  C(1:n,j) = NaN;
end
C(1:n,j+1) = NaN;


text:
disp('Running: TEST_dynQuadRotor3d.m')
s = zeros(12,1) ; % [x, y, z, pitch, roll, yaw, dx, dy, dz, dpitch, droll, dyaw]
u = zeros(4,1) ;

% Coordinate Frame Note:
% World frame is XYZ = [East, North, Up].
% Body frame is XYZ = [port, nose, top], so attitude Euler angles are in
% the order: [pitch, roll, yaw]

% Define environmental and plant model params
% Enviromental params
p.g = -9.81 ; % World Coords is XYZ = [East, North, Up], i.e. gravity is a negative number
p.rho = 1.225 ;

% Inertial params
p.m = 5 ;
p.I = [0.625 0 0; 0 0.625 0; 0 0 1.25] ; % inertia tensor
p.cg = [0 0 0] ; % (m) location of center of gravity

% Propulsion params
qRP.d_prop = 0.305*ones(4,1) ; % propeller diameter (m)
qRP.maxThrust = 25*ones(4,1) ; % thrust at 100% throttle (N)
qRP.maxRPM = 10000*ones(4,1) ; % RPM at 100% throttle (RPM)
qRP.maxTorque = ones(4,1) ;  % torque at 100% throttle (Nm)
qRP.thrustLocations = [0.5 0 0; 0 0.5 0; -0.5 0 0; 0 -0.5 0]; % motor locations (each row one motor in coords: [port, nose, top]
qRP.thrustAxes = repmat([0 0 1],4,1) ; % thrust axes of each motor in coords port, nose, top.
qRP.isSpinDirectionCCW = [1; 0; 1; 0] ; % bool to reverse motor spin direction around 'thrustAxes'.

% Call function that creates the propulsion plant model
[p.propulsion] = definePropulsionModel(qRP); clear variable qRP ;

% Throttle = 0 %
disp('Test 0 - 0% throttle')
[ds] = dynQuadRotor3d(s, u, p) ;

%% Single timestep - full throttle
disp('Test 1 - 100% throttle')
s = zeros(12,1) ;
u = ones(4,1) ;
[ds] = dynQuadRotor3d(s, u, p) ;

%% multiple timesteps
disp('Test 2 - 0% throttle')
u = zeros(4,10) ;
s = zeros(12,10) ;
[ds] = dynQuadRotor3d(s, u, p) ;

%% multiple timesteps
disp('Test 3 - 0% throttle')
u = ones(4,10) ;
s = zeros(12,10) ;
s(4,:) = linspace(0,1,10) ;
[ds] = dynQuadRotor3d(s, u, p) ;

%% Effect of aircraft orientation
% positive pitch; should result in an acceleration 'south' (i.e. negative
% north => ds(8) < 0 ) ;
disp('Test 4 - 100% throttle, pitch = 5 deg')
deflect = 5 ;
s = zeros(12,10) ;
s(4,:) = deg2rad(deflect) ;
u = ones(4,10) ;
[ds] = dynQuadRotor3d(s, u, p) ;

%% negative pitch; should result in an acceleration 'north'
disp('Test 5 - 100% throttle, pitch = -5 deg')
deflect = -5 ;
s = zeros(12,10) ;
s(4,:) = deg2rad(deflect) ;
u = ones(4,10)  ;
[ds] = dynQuadRotor3d(s, u, p) ;

%% positive roll; should result in an acceleration 'east', i.e. ds(7) > 0
disp('Test 6 - 100% throttle, roll = 5 deg')
deflect = 5 ;
s = zeros(12,10) ;
s(5,:) = deg2rad(deflect) ;
u = ones(4,10)  ;
[ds] = dynQuadRotor3d(s, u, p) ;

%% negative roll; should result in an acceleration 'west', i.e. negative element 7
disp('Test 7 - 100% throttle, roll = -5 deg')
deflect = -5 ;
s = zeros(12,1) ;
s(5) = deg2rad(deflect) ;
u = ones(4,1)  ;
[ds] = dynQuadRotor3d(s, u, p) ;

%% positive yaw; should not affect acceleration direction
disp('Test 8 - 100% throttle, yaw = 5 deg')
deflect = 5 ;
s = zeros(12,1) ;
s(6) = deg2rad(deflect) ;
u = ones(4,1)   ;
[ds] = dynQuadRotor3d(s, u, p) ;

%% negative yaw; should not affect acceleration direction
disp('Test 9 - 100% throttle, yaw = -5 deg')
deflect = -5 ;
s = zeros(12,1) ;
s(6) = deg2rad(deflect) ;
u = ones(4,1)  ;
[ds] = dynQuadRotor3d(s, u, p) ;

%% wide vector
disp('Test 10 - ramp all throttles')
s = zeros(12,100) ;
u = repmat(linspace(0,1,100),4,1) ;
[ds] = dynQuadRotor3d(s, u, p) ;

%%
disp('TEST_dynQuadRotor3d.m ran without error')

author: MatthewPeterKelly

repo: OptimTraj


text:
function [F,changed] = convertlogics(F)
%CONVERTLOGICS Internal function to convert logic constraints to mixed integer constraints

changed = 0;
if length(F)>0
    extvariables = yalmip('logicextvariables');
    if ~isempty(extvariables)
        for i = 1:length(F)
            if is(F(i),'elementwise')
                Fi = sdpvar(F(i));
                Fv =getvariables(Fi);
                if length(Fv)==1
                    xb = getbase(Fi);
                    if isequal(xb,[0 1])
                        if ismember(Fv,extvariables)
                            F(i) = (Fi >= 1);
                        end
                    end
                end
            end
        end
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function c = acorr2(x,maxlag)
%ACORR Estimate autocorrelation function of time series
%
%   C = ACORR(X,MAXLAG) returns normalized autocorrelation
%   sequences for each column of X computing correlation via FFT
%
% Copyright (C) 2000-2008 Aki Vehtari
%
% This software is distributed under the GNU General Public
% Licence (version 3 or later); please refer to the file
% Licence.txt, included with the software, for details.

if nargin < 1
  error('Not enough input arguments.');
end
if nargin < 2
  maxlag=length(x)-1;
end
[m,n]=size(x);
c=zeros(maxlag,n);
for i1=1:n
  xn = x(:,i1)-mean(x(:,i1));
  xf = fft(xn,2^nextpow2(2*m-1));
  cf = ifft(abs(xf).^2);
  c(:,i1) = cf(2:(maxlag+1))./cf(1);
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
%% (Internal) Reshape a matrix into a row vector
% Return x as a row vector. This function is useful when a function returns a
% column vector or matrix and you want to immediately reshape it in a functional
% way. Suppose f(a,b,c) returns a column vector, matlab will not let you write
% f(a,b,c)(:)' - you would have to first store the result. With this function you
% can write rowvec(f(a,b,c)) and be assured that the result is a row vector.
%
%    x = colvec(x)
%
% Arguments:
%
%      + x: data
%
% Output:
%
%      + x: columnized data
%
% Example:
%
% See also colvec
%
% This file is from pmtk3.googlecode.com
%
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software Foundation, Inc.,
% 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
%
% Copyright 2008-2015
%
function x = rowvec(x)
    x = x(:)';
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox