text:
function chrpak_test067 ( )

%*****************************************************************************80
%
%% TEST067 tests I4_TO_XXDECODE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST067\n' );
  fprintf ( 1, '  I4_TO_XXDECODE: I => XXDECODE character\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The XXDECODE alphabet\n' );
  fprintf ( 1, '\n' );

  fprintf ( 1, '    "' );
  for i4 = 1 : 64
    c = i4_to_xxdecode ( i4 );
    fprintf ( 1, '%c', c );
  end
  fprintf ( 1, '"\n' );

  return
end


text:
% union_bc - union backward compatible with Matlab versions prior to 2013a

% Copyright (C) 2013 Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [C,IA,IB] = union_bc(A,B,varargin)

errorFlag = error_bc;

v = version;
indp = find(v == '.');
v = str2num(v(1:indp(2)-1));
if v > 7.19, v = floor(v) + rem(v,1)/10; end

if nargin > 2
    ind = strmatch('legacy', varargin);
    if ~isempty(ind)
        varargin(ind) = [];
    end
end

if v >= 7.14
    [C,IA,IB] = union(A,B,varargin{:},'legacy');
    if errorFlag
        [C2,IA2,IB2] = union(A,B,varargin{:});
        if (~isequal(C, C2) || ~isequal(IA, IA2) || ~isequal(IB, IB2))
            warning('backward compatibility issue with call to union function');
        end
    end
else
    [C,IA,IB] = union(A,B,varargin{:});
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
% Impulse response invariant discretization of fractional second
% order filter.
%
% irid_fsof function is prepared to compute a discrete-time finite
% dimensional (z) transfer function to approximate a continuous-time
% fractional second order low-pass filter (LPF) [1/(s^2 + a*s + b)]^r, where "s" is
% the Laplace transform variable; "r" is a real number in the range of
% (0,1); a and b are the time constant of LPF [1/(s^2 + a*s + b)]^r, where a, b >= 0.
%
% The proposed approximation keeps the impulse response "invariant"
%
% IN:
%       a, b: the time constant of (the first order) LPF
%             (a and b are arbitrary positive real numbers)
%       r: the fractional order in (0,1)
%       Ts: the sampling period
%       norder: the finite order of the approximate z-transfer function
%       (the orders of denominator and numerator z-polynomial are the same)
% OUT:
%       sr: returns the LTI object that approximates the [1/(s^2 + a*s + b)]^r
%           in the sense of invariant impulse response.
% TEST CODE
% [sr]=irid_fsof(0.01,3,2,.8,5);
%
% Reference:
%       (1) Yan Li, Hu Sheng and YangQuan Chen.
%           "Analytical Impulse Response of A Fractional Second Order Filter
%           and Its Impulse Response Invariant Discretization".
%       (2) YangQuan Chen.
%           "Impulse-invariant discretization of fractional order low-pass filters".
%           Sept. 2008. CSOIS AFC (Applied Fractional Calculus) Seminar.
%           http://fractionalcalculus.googlepages.com/
% --------------------------------------------------------------------
% Yan Li, Ph.D,
% School of Control Science and Engineering,
% Shandong University, Jinan, Shandong 250061, P. R. China.
% --------------------------------------------------------------------
%
% --------------------------------------------------------------------
% YangQuan Chen, Ph.D, Associate Professor and Graduate Coordinator
% Department of Electrical and Computer Engineering,
% Director, Center for Self-Organizing and Intelligent Systems (CSOIS)
% Utah State University, 4120 Old Main Hill, Logan, UT 84322-4120, USA
% E: yqchen@ece.usu.edu or yqchen@ieee.org, T/F: 1(435)797-0148/3054;
% W: http://www.csois.usu.edu or http://yangquan.chen.googlepages.com
% --------------------------------------------------------------------
%
%--------------------------------------------------------------------
% Hu Sheng, Ph.D,
% Department of Electrical and Computer Engineering,
% Center for Self-Organizing and Intelligent Systems (CSOIS)
% Utah State University, 4120 Old Main Hill, Logan, UT 84322-4120, USA
% --------------------------------------------------------------------
% Note: All the codes were finished under the direction of Dr. YangQuan Chen.
%       The codes were writen by Dr. Yan Li and Hu Sheng,
%
% See also irid_fod.m
%          at http://www.mathworks.com/matlabcentral/files/21342/irid_fod.m
%          irid_doi.m
%          at
%          http://www.mathworks.com/matlabcentral/fileexchange/
%          26380-impulse-response-invariant-discretization-of-distributed-order-integrator

%*************************************************************************
function [sr]=irid_fsof(Ts,a,b,r,norder)

if nargin<5; norder=5; end
if a < 0 | b < 0 , sprintf('%s','a and b constant has to be positive'), return, end
if Ts < 0 , sprintf('%s','Sampling period has to be positive'), return, end
if r>=1 | r<= 0, sprintf('%s','The fractional order should be in (0,1)'), return, end
if norder<2, sprintf('%s','The order of the approximate transfer function has to be greater than 1'), return, end

close all;
wmax0=2*pi/Ts/2;          % rad./sec. Nyquist frequency
wmax=floor(1+ log10(wmax0) );
wmin=wmax-5;
w=logspace(wmin,wmax,1000);
j=sqrt(-1);
L=10/Ts;                  % decides the number of points of the impulse response function h(n)
t=[1:L]*Ts; y=[]; ht=[]; y1=[]; y2=[];

if a^2-4*b<0              % Case 1: a^2-4*b<0
    for k=1:length(t)
        y1(k)=quad(@(tau)realconvolution((-a)/2,sqrt(-a^2+4*b)/2,r,tau,t(k)),0,t(k));
        y2(k)=quad(@(tau)imconvolution((-a)/2,sqrt(-a^2+4*b)/2,r,tau,t(k)),0,t(k));
    end
    ht=y1+y2;
elseif a^2-4*b==0         % Case2: a^2-4*b==0
       ht = exp(-sqrt(b).*t).*t.^(2*r-1)/gamma(2*r);
else                      % Case 3: a^2-4*b>0
    for k=1:length(t)
       ht(k)=quadgk(@(x)integration2(x,a,b,r,t(k)),0,t(k));
    end
    s=(-a+sqrt(abs(a^2-4*b)))/2;
    ht= (exp(s.*t)/gamma(r)/gamma(r)).*ht;
end
h = [ht.*Ts];              % approcimation
q=norder;p=norder;
[B,A]=stmcb((h),q,p);
sprintf('Impulse response invariant Discrete approximated transfer function:')
sr=tf(B,A,Ts)
hht=impulse(sr,t);         % approcimated impulse response
                           % frequency response
srfr=(1./((j*w).^2 +a*j*w+b)).^(r);
srfr1=freqresp(sr,w);      % approcimated frequency response

figure;
subplot(3,1,1)             % comparision of impulse response
plot(t,ht,'b');
hold on;
plot(t,hht./Ts,'r-.')
axis([Ts,Ts.*L,-0.5,1]);
xlabel('Time');ylabel('Impulse response');
grid on;
legend(['impulse response of 1/(s^2 + ',num2str(a), '* s +',num2str(b),' )^{',num2str(abs(r)),'}'],'approximated impulse response');

subplot(3,1,2)            % comparision of magnitude
semilogx(w,20*log10(abs(srfr)),'b');hold on;
semilogx(w,20*log10(abs(reshape(srfr1, 1000, 1))),'r-.');
legend(['mag. Bode of 1/(s^2 + ',num2str(a), '* s +',num2str(b),' )^{',num2str(abs(r)),'}'],'approximated mag. Bode');
xlabel('Frequency (Hz)');ylabel('Magnitude (dB)');grid on;

subplot(3,1,3)            % comparision of phase
semilogx(w,(180/pi) * (angle(srfr)),'b');hold on;
semilogx(w,(180/pi) * (angle(reshape(srfr1, 1000, 1))),'r-.');
grid on
xlabel('Frequency (Hz)');ylabel('Phase (degrees)');
legend(['phase Bode of 1/(s^2 + ',num2str(a), '* s +',num2str(b),' )^{',num2str(abs(r)),'}'],'approximated phase Bode')

end

author: Sable

repo: mcbench-benchmarks


text:
function match(audio_file, gyro_rec, gyro_fs)
    % gyro_fs - Gyro sampling frequency

    % plot_spectrums(gyro_rec, GYRO_FS);

    [audio_samples, audio_fs] = audioread(audio_file);
    [gyro_timestamps, gyro_samples] = read_samples_file(gyro_rec);

    match_samples(audio_samples, audio_fs, gyro_samples, gyro_fs);
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function p29_title ( )

%*****************************************************************************80
%
%% P29_TITLE prints a title for problem 29.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Problem 29\n' );
  fprintf ( 1, '  Name:       Genz #3 / Patterson #8, Corner Peak\n' );
  fprintf ( 1, '  Region:     0 <= X(i) <= 1\n' );
  fprintf ( 1, '  Integrand:  F(X) = 1 / ( 1 + sum( C(i) * X(i) ) )^R\n' );
  fprintf ( 1, '  Parameters:\n' );
  fprintf ( 1, '              R, defaults to 0.3\n' );
  fprintf ( 1, '              C(1:DIM_NUM) defaults to 1/DIM_NUM.\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function SaveEvents(filename,events)

%SaveEvents - Write events to file.
%
%  USAGE
%
%    SaveEvents(filename,events)
%
%    filename            event file name
%    events              event data
%
%  SEE
%
%    See also NewEvents, LoadEvents, SaveRippleEvents.

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

if exist(filename), error('File already exists. Aborting.'); end

file = fopen(filename,'w');
if file == -1,
	error(['Cannot write to ' filename]);
end

for i = 1:length(events.time),
	fprintf(file,'%f\t%s\n',events.time(i)*1000,events.description{i}); % Convert to milliseconds
end

fclose(file);

%

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [hog_data, valid_data, vid_id] = Read_HOG_files_dynamic(users, hog_data_dir)

    hog_data = [];
    vid_id = {};
    valid_data = [];

    for i=1:numel(users)

        hog_files = dir([hog_data_dir, '/' users{i} '*.hog']);

        hog_data_curr_p = [];
        feats_filled = 0;

        for f_num=1:numel(hog_files)
            f = fopen([hog_data_dir,  '/', hog_files(f_num).name], 'r');

            curr_data = [];
            curr_ind = 0;

            while(~feof(f))

                if(curr_ind == 0)
                    num_cols = fread(f, 1, 'int32');
                    if(isempty(num_cols))
                        break;
                    end

                    num_rows = fread(f, 1, 'int32');
                    num_chan = fread(f, 1, 'int32');

                    curr_ind = curr_ind + 1;

                    % preallocate some space
                    if(curr_ind == 1)
                        curr_data = zeros(5000, 1 + num_rows * num_cols * num_chan);
                        num_feats =  1 + num_rows * num_cols * num_chan;
                    end

                    if(curr_ind > size(curr_data,1))
                        curr_data = cat(1, curr_data, zeros(6000, 1 + num_rows * num_cols * num_chan));
                    end
                    feature_vec = fread(f, [1, 1 + num_rows * num_cols * num_chan], 'float32');
                    curr_data(curr_ind, :) = feature_vec;
                else

                    % Reading in batches of 5000

                    feature_vec = fread(f, [4 + num_rows * num_cols * num_chan, 5000], 'float32');
                    feature_vec = feature_vec(4:end,:)';

                    num_rows_read = size(feature_vec,1);

                    curr_data(curr_ind+1:curr_ind+num_rows_read,:) = feature_vec;

                    curr_ind = curr_ind + size(feature_vec,1);

                end

            end

            fclose(f);

            curr_data = curr_data(1:curr_ind,:);
            vid_id_curr = cell(curr_ind,1);
            vid_id_curr(:) = users(i);

            vid_id = cat(1, vid_id, vid_id_curr);

            % Assume same number of frames per video
            if(f_num==1)
                hog_data_curr_p = zeros(curr_ind*numel(users), num_feats);
            end

            if(size(hog_data_curr_p,1) < feats_filled+curr_ind)
               hog_data_curr_p = cat(1, hog_data_curr_p, zeros(size(hog_data_curr_p,1), num_feats));
            end

            hog_data_curr_p(feats_filled+1:feats_filled+curr_ind,:) = curr_data;

            feats_filled = feats_filled + curr_ind;
        end

        hog_data_curr_p(1:feats_filled,2:end) = bsxfun(@plus, hog_data_curr_p(1:feats_filled,2:end), -median(hog_data_curr_p(1:feats_filled,2:end)));
        hog_data = cat(1, hog_data, hog_data_curr_p(1:feats_filled,:));

    end

    if(numel(users) > 0)
        valid_data = hog_data(:,1) > 0;
        hog_data = hog_data(:,2:end);
    end
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function a = bagging(alg,hyper)

% ============================================================================
% BAGGING bagging object
% ============================================================================
% A=bagging(A,H) returns a bagging object initialized with algorithm A and
%                hyperparameters H.
%
% Bagging trains BAGS classifiers on resamplings with replacement of the data
% of size M.  The final classifier is the averaged classifier.
% For classifiers that output a real-valued output in pattern recognition this
% gives you the choice of performing the average before or after taking the
% sign.
%
% You can also use the bagging object to average the results of already
% trained classifiers, simply define BAGGING(A) where A is a GROUP of
% classifiers
%
% Hyperparameters (with defaults)
%   child=svm            -- the algorithm to bag
%   bags=10              -- the number of times to bag
%   m=500                -- the number of points to sample for each bag
%                           (can also be a fraction of the training data)
%
% Example:
%    a1=svm; a1.C=10;
%    a=bagging(a1); a.bags=10; a.m=20;
%    [r a]=train(a,toy2d);
%    loss(test(a,toy2d))
%
% ============================================================================
% Reference : Bagging Predictors
% Author    : Leo Breiman
% Link      : http://citeseer.lcs.mit.edu/breiman96bagging.html
% ============================================================================

a.bags=10; a.m=500; a.child=svm;

p=algorithm('bagging');
a=class(a,'bagging',p);

if nargin>0
	a.child=alg;
end

if strcmp(a.child.algorithm.name,'group')
	%% already defined machines, can choose bags automatically
	a.bags=length(a.child.child);
end

if nargin==2
	eval_hyper;
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
data       = coder.load('GEN_initData.mat');
N          = data.N;
xDim       = data.dim.x;
uDim       = data.dim.u;

% global variable
global ParNMPCGlobalVariable
globalVariable = {'ParNMPCGlobalVariable',coder.Constant(ParNMPCGlobalVariable)};

cfg = coder.config('lib');
cfg.FilePartitionMethod = 'SingleFile';
cfg.EnableOpenMP = true;
cfg.TargetLang = 'C++';
stackUsageMax = (xDim+uDim)*N/360*200000;
cfg.StackUsageMax = stackUsageMax;
cfg.BuildConfiguration = 'Faster Runs';
cfg.EnableMemcpy  = true;
cfg.GenerateExampleMain = 'GenerateCodeOnly';
cfg.GenCodeOnly = true;
cfg.SupportNonFinite = false;
cfg.CustomInitializer = 'iiwa14_init();';      % init function
clear NMPC_Solve

myCCompiler = mex.getCompilerConfigurations(cfg.TargetLang,'Selected');
if ~strcmp(myCCompiler.Manufacturer,'Microsoft')
        cfg.PostCodeGenCommand = 'buildInfo.addLinkFlags(''-fopenmp'')';
end

codegen  -config cfg Simu_Matlab -globals globalVariable

author: deng-haoyang

repo: ParNMPC


text:
function sde_milstein_benchmark(tests,N,Tol)
%SDE_MILSTEIN_BENCHMARK  Performance tests for SDE_MILSTEIN solver function.
%   SDE_MILSTEIN_BENCHMARK will perform all tests. See code (type 'edit
%   SDE_MILSTEIN_BENCHMARK' in the command window) for test details or to add
%   additional tests. The SDETools Toolbox must be on the Matlab path or in the
%   same directory as this function.
%
%   SDE_MILSTEIN_BENCHMARK(TESTS) where TESTS = [N1 N2 ... NM] is a vector of
%   specific test numbers (indices, 1 ... N) to run. Tests can be listed in any
%   order and will be run in the order specified.
%
%   SDE_MILSTEIN_BENCHMARK(TESTS,N)
%   SDE_MILSTEIN_BENCHMARK(TESTS,N,TOL)
%
%   Not part of the the SDETools Toolbox; used only for development.
%
%   See also: SDE_MILSTEIN, SDEARGUMENTS, SDE_MILSTEIN_UNITTEST,
%       SDE_MILSTEIN_VALIDATE, SDE_EULER_BENCHMARK

%   Andrew D. Horchler, horchler @ gmail . com, Created 4-2-11
%   Revision: 1.2, 5-3-13

% Make sure toolbox on path, otherwise ensure we're in right location and add it
if ~isempty(strfind(path,'SDETools'))
    if exist('sde_milstein','file') ~= 2
        error('SDETools:sde_milstein_benchmark:FunctionNotFound',...
             ['The SDETools Toolbox is appears to be on the Matlab path, '...
              'but the SDE_MILSTEIN solver function cannot be found.']);
    end
else
    if exist('../SDETools','dir') ~= 7
        error('SDETools:sde_milstein_benchmark:ToolboxNotFound',...
             ['The SDETools Toolbox is not be on the Matlab path and the '...
              'root directory of the of the toolbox, SDETools, is not in '...
              'the same directory as the ''test'' directory.']);
    end
    addpath('../SDETools');
    Cleanup = onCleanup(@()rmpath('../SDETools'));  % Make sure path is reset
    if exist('sde_milstein','file') ~= 2
        error('SDETools:sde_milstein_benchmark:FunctionNotFoundAfterAdd',...
             ['The SDETools Toolbox was added to the Matlab path, but the '...
              'SDE_MILSTEIN solver function cannot be found.']);
    end
end

% Validate input argument if it exists
if nargin >= 1
    if isempty(tests) || ~isnumeric(tests) || ~all(isfinite(tests))
        error('SDETools:sde_milstein_benchmark:InvalidArg1',...
              'Invalid argument 1.');
    end
    if any(tests < 1)
        error('SDETools:sde_milstein_benchmark:NotAnIndex',...
              'Tests are numbered as indices, from 1 to N.');
    end
    tests = floor(tests);
    RunTests = true;
else
    RunTests = false;
end

if nargin >= 2
    if isempty(N) || ~isnumeric(N) || ~all(isfinite(N))
        error('SDETools:sde_milstein_benchmark:InvalidArg2',...
              'Invalid argument 2.');
    end
    N = floor(N);
else
    N = 1e3;
end

if nargin == 3
    if isempty(Tol) || length(Tol) ~= 1 || ~isnumeric(Tol) || ...
            ~all(isfinite(Tol)) || Tol <= 0
        error('SDETools:sde_milstein_benchmark:InvalidArg2',...
              'Invalid argument 2.');
    end
    Tol = floor(Tol);
else
    Tol = 1e-3;
end

M = 16;
fi = cell(M,1);
gi = cell(M,1);
tspani = cell(M,1);
y0i = cell(M,1);
optsi = cell(M,1);

% no FOR loop cases

% scalar ICs

i = 1;
fi{i} = 1;
gi{i} = 1;
tspani{i} = 0:0.001:1;
y0i{i} = 0;
optsi{i} = [];
oneout{i} = false;
description{i} = 'Scalar, constant drift and diffusion, variable step-size';

i = i+1;
fi{i} = 1;
gi{i} = 1;
tspani{i} = 0:0.5:500;
y0i{i} = 0;
optsi{i} = [];
oneout{i} = false;
description{i} = 'Scalar, constant drift and diffusion, fixed step-size';

i = i+1;
fi{i} = @(t,x)x+1;
gi{i} = @(t,x)x+1;
tspani{i} = 0:0.001:1;
y0i{i} = 0;
optsi{i} = sdeset('ConstFFUN','yes','ConstGFUN','yes');
oneout{i} = false;
description{i} = 'Scalar, constant drift and diffusion functions, variable step-size';

i = i+1;
fi{i} = @(t,x)x+1;
gi{i} = @(t,x)x+1;
tspani{i} = 0:0.5:500;
y0i{i} = 0;
optsi{i} = sdeset('ConstFFUN','yes','ConstGFUN','yes');
oneout{i} = false;
description{i} = 'Scalar, constant drift and diffusion functions, fixed step-size';

i = i+1;
fi{i} = 1;
gi{i} = 1;
tspani{i} = 0:0.001:1;
y0i{i} = 0;
optsi{i} = sdeset('Diagonal','no');
oneout{i} = false;
description{i} = 'Scalar, constant drift, constant non-diagonal diffusion, variable step-size';

i = i+1;
fi{i} = 1;
gi{i} = 1;
tspani{i} = 0:0.5:500;
y0i{i} = 0;
optsi{i} = sdeset('Diagonal','no');
oneout{i} = false;
description{i} = 'Scalar, constant drift, constant non-diagonal diffusion, fixed step-size';

i = i+1;
fi{i} = @(t,x)x+1;
gi{i} = @(t,x)x+1;
tspani{i} = 0:0.001:1;
y0i{i} = 0;
optsi{i} = sdeset('Diagonal','no','ConstFFUN','yes','ConstGFUN','yes');
oneout{i} = false;
description{i} = 'Scalar, constant drift function, constant non-diagonal diffusion function, variable step-size';

i = i+1;
fi{i} = @(t,x)x+1;
gi{i} = @(t,x)x+1;
tspani{i} = 0:0.5:500;
y0i{i} = 0;
optsi{i} = sdeset('Diagonal','no','ConstFFUN','yes','ConstGFUN','yes');
oneout{i} = false;
description{i} = 'Scalar, constant drift function, constant non-diagonal diffusion function, fixed step-size';

% vector ICs

i = i+1;
fi{i} = ones(1000,1);
gi{i} = ones(1000,1);
tspani{i} = 0:0.001:1;
y0i{i}(1,1000) = 0;
optsi{i} = [];
oneout{i} = false;
description{i} = 'Vector, constant drift and diffusion, variable step-size';

i = i+1;
fi{i} = ones(1000,1);
gi{i} = ones(1000,1);
tspani{i} = 0:0.5:500;
y0i{i}(1,1000) = 0;
optsi{i} = [];
oneout{i} = false;
description{i} = 'Vector, constant drift and diffusion, fixed step-size';

i = i+1;
fi{i} = @(t,x)x+1;
gi{i} = @(t,x)x+1;
tspani{i} = 0:0.001:1;
y0i{i} = zeros(1,1000);
optsi{i} = sdeset('ConstFFUN','yes','ConstGFUN','yes');
oneout{i} = false;
description{i} = 'Vector, constant drift and diffusion functions, variable step-size';

i = i+1;
fi{i} = @(t,x)x+1;
gi{i} = @(t,x)x+1;
tspani{i} = 0:0.5:500;
y0i{i}(1,1000) = 0;
optsi{i} = sdeset('ConstFFUN','yes','ConstGFUN','yes');
oneout{i} = false;
description{i} = 'Vector, constant drift and diffusion functions, fixed step-size';

i = i+1;
fi{i} = ones(1000,1);
gi{i} = 1;
tspani{i} = 0:0.001:1;
y0i{i}(1,1000) = 0;
optsi{i} = sdeset('Diagonal','no');
oneout{i} = false;
description{i} = 'Vector, constant drift, constant non-diagonal diffusion, variable step-size';

i = i+1;
fi{i} = ones(1000,1);
gi{i} = 1;
tspani{i} = 0:0.5:500;
y0i{i}(1,1000) = 0;
optsi{i} = sdeset('Diagonal','no');
oneout{i} = false;
description{i} = 'Vector, constant drift, constant non-diagonal diffusion, fixed step-size';

i = i+1;
fi{i} = @(t,x)x+1;
gi{i} = @(t,x)x(1)+1;
tspani{i} = 0:0.001:1;
y0i{i}(1,1000) = 0;
optsi{i} = sdeset('Diagonal','no','ConstFFUN','yes','ConstGFUN','yes');
paramsi{i} = [];
oneout{i} = false;
description{i} = 'Vector, constant drift function, constant non-diagonal diffusion function, variable step-size';

i = i+1;
fi{i} = @(t,x)x+1;
gi{i} = @(t,x)x(1)+1;
tspani{i} = 0:0.5:500;
y0i{i}(1,1000) = 0;
optsi{i} = sdeset('Diagonal','no','ConstFFUN','yes','ConstGFUN','yes');
oneout{i} = false;
description{i} = 'Vector, constant drift function, constant non-diagonal diffusion function, fixed step-size';

i = i+1;
fi{i} = @(t,x)x+1;
gi{i} = @(t,x)x+1;
tspani{i} = 0:0.001:1;
y0i{i} = 0;
optsi{i} = [];
oneout{i} = false;
description{i} = 'Scalar, diagonal diffusion function, variable step-size';

i = i+1;
fi{i} = @(t,x)x+1;
gi{i} = @(t,x)x+1;
tspani{i} = 0:0.5:500;
y0i{i} = 0;
optsi{i} = [];
oneout{i} = false;
description{i} = 'Scalar, diagonal diffusion function, fixed step-size';

i = i+1;
fi{i} = @(t,x)x+1;
gi{i} = @(t,x)x+1;
tspani{i} = 0:0.001:1;
y0i{i} = 0;
optsi{i} = sdeset('SDEType','Ito');
oneout{i} = false;
description{i} = 'Scalar, Ito type, diagonal diffusion function, variable step-size';

i = i+1;
fi{i} = @(t,x)x+1;
gi{i} = @(t,x)x+1;
tspani{i} = 0:0.5:500;
y0i{i} = 0;
optsi{i} = sdeset('SDEType','Ito');
paramsi{i} = [];
oneout{i} = false;
description{i} = 'Scalar, Ito type, diagonal diffusion function, fixed step-size';

% Memory allocation
i = i+1;
fi{i} = ones(500000,1);
gi{i} = ones(500000,1);
tspani{i} = 0:0.5:5;
y0i{i}(1,500000) = 0;
optsi{i} = [];
oneout{i} = false;
description{i} = 'Memory allocation: vector, constant drift and diffusion, fixed step-size';

m = length(fi);
if RunTests
    ts = tests(tests<=m);
else
    ts = 1:m;
end

Stream = RandStream('mt19937ar','Seed',0);

sp = ' ';
t = zeros(1,N);
lts = length(ts);
mt(1,lts) = 0;
st(1,lts) = 0;
ost = ceil(log10((lts+1)));
for k=1:lts
    i = ts(k);

    % Get arguments
    f = fi{i};
    g = gi{i};
    tspan = tspani{i};
    y0 = y0i{i};
    opts = optsi{i};

    % Ensure that all tests start out with same random variates
    try
        RandStream.setGlobalStream(Stream);
    catch                                   	%#ok<CTCH>
        RandStream.setDefaultStream(Stream);	%#ok<SETRS>
    end

    % Warm up function before timing
    if oneout{i}
        if isempty(opts)
            y = sde_milstein(f,g,tspan,y0);     %#ok<*NASGU>
        else
            y = sde_milstein(f,g,tspan,y0,opts);
        end
    else
        if isempty(opts)
            [y,w] = sde_milstein(f,g,tspan,y0);	%#ok<*ASGLU>
        else
            [y,w] = sde_milstein(f,g,tspan,y0,opts);
        end
    end

    j = 0;
    while true  % Timing loop
        j = j+1;

        % Time an iteration
        if oneout{i}
            if isempty(opts)
                tic
                y = sde_milstein(f,g,tspan,y0);
                t(j) = toc;
            else
                tic
                y = sde_milstein(f,g,tspan,y0,opts);
                t(j) = toc;
            end
        else
            if isempty(opts)
                tic
                [y,w] = sde_milstein(f,g,tspan,y0);
                t(j) = toc;
            else
                tic
                [y,w] = sde_milstein(f,g,tspan,y0,opts);
                t(j) = toc;
            end
        end

        % Check if times are sufficient
        if N < 20 && N >= 10
            if (j >= 10 && sum(t(1:j)) > 0.1 && std(t(1:j)) < Tol) || j > N
                break
            end
        elseif N < 10
            if j > N
                break
            end
        else
            if (j >= 10 && sum(t(1:j)) > 0.1 && std(t(1:j)) < Tol) || ...
                    (j >= 20 && sum(t(1:j)) > 2) || j > N || sum(t(1:j)) > 5
                break
            end
        end
    end
    mt(k) = mean(t(1:j));
	st(k) = std(t(1:j));
    fprintf(1,[description{k} '\n']);
    fprintf(1,['Test %s%d took an average of %1.3e sec. (%1.3e '...
               'sec./time-step) over %d runs.\n\n'],...
              sp(ones(1,ost-ceil(log10((i+1))))),i,mt(k),mt(k)/length(tspan),j);
end
fprintf(1,'Total mean time for all %d tests: %4.4f +/- %4.4f sec.\n',...
	lts,sum(mt),mean(st));

author: horchler

repo: SDETools

save_path: github-repos/MATLAB/horchler-SDETools


text:
% Wrapper for vl_nnsplit block
% inputs{1} :   X     : 1 x 1 x 16 x b
% outputs{1}:   y     : 1 x 1 x 10 x b
% outputs{1}:   r     : 1 x 1 x 3 x b
% outputs{1}:   t     : 1 x 1 x 2 x b
% outputs{1}:   s     : 1 x 1 x 1 x b

classdef split < dagnn.Layer

    methods
        function outputs = forward(~, inputs, ~)

            useGPU = isa(inputs{1}, 'gpuArray');
            if useGPU
            [y,r,t,s] = vl_nnsplit(gather(inputs{1}));
            outputs = {gpuArray(y),gpuArray(r),gpuArray(t),gpuArray(s)};
            else
            [y,r,t,s] = vl_nnsplit(inputs{1});
            outputs = {y,r,t,s};
            end

        end

        function [derInputs, derParams] = backward(~, inputs, ~, derOutputs)

            useGPU = isa(inputs{1}, 'gpuArray');
            if useGPU
            derInputs{1} = gpuArray( vl_nnsplit(gather(inputs{1}),gather(derOutputs{1}), gather(derOutputs{2}), gather(derOutputs{3}), gather(derOutputs{4})) );
            else
            derInputs{1} = vl_nnsplit(inputs{1},derOutputs{1},derOutputs{2},derOutputs{3},derOutputs{4});
            end

            derParams = {};
        end

        function outputSizes = getOutputSizes(~, inputSizes)
            outputSizes = {1 1 inputSizes{1}(3) inputSizes{1}(4)} ;
        end

        function obj = split(varargin)
            obj.load(varargin);
        end
    end
end

author: anilbas


text:
function [res] = mne_ex_evoked_grad_amp(inname,bmin,bmax,outname)
%
%   function [res] = mne_ex_evoked_grad_amp(inname,bmin,bmax,outname)
%
%   Compute the magnitude of the tangential gradient at each
%   sensor location using the planar gradiometer data and
%   optionally output the result to a fif file.
%
%   inname      The input file name. All average data sets are
%               read and processed
%   bmin,bmax   Baseline limits in seconds
%   outname     Optional output file name
%
%
%   Function returns the data which was or would have been written
%   to the file
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%

me='MNE:mne_ex_evoked_grad_amp';

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

if nargin == 1
    do_baseline = false;
elseif (nargin == 3 || nargin == 4)
    do_baseline = true;
else
    error(me,'Wrong number of arguments');
end
%
%   Read the data
%
data = fiff_read_evoked_all(inname);
%
%   Figure out the planar gradiometer pairs
%
pairs = zeros(data.info.nchan,2);
npair = 0;
k = 1;
while k < data.info.nchan
    %
    %   First check the coil types
    %
    coil1 = data.info.chs(k).coil_type;
    coil2 = data.info.chs(k+1).coil_type;
    if (coil1 == coil2 && ...
            (coil1 == 2 || coil1 == 3012 || coil1 == 3013))
        one = data.info.ch_names{k};
        two = data.info.ch_names{k+1};
        lastone = one(length(one));
        lasttwo = two(length(two));
        %
        %   Then the channel names
        %
        if (strcmp(one(1:3),'MEG') && strcmp(one(1:3),'MEG'))
            if (strcmp(one(1:(length(one)-1)),two(1:(length(two)-1))) && ...
                    ((lastone == '2' && lasttwo == '3') || ...
                    (lastone == '3' && lasttwo == '2')))
                npair = npair + 1;
                pairs(npair,1) = k;
                pairs(npair,2) = k+1;
                k = k + 1;
            end
        end
    end
    k = k + 1;
end

if npair == 0
    error(me,'No planar gradiometers in these data');
end
%
%   Compute the amplitudes
%
fprintf(1,'Computing the amplitudes');
if do_baseline
    fprintf(1,' (Baseline = %7.1f ... %7.1f ms)',1000*bmin,1000*bmax);
end
fprintf(1,'...');
for k = 1:length(data.evoked)
    epochs = data.evoked(k).epochs;
    %
    %  Setup baseline limits
    %
    if do_baseline
        b1 = double(data.info.sfreq*bmin - data.evoked(k).first);
        b2 = double(data.info.sfreq*bmax - data.evoked(k).first);
        if b1 < 1
            b1 = 1;
        end
        if b2 > size(epochs,2)
            b2 = size(epochs,2)
        end
    else
        b1 = 1;
        b2 = 1;
    end
    %
    %   Go through all pairs
    %
    for p = 1:npair
        one = pairs(p,1);
        two = pairs(p,2);
        if b2 > b1
            base1 = sum(epochs(one,b1:b2))/(b2-b1);
            base2 = sum(epochs(two,b1:b2))/(b2-b1);
            epochs(one,:) = sqrt((epochs(one,:)-base1).*(epochs(one, ...
                :)-base1)+(epochs(two,:)-base2).*(epochs(two,:)-base2));
        else
            epochs(one,:) = sqrt(epochs(one,:).*epochs(one, ...
                :)+epochs(two,:).*epochs(two,:));
        end
    end
    data.evoked(k).epochs = epochs;
    fprintf(1,'.');
end
fprintf(1,'[done]\n');
%
%   Compose the selection name list
%
pairs = pairs(1:npair,:);
for k = 1:npair
    ch_sel_names{k} = data.info.ch_names{pairs(k,1)};
end
%
%   Omit MEG channels but include others
%
for p = 1:data.info.nchan
    if (data.info.chs(p).kind ~= FIFF.FIFFV_MEG_CH)
        k = k + 1;
        ch_sel_names{k} = data.info.ch_names{p};
    end
end
%
%   Modify the bad channel list
%
if ~isempty(data.info.bads)
    nbad = length(data.info.bads);
    for k = 1:npair
        one = data.info.ch_names{pairs(k,1)};
        two = data.info.ch_names{pairs(k,2)};
        %
        %   If one channel of the planar gradiometer is marked bad,
        %   add the other to the bad channel list
        %
        if (~isempty(strmatch(one,data.info.bads)) && ...
                isempty(strmatch(two,data.info.bads)))
            nbad = nbad + 1;
            data.info.bads{nbad} = two;
        elseif (isempty(strmatch(one,data.info.bads)) && ...
                ~isempty(strmatch(two,data.info.bads)))
            nbad = nbad + 1;
            data.info.bads{nbad} = one;
        end
    end
end
%
%   Do the picking
%
res = fiff_pick_channels_evoked(data,ch_sel_names);
%
%   Optionally write an output file
%
if nargin == 4
    fiff_write_evoked(outname,res);
    fprintf(1,'Wrote %s\n',outname);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
filters = haar_filter_bank_2d_spatial;

% Initialize signal and meta
x = uiuc_sample;
U{1}.signal{1} = x;
U{1}.meta.j = zeros(0,1);
options.J = 5;

[A2, W2] = wavelet_layer_2d_pyramid(U{1}, filters, options)
U2 = modulus_layer(W2);

%%
clear
x = uiuc_sample;
filt_opt.type = 'haar';
scat_opt.J = 5;
scat_opt.M = 3;
Wop = wavelet_factory_2d_pyramid(filt_opt, scat_opt);
Sx = scat(x, Wop);

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [ftAllNew,transMdl] = ftTrans_itl(ft,maSrc,target,maLabeled,param)
%Wrapper of the Information-Theoretical Learning (ITL) algorithm
%
%	A feaure-level transfer learning (domain adaptation) algorithm which
% learns a domain-invariant subspace. Application scope:
%	+ two discrete domains (source and target)
%	+ labeled source domain
%	+ unlabeled target domain
%	+ label type: classification
%
% ftAll:	All samples in all domains. n-by-m matrix, n is the number of
%	samples, m is the dimension of features.
% maSrc:	n-by-1 logical vector, maSrc(i)=true if i is from the source
%	domain, 0 if target domain.
% target:	The class labels of the source domain, nsrc-by-1 matrix, nsrc
%	is the number of source	samples. target(i)=j if the i'th sample is from
%	the j'th class.
% maLabeled:	Mask for the labeled samples. n-by-1 matrix, must be the
% same with maSrc.
%
% param: Struct of hyper-parameters, please see the first cell of this
%	program ("default parameters") for details. You can set parameter p to
%	x by setting param.p = x. For parameters that are not set, default
%	values will be used.
%
% ftAllNew:	All samples in the learned subspace.
% transMdl:	A struct containing the model, transMdl.W is the projection
%	matrix.
%
%	The ftProc_pca_tr function is a part of the PRTools toolbox.
%	Notice from the author of the paper: for practice use, you need to
% properly preprocess the data and tune lambda and d.
%
% ref: Y. Shi and F. Sha, "Information-theoretical learning of discriminative
% clusters for unsupervised domain adaptation," in ICML, 2012.
%
% Copyright 2016 Ke YAN, Tsinghua Univ. http://yanke23.com , xjed09@gmail.com

addpath infometric_0.1

%% default parameters
pcaCoef = 0; % see ftProc_pca_tr
lambda = 1; % regularization parameter

defParam

%% sort samples
if any(maSrc~=maLabeled)
	error('maLabeled must be the same with maSrc')
end

ftSrc = ft(maSrc,:);
ftTar = ft(~maSrc,:);

%% compute
[~,pcaModelS] = ftProc_pca_tr(ftSrc,[],struct('pcaCoef',pcaCoef));
[~,pcaModelT] = ftProc_pca_tr(ftTar,[],struct('pcaCoef',pcaCoef));
d = min(size(pcaModelS.W_prj,2),size(pcaModelT.W_prj,2));
W_prjT = pcaModelT.W_prj(:,1:d);

L = infometric(W_prjT, ftSrc, target, ftTar, lambda);

%% project the samples
ftAllNew = zeros(size(ft,1),d);
ftAllNew(maSrc,:) = ftSrc*L;
ftAllNew(~maSrc,:) = ftTar*L;
transMdl.L = L;

end

author: viggin

repo: domain-adaptation-toolbox


text:
function coeff = calc_energy_coeff(dim,N,s,energy)
%CALC_ENERGY_COEFF Coefficient of second term in expansion of energy
%
%Syntax
% coeff = calc_energy_coeff(d,N,s,energy);
%
%Description
% COEFF = CALC_ENERGY_COEFF(dim,N,s,ENERGY) sets COEFF to be the coefficient of
% the second term of an expansion of ENERGY with the same form as the expansion
% of E(dim,N,s), the minimum r^(-s) energy of a set of N points on S^dim.
%
% Specifically, for s not equal to 0, COEFF is the solution to
%
% ENERGY == (SPHERE_INT_ENERGY(dim,s)/2) N^2 + COEFF N^(1+s/dim),
%
% and for s == 0 (the logarithmic potential), COEFF is the solution to
%
% ENERGY == (SPHERE_INT_ENERGY(dim,0)/2) N^2 + COEFF N LOG(N).
%
% The argument dim must be a positive integer.
% The argument N must be a positive integer or an array of positive integers.
% The argument ENERGY must an array of real numbers of the same array size as N.
% The result COEFF will be an array of the same size as N.
%
%Notes
% 1) The energy expansion is not valid for N == 1, and in particular,
%
% EQ_ENERGY_COEFF(dim,N,0,energy) := 0.
%
% 2) For s > 0, [KuiS98 (1.6) p524] has
%
% E(dim,N,s) == (SPHERE_INT_ENERGY(dim,s)/2) N^2 + COEFF N^(1+s/dim) + ...
%
% where SPHERE_INT_ENERGY(dim,s) is the energy integral of the r^(-s) potential
% on S^dim.
%
% The case s == 0 (logarithmic potential) can be split into subcases.
% For s == 0 and dim == 1, E(1,N,0) is obtained by equally spaced points on S^1,
% and the formula for the log potential for N equally spaced points on a circle
% gives
%
% E(1,N,0) == (-1/2) N LOG(N) exactly.
%
% For s == 0 and dim == 2, [SafK97 (4) p7] has
%
% E(2,N,0) == (SPHERE_INT_ENERGY(2,0)/2) N^2 + COEFF N LOG(N) + o(N LOG(N)).
%
% In general, for s == 0,
%
% E(dim,N,0) == (SPHERE_INT_ENERGY(dim,0)/2) N^2 + COEFF N LOG(N) + ...
%
% with sphere_int_energy(1,0) == 0.
%
% CALC_ENERGY_COEFF just uses this general formula for s == 0, so for s == 0 and
% dim == 1, the coefficient returned is actually the coefficient of the first
% non-zero term.
%
%Examples
% > N=2:6
%  N =
%       2     3     4     5     6
%
% > energy=eq_energy_dist(2,N,0)
%  energy =
%     -0.6931   -1.3863   -2.7726   -4.4205   -6.2383
%
% > calc_energy_coeff(2,N,0,energy)
%  ans =
%     -0.2213   -0.1569   -0.2213   -0.2493   -0.2569
%
%See also
% EQ_ENERGY_DIST, EQ_ENERGY_COEFF

% Copyright 2004-2005 Paul Leopardi for the University of New South Wales.
% $Revision 1.10 $ $Date 2005-06-01 $
% Documentation files renamed
% $Revision 1.00 $ $Date 2005-02-12 $
%
% For licensing, see COPYING.
% For references, see AUTHORS.
% For revision history, see CHANGELOG.

%
% Compute the energy coefficient: subtract the first term in the expansion of
% the minimum energy and divide by the power of N in the second term.
%
if s > 0
    %
    % The first term in the expansion of the minimum energy.
    % Ref: [KuiS98 (1.6) p524]
    %
    first_term = (sphere_int_energy(dim,s)/2) * N.^2;
    coeff = (energy-first_term) ./ (N.^(1+s/dim));
else
    %
    % Flatten N into a row vector.
    %
    shape = size(N);
    n_partitions = prod(shape);
    N = reshape(N,1,n_partitions);
    %
    % Refs for s==0, dim == 2:
    % [SafK97 (4) p. 7] [Zho95 (5.6) p. 68, 3.11 - corrected) p. 42]
    %
    first_term = (sphere_int_energy(dim,s)/2) * N.^2;
    %
    % Avoid division by zero.
    %
    coeff = zeros(size(N));
    neq1 = (N ~= 1);
    coeff(neq1) = (energy(neq1)-first_term(neq1)) ./ (N(neq1) .* log(N(neq1)));
    %
    % Reshape output to same array size as original N.
    %
    coeff = reshape(coeff,shape);
    %
end
%
% end function

function energy = sphere_int_energy(dim,s)
%SPHERE_INT_ENERGY Energy integral of r^(-s) potential
%
%Syntax
% energy = sphere_int_energy(d,s);
%
%Description
% ENERGY = SPHERE_INT_ENERGY(dim,s) sets ENERGY to be the energy integral
% on S^dim of the r^(-s) potential, defined using normalized Lebesgue measure.
%
% Ref for s > 0: [KuiS98 (1.6) p524]
% Ref for s == 0 and dim == 2: SafK97 (4) p. 7]
% For s == 0 and dim >= 2, integral was obtained using Maple:
%     energy = (1/2)*(omega(dim)/omega(dim+1)* ...
%          int(-log(2*sin(theta/2)*(sin(theta))^(dim-1),theta=0..Pi),
%     where omega(dim+1) == area_of_sphere(dim).

if s ~= 0
    energy = (gamma((dim+1)/2)*gamma(dim-s)/(gamma((dim-s+1)/2)*gamma(dim-s/2)));
elseif dim ~= 1
    energy = (psi(dim)-psi(dim/2)-log(4))/2;
else
    energy = 0;
end
%
% end function

repo: gspbox


text:
function [trl, event] = ft_trialfun_edf(cfg)

% FT_TRIALFUN_EDF is an example trial function for EDF data. It searches for events
% of type "up" in an analog data channel, as indentified by thresholding. This
% threshold can be a hard threshold, i.e. a numeric, or can flexibly be defined
% depending on the data, for example calculating the 'median' of an analog signal.
%
% You can use this as a template for your own conditial trial definitions.
%
% Use this function by calling
%   [cfg] = ft_definetrial(cfg)
% where the configuration structure should contain
%   cfg.dataset  = string with the filename
%   cfg.trialfun = 'ft_trialfun_edf'
%
% See also FT_DEFINETRIAL, FT_TRIALFUN_GENERAL

ft_warning('this trial function is only an example, please copy it and adapt it to your specific EDF situation - see http://www.fieldtriptoolbox.org/getting_started/edf');

% read the header information
hdr           = ft_read_header(cfg.dataset);

% read the events from the data
chanindx      = 1; % this should be adapted to your data
detectflank   = 'up';
threshold     = '(3/2)*nanmedian'; % or, e.g., 1/2 times the median for down flanks
event         = ft_read_event(cfg.dataset, 'trigindx', trigindx, 'detectflank', detectflank, 'threshold', threshold);

% define trials around the events
trl           = [];
pretrig       = 1 * hdr.Fs; % e.g., 1 sec before trigger
posttrig      = 2 * hdr.Fs; % e.g., 2 sec after trigger
for i = 1:numel(event)
  offset    = -hdr.nSamplesPre;  % number of samples prior to the trigger
  trlbegin  = event(i).sample - pretrig;
  trlend    = event(i).sample + posttrig;
  newtrl    = [trlbegin trlend offset];
  trl       = [trl; newtrl]; % store in the trl matrix
end


text:
function frame = getframe(vr)
%FRAME=GETFRAME(VR)
%  Returns the current frame that has been read in from video VR.  Calling
%  GETFRAME does not automatically advance the current frame so the user
%  must first call NEXT, STEP, or SEEK before calling GETFRAME the
%  first time.  This makes it easy to write loops such as:
%    vr = videoReader(...);
%    while (next(vr))
%      img = getframe(vr);
%      ...do something...
%    end
%    vr = close(vr);
%
%SEE ALSO
%  videoReader
%
%Copyright (c) 2006 Gerald Dalley
%See "MIT.txt" in the installation directory for licensing details (especially
%when using this library on GNU/Linux).

frame = feval(vr.plugin, 'getframe', vr.handle);


text:
function [Csig,Cbak,Covl]= composite(cleanFile, enhancedFile);
% ----------------------------------------------------------------------
%          Composite Objective Speech Quality Measure
%
%   This function implements the composite objective measure proposed in
%   [1].
%
%   Usage:  [sig,bak,ovl]=composite(cleanFile.wav, enhancedFile.wav)
%
%         cleanFile.wav - clean input file in .wav format
%         enhancedFile  - enhanced output file in .wav format
%         sig           - predicted rating [1-5] of speech distortion
%         bak           - predicted rating [1-5] of noise distortion
%         ovl           - predicted rating [1-5] of overall quality
%
%       In addition to the above ratings (sig, bak, & ovl) it returns
%       the individual values of the LLR, SNRseg, WSS and PESQ measures.
%
%  Example call:  [sig,bak,ovl] =composite('sp04.wav','enhanced.wav')
%
%
%  References:
%
%     [1]   Hu, Y. and Loizou, P. (2006). Evaluation of objective measures
%           for speech enhancement. Proc. Interspeech, Pittsburg, PA.
%
%   Authors: Yi Hu and Philipos C. Loizou
%   (the LLR, SNRseg and WSS measures were based on Bryan Pellom and John
%     Hansen's implementations)
%
% Copyright (c) 2006 by Philipos C. Loizou
% $Revision: 0.0 $  $Date: 10/09/2006 $

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

if nargin~=2
    fprintf('USAGE: [sig,bak,ovl]=composite(cleanFile.wav, enhancedFile.wav)\n');
    fprintf('For more help, type: help composite\n\n');
    return;
end

alpha= 0.95;

[data1, Srate1, Nbits1]= wavread(cleanFile);
[data2, Srate2, Nbits2]= wavread(enhancedFile);
if ( Srate1~= Srate2) | ( Nbits1~= Nbits2)
    error( 'The two files do not match!\n');
end

len= min( length( data1), length( data2));
data1= data1( 1: len)+eps;
data2= data2( 1: len)+eps;

% -- compute the WSS measure ---
%
wss_dist_vec= wss( data1, data2,Srate1);
wss_dist_vec= sort( wss_dist_vec);
wss_dist= mean( wss_dist_vec( 1: round( length( wss_dist_vec)*alpha)));

% --- compute the LLR measure ---------
%
LLR_dist= llr( data1, data2,Srate1);
LLRs= sort(LLR_dist);
LLR_len= round( length(LLR_dist)* alpha);
llr_mean= mean( LLRs( 1: LLR_len));

% --- compute the SNRseg ----------------
%
[snr_dist, segsnr_dist]= snr( data1, data2,Srate1);
snr_mean= snr_dist;
segSNR= mean( segsnr_dist);

% -- compute the pesq ----
%
% if     Srate1==8000,    mode='nb';
% elseif Srate1 == 16000, mode='wb';
% else,
%      error ('Sampling freq in PESQ needs to be 8 kHz or 16 kHz');
% end

 [pesq_mos_scores]= pesq(cleanFile, enhancedFile);

 if length(pesq_mos_scores)==2
     pesq_mos=pesq_mos_scores(1); % take the raw PESQ value instead of the
                                  % MOS-mapped value (this composite
                                  % measure was only validated with the raw
                                  % PESQ value)
 else
     pesq_mos=pesq_mos_scores;
 end

% --- now compute the composite measures ------------------
%
Csig = 3.093 - 1.029*llr_mean + 0.603*pesq_mos-0.009*wss_dist;
  Csig = max(1,Csig);  Csig=min(5, Csig); % limit values to [1, 5]
Cbak = 1.634 + 0.478 *pesq_mos - 0.007*wss_dist + 0.063*segSNR;
  Cbak = max(1, Cbak); Cbak=min(5,Cbak); % limit values to [1, 5]
Covl = 1.594 + 0.805*pesq_mos - 0.512*llr_mean - 0.007*wss_dist;
  Covl = max(1, Covl); Covl=min(5, Covl); % limit values to [1, 5]

fprintf('\n LLR=%f   SNRseg=%f   WSS=%f   PESQ=%f\n',llr_mean,segSNR,wss_dist,pesq_mos);

return; %=================================================================

function distortion = wss(clean_speech, processed_speech,sample_rate)

% ----------------------------------------------------------------------
% Check the length of the clean and processed speech.  Must be the same.
% ----------------------------------------------------------------------

clean_length      = length(clean_speech);
processed_length  = length(processed_speech);

if (clean_length ~= processed_length)
  disp('Error: Files  musthave same length.');
  return
end

% ----------------------------------------------------------------------
% Global Variables
% ----------------------------------------------------------------------

winlength   = round(30*sample_rate/1000); %240;		   % window length in samples
skiprate    = floor(winlength/4);		   % window skip in samples
max_freq    = sample_rate/2;	   % maximum bandwidth
num_crit    = 25;		   % number of critical bands

USE_FFT_SPECTRUM = 1;		   % defaults to 10th order LP spectrum
n_fft       = 2^nextpow2(2*winlength);
n_fftby2    = n_fft/2;		   % FFT size/2
Kmax        = 20;		   % value suggested by Klatt, pg 1280
Klocmax     = 1;		   % value suggested by Klatt, pg 1280

% ----------------------------------------------------------------------
% Critical Band Filter Definitions (Center Frequency and Bandwidths in Hz)
% ----------------------------------------------------------------------

cent_freq(1)  = 50.0000;   bandwidth(1)  = 70.0000;
cent_freq(2)  = 120.000;   bandwidth(2)  = 70.0000;
cent_freq(3)  = 190.000;   bandwidth(3)  = 70.0000;
cent_freq(4)  = 260.000;   bandwidth(4)  = 70.0000;
cent_freq(5)  = 330.000;   bandwidth(5)  = 70.0000;
cent_freq(6)  = 400.000;   bandwidth(6)  = 70.0000;
cent_freq(7)  = 470.000;   bandwidth(7)  = 70.0000;
cent_freq(8)  = 540.000;   bandwidth(8)  = 77.3724;
cent_freq(9)  = 617.372;   bandwidth(9)  = 86.0056;
cent_freq(10) = 703.378;   bandwidth(10) = 95.3398;
cent_freq(11) = 798.717;   bandwidth(11) = 105.411;
cent_freq(12) = 904.128;   bandwidth(12) = 116.256;
cent_freq(13) = 1020.38;   bandwidth(13) = 127.914;
cent_freq(14) = 1148.30;   bandwidth(14) = 140.423;
cent_freq(15) = 1288.72;   bandwidth(15) = 153.823;
cent_freq(16) = 1442.54;   bandwidth(16) = 168.154;
cent_freq(17) = 1610.70;   bandwidth(17) = 183.457;
cent_freq(18) = 1794.16;   bandwidth(18) = 199.776;
cent_freq(19) = 1993.93;   bandwidth(19) = 217.153;
cent_freq(20) = 2211.08;   bandwidth(20) = 235.631;
cent_freq(21) = 2446.71;   bandwidth(21) = 255.255;
cent_freq(22) = 2701.97;   bandwidth(22) = 276.072;
cent_freq(23) = 2978.04;   bandwidth(23) = 298.126;
cent_freq(24) = 3276.17;   bandwidth(24) = 321.465;
cent_freq(25) = 3597.63;   bandwidth(25) = 346.136;

bw_min      = bandwidth (1);	   % minimum critical bandwidth

% ----------------------------------------------------------------------
% Set up the critical band filters.  Note here that Gaussianly shaped
% filters are used.  Also, the sum of the filter weights are equivalent
% for each critical band filter.  Filter less than -30 dB and set to
% zero.
% ----------------------------------------------------------------------

min_factor = exp (-30.0 / (2.0 * 2.303));       % -30 dB point of filter

for i = 1:num_crit
  f0 = (cent_freq (i) / max_freq) * (n_fftby2);
  all_f0(i) = floor(f0);
  bw = (bandwidth (i) / max_freq) * (n_fftby2);
  norm_factor = log(bw_min) - log(bandwidth(i));
  j = 0:1:n_fftby2-1;
  crit_filter(i,:) = exp (-11 *(((j - floor(f0)) ./bw).^2) + norm_factor);
  crit_filter(i,:) = crit_filter(i,:).*(crit_filter(i,:) > min_factor);
end

% ----------------------------------------------------------------------
% For each frame of input speech, calculate the Weighted Spectral
% Slope Measure
% ----------------------------------------------------------------------

num_frames = clean_length/skiprate-(winlength/skiprate); % number of frames
start      = 1;					% starting sample
window     = 0.5*(1 - cos(2*pi*(1:winlength)'/(winlength+1)));

for frame_count = 1:num_frames

   % ----------------------------------------------------------
   % (1) Get the Frames for the test and reference speech.
   %     Multiply by Hanning Window.
   % ----------------------------------------------------------

   clean_frame = clean_speech(start:start+winlength-1);
   processed_frame = processed_speech(start:start+winlength-1);
   clean_frame = clean_frame.*window;
   processed_frame = processed_frame.*window;

   % ----------------------------------------------------------
   % (2) Compute the Power Spectrum of Clean and Processed
   % ----------------------------------------------------------

    if (USE_FFT_SPECTRUM)
       clean_spec     = (abs(fft(clean_frame,n_fft)).^2);
       processed_spec = (abs(fft(processed_frame,n_fft)).^2);
    else
       a_vec = zeros(1,n_fft);
       a_vec(1:11) = lpc(clean_frame,10);
       clean_spec     = 1.0/(abs(fft(a_vec,n_fft)).^2)';

       a_vec = zeros(1,n_fft);
       a_vec(1:11) = lpc(processed_frame,10);
       processed_spec = 1.0/(abs(fft(a_vec,n_fft)).^2)';
    end

   % ----------------------------------------------------------
   % (3) Compute Filterbank Output Energies (in dB scale)
   % ----------------------------------------------------------

   for i = 1:num_crit
      clean_energy(i) = sum(clean_spec(1:n_fftby2) ...
		            .*crit_filter(i,:)');
      processed_energy(i) = sum(processed_spec(1:n_fftby2) ...
			        .*crit_filter(i,:)');
   end
   clean_energy = 10*log10(max(clean_energy,1E-10));
   processed_energy = 10*log10(max(processed_energy,1E-10));

   % ----------------------------------------------------------
   % (4) Compute Spectral Slope (dB[i+1]-dB[i])
   % ----------------------------------------------------------

   clean_slope     = clean_energy(2:num_crit) - ...
		     clean_energy(1:num_crit-1);
   processed_slope = processed_energy(2:num_crit) - ...
		     processed_energy(1:num_crit-1);

   % ----------------------------------------------------------
   % (5) Find the nearest peak locations in the spectra to
   %     each critical band.  If the slope is negative, we
   %     search to the left.  If positive, we search to the
   %     right.
   % ----------------------------------------------------------

   for i = 1:num_crit-1

       % find the peaks in the clean speech signal

       if (clean_slope(i)>0) 		% search to the right
	  n = i;
          while ((n<num_crit) & (clean_slope(n) > 0))
	     n = n+1;
 	  end
	  clean_loc_peak(i) = clean_energy(n-1);
       else				% search to the left
          n = i;
	  while ((n>0) & (clean_slope(n) <= 0))
	     n = n-1;
 	  end
	  clean_loc_peak(i) = clean_energy(n+1);
       end

       % find the peaks in the processed speech signal

       if (processed_slope(i)>0) 	% search to the right
	  n = i;
          while ((n<num_crit) & (processed_slope(n) > 0))
	     n = n+1;
	  end
	  processed_loc_peak(i) = processed_energy(n-1);
       else				% search to the left
          n = i;
	  while ((n>0) & (processed_slope(n) <= 0))
	     n = n-1;
 	  end
	  processed_loc_peak(i) = processed_energy(n+1);
       end

   end

   % ----------------------------------------------------------
   %  (6) Compute the WSS Measure for this frame.  This
   %      includes determination of the weighting function.
   % ----------------------------------------------------------

   dBMax_clean       = max(clean_energy);
   dBMax_processed   = max(processed_energy);

   % The weights are calculated by averaging individual
   % weighting factors from the clean and processed frame.
   % These weights W_clean and W_processed should range
   % from 0 to 1 and place more emphasis on spectral
   % peaks and less emphasis on slope differences in spectral
   % valleys.  This procedure is described on page 1280 of
   % Klatt's 1982 ICASSP paper.

   Wmax_clean        = Kmax ./ (Kmax + dBMax_clean - ...
		 	    clean_energy(1:num_crit-1));
   Wlocmax_clean     = Klocmax ./ ( Klocmax + clean_loc_peak - ...
				clean_energy(1:num_crit-1));
   W_clean           = Wmax_clean .* Wlocmax_clean;

   Wmax_processed    = Kmax ./ (Kmax + dBMax_processed - ...
			        processed_energy(1:num_crit-1));
   Wlocmax_processed = Klocmax ./ ( Klocmax + processed_loc_peak - ...
			            processed_energy(1:num_crit-1));
   W_processed       = Wmax_processed .* Wlocmax_processed;

   W = (W_clean + W_processed)./2.0;

   distortion(frame_count) = sum(W.*(clean_slope(1:num_crit-1) - ...
		       processed_slope(1:num_crit-1)).^2);

   % this normalization is not part of Klatt's paper, but helps
   % to normalize the measure.  Here we scale the measure by the
   % sum of the weights.

   distortion(frame_count) = distortion(frame_count)/sum(W);

   start = start + skiprate;

end

%-----------------------------------------------
function distortion = llr(clean_speech, processed_speech,sample_rate)

% ----------------------------------------------------------------------
% Check the length of the clean and processed speech.  Must be the same.
% ----------------------------------------------------------------------

clean_length      = length(clean_speech);
processed_length  = length(processed_speech);

if (clean_length ~= processed_length)
  disp('Error: Both Speech Files must be same length.');
  return
end

% ----------------------------------------------------------------------
% Global Variables
% ----------------------------------------------------------------------

winlength   = round(30*sample_rate/1000); %  window length in samples
skiprate    = floor(winlength/4);		   % window skip in samples
if sample_rate<10000
   P           = 10;		   % LPC Analysis Order
else
    P=16;     % this could vary depending on sampling frequency.
end

% ----------------------------------------------------------------------
% For each frame of input speech, calculate the Log Likelihood Ratio
% ----------------------------------------------------------------------

num_frames = clean_length/skiprate-(winlength/skiprate); % number of frames
start      = 1;					% starting sample
window     = 0.5*(1 - cos(2*pi*(1:winlength)'/(winlength+1)));

for frame_count = 1:num_frames

   % ----------------------------------------------------------
   % (1) Get the Frames for the test and reference speech.
   %     Multiply by Hanning Window.
   % ----------------------------------------------------------

   clean_frame = clean_speech(start:start+winlength-1);
   processed_frame = processed_speech(start:start+winlength-1);
   clean_frame = clean_frame.*window;
   processed_frame = processed_frame.*window;

   % ----------------------------------------------------------
   % (2) Get the autocorrelation lags and LPC parameters used
   %     to compute the LLR measure.
   % ----------------------------------------------------------

   [R_clean, Ref_clean, A_clean] = ...
      lpcoeff(clean_frame, P);
   [R_processed, Ref_processed, A_processed] = ...
      lpcoeff(processed_frame, P);

   % ----------------------------------------------------------
   % (3) Compute the LLR measure
   % ----------------------------------------------------------

   numerator   = A_processed*toeplitz(R_clean)*A_processed';
   denominator = A_clean*toeplitz(R_clean)*A_clean';
   distortion(frame_count) = log(numerator/denominator);
   start = start + skiprate;

end

%---------------------------------------------
function [acorr, refcoeff, lpparams] = lpcoeff(speech_frame, model_order)

   % ----------------------------------------------------------
   % (1) Compute Autocorrelation Lags
   % ----------------------------------------------------------

   winlength = max(size(speech_frame));
   for k=1:model_order+1
      R(k) = sum(speech_frame(1:winlength-k+1) ...
		     .*speech_frame(k:winlength));
   end

   % ----------------------------------------------------------
   % (2) Levinson-Durbin
   % ----------------------------------------------------------

   a = ones(1,model_order);
   E(1)=R(1);
   for i=1:model_order
      a_past(1:i-1) = a(1:i-1);
      sum_term = sum(a_past(1:i-1).*R(i:-1:2));
      rcoeff(i)=(R(i+1) - sum_term) / E(i);
      a(i)=rcoeff(i);
      a(1:i-1) = a_past(1:i-1) - rcoeff(i).*a_past(i-1:-1:1);
      E(i+1)=(1-rcoeff(i)*rcoeff(i))*E(i);
   end

   acorr    = R;
   refcoeff = rcoeff;
   lpparams = [1 -a];

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

function [overall_snr, segmental_snr] = snr(clean_speech, processed_speech,sample_rate)

% ----------------------------------------------------------------------
% Check the length of the clean and processed speech.  Must be the same.
% ----------------------------------------------------------------------

clean_length      = length(clean_speech);
processed_length  = length(processed_speech);

if (clean_length ~= processed_length)
  disp('Error: Both Speech Files must be same length.');
  return
end

% ----------------------------------------------------------------------
% Scale both clean speech and processed speech to have same dynamic
% range.  Also remove DC component from each signal
% ----------------------------------------------------------------------

%clean_speech     = clean_speech     - mean(clean_speech);
%processed_speech = processed_speech - mean(processed_speech);

%processed_speech = processed_speech.*(max(abs(clean_speech))/ max(abs(processed_speech)));

overall_snr = 10* log10( sum(clean_speech.^2)/sum((clean_speech-processed_speech).^2));

% ----------------------------------------------------------------------
% Global Variables
% ----------------------------------------------------------------------

winlength   = round(30*sample_rate/1000); %240;		   % window length in samples
skiprate    = floor(winlength/4);		   % window skip in samples
MIN_SNR     = -10;		   % minimum SNR in dB
MAX_SNR     =  35;		   % maximum SNR in dB

% ----------------------------------------------------------------------
% For each frame of input speech, calculate the Segmental SNR
% ----------------------------------------------------------------------

num_frames = clean_length/skiprate-(winlength/skiprate); % number of frames
start      = 1;					% starting sample
window     = 0.5*(1 - cos(2*pi*(1:winlength)'/(winlength+1)));

for frame_count = 1: num_frames

   % ----------------------------------------------------------
   % (1) Get the Frames for the test and reference speech.
   %     Multiply by Hanning Window.
   % ----------------------------------------------------------

   clean_frame = clean_speech(start:start+winlength-1);
   processed_frame = processed_speech(start:start+winlength-1);
   clean_frame = clean_frame.*window;
   processed_frame = processed_frame.*window;

   % ----------------------------------------------------------
   % (2) Compute the Segmental SNR
   % ----------------------------------------------------------

   signal_energy = sum(clean_frame.^2);
   noise_energy  = sum((clean_frame-processed_frame).^2);
   segmental_snr(frame_count) = 10*log10(signal_energy/(noise_energy+eps)+eps);
   segmental_snr(frame_count) = max(segmental_snr(frame_count),MIN_SNR);
   segmental_snr(frame_count) = min(segmental_snr(frame_count),MAX_SNR);

   start = start + skiprate;

end

author: jtkim-kaist

repo: Speech-enhancement


text:
% Modelling Nonlinear Wave Propagation Example
%
% This example describes the characteristics of the nonlinearity
% encapsulated by the first-order k-Wave simulation functions.
%
% author: Bradley Treeby
% date: 8th December 2011
% last update: 6th September 2013
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

clear all;

% =========================================================================
% DEFINE SIMULATION PROPERTIES
% =========================================================================

% define the properties used in the simulation
p0 = 5e6;                   % source pressure [Pa]
c0 = 1500;                  % sound speed [m/s]
rho0 = 1000;                % density [kg/m^3]
alpha_0 = 0.25;             % absorption coefficient [dB/(MHz^2 cm)]
sigma = 1;                  % shock parameter
source_freq = 1e6;          % frequency [Hz]
points_per_wavelength = 50; % number of grid points per wavelength at f0
wavelength_separation = 15; % separation between the source and detector
pml_size = 64;              % PML size
CFL = 0.25;                 %

% compute corresponding grid spacing
dx = c0/(points_per_wavelength*source_freq);  % [m]

% compute corresponding grid size
Nx = wavelength_separation*points_per_wavelength + 20;

% =========================================================================
% RUN SIMULATION
% =========================================================================

% create the computational grid
kgrid = makeGrid(Nx, dx);

% assign the properties of the propagation medium
medium.sound_speed = c0;
medium.density = rho0;
medium.alpha_power = 2;
medium.alpha_coeff = alpha_0;
medium.alpha_mode = 'no_dispersion';

% extract the maximum frequency supported by the grid
f_max = min(medium.sound_speed)/(2*dx);     % [Hz]

% define a single source element
source.p_mask = zeros(Nx, 1);
source.p_mask(10) = 1;

% define a single sensor position an integer number of wavelengths away
sensor.mask = zeros(Nx, 1);
x_px = wavelength_separation*points_per_wavelength;
sensor.mask(10 + x_px) = 1;
x = x_px*dx;

% compute the nonlinearity coefficient required to give the correct shock
% parameter
mach_num = p0/(rho0*c0.^2);
k = 2*pi*source_freq/c0;
BonA = 2*(sigma/(mach_num*k*x) - 1);
medium.BonA = BonA;

% set the simulation options
input_args = {'PlotFreq', 20, 'PlotScale', [-p0*1.05, p0*1.05],...
    'PMLInside', false, 'PMLSize', pml_size, 'PMLAlpha', 1.5};

% compute points per temporal period
points_per_period = round(points_per_wavelength / CFL);

% compute corresponding time spacing
dt = 1/(points_per_period*source_freq);

% create the time array using an integer number of points per period
t_end = 25e-6;
kgrid.t_array = 0:dt:t_end;

% create the source term, offset by dt/2 so there is a point that lands on the axis
source.p = p0*sin(2*pi*source_freq*(kgrid.t_array + dt/2));

% run the simulation
sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor, input_args{:});

% extract a single wavelength
sensor_data = sensor_data((wavelength_separation + 4)*points_per_period:(wavelength_separation + 5)*points_per_period);

% create time axis for mendousse solution
t_axis = (0:dt:dt*(length(sensor_data) - 1));

% compute mendousse solution for comparison
p_mendousse = mendousse(x*ones(size(t_axis)), t_axis, source_freq, p0, c0, rho0, BonA, alpha_0);

% =========================================================================
% VISUALISATION
% =========================================================================

% plot the time series
figure;
subplot(2, 1, 1), plot(t_axis*1e6, p_mendousse/1e6, 'k-');
hold on;
plot(t_axis*1e6, sensor_data/1e6, 'kx')
xlabel('Time [\mus]');
ylabel('Pressure [MPa]');
legend('Mendousse', 'k-Wave');

% get the amplitude spectra
[f, as_mendousse] = spect(p_mendousse, 1/dt);
f = f./1e6;
as_kspace = spect(sensor_data, 1/dt);

% extract and plot the data at the harmonics
subplot(2, 1, 2);
for harm = 1:10

    % find index of frequency
    [f_val, f_index] = findClosest(f, harm);

    % plot reference
    stem(f_val, as_mendousse(f_index)/1e6, 'Color', 'k', 'Marker', 'o');
    hold on;

    % plot simulation
    plot(f_val, as_kspace(f_index)/1e6, 'kx');

end

% annotate the plot
xlabel('Frequency [MHz]');
ylabel('Amplitude [MPa]');
legend('Mendousse', 'k-Wave');
set(gca, 'XLim', [0.5, 10.5]);

repo: Aero-matlab


text:
function varargout = createRhombododecahedron()
%CREATERHOMBODODECAHEDRON Create a 3D mesh representing a rhombododecahedron
%
%   [V, E, F] = createRhombododecahedron
%   V is a 14-by-3 array with vertex coordinate,
%   E is a 12-by-2 array containing indices of neighbour vertices,
%   F is a 8-by-3 array containing vertices array of each face.
%
%   [V, F] = createRhombododecahedron;
%   Returns only the vertices and the face vertex indices.
%
%   MESH = createRhombododecahedron;
%   Returns the data as a mesh structure, with fields 'vertices', 'edges'
%   and 'faces'.
%
%   Example
%   [v, e, f] = createRhombododecahedron;
%   drawMesh(v, f);
%
%
%   See also
%   meshes3d, drawMesh

%   ---------
%   author : David Legland
%   e-mail: david.legland@inra.fr
%   INRA - TPV URPOI - BIA IMASTE
%   created the 10/02/2005.
%

%   HISTORY
%   04/01/2007: remove unused variables

nodes = [0 0 2;...
    1 -1 1;1 1 1;-1 1 1;-1 -1 1;...
    2 0 0;0 2 0;-2 0 0;0 -2 0;...
    1 -1 -1;1 1 -1;-1 1 -1;-1 -1 -1;...
    0 0 -2];

edges = [...
    1 2;1 3;1 4;1 5;...
    2 6;2 9;3 6;3 7;4 7;4 8;5 8;5 9;...
    6 10;6 11;7 11;7 12;8 12;8 13;9 10;9 13; ...
    10 14;11 14;12 14;13 14];

faces = [...
    1 2 6 3;...
    1 3 7 4;...
    1 4 8 5;...
    1 5 9 2;...
    2 9 10 6;...
    3 6 11 7;...
    4 7 12 8;...
    5 8 13 9;...
    6 10 14 11;...
    7 11 14 12;...
    8 12 14 13;...
    9 13 14 10];

% format output
varargout = formatMeshOutput(nargout, nodes, edges, faces);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function h = surfo(varargin)

% surfo surface plot with workaround surface normals for OpenGL renderer
%
% Syntax:
%     surfo(Z)
%     surfo(Z,C)
%     surfo(X,Y,Z)
%     surfo(X,Y,Z,C)
%     surfo(...,'PropertyName',PropertyValue)
%     surfo(axes_handles,...)
%
% Only use surfo as a subsitute for surf if you need both transparency and
% flat lighting and the OpenGL renderer produces wrong results (see
% Example).
%
% surfo replaces the surface's standard 'VertexNormals' data with the face
% normal data computed by facenorm.
% surfo always switches the figure's renderer to OpenGL as the new normal
% data is not suited for lighting with the zbuffer renderer.
%
% Example:
%
%   % generate some demo data
%   [Y,Z,X] = cylinder([0.8 1 1 0.8],144);
%   idxShift = mod(1:size(X,2),4)>1;
%   X(2,idxShift) = X(2,idxShift)-0.2;
%   X(3,idxShift) = X(3,idxShift)+0.2;
%
%   % standard surf plot
%   figure(1)
%   surf(X,Y,Z,'EdgeAlpha',0.2,'FaceColor','c')
%   axis equal
%   light('Position',[0 0 1],'Style','infinite')
%
%   % same plot with surfo
%   % (giving symmetrical lighting, as it should)
%   figure(10)
%   surfo(X,Y,Z,'EdgeAlpha',0.2,'FaceColor','c')
%   axis equal
%   light('Position',[0 0 1],'Style','infinite')
%
% See also surf, facenorm

% Andres, 2009-07-08

args = varargin;

% plot with surf
hs = surf(args{:});

% switch to OpenGL renderer
fig = ancestor(hs,'figure');
set(fig,'Renderer','OpenGL')

% extract surface coordinate data
X = get(hs,'XData');
Y = get(hs,'YData');
Z = get(hs,'ZData');

% use workaround normals
set(hs,'VertexNormals',facenorm(X,Y,Z));

if nargout > 0
    h = hs;
end


text:
function ob = mridft(kx, ky, t, we, xval, yval, fov, N)
%function ob = mridft(kx, ky, t, we, xval, yval, fov, N)
%	Construct MRI object, which can do Ax and A'y operations
%       Brad Sutton    University of Michigan

%	default object
ob.kx = 0;	% should these be []'s
ob.ky = 0;
ob.t = 0;
ob.we = 0;
ob.xval = 0;
ob.yval = 0;
ob.is.empty	= true;
ob.is.transpose = false;
ob.fov = 0;  %FOV in cm's
ob.n = 0;
%ob.version = 1.0;

if nargin == 0
	ob = class(ob, 'mridft');
	return
end

if isa(kx, 'mridft')
	ob = kx;
	return
end

if nargin ~= 8
	help mridft
	error nargin
end

	%	fill object
	ob.kx = kx;
	ob.ky = ky;
	ob.t = t;
	ob.we = we;
	ob.xval = xval;
	ob.yval = yval;
        ob.fov = fov;
        ob.n = N;

	ob.is.empty	= false;

%	ob.m = size(we);	% image size
%	ob.n = size(we);

	ob = class(ob, 'mridft');

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% This demo is for classification problem that has parameter(s) to
% optimize. If more parameters to be optimized, more inner LOOCV runs
% (additional for-loop) should be used. All networks corresponding to
% different combinations of parameters should be prepared beforehand.
% Note that this function uses the LOOCV.

%
% Input:
%     result_dir: the directory you want to store all the results in;
%     meth_Net: the brain network construction method;
%     BOLD: time courses extracted using a template;
%     label: the label for each subject; e.g., -1 for normal controls and 1 for patients;
%     para_test_flag: whether or not to perform parameter sensitivity test.
%     varargin: the different parameters for different brain network construction method;
%
% Output:
%     opt_paramt: the best parameter combination in the parameter sensitivity test;
%     opt_t: the selected parameter combination in each parameter optimization fold;
%     AUC,SEN,SPE,F1,Acc,Youden,BalanceAccuracy: model performance;
%     w: weight of each selected features;
%     varargout: the feature selection index for the following finding features section.

% Written by Zhen Zhou, zzstefan@email.unc.edu
% IDEA lab, https://www.med.unc.edu/bric/ideagroup
% Department of Radiology and BRIC, University of North Carolina at Chapel Hill

function [opt_paramt,AUC,SEN,SPE,F1,Acc,opt_t,ttest_p,midw_lasso,w,plot_ROC]=param_select_demo(result_dir,meth_Net,BOLD,label,para_test_flag,varargin)
switch meth_Net
    case {'SR','WSR','GSR'}
        lambda_1=varargin{1};
        lambda_lasso=varargin{2};
    case {'SLR','SGR','WSGR','SSGSR'}
        lambda_1=varargin{1};
        lambda_2=varargin{2};
        lambda_lasso=varargin{3};
    case 'dHOFC'
        W=varargin{1};
        s=varargin{2};
        C=varargin{3};
        lambda_lasso=varargin{4};
end
%clear; close all;
%meth_Net='WSR'
root=pwd;
addpath(genpath([root '/Function']));
addpath(genpath([root '/Toolbox']));

%% Load BOLD signals
%load('demo_param_select_toydata.mat');    % This is just an example data including 10 subjects only
%load('demo_framwk_toydata.mat');
% You can replace this by your own dataset
% BOLD: A cell array, each matrix includes BOLD signals (137 volumes x 116 ROIs) of one subject
% label: class labels, 1 for patient and -1 for healthy controls

[~,nROI]=size(BOLD{1});
nSubj=length(BOLD);

% Construct FC networks under different parameters
% mkdir ./Generated_BrainNet_test
% dir=['./Generated_BrainNet_test/'];

fprintf('Begin network construction\n');
switch meth_Net
    case 'SR'       % Sparse representation
        lambda=lambda_1; % parameter for sparsity
        parfor i=1:length(lambda)
            BrainNet{i}=SR(BOLD,lambda(i));
        end
        %save (char(strcat(dir,meth_Net,'_net.mat')),'BrainNet','-v7.3');
    case 'WSR'      % PC weighted SR
        lambda=lambda_1; % parameter for sparsity
        parfor i=1:length(lambda)
            BrainNet{i}=WSR(BOLD,lambda(i));
        end
        %save (char(strcat(dir,meth_Net,'_net.mat')),'BrainNet','-v7.3');

    case 'SLR'      % Sparse low-rank representation
        lambda1=lambda_1;
        lambda2=lambda_2;
        num_lambda1=length(lambda1);
        num_lambda2=length(lambda2);
        parfor i=1:num_lambda1
            for j=1:num_lambda2
                BrainNet{i,j}=SLR(BOLD,lambda1(i),lambda2(j));
            end
        end
        BrainNet=reshape(BrainNet,1,num_lambda1*num_lambda2);
        %save (char(strcat(dir,meth_Net,'_net.mat')),'BrainNet','-v7.3');

    case 'SGR'      % Sparse group representation
        lambda1=lambda_1; % parameter for sparsity
        lambda2=lambda_2; % parameter for group sparsity
        num_lambda1=length(lambda1);
        num_lambda2=length(lambda2);
        parfor i=1:num_lambda1
            for j=1:num_lambda2
                BrainNet{i,j}=SGR(BOLD,lambda1(i),lambda2(j));
            end
        end
        BrainNet=reshape(BrainNet,1,num_lambda1*num_lambda2);
        save (char(strcat(result_dir,meth_Net,'_net.mat')),'BrainNet','-v7.3');

    case 'WSGR'     % PC weighted SGR
        lambda1=lambda_1; % parameter for sparsity
        lambda2=lambda_2; % parameter for group sparsity
        num_lambda1=length(lambda1);
        num_lambda2=length(lambda2);
        parfor i=1:num_lambda1
            for j=1:num_lambda2
                BrainNet{i,j}=WSGR(BOLD,lambda1(i),lambda2(j));
            end
        end
        BrainNet=reshape(BrainNet,1,num_lambda1*num_lambda2);
        %save (char(strcat(dir,meth_Net,'_net.mat')),'BrainNet','-v7.3');

    case 'GSR'      % Group sparse representation
        lambda=lambda_1;
        parfor i=1:length(lambda)
            BrainNet{i}=GSR(BOLD,lambda(i));
        end
        %save (char(strcat(dir,meth_Net,'_net.mat')),'BrainNet','-v7.3');

    case 'SSGSR'    % Strength and Similarity guided GSR
        lambda1=lambda_1; % parameter for group sparsity
        lambda2=lambda_2; % parameter for inter-subject LOFC-pattern similarity
        num_lambda1=length(lambda1);
        num_lambda2=length(lambda2);
        parfor i=1:num_lambda1
            for j=1:num_lambda2
                BrainNet{i,j}=SSGSR(BOLD,lambda1(i),lambda2(j));
            end
        end
        BrainNet=reshape(BrainNet,1,num_lambda1*num_lambda2);
        %save (char(strcat(dir,meth_Net,'_net.mat')),'BrainNet','-v7.3');

    case 'dHOFC'    % Dynamic high-order FC
        num_C=length(C);
        num_W=length(W);
        parfor i=1:num_W % number of clusters
            for j=1:num_C
                [BrainNet{i,j},IDX{i,j}]=dHOFC(BOLD,W(i),s,C(j));
            end
        end
        BrainNet=reshape(BrainNet,1,num_W*num_C);
        %save (char(strcat(dir,meth_Net,'_net.mat')),'BrainNet','-v7.3');
end
fprintf('Network construction finished\n');

% Note, usually, multiple networks with different parameters will be
% established beforehand so that the results can be stored and called
% during parameter selection to save computational time

% Extract features from FC networks under different parameters
if ~strcmpi(meth_Net,'dHOFC')
    idxtu=triu(ones(nROI,nROI));
    for i=1:length(BrainNet)
        Feat=zeros(nSubj,nROI*(nROI-1)/2);  % choose FC strengths as features
        for j=1:nSubj
            tempNet=BrainNet{i}(:,:,j);
            %tempA=tempNet(idxtu~=0);
            Feat(j,:)=tempNet(idxtu==0); % only use the upper trangle coefficients since
            % brain network has been symmetric
        end
        All_Feat{i} = Feat;
    end

else
    %idxtu=triu(ones(nROI,nROI),1);
    for i=1:length(BrainNet)
        temp=ceil(i/length(W));
        Feat=zeros(nSubj,C(temp));
        flag=2;
        for j=1:nSubj
            Feat(j,:)=wlcc(BrainNet{i}(:,:,j),flag);
        end
        All_Feat{i} = Feat;
    end
end

e = 1:nSubj;
indx = floor(nSubj*[1:10]/10);
cpred = zeros(nSubj,1);
acc = zeros(nSubj,1);
score = zeros(nSubj,1);

% LOOCV for testing
for i=1:nSubj
    ind = find(indx == i);
    if ~isempty(ind)
        fprintf(1,'Begin process %d%%...\n',ind*10);
    end

    Tst_ind = i;
    telabel = label(i);
    Trn_ind = e;
    Trn_ind(i) = [];
    trlabel = label;
    trlabel(i) = [];

    % Nested LOOCV on Train data for Model selection
    max_acc = 0;
    for t = 1:length(All_Feat)
        % Feature generation
        Feat = All_Feat{t};
        tmpTestCorr = zeros(length(trlabel),1);
        for j =1:length(trlabel)
            % LOOCV on Training Set
            Tst1_ind = Trn_ind(j);
            telabel1 = trlabel(j);
            Trn1_ind = Trn_ind;
            Trn1_ind(j) = [];
            trlabel1 = trlabel;
            trlabel1(j) = [];

            trFe = Feat(Trn1_ind,:);
            teFe = Feat(Tst1_ind,:);

            % Feature selection using t-test
            if ~strcmpi(meth_Net,'dHOFC')
                pval=0.05;  % generally use pval<0.05 as a threshold for feature selection
                [~,p]=ttest2(trFe(trlabel1==-1,:),trFe(trlabel1==1,:));
                trFe=trFe(:,p<pval);
                teFe=teFe(:,p<pval);

                midw=lasso(trFe,trlabel1,'Lambda',lambda_lasso);  % parameter lambda for sparsity
                trFe=trFe(:,midw~=0);
                teFe=teFe(:,midw~=0);
            else
                midw=lasso(trFe,trlabel1,'Lambda',lambda_lasso);  % parameter lambda for sparsity
                trFe=trFe(:,midw~=0);
                teFe=teFe(:,midw~=0);
            end

            % Feature normalization
            Mtr=mean(trFe);
            Str=std(trFe);
            trFe=trFe-repmat(Mtr,size(trFe,1),1);
            trFe=trFe./repmat(Str,size(trFe,1),1);
            teFe=teFe-Mtr;
            teFe=teFe./Str;

            % train SVM model
            classmodel=svmtrain(trlabel1,trFe,'-t 0 -c 1 -q'); % linear SVM (require LIBSVM toolbox)
            % classify
            [~,acc,~]=svmpredict(telabel1,teFe,classmodel,'-q');
            tmpTestCorr(j,1) = acc(1);
        end
        mTestCorr(t) = sum(tmpTestCorr)/length(tmpTestCorr);
        if mTestCorr(t)>max_acc
            max_acc = mTestCorr(t);
            opt_t(i) = t;
        end
    end
    % Feature generation
    Feat = All_Feat{opt_t(i)};
    trFe = Feat(Trn_ind,:);
    teFe = Feat(Tst_ind,:);

    % Feature selection ag
    if ~strcmpi(meth_Net,'dHOFC')
        pval=0.05;  % generally use pval<0.05 as a threshold for feature selection
        [~,p]=ttest2(trFe(trlabel==-1,:),trFe(trlabel==1,:));
        trFe=trFe(:,p<pval);
        teFe=teFe(:,p<pval);
        ttest_p{i}=p;

        midw=lasso(trFe,trlabel,'Lambda',lambda_lasso);  % parameter lambda for sparsity
        trFe=trFe(:,midw~=0);
        teFe=teFe(:,midw~=0);
        midw_lasso{i}=midw;
    else
        midw=lasso(trFe,trlabel,'Lambda',lambda_lasso);  % parameter lambda for sparsity
        trFe=trFe(:,midw~=0);
        teFe=teFe(:,midw~=0);
        midw_lasso{i}=midw;
        ttest_p{i}=midw;

    end

    % Feature normalization ag
    Mtr=mean(trFe);
    Str=std(trFe);
    trFe=trFe-repmat(Mtr,size(trFe,1),1);
    trFe=trFe./repmat(Str,size(trFe,1),1);
    teFe=teFe-Mtr;
    teFe=teFe./Str;

    % train SVM model ag
    classmodel=svmtrain(trlabel,trFe,'-t 0 -c 1 -q'); % linear SVM (require LIBSVM toolbox)
    % classify ag
    w{i}=classmodel.SVs'*classmodel.sv_coef;
    [cpred(i),~,score(i)]=svmpredict(telabel,teFe,classmodel,'-q');
end

Acc=100*sum(cpred==label)/nSubj;
[AUC,SEN,SPE,F1,Youden,BalanceAccuracy,plot_ROC]=perfeval(label,cpred,score,result_dir);

if para_test_flag==1
    fprintf('Begin parameter sensitivity test\n');
    for i=1:length(All_Feat)
        Acc_para(i)=parameter_sensitivity_test_loocv(All_Feat{i},nSubj,label,meth_Net,lambda_lasso);
    end
    save (char(strcat(result_dir,'/Acc_para.mat')),'Acc_para');
      switch meth_Net
            case {'SR','WSR','GSR'}
                [opt_paramt]=select_para(meth_Net,Acc_para,lambda);
                save_model(BrainNet,meth_Net,label,result_dir,opt_paramt,Acc_para,lambda_lasso);
                x=lambda;
                for j=1:length(x)
                    x_label{j}=num2str(x(j));
                end
                ind=1:length(x);
                y=Acc_para;
                figure('visible','off');
                bar(ind,y);

                ylim([0,100]);
                set(gca,'XTickLabel',x_label);
                xlabel(sprintf('\x3bb_1'));
                ylabel('Accuracy');
                title('Parameter sensitivity test');
                print(gcf,'-r1000','-dtiff',char(strcat(result_dir,'/para_sensitivity.tiff')));
                %print(gcf,'-depsc',char(strcat(result_dir,'/para_sensitivity.eps')));
             case {'SLR','SGR','WSGR','SSGSR'}
                 [opt_paramt]=select_para(meth_Net,Acc_para,lambda1,lambda2);
                 save_model(BrainNet,meth_Net,label,result_dir,opt_paramt,Acc_para,lambda_lasso);
                 x=lambda2;
                 y=lambda1;
                 for i=1:length(x)
                    x_label{i}=num2str(x(i));
                 end
                 for j=1:length(y)
                    y_label{j}=num2str(y(j));
                 end
                 z=reshape(Acc_para,num_lambda1,num_lambda2);
                 figure('visible','off');
                 Bar1=bar3(z);
                 for Element = 1:length(Bar1)
                    ZData = get(Bar1(Element),'ZData');
                    set(Bar1(Element), 'CData', ZData,...
                        'FaceColor', 'interp');
                 end
                 colorbar
                 zlim([0,100]);
                 xlabel(sprintf('\x3bb_2'));
                 ylabel(sprintf('\x3bb_1'));
                 zlabel(sprintf('Accuracy'));
                 if length(x)>=8
                     ind_x=1:2:length(x);
                     set(gca,'XTick',ind_x,'XTickLabel',x_label(1:2:end));
                 else
                     ind_x=1:length(x);
                     set(gca,'XTick',ind_x,'XTickLabel',x_label);
                 end
                 if length(y)>=8
                     ind_y=1:2:length(y);
                     set(gca,'YTick',ind_y,'YTickLabel',y_label(1:2:end));
                 else
                     ind_y=1:length(y);
                     set(gca,'YTick',ind_y,'YTickLabel',y_label);
                 end
                 title('Parameter sensitivity test');
                 print(gcf,'-r1000','-dtiff',char(strcat(result_dir,'/para_sensitivity.tiff')));
                 %print(gcf,'-depsc',char(strcat(result_dir,'/para_sensitivity.eps')));
             case 'dHOFC'
                 [opt_paramt]=select_para(meth_Net,Acc_para,W,C);
                 save_model(BrainNet,meth_Net,label,result_dir,opt_paramt,Acc_para,lambda_lasso);
                 x=C;
                 y=W;
                 for i=1:length(x)
                    x_label{i}=num2str(x(i));
                 end
                 for j=1:length(y)
                    y_label{j}=num2str(y(j));
                 end
                 z=reshape(Acc_para,length(W),length(C));
                 figure('visible','off');
                 Bar1=bar3(z);
                 for Element = 1:length(Bar1)
                    ZData = get(Bar1(Element),'ZData');
                    set(Bar1(Element), 'CData', ZData,...
                        'FaceColor', 'interp');
                 end
                 colorbar
                 zlim([0,100]);
                 xlabel(sprintf('Number of Clusters'));
                 ylabel(sprintf('Window Length'));
                 zlabel(sprintf('Accuracy'));
                 if length(x)>=8
                     ind_x=1:2:length(x);
                     set(gca,'XTick',ind_x,'XTickLabel',x_label(1:2:end));
                 else
                     ind_x=1:length(x);
                     set(gca,'XTick',ind_x,'XTickLabel',x_label);
                 end
                 if length(y)>=8
                     ind_y=1:2:length(y);
                     set(gca,'YTick',ind_y,'YTickLabel',y_label(1:2:end));
                 else
                     ind_y=1:length(y);
                     set(gca,'YTick',ind_y,'YTickLabel',y_label);
                 end
                 title('Parameter sensitivity test');
                 print(gcf,'-r1000','-dtiff',char(strcat(result_dir,'/para_sensitivity.tiff')));
      end
      fprintf('End parameter sensitivity test\n');
end

k_times=10;
switch meth_Net
        case {'SR','WSR','GSR'}
            [result_features]=back_find_low_node_para(result_dir,nSubj,k_times,nROI,w,'loocv',ttest_p,midw_lasso);
            write_log(result_dir,meth_Net,'loocv',AUC,SEN,SPE,F1,Acc,Youden,BalanceAccuracy,lambda_1,lambda_lasso,opt_paramt,k_times,opt_t);
        case {'SGR','WSGR','SSGSR','SLR'}
            [result_features]=back_find_low_node_para(result_dir,nSubj,k_times,nROI,w,'loocv',ttest_p,midw_lasso);
            write_log(result_dir,meth_Net,'loocv',AUC,SEN,SPE,F1,Acc,Youden,BalanceAccuracy,lambda_1,lambda_2,lambda_lasso,opt_paramt,k_times,opt_t);
        case 'dHOFC'
            [result_features]=back_find_high_node(W,C,nROI,w,midw_lasso,IDX,opt_t);
            write_log(result_dir,meth_Net,'loocv',AUC,SEN,SPE,F1,Acc,Youden,BalanceAccuracy,W,s,C,lambda_lasso,opt_paramt,k_times,opt_t);
end

save (char(strcat(result_dir,'/result_features.mat')),'result_features');

author: zzstefan

repo: BrainNetClass


text:
function ma_handlePlayButton(handles, frameRate)
%ma_handlePlayButton Summary of this function goes here
%   Detailed explanation goes here
    global GLOBAL_VideoWriter GLOBAL_isRecording;
    persistent mTimer;

    if(get(handles.playButton,'Value')==1)
        set(handles.playButton,'String','<html>■');
        delete(mTimer);

        maData = getappdata(handles.ma_MainGUI,'ma_data');
        celBodyData = getappdata(handles.ma_MainGUI,'celBodyData');
        stateLog = maData.stateLog;

        time = get(handles.movieFrameSlider,'Value');
        [orbitToPlot] = ma_getOrbitToPlot(stateLog, time);

        hold(handles.movieAxes,'on');
        period = 1/frameRate;
        period = round(period*1000)/1000;
        tFunc = @(obj,event) ma_playMovieTimerCallback(obj, event, handles.playButton, handles.ma_MissionAnimatorGUI, handles.movieAxes, maData, frameRate, celBodyData, handles);
        mTimer = timer('TimerFcn',tFunc,'Period',period,'ExecutionMode','fixedRate','UserData',[NaN, NaN, orbitToPlot]);
        hold(handles.movieAxes,'off');

        if(GLOBAL_isRecording)
            GLOBAL_VideoWriter.FrameRate = frameRate;
            open(GLOBAL_VideoWriter);
        end

        drawnow;
        start(mTimer);
    else
        set(handles.playButton,'String','<html>&#9658');

        if(GLOBAL_isRecording)
            try
                close(GLOBAL_VideoWriter);
            catch
                %no action
            end
            GLOBAL_isRecording = false;
            GLOBAL_VideoWriter = [];
        end
        drawnow;
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot