text:
function song = sing(str, bpm)
%SING   A basic keyboard for MATLAB using CHEBFUNs.
%   S = SING('STR1 STR2 STR3 ...') creates a CHEBFUN F corresponding to the
%   musical notes in each of the STR which vaguely correspond to Helmholtz pitch
%   notation (http://en.wikipedia.org/wiki/Helmholtz_pitch_notation).
%
%   S = SING('STR1 ...', BPM) alters the temp of the tune (in beats per minute).
%   The default value is 60bpm.
%
%   For example, S = SING('A'), produces a function corresponding to the note
%   A, S = SING('B') produces a B, and so on.
%
%   Basic chords are provided, for example S = SING('CEG') will produce a major
%   triad (http://en.wikipedia.org/wiki/Major_chord). The convention is for the
%   kth note to have 1/k times the amplitude of the first.
%
%   Empty spaces correspond to pause of a sixteenth note (semiquaver), and
%   commas separate notes without a pause.
%
%   If STR is in uppercase, it will last for a quarter note (crotchet), and
%   lower case notes for a sixteenth note (semiquaver).
%
%   The SING keyboard has three octaves (in A to G). S = SING('C-') produces a
%   low C, S = SING('C') produces a middle C, and S = SING('C+') produces a
%   high C.
%
%   Sharps are supported by using '#', but there is currently no notation for
%   flats.
%
%   Examples:
%    1. Ode to Joy:
%       S = SING('BG- BG- CA DB DB CA BG- AD- G-B G-B AD- BG- BG- AD- AD-');
%    2. God Save The Queen
%       S = SING('CG CG DG BG CG DG EC EC FC EC DG CG DG CG BG CG',30);
%
%   CHEBFUNs for these examples may be generated with SING(1) and SING(2),
%   respectively.
%
% See also SOUND, CHEBTUNE.

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

% Set up basic frequencies.
s0 = 2^(1/12);
f0 = 440*pi;
fp1 = 2*f0;
fm1 = .5*f0;

% Set tempo.
if ( nargin < 2 )
    bpm = 60;
    t0 = 1/4;
else
    t0 = 60/bpm/4;
end

q = chebfun('x', [0 t0]);   % Quarter note.
s = chebfun('x', [0 t0/4]); % Sixteenth note.

%% Middle

% Quarter notes
A = sin(f0*q);
As = sin(f0*s0*q);
B = sin(f0*s0^2*q);
C = sin(f0*s0^3*q);
Cs = sin(f0*s0^4*q);
D = sin(f0*s0^5*q);
Ds = sin(f0*s0^6*q);
E = sin(f0*s0^7*q);
F = sin(f0*s0^8*q);
Fs = sin(f0*s0^9*q);
G = sin(f0*s0^10*q);
Gs = sin(f0*s0^11*q);

% Sixteenth notes
a = sin(f0*s);
as = sin(f0*s0*s);
b = sin(f0*s0^2*s);
c = sin(f0*s0^3*s);
cs = sin(f0*s0^4*s);
d = sin(f0*s0^5*s);
ds = sin(f0*s0^6*s);
e = sin(f0*s0^7*s);
f = sin(f0*s0^8*s);
fs = sin(f0*s0^9*s);
g = sin(f0*s0^10*s);
gs = sin(f0*s0^11*s);

%% -

% Quarter notes
Am = sin(fm1*q);
Ams = sin(fm1*s0*q);
Bm = sin(fm1*s0^2*q);
Cm = sin(fm1*s0^3*q);
Cms = sin(fm1*s0^4*q);
Dm = sin(fm1*s0^5*q);
Dms = sin(fm1*s0^6*q);
Em = sin(fm1*s0^7*q);
Fm = sin(fm1*s0^8*q);
Fms = sin(fm1*s0^9*q);
Gm = sin(fm1*s0^10*q);
Gms = sin(fm1*s0^11*q);

% Sixteenth notes
am = sin(fm1*s);
ams = sin(fm1*s0*s);
bm = sin(fm1*s0^2*s);
cm = sin(fm1*s0^3*s);
cms = sin(fm1*s0^4*s);
dm = sin(fm1*s0^5*s);
dms = sin(fm1*s0^6*s);
em = sin(fm1*s0^7*s);
fm = sin(fm1*s0^8*s);
fms = sin(fm1*s0^9*s);
gm = sin(fm1*s0^10*s);
gms = sin(fm1*s0^11*s);

%% +

% Quarter notes
Ap = sin(fp1*q);
Aps = sin(fp1*s0*q);
Bp = sin(fp1*s0^2*q);
Cp = sin(fp1*s0^3*q);
Cps = sin(fp1*s0^4*q);
Dp = sin(fp1*s0^5*q);
Dps = sin(fp1*s0^6*q);
Ep = sin(fp1*s0^7*q);
Fp = sin(fp1*s0^8*q);
Fps = sin(fp1*s0^9*q);
Gp = sin(fp1*s0^10*q);
Gps = sin(fp1*s0^11*q);

% Sixteenth notes
ap = sin(fp1*s);
aps = sin(fp1*s0*s);
bp = sin(fp1*s0^2*s);
cp = sin(fp1*s0^3*s);
cps = sin(fp1*s0^4*s);
dp = sin(fp1*s0^5*s);
dps = sin(fp1*s0^6*s);
ep = sin(fp1*s0^7*s);
fp = sin(fp1*s0^8*s);
fps = sin(fp1*s0^9*s);
gp = sin(fp1*s0^10*s);
gps = sin(fp1*s0^11*s);

%%

% Easy access to example tunes.
if ( (nargin == 0) || isnumeric(str) )
    if ( nargin == 0 )
        str = 1;
    end

    switch str
        case 1 % Ode to Joy
            str = 'BG- BG- CA DB DB CA BG- AD- G-B G-B AD- BG- BG- AD- AD-';
        case 2 % God Save The Queen
            str = 'CG CG DG BG CG DG EC EC FC EC DG CG DG CG BG CG';
    end
end

% Build the CHEBFUN from the string of notes for the tune.
song = chebfun;
l = 0;
while ( ~isempty(str) )
    l = l + 1;
    k = 0;

    % Get the string for the next chord.
    for j = 1:numel(str)
        k = k + 1;
        if ( strcmp(str(k), ' ') || strcmp(str(k), ',') )
            break
        end
    end

    strk = str(1:k);
    str(1:k) = [];

    % Quarter note or sixteenth note?
    if ( isstrprop(strk(1), 'upper') )
        strk = upper(strk);
        songtmp = 0*A;
    else
        strk = lower(strk);
        songtmp = 0*a;
    end

    % Assemble the chord from its constituent notes.
    j = 0;
    while ( ~isempty(strk) )
        % Parse out one note from the chord.
        if ( (length(strk) > 1) && any(strcmpi(strk(2), {'+' '-' '#' '*'})) )
            if ( (length(strk) > 2) && any(strcmpi(strk(3), {'#' '*'})) )
                strjk = strk(1:3);
                strk(1:3) = [];
            else
                strjk = strk(1:2);
                strk(1:2) = [];
            end
        else
            strjk = strk(1);
            strk(1) = [];
        end

        % Set the amplitude of this note in the chord.
        j = j + 1;
        amp = .4/j;

        % Decode and add in the note.
        if ( strcmp(strjk, 'r') || strcmp(strjk, ' ') )  % Sixteenth rest.
            songtmp = join(songtmp, 0*s);
            break
        elseif ( strcmp(strjk, 'R') )                    % Quarter rest.
            songtmp = join(songtmp, 0*q);
            break
        else                                             % Actual note.
            songtmp = songtmp + amp*getNote(strjk);
        end
    end

    % Tack the chord onto the end of the song.
    song = join(song, songtmp);
end

    % Nested function for mapping note strings to note CHEBFUNs defined above.
    function note = getNote(s)
        switch (s)
            case 'A', note = A;
            case 'B', note = B;
            case 'C', note = C;
            case 'D', note = D;
            case 'E', note = E;
            case 'F', note = F;
            case 'G', note = G;

            case 'a', note = a;
            case 'b', note = b;
            case 'c', note = c;
            case 'd', note = d;
            case 'e', note = e;
            case 'f', note = f;
            case 'g', note = g;

            case 'A-', note = Am;
            case 'B-', note = Bm;
            case 'C-', note = Cm;
            case 'D-', note = Dm;
            case 'E-', note = Em;
            case 'F-', note = Fm;
            case 'G-', note = Gm;

            case 'a-', note = am;
            case 'b-', note = bm;
            case 'c-', note = cm;
            case 'd-', note = dm;
            case 'e-', note = em;
            case 'f-', note = fm;
            case 'g-', note = gm;

            case 'A+', note = Ap;
            case 'B+', note = Bp;
            case 'C+', note = Cp;
            case 'D+', note = Dp;
            case 'E+', note = Ep;
            case 'F+', note = Fp;
            case 'G+', note = Gp;

            case 'a+', note = ap;
            case 'b+', note = bp;
            case 'c+', note = cp;
            case 'd+', note = dp;
            case 'e+', note = ep;
            case 'f+', note = fp;
            case 'g+', note = gp;

            case 'a#', note = as;
            case 'c#', note = cs;
            case 'd#', note = ds;
            case 'f#', note = fs;
            case 'g#', note = gs;

            case 'A#', note = As;
            case 'C#', note = Cs;
            case 'D#', note = Ds;
            case 'F#', note = Fs;
            case 'G#', note = Gs;

            case 'a+#', note = aps;
            case 'c+#', note = cps;
            case 'd+#', note = dps;
            case 'f+#', note = fps;
            case 'g+#', note = gps;

            case 'A+#', note = Aps;
            case 'C+#', note = Cps;
            case 'D+#', note = Dps;
            case 'F+#', note = Fps;
            case 'G+#', note = Gps;

            case 'a-#', note = ams;
            case 'c-#', note = cms;
            case 'd-#', note = dms;
            case 'f-#', note = fms;
            case 'g-#', note = gms;

            case 'A-#', note = Ams;
            case 'C-#', note = Cms;
            case 'D-#', note = Dms;
            case 'F-#', note = Fms;
            case 'G-#', note = Gms;
        end
    end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function varargout = plot_spline(P,C)
  % PLOT_SPLINE Plot a cubic Bezier spline.
  %
  % [pe,p] = plot_spline(P,C)
  %
  % Inputs:
  %   P  #P by dim list of control point locations
  %   C  #C by 4 list of indices into P of cubic Bezier curves
  % Outputs:
  %   pe  plot handles for UI
  %   p  plot handle for curves
  % Example:
  %   cellfun(@(pe) arrayfun(@(p) set(p,'Color','r'),pe),plot_spline(P,C))
  p = {};
  pe = {};
  ish = ishold;
  for c = 1:size(C,1)
    [pe{c},p{c}] = plot_cubic(P(C(c,:),:));
    hold on;
  end
  hold off;
  if ish
    hold on
  end
  if nargout>=1
    varargout{1} = pe;
    if nargout >= 2
      varargout{2} = p;
    end
  end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function value = r8_gamma_sample ( a, r )

%*****************************************************************************80
%
%% R8_GAMMA_SAMPLE generates a Gamma random deviate.
%
%  Discussion:
%
%    This procedure generates random deviates from the gamma distribution whose
%    density is (A^R)/Gamma(R) * X^(R-1) * Exp(-A*X)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2013
%
%  Author:
%
%    Original FORTRAN77 version by Barry Brown, James Lovato.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Joachim Ahrens, Ulrich Dieter,
%    Generating Gamma Variates by a Modified Rejection Technique,
%    Communications of the ACM,
%    Volume 25, Number 1, January 1982, pages 47-54.
%
%    Joachim Ahrens, Ulrich Dieter,
%    Computer Methods for Sampling from Gamma, Beta, Poisson and
%    Binomial Distributions,
%    Computing,
%    Volume 12, Number 3, September 1974, pages 223-246.
%
%  Parameters:
%
%    Input, real A, the rate parameter.
%
%    Input, real R, the shape parameter.
%
%    Output, real VALUE, a random deviate
%    from the distribution.
%
  value = r8_gamma_01_sample ( r ) / a;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%
% surf_registration_table.m
%
% Original Author: Laurence Wastiaux
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:13 $
%    $Revision: 1.3 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

[Dsubj, nf, r90, thdemean, M] = surf_registration_stats('/space/neo/2/recon/buckner')
%load('/space/okapi/3/data/laurence/ADF/surf_registration/residual_var90.mat') %loads r90
%load('/space/okapi/3/data/laurence/ADF/surf_registration/Subj_Names.mat') %loads nf
%files=dir('/space/neo/2/recon/buckner');
fp=fopen('~/surfreg_table.txt', 'w');
fprintf(fp,'# FSAFD SurfRegistrationCheck 1\n# date 20050516\n# $Id: surf_registration_table.m,v 1.3 2011/03/02 00:04:13 nicks Exp $\n# Buckner data set\n# hemi lh\n# PercentRemovedVar 90\n# ncols 1\n# nrows %d\n# label_col 1 PercentResidualVariance\n', length(r90));
nn=1;
for i=1:length(r90)
    %nb=I(i);
    if(r90(i)>0)
        fprintf(fp, '# label_row %d %s\n', nn, char(nf(i).name));
        nn=nn+1;
    end

end
for j=1:length(r90)
    if(r90(j)>0)
        fprintf(fp, '%f\n', r90(j));
    end
end
fclose(fp)

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef testNotCommutingHomogPlaneStress < ...
            testCommutingHomogPlaneStress

    methods (Access = protected, Static)

        function nu = createPoissonValue()
            nu = 1/3;
        end
    end

    methods (Access = public)

        function hasPassed = hasPassed(obj)
            c1 = obj.vhpTensor.getValue();
            c2 = obj.vphTensor.getValue();
            error = norm(c2-c1)/norm(c1);
            hasPassed = error > 1e-6;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [err,Loc] = FindChar (C, cfnd)
%
%   [err,Loc] = FindChar (C, cfnd)
%
%Purpose:
%   find certain characters in an array of characters, or a string
%
%
%Input Parameters:
%   C : an array of characters or a string
%   cfnd : the character you'r looking for like 'f'
%   	pass the following for special characters
%		'NewLine', 'Tab', 'Space', 'Comma'
%
%Output Parameters:
%   err : 0 No Problem
%       : 1 Mucho Problems
%   Loc : a vector of the location in C containing cfnd
%      ie any element Loc of C is cnfd. C(Loc) = cfnd
%
%
%More Info :
%
%   example :
%	ks2 = sprintf('d, v h\t\nc c , c d\n vf');
%	[err, Loc] = FindChar(ks2, ',')
% or in the previous case, you could use
%  [err, Loc] = FindChar(ks2, 'Comma')
% or
%  [err, Loc] = FindChar(ks2, 'Tab')
%
% ks2 could be a whole ascii file or a chunk of one, like
% fid = fopen('thisfile','r'); ks2 = fscanf(fid,'%c',100); fclose (fid);
% the previous line reads the 1st 100 characters
% carefull, if you use
% ks2 = fscanf(fid,'%s',100);
% you'll read the first 100 strings (space, tab and newline delimited)
%
% see also
%  SkipMatlabHelp
%  PurgeComments
%	NextString
%
%     Author : Ziad Saad
%     Date : Sat Mar 27 13:42:26 CST 1999

%Define the function name for easy referencing
FuncName = 'FindChar';

%Debug Flag
DBG = 1;

%initailize return variables
err = 1;
Loc = [];

switch cfnd,
	%to find the ascii number for new characters, try
	%	ks2 = sprintf('d, v h\t\n');
	%	ks2num = str2num(sprintf('%d\n',ks2))
	case 'NewLine',
		cfndInt = 10;
	case '\n',
		cfndInt = 10;
	case 'Tab',
		cfndInt = 9;
	case '\t',
		cfndInt = 9;
	case 'Space',
		cfndInt = 32;
	case 'Comma',
		cfndInt = 44;
	otherwise,
		if (length(cfnd) > 1),
			err = ErrEval(FuncName,'Err_Special Character not supported');
			return;
		else
			cfndInt = str2num(sprintf('%d\n',cfnd));
		end
end

%Now find such characters

	Loc = find (C == cfndInt(1));

%that's it, get back

err = 0;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function days = month_length_greek ( y, m )

%*****************************************************************************80
%
%% MONTH_LENGTH_GREEK returns the number of days in a Greek month.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year in which the month occurred.
%
%    Input, integer M, the number of the month.
%
%    Output, integer DAYS, the number of days
%    in the month.
%
  mdays = [ 30, 29, 30, 29, 30, 29, 29, 30, 29, 30, 29, 30, 29 ];
%
%  Copy the input.
%
  m2 = m;
  y2 = y;

  if ( m2 < 1 )
    days = 0;
    return
  end
%
%  A 13-month year.
%
  if ( year_is_embolismic_greek ( y2 ) )

    if ( 13 < m2 )
      days = 0;
      return
    end

    days = mdays(m2);

    if ( m2 == 7 && year_is_leap_greek ( y2 ) )
      days = days + 1;
    end
%
%  A 12 month year.
%
  else

    if ( m2 <= 6 )
      days = mdays(m2);
    elseif ( m2 <= 12 )
      days = mdays(m2+1);
    else
      days = 0;
    end

  end

  return
end


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function [c,fc] = FInvSABR2(x, a, b, r, n, f, t)
% This function computes the inverse of a SABR cumulative distribution
% Since it is based on FSABR2 it is very fast!
% x is the quantile for which we have to find the value to use
f = @(k) psabr_4(a, b, r, n, f, k, t); % handle
f2 = @(t) (x - FSABR2(t,f));           % objective function

% Use simple bisection
a = 0.001;  b = 1;                     % lower and upper bound
eps = 1e-3;         % threshold
k = 0;              % index

fc = zeros(1,10);   % initialize

while(k <= 100)
    k = k+1;
    c=0.5*(b+a);
    fc = f2(c);
    if( (-eps <= fc) && (fc <= eps))
        break;
    elseif (f2(c) < 0)
        b = c;
    else
        a = c;
    end
end

end

author: Sable

repo: mcbench-benchmarks


text:
function z = isequal(x,y)
%ISEQUAL for sptensors.
%
%   ISEQUAL(A,B) compares the sparse tensors A and B for equality.
%
%   See also SPTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

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

%% Observations for sparse matrix case.
% The result of isequal(a,full(a)) is true!

%%
if ~isequal(x.size,y.size)
    z = false;
elseif isa(x,'sptensor') && isa(y,'sptensor')
    z = (nnz(x-y) == 0);
elseif isa(y,'tensor')
    z = isequal(full(x),y);
else
    z = false;
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
classdef GlobalPooling < dagnn.Filter
  properties
    method = 'avg';
    poolSize = [1 1];
  end

  methods
    function outputs = forward(self, inputs, params)
      outputs{1} = vl_nnglobalpool(inputs{1}, 'method', self.method) ;
    end

    function [derInputs, derParams] = backward(self, inputs, params, derOutputs)
      derInputs{1} = vl_nnglobalpool(inputs{1}, derOutputs{1}, 'method', self.method) ;
      derParams = {} ;
    end

    function kernelSize = getKernelSize(obj)
      kernelSize = obj.poolSize ;
    end

    function outputSizes = getOutputSizes(obj, inputSizes)
      %outputSizes = getOutputSizes@dagnn.Filter(obj, inputSizes) ;
      %outputSizes{1}(3) = inputSizes{1}(3) ;
      poolSize = [inputSizes{1}(1)  inputSizes{1}(2)];
      outputSizes{1}=[1,1,inputSizes{1}(3),inputSizes{1}(4)];
    end

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

repo: MFT


text:
function output=SSBerouti79(signal,fs,IS)

% OUTPUT=SSBEROUTI79(S,FS,IS)
% Nonlinear Spectral Subtraction based on Berouti 79. Power spectral
% subtraction with adjusting subtraction factor. the adjustment is
% according to local a postriori SNR.
% S is the noisy signal, FS is the sampling frequency and IS is the initial
% silence (noise only) length in seconds (default value is .25 sec)
%
% Required functions:
% SEGMENT
% VAD
% Sep-04
% Esfandiar Zavarehei

if (nargin<3 | isstruct(IS))
    IS=.25; %seconds
end
W=fix(.025*fs); %Window length is 25 ms
nfft=W;
SP=.4; %Shift percentage is 40% (10ms) %Overlap-Add method works good with this value(.4)
wnd=hamming(W);

% IGNORE THIS SECTION FOR CAMPATIBALITY WITH ANOTHER PROGRAM FROM HERE.....
if (nargin>=3 & isstruct(IS))%This option is for compatibility with another programme
    W=IS.windowsize
    SP=IS.shiftsize/W;
    nfft=IS.nfft;
    wnd=IS.window;
    if isfield(IS,'IS')
        IS=IS.IS;
    else
        IS=.25;
    end
end
% .......IGNORE THIS SECTION FOR CAMPATIBALITY WITH ANOTHER PROGRAM T0 HERE

NIS=fix((IS*fs-W)/(SP*W) +1);%number of initial silence segments
Gamma=2;%Magnitude Power (1 for magnitude spectral subtraction 2 for power spectrum subtraction)
%Change Gamma to 1 to get a completely different performance
y=segment(signal,W,SP,wnd);
Y=fft(y,nfft);
YPhase=angle(Y(1:fix(end/2)+1,:)); %Noisy Speech Phase
Y=abs(Y(1:fix(end/2)+1,:)).^Gamma;%Specrogram
numberOfFrames=size(Y,2);
FreqResol=size(Y,1);

N=mean(Y(:,1:NIS)')'; %initial Noise Power Spectrum mean

NoiseCounter=0;
NoiseLength=9;%This is a smoothing factor for the noise updating

Beta=.03;
minalpha=1;
maxalpha=3;
minSNR=-5;
maxSNR=20;
alphaSlope=(minalpha-maxalpha)/(maxSNR-minSNR);
alphaShift=maxalpha-alphaSlope*minSNR;

BN=Beta*N;

for i=1:numberOfFrames
    [NoiseFlag, SpeechFlag, NoiseCounter, Dist]=vad(Y(:,i).^(1/Gamma),N.^(1/Gamma),NoiseCounter); %Magnitude Spectrum Distance VAD
    if SpeechFlag==0
        N=(NoiseLength*N+Y(:,i))/(NoiseLength+1); %Update and smooth noise
        BN=Beta*N;
    end

    SNR=10*log(Y(:,i)./N);
    alpha=alphaSlope*SNR+alphaShift;
    alpha=max(min(alpha,maxalpha),minalpha);

    D=Y(:,i)-alpha.*N; %Nonlinear (Non-uniform) Power Specrum Subtraction

    X(:,i)=max(D,BN); %if BY>D X=BY else X=D which sets very small values of subtraction result to an attenuated
                      %version of the input power spectrum.
end

output=OverlapAdd2(X.^(1/Gamma),YPhase,W,SP*W);

function ReconstructedSignal=OverlapAdd2(XNEW,yphase,windowLen,ShiftLen);

%Y=OverlapAdd(X,A,W,S);
%Y is the signal reconstructed signal from its spectrogram. X is a matrix
%with each column being the fft of a segment of signal. A is the phase
%angle of the spectrum which should have the same dimension as X. if it is
%not given the phase angle of X is used which in the case of real values is
%zero (assuming that its the magnitude). W is the window length of time
%domain segments if not given the length is assumed to be twice as long as
%fft window length. S is the shift length of the segmentation process ( for
%example in the case of non overlapping signals it is equal to W and in the
%case of %50 overlap is equal to W/2. if not givven W/2 is used. Y is the
%reconstructed time domain signal.
%Sep-04
%Esfandiar Zavarehei

if nargin<2
    yphase=angle(XNEW);
end
if nargin<3
    windowLen=size(XNEW,1)*2;
end
if nargin<4
    ShiftLen=windowLen/2;
end
if fix(ShiftLen)~=ShiftLen
    ShiftLen=fix(ShiftLen);
    disp('The shift length have to be an integer as it is the number of samples.')
    disp(['shift length is fixed to ' num2str(ShiftLen)])
end

[FreqRes FrameNum]=size(XNEW);

Spec=XNEW.*exp(j*yphase);

if mod(windowLen,2) %if FreqResol is odd
    Spec=[Spec;flipud(conj(Spec(2:end,:)))];
else
    Spec=[Spec;flipud(conj(Spec(2:end-1,:)))];
end
sig=zeros((FrameNum-1)*ShiftLen+windowLen,1);
weight=sig;
for i=1:FrameNum
    start=(i-1)*ShiftLen+1;
    spec=Spec(:,i);
    sig(start:start+windowLen-1)=sig(start:start+windowLen-1)+real(ifft(spec,windowLen));
end
ReconstructedSignal=sig;

function [NoiseFlag, SpeechFlag, NoiseCounter, Dist]=vad(signal,noise,NoiseCounter,NoiseMargin,Hangover)

%[NOISEFLAG, SPEECHFLAG, NOISECOUNTER, DIST]=vad(SIGNAL,NOISE,NOISECOUNTER,NOISEMARGIN,HANGOVER)
%Spectral Distance Voice Activity Detector
%SIGNAL is the the current frames magnitude spectrum which is to labeld as
%noise or speech, NOISE is noise magnitude spectrum template (estimation),
%NOISECOUNTER is the number of imediate previous noise frames, NOISEMARGIN
%(default 3)is the spectral distance threshold. HANGOVER ( default 8 )is
%the number of noise segments after which the SPEECHFLAG is reset (goes to
%zero). NOISEFLAG is set to one if the the segment is labeld as noise
%NOISECOUNTER returns the number of previous noise segments, this value is
%reset (to zero) whenever a speech segment is detected. DIST is the
%spectral distance.
%Saeed Vaseghi
%edited by Esfandiar Zavarehei
%Sep-04

if nargin<4
    NoiseMargin=3;
end
if nargin<5
    Hangover=8;
end
if nargin<3
    NoiseCounter=0;
end

FreqResol=length(signal);

SpectralDist= 20*(log10(signal)-log10(noise));
SpectralDist(find(SpectralDist<0))=0;

Dist=mean(SpectralDist);
if (Dist < NoiseMargin)
    NoiseFlag=1;
    NoiseCounter=NoiseCounter+1;
else
    NoiseFlag=0;
    NoiseCounter=0;
end

% Detect noise only periods and attenuate the signal
if (NoiseCounter > Hangover)
    SpeechFlag=0;
else
    SpeechFlag=1;
end

function Seg=segment(signal,W,SP,Window)

% SEGMENT chops a signal to overlapping windowed segments
% A= SEGMENT(X,W,SP,WIN) returns a matrix which its columns are segmented
% and windowed frames of the input one dimentional signal, X. W is the
% number of samples per window, default value W=256. SP is the shift
% percentage, default value SP=0.4. WIN is the window that is multiplied by
% each segment and its length should be W. the default window is hamming
% window.
% 06-Sep-04
% Esfandiar Zavarehei

if nargin<3
    SP=.4;
end
if nargin<2
    W=256;
end
if nargin<4
    Window=hamming(W);
end
Window=Window(:); %make it a column vector

L=length(signal);
SP=fix(W.*SP);
N=fix((L-W)/SP +1); %number of segments

Index=(repmat(1:W,N,1)+repmat((0:(N-1))'*SP,1,W))';
hw=repmat(Window,1,N);
Seg=signal(Index).*hw;


text:
function varargout = entropy(varargin)
%ENTROPY
%
% y = ENTROPY(x)
%
% Computes/declares concave entropy -sum(x.*log(x))
%
% Implemented as evalutation based nonlinear operator. Hence, the concavity
% of this function is exploited to perform convexity analysis and rigorous
% modelling.
%
% See also CROSSENTROPY, KULLBACKLEIBLER.

switch class(varargin{1})

    case 'double'
        x = varargin{1};
        % Safe version with defined negative values (helps fmincon when
        % outside feasible region)
        if any(x<=0)
            z = abs(x);z(z==0)=1;
            y = z.*log(z);
            y(x==0) = 0;
            y(x<0) = 3*(x(x<0)-1).^2-3;
            varargout{1} = -sum(y);
        else
            varargout{1} = -sum(x.*log(x));
        end

    case {'sdpvar','ndsdpvar'}

        varargin{1} = reshape(varargin{1},[],1);
        varargout{1} = yalmip('define',mfilename,varargin{1});

    case 'char'

        X = varargin{3};
        F = (X >= 0);

        operator = struct('convexity','concave','monotonicity','none','definiteness','none','model','callback');
        operator.range = [-inf exp(-1)*length(X)];
        operator.domain = [0 inf];
        operator.bounds = @bounds;
        operator.convexhull = @convexhull;
        operator.derivative = @derivative;
        varargout{1} = F;
        varargout{2} = operator;
        varargout{3} = X;

    otherwise
        error('SDPVAR/ENTROPY called with CHAR argument?');
end

function df = derivative(x)
x(x<=0)=eps;
df = (-1-log(x));

function [L, U] = bounds(xL,xU)

t = find(xL==0);
u = find(xL<0);
xL(t)=1;
LU = [-xL.*log(xL) -xU.*log(xU)];
LU(t,1)=0;
xL(t) = 0;
L = min(LU,[],2);
U = max(LU,[],2);
U((xL < exp(-1)) & (xU > exp(-1))) = exp(-1);
L = sum(L);
U = sum(U);

function [Ax, Ay, b] = convexhull(xL,xU)

if length(xL)==1
    xM = (xU+xL)/2;
    f1 = entropy(xL);
    f2 = entropy(xM);
    f3 = entropy(xU);
    df1 = derivative(xL);
    df2 = derivative(xM);
    df3 = derivative(xU);
    [Ax,Ay,b] = convexhullConcave(xL,xM,xU,f1,f2,f3,df1,df2,df3);

elseif length(xL)==2
    x1 = [xL(1);xL(2)];
    x2 = [xU(1);xL(2)];
    x3 = [xL(1);xU(2)];
    x4 = [xU(1);xU(2)];
    x5 = (xL+xU)/2;

    f1 = entropy(x1);
    f2 = entropy(x2);
    f3 = entropy(x3);
    f4 = entropy(x4);
    f5 = entropy(x5);

    df1 = derivative(x1);
    df2 = derivative(x2);
    df3 = derivative(x3);
    df4 = derivative(x4);
    df5 = derivative(x5);

    [Ax,Ay,b] = convexhullConcave2D(x1,f1,df1,x2,f2,df2,x3,f3,df3,x4,f4,df4,x5,f5,df5);
else
    Ax = [];
    Ay = [];
    b = [];
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% POP_NEWTIMEF - Returns estimates and plots of event-related (log) spectral
%           perturbation (ERSP) and inter-trial coherence (ITC) phenomena
%           timelocked to a set of single-channel input epochs
%
% Usage:
%   >> pop_newtimef(EEG, typeplot);          % pop_up window
%   >> pop_newtimef(EEG, typeplot, lastcom); % pop_up window
%   >> pop_newtimef(EEG, typeplot, channel); % do not pop-up window
%   >> pop_newtimef(EEG, typeproc, num, tlimits,cycles,
%                        'key1',value1,'key2',value2, ... );
%
% Graphical interface:
%   "Channel/component number" - [edit box] this is the index of the data
%              channel or the index of the component for which to plot the
%              time-frequency decomposition.
%   "Sub-epoch time limits" - [edit box] sub epochs may be extracted (note that
%              this function aims at plotting data epochs not continuous data).
%              You may select the new epoch limits in this edit box.
%   "Use n time points" - [multiple choice list] this is the number of time
%              points to use for the time-frequency decomposition. The more
%              time points, the longer the time-frequency decomposition
%              takes to compute.
%   "Frequency limits" - [edit box] these are the lower and upper
%              frequency limit of the time-frequency decomposition. Instead
%              of limits, you may also enter a sequence of frequencies. For
%              example to compute the time-frequency decomposition at all
%              frequency between 5 and 50 hertz with 1 Hz increment, enter "1:50"
%   "Use limits, padding n" - [multiple choice list] "using limits" means
%              to use the upper and lower limits in "Frequency limits" with
%              a specific padding ratio (padratio argument of newtimef).
%              The last option "use actual frequencies" forces newtimef to
%              ignore the padratio argument and use the vector of frequencies
%              given as input in the "Frequency limits" edit box.
%   "Log spaced" - [checkbox] you may check this box to compute log-spaced
%              frequencies. Note that this is only relevant if you specify
%              frequency limits (in case you specify actual frequencies,
%              this parameter is ignored).
%   "Use divisive baseline" - [multiple choice list] there are two types of
%              baseline correction, additive (the baseline is subtracted)
%              or divisive (the data is divided by the baseline values).
%              The choice is yours. There is also the option to perform
%              baseline correction in single trials. See the 'trialbase' "full"
%              option in the newtimef.m documentation for more information.
%   "No baseline" - [checkbox] check this box to compute the raw time-frequency
%              decomposition with no baseline removal.
%   "Wavelet cycles" - [edit box] specify the number of cycle at the lowest
%              and highest frequency. Instead of specifying the number of cycle
%              at the highest frequency, you may also specify a wavelet
%              "factor" (see newtimef help message). In addition, it is
%              possible to specify actual wavelet cycles for each frequency
%              by entering a sequence of numbers.
%   "Use FFT" - [checkbox] check this checkbox to use FFT instead of
%              wavelet decomposition.
%   "ERSP color limits" - [edit box] set the upper and lower limit for the
%              ERSP image.
%   "see log power" - [checkbox] the log power values (in dB) are plotted.
%              Uncheck this box to plot the absolute power values.
%   "ITC color limits" - [edit box] set the upper and lower limit for the
%              ITC image.
%   "plot ITC phase" - [checkbox] check this box plot plot (overlaid on
%              the ITC amplitude) the polarity of the ITC complex value.
%   "Bootstrap significance level" - [edit box] use this edit box to enter
%              the p-value threshold for masking both the ERSP and the ITC
%              image for significance (masked values appear as light green)
%   "FDR correct" - [checkbox] this correct the p-value for multiple comparisons
%              (across all time and frequencies) using the False Discovery
%              Rate method. See the fdr.m function for more details.
%   "Optional newtimef arguments" - [edit box] addition argument for the
%              newtimef function may be entered here in the 'key', value
%              format.
%   "Plot Event Related Spectral Power" - [checkbox] plot the ERSP image
%              showing event related spectral stimulus induced changes
%   "Plot Inter Trial Coherence" - [checkbox] plot the ITC image.
%   "Plot Curve at each frequency" - [checkbox] instead of plotting images,
%              it is also possible to display curves at each frequency.
%              This functionality is beta and might not work in all cases.
%
% Inputs:
%   INEEG    - input EEG dataset
%   typeproc - type of processing: 1 process the raw channel data
%                                  0 process the ICA component data
%   num      - component or channel number
%   tlimits  - [mintime maxtime] (ms) sub-epoch time limits to plot
%   cycles   -  > 0 --> Number of cycles in each analysis wavelet
%               = 0 --> Use FFTs (with constant window length
%                       at all frequencies)
%
% Optional inputs:
%    See the NEWTIMEF function.
%
% Outputs: Same as NEWTIMEF; no outputs are returned when a
%          window pops-up to ask for additional arguments
%
% Saving the ERSP and ITC output values:
%    Simply look up the history using the eegh function (type eegh).
%    Then copy and paste the POP_NEWTIMEF command and add output args.
%    See the NEWTIMEF function for a list of outputs. For instance,
% >> [ersp itc powbase times frequencies] = pop_newtimef( EEG, ....);
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: NEWTIMEF, EEGLAB

% Copyright (C) 2002 University of California San Diego
%
% 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.

% 01-25-02 reformated help & license -ad
% 03-08-02 add eeglab option & optimize variable sizes -ad
% 03-10-02 change newtimef call -ad
% 03-18-02 added title -ad & sm
% 04-04-02 added outputs -ad & sm

function varargout = pop_newtimef(EEG, typeproc, num, tlimits, cycles, varargin );

varargout{1} = '';
% display help if not enough arguments
% ------------------------------------
if nargin < 2
	help pop_newtimef;
	return;
end;
lastcom = [];
if nargin < 3
	popup = 1;
else
	popup = ischar(num) | isempty(num);
	if ischar(num)
		lastcom = num;
	end
end

% pop up window
% -------------
if popup
	[txt, vars] = gethelpvar('newtimef.m');
    commandchan = 'tmpEEG = get(gcbf, ''userdata''); tmpchanlocs = tmpEEG(1).chanlocs; [tmp tmpval] = pop_chansel({tmpchanlocs.labels}, ''withindex'', ''on'', ''selectionmode'', ''single''); set(findobj(gcbf, ''tag'', ''chan''), ''string'',tmpval); clear tmpEEG tmp tmpchanlocs tmpval';

    g = [1 0.3 0.6 0.4];
    g2 = [1 0.3 0.25 0.75];
	geometry = { g2 g g g g g g g [0.975 1.27] [1] [1.2 1 1.2]};
    uilist = { ...
               { 'Style', 'text', 'string', fastif(typeproc, 'Channel number', 'Component number'), 'fontweight', 'bold'  } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,3,[],'1') 'tag' 'chan'} ...
               { 'style' 'pushbutton' 'string'  '...', 'enable' fastif(~isempty(EEG(1).chanlocs) && typeproc, 'on', 'off') 'callback' commandchan } ...
               {} ...
               ...
			   { 'Style', 'text', 'string', 'Sub epoch time limits [min max] (msec)', 'fontweight', 'bold' } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,4,[],[ int2str(EEG.xmin*1000) ' ' int2str(EEG.xmax*1000) ]) 'tag' 'tlimits' } ...
               { 'Style', 'popupmenu', 'string', 'Use 50 time points|Use 100 time points|Use 150 time points|Use 200 time points|Use 300 time points|Use 400 time points' 'tag' 'ntimesout' 'value' 4} { } ...
			   ...
			   { 'Style', 'text', 'string', 'Frequency limits [min max] (Hz) or sequence', 'fontweight', 'bold' } ...
			   { 'Style', 'edit', 'string', '' 'tag' 'freqs'  } ...
               { 'Style', 'popupmenu', 'string', 'Use limits, padding 1|Use limits, padding 2|Use limits, padding 4|Use actual freqs.' 'tag' 'nfreqs' }  ...
               { 'Style', 'checkbox', 'string' 'Log spaced' 'value' 0 'tag' 'freqscale' } ...
			   ...
			   { 'Style', 'text', 'string', 'Baseline limits [min max] (msec) (0->pre-stim.)', 'fontweight', 'bold' } ...
			   { 'Style', 'edit', 'string', '0' 'tag' 'baseline' } ...
			   { 'Style', 'popupmenu',  'string', 'Use divisive baseline (DIV)|Use standard deviation (STD)|Use single trial DIV baseline|Use single trial STD baseline' 'tag' 'basenorm' } ...
               { 'Style', 'checkbox', 'string' 'No baseline' 'tag' 'nobase' } ...
               ...
               { 'Style', 'text', 'string', 'Wavelet cycles [min max/fact] or sequence', 'fontweight', 'bold' } ...
               { 'Style', 'edit', 'string', getkeyval(lastcom,5,[],'3 0.8') 'tag' 'cycle' } ...
               { 'Style', 'checkbox', 'string' 'Use FFT' 'value' 0 'tag' 'fft' } ...
               { } ...
			   ...
			   { 'Style', 'text', 'string', 'ERSP color limits [max] (min=-max)', 'fontweight', 'bold' } ...
               { 'Style', 'edit', 'string', '' 'tag' 'erspmax'} ...
               { 'Style', 'checkbox', 'string' 'see log power (set)' 'tag' 'scale' 'value' 1} {} ...
			   ...
			   { 'Style', 'text', 'string', 'ITC color limits [max]', 'fontweight', 'bold' } ...
               { 'Style', 'edit', 'string', '' 'tag' 'itcmax'} ...
               { 'Style', 'checkbox', 'string' 'plot ITC phase (set)' 'tag' 'plotphase' } {} ...
			   ...
			   { 'Style', 'text', 'string', 'Bootstrap significance level (Ex: 0.01 -> 1%)', 'fontweight', 'bold' } ...
               { 'Style', 'edit', 'string', getkeyval(lastcom,'alpha') 'tag' 'alpha'} ...
               { 'Style', 'checkbox', 'string' 'FDR correct (set)' 'tag' 'fdr' } {} ...
			   ...
			   { 'Style', 'text', 'string', 'Optional newtimef() arguments (see Help)', 'fontweight', 'bold', ...
				 'tooltipstring', 'See newtimef() help via the Help button on the right...' } ...
			   { 'Style', 'edit', 'string', '' 'tag' 'options' } ...
			   {} ...
			   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'plotersp','present',0), 'string', ...
				 'Plot Event Related Spectral Power', 'tooltipstring', ...
				 'Plot log spectral perturbation image in the upper panel' 'tag' 'plotersp' } ...
			   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'plotitc','present',0), 'string', ...
				 'Plot Inter Trial Coherence', 'tooltipstring', ...
				 'Plot the inter-trial coherence image in the lower panel' 'tag' 'plotitc' } ...
			   { 'Style', 'checkbox', 'value', 0, 'string', ...
				 'Plot curve at each frequency' 'tag' 'plotcurve' } ...
			 };
      % { 'Style', 'edit', 'string', '''padratio'', 4, ''plotphase'', ''off''' } ...
			   %{ 'Style', 'text', 'string',  '[set] -> Plot ITC phase sign', 'fontweight', 'bold', ...
			%	 'tooltipstring', ['Plot the sign (+/-) of inter-trial coherence phase' 10 ...
			%		'as red (+) or blue (-)'] } ...
			%   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'plotphase','present',1) } { } ...

	[ tmp1, tmp2, strhalt, result ] = inputgui( 'geometry', geometry, 'uilist', uilist, 'helpcom', 'pophelp(''pop_newtimef'');', ...
					   'title', fastif(typeproc, 'Plot channel time frequency -- pop_newtimef()', ...
							  'Plot component time frequency -- pop_newtimef()'), 'userdata', EEG);
	if length( tmp1 ) == 0 return; end

	if result.fft,      result.cycle = '0'; end
	if result.nobase,   result.baseline = 'NaN'; end

    num   = eeg_decodechan(EEG.chanlocs, result.chan );
	tlimits	 = eval( [ '[' result.tlimits ']' ] );
	cycles	 = eval( [ '[' result.cycle   ']' ] );
    freqs    = eval( [ '[' result.freqs   ']' ] );
    %result.ncycles == 2 is ignored

    % add topoplot
    % ------------
    options = [];
    if isfield(EEG.chanlocs, 'theta') && ~isempty(EEG.chanlocs(num).theta)
        if ~isfield(EEG, 'chaninfo'), EEG.chaninfo = []; end
        if typeproc == 1
            if isempty(EEG.chanlocs), caption = [ 'Channel ' int2str(num) ]; else caption = EEG.chanlocs(num).labels; end
            options = [options ', ''topovec'', ' int2str(num) ...
                        ', ''elocs'', EEG.chanlocs, ''chaninfo'', EEG.chaninfo, ''caption'', ''' caption '''' ];
        else
            options = [options ', ''topovec'', EEG.icawinv(:,' int2str(num) ...
                       '), ''elocs'', EEG.chanlocs, ''chaninfo'', EEG.chaninfo, ''caption'', [''IC ' num2str(num) ''']' ];
      end
    end

	if ~isempty( result.baseline ),  options = [ options ', ''baseline'',[' result.baseline ']' ]; end
    if ~isempty( result.alpha ),     options = [ options ', ''alpha'',' result.alpha ];   end
	if ~isempty( result.options ),   options = [ options ',' result.options ];            end
	if ~isempty( result.freqs ),     options = [ options ', ''freqs'', [' result.freqs ']'   ]; end
	if ~isempty( result.erspmax ),   options = [ options ', ''erspmax'', [' result.erspmax ']' ]; end
	if ~isempty( result.itcmax ),    options = [ options ', ''itcmax'','  result.itcmax ];      end
	if ~result.plotersp,             options = [ options ', ''plotersp'', ''off''' ];     end
	if ~result.plotitc,              options = [ options ', ''plotitc'' , ''off''' ];     end
	if result.plotcurve,             options = [ options ', ''plottype'', ''curve''' ];   end
	if result.fdr,                   options = [ options ', ''mcorrect'', ''fdr''' ];     end
	if result.freqscale,             options = [ options ', ''freqscale'', ''log''' ];    end
	if ~result.plotphase,            options = [ options ', ''plotphase'', ''off''' ];    end
	if ~result.scale,                options = [ options ', ''scale'', ''abs''' ];        end
    if result.ntimesout == 1,        options = [ options ', ''ntimesout'', 50' ];         end
    if result.ntimesout == 2,        options = [ options ', ''ntimesout'', 100' ];        end
    if result.ntimesout == 3,        options = [ options ', ''ntimesout'', 150' ];        end
    if result.ntimesout == 5,        options = [ options ', ''ntimesout'', 300' ];        end
    if result.ntimesout == 6,        options = [ options ', ''ntimesout'', 400' ];        end
    if result.nfreqs == 1,           options = [ options ', ''padratio'', 1' ];           end
    if result.nfreqs == 2,           options = [ options ', ''padratio'', 2' ];           end
    if result.nfreqs == 3,           options = [ options ', ''padratio'', 4' ];           end
    if result.nfreqs == 4,           options = [ options ', ''nfreqs'', ' int2str(length(freqs)) ]; end
    if result.basenorm == 2,         options = [ options ', ''basenorm'', ''on''' ];      end
    if result.basenorm == 4,         options = [ options ', ''basenorm'', ''on''' ];      end
    if result.basenorm >= 3,         options = [ options ', ''trialbase'', ''full''' ];      end

    % add title
    % ---------
	if isempty( findstr(  '''title''', result.options))
        if ~isempty(EEG.chanlocs) && typeproc
            chanlabel = EEG.chanlocs(num).labels;
        else
            chanlabel = int2str(num);
        end
	end

    % compute default winsize
    % -----------------------
    if EEG.xmin < 0 && isempty(findstr( '''winsize''', result.options)) && isempty( result.freqs )
        fprintf('Computing window size in pop_newtimef based on half of the length of the baseline period');
        options = [ options ', ''winsize'', ' int2str(-EEG.xmin*EEG.srate) ];
    end

	figure; try, icadefs; set(gcf, 'color', BACKCOLOR); catch, end
else
    options = [ ',' vararg2str(varargin) ];
end

% compute epoch limits
% --------------------
if isempty(tlimits)
	tlimits = [EEG.xmin, EEG.xmax]*1000;
end
pointrange1 = round(max((tlimits(1)/1000-EEG.xmin)*EEG.srate, 1));
pointrange2 = round(min((tlimits(2)/1000-EEG.xmin)*EEG.srate+1, EEG.pnts));
pointrange = [pointrange1:pointrange2];

% call function sample either on raw data or ICA data
% ---------------------------------------------------
if typeproc == 1
	tmpsig = EEG.data(num,pointrange,:);
else
	if ~isempty( EEG.icasphere )
        if ~isempty(EEG.icaact)
    		tmpsig = EEG.icaact(num,pointrange,:);
 	    else
            tmpsig = (EEG.icaweights(num,:)*EEG.icasphere)*reshape(EEG.data(:,pointrange,:), EEG.nbchan, EEG.trials*length(pointrange));
        end
	else
		error('You must run ICA first');
	end;
end;
tmpsig = reshape( tmpsig, length(num), size(tmpsig,2)*size(tmpsig,3));

% outputs
% -------
outstr = '';
if ~popup
    for io = 1:nargout, outstr = [outstr 'varargout{' int2str(io) '},' ]; end
    if ~isempty(outstr), outstr = [ '[' outstr(1:end-1) '] =' ]; end
end

% plot the datas and generate output command
% --------------------------------------------
if length( options ) < 2
    options = '';
end
if nargin < 4
    varargout{1} = sprintf('figure; pop_newtimef( EEG, %d, %d, [%s], [%s] %s);', typeproc, num, ...
			int2str(tlimits), num2str(cycles), options);
end
com = sprintf('%s newtimef( tmpsig(:, :), length(pointrange), [tlimits(1) tlimits(2)], EEG.srate, cycles %s);', outstr, options);
eval(com)

return;

% get contextual help
% -------------------
function txt = context(var, allvars, alltext);
	loc = strmatch( var, allvars);
	if ~isempty(loc)
		txt= alltext{loc(1)};
	else
		disp([ 'warning: variable ''' var ''' not found']);
		txt = '';
	end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function varargout = atan(varargin)

switch class(varargin{1})

    case 'sdpvar'
        varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char'

        operator = CreateBasicOperator('increasing','callback');
        operator.derivative = @(x)((1+x.^2).^-1);
        operator.inverse = @(x)(tan(x));
        operator.inflection = [-inf 1 0 -1];
        operator.range = [-pi/2 pi/2];

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error(['SDPVAR/' upper(mfilename) ' called with weird argument']);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function planC = createIMuids(planC)
% function planC = createIMuids(planC)
%
% APA, 10/21/2014

if ~exist('planC','var')
    global planC
end

indexS = planC{end};

nIM = length(planC{indexS.IM});
for i = 1:nIM
    % Copy beams, goals, solution, params fields from IMsetup to IMDosimetry
    if ~isfield(planC{indexS.IM}(i),'IMDosimetry') && isfield(planC{indexS.IM}(i),'IMSetup')
        planC{indexS.IM}(i).IMDosimetry = planC{indexS.IM}(i).IMSetup;
    else
        %do nothing asuming IMDosimerty was properly populated
    end

    % check for empty dosimetry and assign default name
    if ~isempty(planC{indexS.IM}(i).IMDosimetry)
        % Assign a default name 'IM doseSet'
        planC{indexS.IM}(i).IMDosimetry.name = 'IM doseSet';
        % Assume Dosimetry (beamlets) is fresh (1=fresh, 2=stale)
        planC{indexS.IM}(i).IMDosimetry.isFresh = 1;
    end

    % check for empty IMUID field assign default IMUID
    if isfield(planC{indexS.IM}(i),'IMUID') && isempty(planC{indexS.IM}(i).IMUID)
        planC{indexS.IM}(i).IMUID = createUID('IM');
    end

    % check for empty beamUID field assign default beamUID
    if (isfield(planC{indexS.IM}(i).IMDosimetry,'beams') && ~isempty(planC{indexS.IM}(i).IMDosimetry.beams)) ...
            && ((isfield(planC{indexS.IM}(i).IMDosimetry.beams(1),'beamUID') && isempty(planC{indexS.IM}(i).IMDosimetry.beams(1).beamUID)) || ~isfield(planC{indexS.IM}(i).IMDosimetry.beams(1),'beamUID'))
        for iBeam = 1:length(planC{indexS.IM}(i).IMDosimetry.beams)
            planC{indexS.IM}(i).IMDosimetry.beams(iBeam).beamUID = createUID('BEAM');
        end
    end

    % Link structure UID to goals
    if isfield(planC{indexS.IM}(i).IMDosimetry,'goals')
        for j = 1:length(planC{indexS.IM}(i).IMDosimetry.goals)
            structNum = planC{indexS.IM}(i).IMDosimetry.goals(j).structNum;
            if isempty(planC{indexS.IM}(i).IMDosimetry.goals(j).strUID)
                planC{indexS.IM}(i).IMDosimetry.goals(j).strUID = planC{indexS.structures}(structNum).strUID;
            end
        end

        % Link scan UID to IM
        if isempty(planC{indexS.IM}(i).IMDosimetry.assocScanUID)
            planC{indexS.IM}(i).IMDosimetry.assocScanUID = planC{indexS.scan}(1).scanUID;
        end

        % put beamlets under its parent beam
        if isfield(planC{indexS.IM}(i).IMDosimetry,'beamlets')
            structPresentV = ~cellfun('isempty',squeeze({planC{indexS.IM}(i).IMDosimetry.beamlets.beamNum}));
            structPresentM = reshape(structPresentV,size(planC{indexS.IM}(i).IMDosimetry.beamlets));
            beamletTotal = 0;
            for beamNum = 1:length(planC{indexS.IM}(i).IMDosimetry.beams)
                numBeamlets = length(planC{indexS.IM}(i).IMDosimetry.beams(beamNum).xPBPosV);
                planC{indexS.IM}(i).IMDosimetry.beams(beamNum).beamlets = planC{indexS.IM}(i).IMDosimetry.beamlets(structPresentM(:,1),beamletTotal+1:beamletTotal+numBeamlets);
                beamletTotal = beamletTotal + numBeamlets;
            end

            % assign strUID
            count = 0;
            strPresentV = find(structPresentM(:,1));
            for strNum = 1:length(strPresentV)
                count = count + 1;
                for beamNum = 1:length(planC{indexS.IM}(i).IMDosimetry.beams)
                    [planC{indexS.IM}(i).IMDosimetry.beams(beamNum).beamlets(count,:).strUID] = deal(planC{indexS.structures}(strPresentV(strNum)).strUID);
                end
            end

            % remove beamlets field from IMDosimetry since it is a part of beams now
            planC{indexS.IM}(i).IMDosimetry = rmfield(planC{indexS.IM}(i).IMDosimetry,'beamlets');

        end

        %rename solution to solutions
        if isfield(planC{indexS.IM}(i).IMDosimetry,'solution')
            planC{indexS.IM}(i).IMDosimetry.solutions =  planC{indexS.IM}(i).IMDosimetry.solution;
            planC{indexS.IM}(i).IMDosimetry = rmfield(planC{indexS.IM}(i).IMDosimetry,'solution');
        end

    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% STD_GETDATASET - Constructs and returns EEG dataset from STUDY design.
%
% Usage:
%     >> EEG = std_getdataset(STUDY, ALLEEG, 'key', 'val', ...);
%
% Inputs:
%   STUDY      - EEGLAB STUDY set
%   ALLEEG     - vector of the EEG datasets included in the STUDY structure
%
% Optional inputs:
%   'design'     - [numeric vector] STUDY design index. Default is to use
%                  the current design.
%   'cell'       - [integer] index of the STUDY design cell to convert to
%                  an EEG dataset. Default is the first cell.
%   'dataset'    - [integer] bypass cell and design above. Simply get
%                  dataset with a given index.
%   'rmcomps'    - [integer array] remove artifactual components (this entry
%                  is ignored when plotting components). This entry contains
%                  the indices of the components to be removed. Default is none.
%   'rmclust'    - [integer array] which cluster(s) to remove from the data
%                  Default is none.
%   'interp'     - [struct] channel location structure containing electrode
%                  to interpolate ((this entry is ignored when plotting
%                  components). Default is no interpolation.
%   'onecomppercluster' - ['on'|'off'] when 'on' enforces one component per
%                  cluster. Default is 'off'.
%   'interpcomponent' - ['on'|'off'] when 'on' interpolating component
%                  scalp maps. Default is 'off'.
%   'cluster'    - [integer] which cluster(s). When this option is being
%                  used, only the component contained in the selected
%                  clusters are being loaded in the dataset.
%   'checkonly'  - ['on'|'off'] use in conjunction with the option above.
%                  When 'on', no dataset is returned. Default is 'off'.
%
% Outputs:
%   EEG          - EEG dataset structure corresponding to the selected
%                  STUDY design cell (element)
%   complist     - list of components selected
%
% Example to build the dataset corresponding to the first cell of the first
% design:
%   EEG = std_getdataset(STUDY, ALLEEG, 'design', 1, 'cell', 1);
%
% Example to check that all datasets in the design have exactly one
% component per cluster for cluster 2, 3, 4 and 5.
%   std_getdataset(STUDY, ALLEEG, 'design', 1, 'cell',
%   [1:length(STUDY.design(1).cell)], 'cluster', [2 3 4 5], 'checkonly', 'on');
%
% Authors: Arnaud Delorme, SCCN, INC, UCSD, January, 2012

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, October 12, 2004, arno@sccn.ucsd.edu
%
% 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 [EEGOUT listcomp rmlistcomp] = std_getdataset(STUDY, ALLEEG, varargin);

if nargin < 2
    help std_getdataset;
    return;
end

opt = finputcheck( varargin, { 'design'  'integer'   []    STUDY.currentdesign;
                               'interp'  'struct'    { }   struct([]);
                               'rmcomps' 'integer'	 []    [];
                               'cluster' 'integer'   []    [];
                               'rmclust' 'integer'   []    [];
                               'dataset' 'integer'   []    [];
                               'onecomppercluster' 'string'   {'on' 'off'}    'off';
                               'interpcomponent'   'string'   {'on' 'off'}    'off';
                               'checkonly' 'string'   {'on' 'off'}    'off';
                               'cell'    'integer'   []    1                    }, 'std_getdataset');
%                               'mode'   'string'   { 'channels' 'components' }  'channels';
if ischar(opt), error(opt); end

if length(opt.cell) > 1
    % recursive call if more than one dataset
    % ---------------------------------------
    indcell = strmatch('cell', varargin(1:2:end));
    for index = 1:length(opt.cell)
        varargin{2*indcell} = index;
        EEGOUT(index) = std_getdataset(STUDY, ALLEEG, varargin{:});
    end
else
    if ~isempty(opt.dataset)
        mycell.dataset = opt.dataset;
        mycell.trials  = { [1:ALLEEG(opt.dataset).trials] };
    else
        mycell = STUDY.design(opt.design).cell(opt.cell);
    end

    % find components in non-artifactual cluster
    if ~isempty(opt.cluster)
        listcomp = [];
        for index = 1:length(opt.cluster)
            clsset  = STUDY.cluster(opt.cluster(index)).sets;
            clscomp = STUDY.cluster(opt.cluster(index)).comps;

            [indrow indcomp] = find(clsset == mycell.dataset(1));
            if length(indcomp) ~= 1 && strcmpi(opt.onecomppercluster, 'on')
                error(sprintf('Dataset %d must have exactly 1 components in cluster %d', mycell.dataset(1), opt.cluster(index)));
            end
            listcomp = [listcomp clscomp(indcomp')];
        end
    end

    % find components in artifactual clusters
    if ~isempty(opt.rmclust)
        rmlistcomp = [];
        for index = 1:length(opt.rmclust)
            clsset  = STUDY.cluster(opt.rmclust(index)).sets;
            clscomp = STUDY.cluster(opt.rmclust(index)).comps;
            [indrow indcomp] = find(clsset == mycell.dataset(1));
            rmlistcomp = [rmlistcomp clscomp(indcomp')];
        end
        if ~isempty(opt.rmcomps)
            disp('Both ''rmclust'' and ''rmcomps'' are being set. Artifact components will be merged');
        end
        opt.rmcomps = [ opt.rmcomps rmlistcomp ];
    end;
    rmlistcomp  = opt.rmcomps;
    if strcmpi(opt.checkonly, 'on'), EEGOUT = 0; return; end

    % get data
    EEG = ALLEEG(mycell.dataset);
    EEGOUT = EEG(1);
    EEGOUT.data = eeg_getdatact(EEG, 'channel', [], 'trialindices', mycell.trials, 'rmcomps', opt.rmcomps, 'interp', opt.interp);
    EEGOUT.trials   = size(EEGOUT.data,3);
    EEGOUT.nbchan   = size(EEGOUT.data,1);
    if ~isempty(opt.interp)
        EEGOUT.chanlocs = opt.interp;
    end
    EEGOUT.event  = [];
    EEGOUT.epoch  = [];
    if isfield(mycell, 'filebase')
        EEGOUT.filename = mycell.filebase;
        EEGOUT.condition = mycell.value{1};
        EEGOUT.group     = mycell.value{2};
        EEGOUT.subject   = mycell.case;
    end
    if ~isempty(opt.cluster)
        if ~isempty(opt.interp) && strcmpi(opt.interpcomponent, 'on')
            TMPEEG = EEGOUT;
            TMPEEG.chanlocs = EEG(1).chanlocs;
            TMPEEG.data = EEG(1).icawinv;
            TMPEEG.nbchan = size(TMPEEG.data,1);
            TMPEEG.pnts   = size(TMPEEG.data,2);
            TMPEEG.trials = 1;
            TMPEEG = eeg_interp(TMPEEG, opt.interp, 'spherical');
            EEGOUT.icawinv = TMPEEG.data(:, listcomp);
        else
            EEGOUT.icawinv = EEGOUT.icawinv(:, listcomp);
        end
        EEGOUT.icaact         = eeg_getdatact(EEG, 'component', listcomp, 'trialindices', mycell.trials );
        EEGOUT.icaweights     = EEGOUT.icaweights(listcomp,:);
        EEGOUT.etc.clustid    = { STUDY.cluster(opt.cluster).name };  % name of each cluster
        EEGOUT.etc.clustcmpid = listcomp;                             % index of each component in original ICA matrix
    else
        EEGOUT = eeg_checkset(EEGOUT);
    end
end


text:
function [yn,en,S] = FXNLMSadapt(un,dn,S)

% FXNLMSadapt       Normalized FXLMS Algorithm
%
%                   Perform over the entire length of the input sequence. The
%                   history of output, square error and coefficients of FIR
%                   filters are passed out to extenal
% Arguments:
% un                Input signal
% dn                Desired signal
% S                 Adptive filter parameters as defined in NLMSinit.m
% yn                History of output signal
% en                History of error signal
%
% by Lee, Gan, and Kuo, 2008
% Subband Adaptive Filtering: Theory and Implementation
% Publisher: John Wiley and Sons, Ltd

M = length(S.coeffs);                     % Length of input sequence
mu = S.step;                              % Step size of NLMS algorithm (between 0 and 2)
leak = S.leakage;                         % Leaky factor for leaky LMS algorithm
alpha = S.alpha;                          % Small constant
AdaptStart = S.AdaptStart;
w = S.coeffs;                             % Coefficients of FIR filter
u = zeros(M,1);
fu = zeros(M,1);

ITER = length(un);                        % Length of input sequence
yn = zeros(1,ITER);                       % Initialize output sequence to zero
en = zeros(1,ITER);                       % Initialize error sequence to zero
%filt_un = zeros(1,ITER);                 % Initialize filtered input to zero
%est_filt_un = zeros(1,ITER);             % Initialize filtered input to zero (est)
filt_un = filter(S.sec_num,S.sec_den,un); % Filtered input
est_filt_un = filter(S.estsec,1,un);      % Filtered input (est)

if isfield(S,'unknownsys')
    b = S.unknownsys;
    norm_b = norm(b);
    eml = zeros(1,ITER);
    ComputeEML = 1;
else
    ComputeEML = 0;
end

for n = 1:ITER
   u = [filt_un(n); u(1:end-1)];          % Input signal vector (through the actual
                                          %   secondary path)
   fu = [est_filt_un(n); fu(1:end-1)];    % Filtered input signal vector (through
                                          %   the estimated secondary path)
   yn(n) = w'*u;                          % Inner product of filter coefficients and
                                          %   tappd delay line
   en(n) = dn(n)-yn(n);                   % Compute error signal
   if ComputeEML == 1;
        eml(n) = norm(b-w)/norm_b;        % System error norm (normalized)
   end
   if n >= AdaptStart
        w = (1-mu*leak)*w + ((mu*en(n))/(fu'*fu+alpha))*fu;
                                          % LMS algorithm in leaky mode
        S.iter = S.iter + 1;
   end
end

S.coeffs = w;                             % Coefficient values at the final iteration
if ComputeEML == 1;
    S.eml = eml;
end

author: CharlesThaCat

repo: acoustic-interference-cancellation


text:
classdef BriefDescriptorExtractor < handle
    %BRIEFDESCRIPTOREXTRACTOR  Class for computing BRIEF descriptors
    %
    % As described in [calon2010].
    %
    % ## References
    % [calon2010]:
    % > Michael Calonder, Vincent Lepetit, Christoph Strecha, and Pascal Fua.
    % > "Brief: Binary robust independent elementary features".
    % > In Computer Vision-ECCV 2010, pages 778-792. Springer, 2010.
    %
    % See also: cv.BriefDescriptorExtractor.BriefDescriptorExtractor,
    %  cv.DescriptorExtractor
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    methods
        function this = BriefDescriptorExtractor(varargin)
            %BRIEFDESCRIPTOREXTRACTOR  Constructor
            %
            %     obj = cv.BriefDescriptorExtractor()
            %     obj = cv.BriefDescriptorExtractor(..., 'OptionName',optionValue, ...)
            %
            % ## Options
            % * __Bytes__ legth of the descriptor in bytes, valid values are:
            %   16, 32 (default) or 64.
            % * __UseOrientation__ sample patterns using keypoints orientation,
            %   disabled by default. default false.
            %
            % See also: cv.BriefDescriptorExtractor.compute
            %
            this.id = BriefDescriptorExtractor_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.BriefDescriptorExtractor
            %
            if isempty(this.id), return; end
            BriefDescriptorExtractor_(this.id, 'delete');
        end

        function typename = typeid(this)
            %TYPEID  Name of the C++ type (RTTI)
            %
            %     typename = obj.typeid()
            %
            % ## Output
            % * __typename__ Name of C++ type
            %
            typename = BriefDescriptorExtractor_(this.id, 'typeid');
        end
    end

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.BriefDescriptorExtractor.empty,
            %  cv.BriefDescriptorExtractor.load
            %
            BriefDescriptorExtractor_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Checks if detector object is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the detector object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.BriefDescriptorExtractor.clear,
            %  cv.BriefDescriptorExtractor.load
            %
            b = BriefDescriptorExtractor_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.BriefDescriptorExtractor.load
            %
            BriefDescriptorExtractor_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.BriefDescriptorExtractor.save
            %
            BriefDescriptorExtractor_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.BriefDescriptorExtractor.save,
            %  cv.BriefDescriptorExtractor.load
            %
            name = BriefDescriptorExtractor_(this.id, 'getDefaultName');
        end
    end

    %% Features2D: DescriptorExtractor
    methods
        function ntype = defaultNorm(this)
            %DEFAULTNORM  Returns the default norm type
            %
            %     ntype = obj.defaultNorm()
            %
            % ## Output
            % * __ntype__ Norm type. One of `cv::NormTypes`:
            %   * __Inf__
            %   * __L1__
            %   * __L2__
            %   * __L2Sqr__
            %   * __Hamming__
            %   * __Hamming2__
            %
            % Always `Hamming` for BRIEF.
            %
            % See also: cv.BriefDescriptorExtractor.compute,
            %  cv.DescriptorMatcher
            %
            ntype = BriefDescriptorExtractor_(this.id, 'defaultNorm');
        end

        function sz = descriptorSize(this)
            %DESCRIPTORSIZE  Returns the descriptor size in bytes
            %
            %     sz = obj.descriptorSize()
            %
            % ## Output
            % * __sz__ Descriptor size, as specified in `Bytes` argument in
            %   constructor.
            %
            % See also: cv.BriefDescriptorExtractor.descriptorType,
            %  cv.BriefDescriptorExtractor.compute
            %
            sz = BriefDescriptorExtractor_(this.id, 'descriptorSize');
        end

        function dtype = descriptorType(this)
            %DESCRIPTORTYPE  Returns the descriptor type
            %
            %     dtype = obj.descriptorType()
            %
            % ## Output
            % * __dtype__ Descriptor type, one of numeric MATLAB class names.
            %
            % Always `uint8` for BRIEF.
            %
            % See also: cv.BriefDescriptorExtractor.descriptorSize,
            %  cv.BriefDescriptorExtractor.compute
            %
            dtype = BriefDescriptorExtractor_(this.id, 'descriptorType');
        end

        function [descriptors, keypoints] = compute(this, img, keypoints)
            %COMPUTE  Computes the descriptors for a set of keypoints detected in an image or image set
            %
            %     [descriptors, keypoints] = obj.compute(img, keypoints)
            %     [descriptors, keypoints] = obj.compute(imgs, keypoints)
            %
            % ## Input
            % * __img__ Image (first variant), 8-bit grayscale image.
            % * __imgs__ Image set (second variant), cell array of images.
            % * __keypoints__ Input collection of keypoints. Keypoints for
            %   which a descriptor cannot be computed are removed. Sometimes
            %   new keypoints can be added, for example: cv.SIFT duplicates
            %   keypoint with several dominant orientations (for each
            %   orientation). In the first variant, this is a struct-array of
            %   detected keypoints. In the second variant, it is a cell-array,
            %   where `keypoints{i}` is a set of keypoints detected in
            %   `imgs{i}` (a struct-array like before).
            %
            % ## Output
            % * __descriptors__ Computed descriptors. In the second variant of
            %   the method `descriptors{i}` are descriptors computed for a
            %   `keypoints{i}`. Row `j` in `descriptors` (or `descriptors{i}`)
            %   is the descriptor for `j`-th keypoint.
            % * __keypoints__ Optional output with possibly updated keypoints.
            %
            % See also: cv.BriefDescriptorExtractor.BriefDescriptorExtractor
            %
            [descriptors, keypoints] = BriefDescriptorExtractor_(this.id, 'compute', img, keypoints);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [T]=affineTransformationMatrixDirect(V1,V2)

%%
%Force input to 3D
if size(V1,2)==2
    V1(:,3)=0;
end

if size(V2,2)==2
    V2(:,3)=0;
end

%Expand to nx4
V1_M=V1;
V1_M(:,4)=1;
V2_M=V2;
V2_M(:,4)=1;

%Get transformation using left devide
T=(V1_M\V2_M)';

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function[upa,upb,upc,upd,upe,upf]=ellband(varargin)
%ELLBAND  Bandwidth of modulated elliptical signals in two or three dimensions.
%
%   [A,B,C]=ELLBAND(KAPPA,LAMBDA,THETA,PHI) computes the instantaneous
%   bandwidth of the elliptical signal characterized by RMS amplitude
%   KAPPA, linearity LAMBDA, orientation THETA, and orbital phase PHI.
%
%   The three output arguments are
%
%          A  - Amplitude modulation bandwidth
%          B  - Deformation bandwidth
%          C  - Precession bandwidth.
%
%   and these satisfy UPSILON^2=A^2+B^2+C^2 where UPSILON is the joint
%   instantaneous bandwidth of the bivariate signal.
%
%   The form of these terms is as follows:
%
%         A = 1/KAPPA d/dt KAPPA
%         B = 1/2 * 1/SQRT(1-LAMBDA^2) *  d/dt LAMBDA
%         C = LAMBDA d/dt THETA
%
%   [A,B,C,UPSILON]=ELLBAND(KAPPA,LAMBDA,THETA,PHI) also returns the total
%   instantaneous bandwith UPSILON=SQRT(A^2+B^2+C^2).
%
%   ELLBAND(...,DIM) performs the analysis with time running along
%   dimension DIM, as opposed to the default behavior of DIM=1.
%
%   For details see Lilly and Olhede (2010).
%
%   ELLBAND also works if the input arguments are cell arrays of numerical
%   arrays, in which case the output will be similarly sized cell arrays.
%   __________________________________________________________________
%
%   Three dimensions
%
%   ELLBAND can also compute the instantaneous bandwidth of modulated
%   elliptical signals in three dimensions.
%
%   [A,B,C,D,E]=ELLBAND(KAPPA,LAMBDA,THETA,PHI,ALPHA,BETA) returns the
%   terms in the bandwidth from a modulated ellipical signal in a plane
%   with a normal vector having azimuth angle ALPHA and zenith angle BETA.
%
%   The five output arguments are
%
%          A   - Amplitude modulation bandwidth, as in 2D
%          B   - Deformation bandwidth, as in 2D
%          C   - Precession bandwidth, as in 2D
%          D   - Precession bandwidth with full 3D effects
%          E   - Bandwidth due to motion of the normal to the plane
%
%   and these, in principle, satisfy UPSILON^2=A^2+B^2+C^2+D^2+|E|^2 where
%   UPSILON is the joint instantaneous bandwidth of the trivariate signal.
%   See below for a caveat on this statement.
%
%   Terms A--C are just as in the bivariate case.  The new terms are:
%
%         D = LAMBDA [d/dt THETA + COS(BETA) * d/dt ALPHA]
%         E = N^T X_+ / |X_+^H X_+|
%
%   where N is the trivariate normal vector, X_+ is the trivariate analytic
%   signal vector, and "T" denotes the matrix transpose, and "H" the
%   Hermitian transpose.  Note that term E may be complex-valued.
%
%   Note that term C does not contribute to the full bandwidth, but is
%   output in order to compare the two-dimensional and three-dimensional
%   effects in the full precession bandwidth, term D.
%
%   An important point is that the trivariate ellipse parameters can be
%   ill-defined for a nearly linear signal, and the elliptical bandwidth
%   terms can give erroneously large values at isolated points.  To check
%   for this, compare with the joint bandwidth from INSTMOM.
%
%   [A,B,C,D,E,UPSILON]=ELLBAND(KAPPA,LAMBDA,THETA,PHI,ALPHA,BETA) also
%   returns the total bandwith UPSILON=SQRT(A^2+B^2+C^2+D^2+|E|^2).
%
%   For details see Lilly (2011).
%   __________________________________________________________________
%
%   ELLBAND(DT,...) sets the sample interval DT, which defaults to DT=1.
%   DT may be a scalar, or if the input fields are cell arrays having
%   length N, DT may be a numerical array of length N.
%
%   See also ANATRANS, WAVETRANS, INSTMOM.
%
%   'ellband --t' runs a test.
%   'ellband --f' generates a figure from Lilly and Olhede (2010).
%
%   Usage:  [a,b,c]=ellband(kappa,lambda,theta,phi);
%           [a,b,c]=ellband(dt,kappa,lambda,theta,phi);
%           [a,b,c]=ellband(dt,kappa,lambda,theta,phi,dim);
%           [a,b,c,upsilon]=ellband(dt,kappa,lambda,theta,phi,dim);
%           [a,b,c,d,e]=ellband(kappa,lambda,theta,phi,alpha,beta);
%           [a,b,c,d,e]=ellband(dt,kappa,lambda,theta,phi,alpha,beta);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2006--2020 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--t')
    ellband_test,return
end

if strcmpi(varargin{1}, '--f')
    type makefigs_ellband
    makefigs_ellband;
    return
end

if ~iscell(varargin{end})&&length(varargin{end})==1
    dim=varargin{end};
    varargin=varargin(1:end-1);
else
    dim=1;
end

if length(varargin)==5||length(varargin)==7
    dt=varargin{1};
    varargin=varargin(2:end);
else
    dt=1;
end

kappa=varargin{1};
lambda=varargin{2};
theta=varargin{3};
phi=varargin{4};
if length(varargin)>4
    alpha=varargin{5};
    beta=varargin{6};
    dims=3;
else
    alpha=[];
    beta=[];
    dims=2;
end

%Redo some initializing in case of cell arrays
if iscell(kappa)
    if isempty(alpha)
        for i=1:length(kappa)
            alpha{i,1}=[];
            beta{i,1}=[];
        end
    end
    if length(dt)==1
        dto=dt;
        for i=1:length(kappa)
            dt(i)=dto;
        end
    end
end

if ~isempty(kappa)
    if iscell(kappa)
        [upa,upb,upc,upd,upe,upf]=vempty;
        for i=1:length(kappa)
            [upa{i,1},upb{i,1},upc{i,1},upd{i,1},upe{i,1},upf{i,1}]=...
                ellband_one(dt(i),kappa{i,1},lambda{i,1},theta{i,1},phi{i,1},alpha{i,1},beta{i,1},dim,dims);
        end
    else
        [upa,upb,upc,upd,upe,upf]=ellband_one(dt,kappa,lambda,theta,phi,alpha,beta,dim,dims);
    end
else
    upa=kappa;upb=kappa;upc=kappa;
    upd=kappa;upe=kappa;upf=kappa;
end

function[upa,upb,upc,upd,upe,upf]=ellband_one(dt,kappa,lambda,theta,phi,alpha,beta,dim,dims)

omtheta=frac(vdiff(unwrap(theta,[],dim),dim),dt);
zeta=sign(lambda).*sqrt(1-lambda.^2);

%These do not seem much different
%upa=frac((vdiff(kappa,dim)),kappa*dt);
upa=frac((vdiff(log(kappa),dim)),dt);
upb=frac(sqrt(frac(1,1-lambda.^2)).*(frac(1,2)*vdiff(lambda,dim)),dt);

upb2=frac(sqrt(frac(1,1-zeta.^2)).*(frac(1,2)*vdiff(zeta,dim)),dt);
upb(abs(lambda)>sqrt(1/2))=-upb2(abs(lambda)>sqrt(1/2));
%To correct numerical instability when lambda is close to unity

%This is the same
%[a,b]=kl2ab(kappa,lambda);
%upc=abs(frac(a.*b,a.^2+b.^2).*vdiff(log(abs(b./a)),1))./dt;
upc=(lambda.*omtheta);

%Replace isolated nans that may have been lost
bool=isnan(kappa);
upa(bool)=nan;
upb(bool)=nan;
upc(bool)=nan;

if dims==3
     omalpha=frac(vdiff(unwrap(alpha,[],dim),dim),dt);
     ombeta=frac(vdiff(unwrap(beta,[],dim),dim),dt);

     upd=lambda.*(omtheta+omalpha.*cos(beta));

     [x,y,z]=ellsig(kappa,lambda,theta,phi,alpha,beta);

     %Now find x-tilde, 2-d vector
     [x,y,z]=vectmult(jmat3(-alpha,3),x,y,z);
     [x,y,z]=vectmult(jmat3(-beta,1),x,y,z);

     numer=-omalpha.*sin(beta).*x+ombeta.*y;
     denom=sqrt(abs(x).^2+abs(y).^2);
     upe=frac(numer,denom);
     upf=sqrt(upa.^2+upb.^2+upc.^2+upd.^2+abs(upe).^2);

     %Replace isolated nans that may have been lost
     upd(bool)=nan;
     upe(bool)=nan;
     upf(bool)=nan;

else
    upd=sqrt(upa.^2+upb.^2+upc.^2);
    upe=[];
    upf=[];
end

function[]=ellband_test
ellband_test1;
ellband_test2;
ellband_test3;
ellband_test4;

function[]=ellband_test1

load npg2006
use npg2006

%Decide on frequencies
fs=2*pi./(logspace(log10(10),log10(100),50)');

%Compute wavelet transforms using generalized Morse wavelets
[wx,wy]=wavetrans(real(cx),imag(cx),{1,2,4,fs,'bandpass'},'mirror');
[wxr,wyr,ir,jr]=ridgewalk(dt,wx,wy,fs,sqrt(2*4),1);

[kappa,lambda,theta,phi]=ellparams(wxr,wyr);
[ba,bd,bp]=ellband(kappa,lambda,theta,phi);
[a,om,upbar]=instmom([wxr wyr],1,2);

wxr2=permute(wxr,[3 2 1]);
wyr2=permute(wyr,[3 2 1]);
clear wr2
wr2(1,1,:)=wxr2;wr2(1,2,:)=wyr2;
[kappa2,lambda2,theta2,phi2]=ellparams(wxr2,wyr2,3);
[ba2,bd2,bp2]=ellband(kappa2,lambda2,theta2,phi2,3);
[a2,om2,upbar2]=instmom(wr2,3,2);
vcolon(kappa2,lambda2,theta2,phi2,ba2,bd2,bp2,a2,om2,upbar2);
bool=aresame([kappa2 lambda2 theta2 phi2 ba2 bd2 bp2 a2 om2 upbar2],...
    [kappa lambda theta phi ba bd bp a om upbar]);

reporttest('ELLBAND time running in pages', bool);

vindex(upbar,ba,bd,bp,2:length(upbar)-1,1);
err=vsum((upbar-sqrt(ba.^2+bd.^2+bp.^2)).^2,1)./vsum(upbar.^2,1);
reporttest('ELLBAND terms sum to bandwidth from INSTMOM', err<1e-5);

function[]=ellband_test2
t=(0:1:925)';
cxe=zeros(length(t),3);

kappa=3*exp(2*0.393*(t/1000-1));
lambda=0.4+0*t;
phi=(t/1000*5)*2*pi;
theta=pi/4+0*t;

om=vdiff(phi,1);  %Since theta is constant

[x,y]=ellsig(kappa,lambda,theta,phi);
[a,om,up]=instmom([x y],1,2);

ups=sqrt(vsum(abs(up).^2,2));
[upsa,upsb,upsc]=ellband(kappa,lambda,theta,phi);
b1=aresame(vmean([ups upsa upsb upsc]./[om om om om],1),[1 1 0 0]*0.025,1e-4);

om=vdiff(phi,1);  %Since theta is constant

kappa=2.5+0*t;
lambda=zeros(size(t));
for i=2:length(t)
    lambda(i)=real(lambda(i-1)+2*sqrt(1-lambda(i-1).^2)*0.025.*om(i));
end
lambda(1)=nan;
lambda(lambda>1)=1;

[x,y]=ellsig(kappa,lambda,theta,phi);
[a,om,ups]=instmom([x y],1,2);

[upsa,upsb,upsc]=ellband(kappa,lambda,theta,phi);
b2=aresame(vmean([ups upsa upsb upsc]./[om om om om],1),[1 0 1 0]*0.025,1e-4);

[kappa,lambda]=ab2kl(3+zeros(size(t)),2+zeros(size(t)));

theta=phi/14.45;

[x,y]=ellsig(kappa,lambda,theta,phi);
[a,om,ups]=instmom([x y],1,2);

[upsa,upsb,upsc]=ellband(kappa,lambda,theta,phi);
b3=aresame(vmean([ups upsa upsb upsc]./[om om om om],1),[1 0 0 1]*0.025,1e-4);

reporttest('ELLBAND three panels of figure each have bandwidth 0.025',b1&&b2&&b3)

function[]=ellband_test3

kappao=10;
lambdao=2/3;
omtheta=linspace(-2,2,100);
omphi=1;

dt=0.1;
t=[0:dt:100]';

zo=ellsig(kappao,lambdao,0,omphi.*t);
psi=sleptap(length(zo));

z=vzeros(length(zo),length(omtheta));
for i=1:length(omtheta)
    z(:,i)=rot(omtheta(i).*t).*zo;
end

%Taper for better computation
z=z.*vrep(psi(:,1),size(z,2),2);
[zp,zn]=anatrans(z,conj(z));

[ap,omp,upp]=instmom(dt,zp);
[an,omn,upn]=instmom(dt,zn);

kappa=sqrt(abs(zp).^2+abs(zn).^2);
om=frac(omp.*abs(zp).^2+omn.*abs(zn).^2,kappa.^2);
ombar=vmean(om,1,squared(kappa));

sig=sqrt(frac((upp.^2+(omp-vrep(ombar,size(om,1),1)).^2).*abs(zp).^2+...
    (upn.^2+(omn-vrep(ombar,size(om,1),1)).^2).*abs(zn).^2,kappa.^2));
sigbar=sqrt(vmean(sig.^2,1,squared(kappa)));

[x,y]=vectmult(sqrt(2)*tmat',zp,zn);
[kappa,lambda,theta,phi]=ellparams(x,y);
[ba,bd,bp]=ellband(dt,kappa,lambda,theta,phi);
bom=om-vrep(ombar,size(om,1),1);

sig2=sqrt(ba.^2+bd.^2+bp.^2+bom.^2);
sigbar2=sqrt(vmean(sig2.^2,1,squared(kappa)));

%plot(sig,'b'),hold on,plot(sig2,'r')  %Visually identical
err=sum(squared(sigbar-sigbar2))./sum(squared(sigbar));
reporttest('ELLBAND rotary and elliptical bandwidths match for shifted ellipse, non-unit sample rate',err<1e-5)

function[]=ellband_test4

load solomon
use solomon

%Choose central portion where polarization is elliptical
vindex(x,y,z,420:580,1);
vfilt(x,y,z,10,'zeros');
[x,y,z]=anatrans(x,y,z,'mirror');

[kappa,lambda,theta,phi,alpha,beta]=ellparams(x,y,z);
[a,ombar,upbar]=instmom([x,y,z],1,2);
[a,b,c,d,e]=ellband(kappa,lambda,theta,phi,alpha,beta);
upbar1=sqrt(a.^2+b.^2+d.^2+abs(e).^2);

err=abs(upbar-upbar1).^2./abs(upbar).^2;
err=sort(err,'descend');
err=err(20:end);
%figure, plot([upbar upbar1])
reporttest('ELLBAND trivariate case, sum of bandwith terms for Solomon Islands (removing worst outliers)',allall(err<0.05))

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
% RANSACFITFUNDMATRIX - fits fundamental matrix using RANSAC
%
% Usage:   [F, inliers] = ransacfitfundmatrix(x1, x2, t)
%
% Arguments:
%          x1  - 2xN or 3xN set of homogeneous points.  If the data is
%                2xN it is assumed the homogeneous scale factor is 1.
%          x2  - 2xN or 3xN set of homogeneous points such that x1<->x2.
%          t   - The distance threshold between data point and the model
%                used to decide whether a point is an inlier or not.
%                Note that point coordinates are normalised to that their
%                mean distance from the origin is sqrt(2).  The value of
%                t should be set relative to this, say in the range
%                0.001 - 0.01
%
% Note that it is assumed that the matching of x1 and x2 are putative and it
% is expected that a percentage of matches will be wrong.
%
% Returns:
%          F       - The 3x3 fundamental matrix such that x2'Fx1 = 0.
%          inliers - An array of indices of the elements of x1, x2 that were
%                    the inliers for the best model.
%
% See Also: RANSAC, FUNDMATRIX

% Copyright (c) 2004-2005 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.

% February 2004  Original version
% August   2005  Distance error function changed to match changes in RANSAC

function [F, inliers] = ransacfitfundmatrix(x1, x2, t, feedback)

    if ~all(size(x1)==size(x2))
        error('Data sets x1 and x2 must have the same dimension');
    end

    if nargin == 3
	feedback = 0;
    end

    [rows,npts] = size(x1);
    if rows~=2 & rows~=3
        error('x1 and x2 must have 2 or 3 rows');
    end

    if rows == 2    % Pad data with homogeneous scale factor of 1
        x1 = [x1; ones(1,npts)];
        x2 = [x2; ones(1,npts)];
    end

    % Normalise each set of points so that the origin is at centroid and
    % mean distance from origin is sqrt(2).  normalise2dpts also ensures the
    % scale parameter is 1.  Note that 'fundmatrix' will also call
    % 'normalise2dpts' but the code in 'ransac' that calls the distance
    % function will not - so it is best that we normalise beforehand.
    [x1, T1] = normalise2dpts(x1);
    [x2, T2] = normalise2dpts(x2);

    s = 8;  % Number of points needed to fit a fundamental matrix. Note that
            % only 7 are needed but the function 'fundmatrix' only
            % implements the 8-point solution.

    fittingfn = @fundmatrix;
    distfn    = @funddist;
    degenfn   = @isdegenerate;
    % x1 and x2 are 'stacked' to create a 6xN array for ransac
    [F, inliers] = ransac([x1; x2], fittingfn, distfn, degenfn, s, t, feedback);

    % Now do a final least squares fit on the data points considered to
    % be inliers.
    F = fundmatrix(x1(:,inliers), x2(:,inliers));

    % Denormalise
    F = T2'*F*T1;

%--------------------------------------------------------------------------
% Function to evaluate the first order approximation of the geometric error
% (Sampson distance) of the fit of a fundamental matrix with respect to a
% set of matched points as needed by RANSAC.  See: Hartley and Zisserman,
% 'Multiple View Geometry in Computer Vision', page 270.
%
% Note that this code allows for F being a cell array of fundamental matrices of
% which we have to pick the best one. (A 7 point solution can return up to 3
% solutions)

function [bestInliers, bestF] = funddist(F, x, t);

    x1 = x(1:3,:);    % Extract x1 and x2 from x
    x2 = x(4:6,:);

    if iscell(F)  % We have several solutions each of which must be tested

	nF = length(F);   % Number of solutions to test
	bestF = F{1};     % Initial allocation of best solution
	ninliers = 0;     % Number of inliers

	for k = 1:nF
	    x2tFx1 = zeros(1,length(x1));
	    for n = 1:length(x1)
		x2tFx1(n) = x2(:,n)'*F{k}*x1(:,n);
	    end

	    Fx1 = F{k}*x1;
	    Ftx2 = F{k}'*x2;

	    % Evaluate distances
	    d =  x2tFx1.^2 ./ ...
		 (Fx1(1,:).^2 + Fx1(2,:).^2 + Ftx2(1,:).^2 + Ftx2(2,:).^2);

	    inliers = find(abs(d) < t);     % Indices of inlying points

	    if length(inliers) > ninliers   % Record best solution
		ninliers = length(inliers);
		bestF = F{k};
		bestInliers = inliers;
	    end
	end

    else     % We just have one solution
	x2tFx1 = zeros(1,length(x1));
	for n = 1:length(x1)
	    x2tFx1(n) = x2(:,n)'*F*x1(:,n);
	end

	Fx1 = F*x1;
	Ftx2 = F'*x2;

	% Evaluate distances
	d =  x2tFx1.^2 ./ ...
	     (Fx1(1,:).^2 + Fx1(2,:).^2 + Ftx2(1,:).^2 + Ftx2(2,:).^2);

	bestInliers = find(abs(d) < t);     % Indices of inlying points
	bestF = F;                          % Copy F directly to bestF

    end

%----------------------------------------------------------------------
% (Degenerate!) function to determine if a set of matched points will result
% in a degeneracy in the calculation of a fundamental matrix as needed by
% RANSAC.  This function assumes this cannot happen...

function r = isdegenerate(x)
    r = 0;

author: jianxiongxiao

repo: ProfXkit


text:
classdef ZCR < Algorithm

    methods (Access = public)

        function obj = ZCR()
            obj.name = 'ZCR';
            obj.inputPort = DataType.kSignal;
            obj.outputPort = DataType.kFeature;
        end

        function result = compute(~,signal)
            N = length(signal);
            result = sum(abs(diff(signal>0))) / N;
        end

        function metrics = computeMetrics(~,input)
            n = size(input,1);
            flops = 5 * n;
            memory = 1;
            outputSize = Constants.kFeatureBytes;
            metrics = Metric(flops,memory,outputSize);
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function isargdir(varargin)
%ISARGDIR throws an error if not a directory
%
%   Usage: isargdirectory(arg1,arg2,...)
%
%   Input parameters:
%       args        - variable number of args
%
%   See also: isargfile

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking for directory ==========================================
for ii = 1:nargin
    if ~ischar(varargin{ii}) || ~exist(varargin{ii},'dir')
        error('%s need to be a directory.',inputname(ii));
    end
end

author: sfstoolbox

repo: sfs-matlab


text:
function [goal]=Target_Allocation(goal,quad_init_x,quad_init_y,num,target_num,goal_series,temp_goal,mirror_dis)
mirror_count=0;
row=0;
line=0;
temp_sum_dis=0;
sum_dis=[0;0;0;0];
limited_init_distance=0;%?????????
quad_init=[quad_init_x quad_init_y];

distance=[sqrt((quad_init(1,1)-goal(1,1))^2+(quad_init(1,2)-goal(1,2))^2) sqrt((quad_init(1,1)-goal(2,1))^2+(quad_init(1,2)-goal(2,2))^2) sqrt((quad_init(1,1)-goal(3,1))^2+(quad_init(1,2)-goal(3,2))^2) sqrt((quad_init(1,1)-goal(4,1))^2+(quad_init(1,2)-goal(4,2))^2);
          sqrt((quad_init(2,1)-goal(1,1))^2+(quad_init(2,2)-goal(1,2))^2) sqrt((quad_init(2,1)-goal(2,1))^2+(quad_init(2,2)-goal(2,2))^2) sqrt((quad_init(2,1)-goal(3,1))^2+(quad_init(2,2)-goal(3,2))^2) sqrt((quad_init(2,1)-goal(4,1))^2+(quad_init(2,2)-goal(4,2))^2);
          sqrt((quad_init(3,1)-goal(1,1))^2+(quad_init(3,2)-goal(1,2))^2) sqrt((quad_init(3,1)-goal(2,1))^2+(quad_init(3,2)-goal(2,2))^2) sqrt((quad_init(3,1)-goal(3,1))^2+(quad_init(3,2)-goal(3,2))^2) sqrt((quad_init(3,1)-goal(4,1))^2+(quad_init(3,2)-goal(4,2))^2);
          sqrt((quad_init(4,1)-goal(1,1))^2+(quad_init(4,2)-goal(1,2))^2) sqrt((quad_init(4,1)-goal(2,1))^2+(quad_init(4,2)-goal(2,2))^2) sqrt((quad_init(4,1)-goal(3,1))^2+(quad_init(4,2)-goal(3,2))^2) sqrt((quad_init(4,1)-goal(4,1))^2+(quad_init(4,2)-goal(4,2))^2)];

distance_list=[sqrt((quad_init(1,1)-goal(1,1))^2+(quad_init(1,2)-goal(1,2))^2) sqrt((quad_init(1,1)-goal(2,1))^2+(quad_init(1,2)-goal(2,2))^2) sqrt((quad_init(1,1)-goal(3,1))^2+(quad_init(1,2)-goal(3,2))^2) sqrt((quad_init(1,1)-goal(4,1))^2+(quad_init(1,2)-goal(4,2))^2) ...,
              sqrt((quad_init(2,1)-goal(1,1))^2+(quad_init(2,2)-goal(1,2))^2) sqrt((quad_init(2,1)-goal(2,1))^2+(quad_init(2,2)-goal(2,2))^2) sqrt((quad_init(2,1)-goal(3,1))^2+(quad_init(2,2)-goal(3,2))^2) sqrt((quad_init(2,1)-goal(4,1))^2+(quad_init(2,2)-goal(4,2))^2) ...,
              sqrt((quad_init(3,1)-goal(1,1))^2+(quad_init(3,2)-goal(1,2))^2) sqrt((quad_init(3,1)-goal(2,1))^2+(quad_init(3,2)-goal(2,2))^2) sqrt((quad_init(3,1)-goal(3,1))^2+(quad_init(3,2)-goal(3,2))^2) sqrt((quad_init(3,1)-goal(4,1))^2+(quad_init(3,2)-goal(4,2))^2) ...,
              sqrt((quad_init(4,1)-goal(1,1))^2+(quad_init(4,2)-goal(1,2))^2) sqrt((quad_init(4,1)-goal(2,1))^2+(quad_init(4,2)-goal(2,2))^2) sqrt((quad_init(4,1)-goal(3,1))^2+(quad_init(4,2)-goal(3,2))^2) sqrt((quad_init(4,1)-goal(4,1))^2+(quad_init(4,2)-goal(4,2))^2);
              1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4;
              1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4];
distance_sort=distance_list;
for k=1:(num*target_num)
    for j=2:(num*target_num-k+1)
        if distance_sort(1,j-1)>distance_sort(1,j)
            Q=distance_sort(:,j-1);
            distance_sort(:,j-1)=distance_sort(:,j);
            distance_sort(:,j)=Q;
        end
    end
end
for k=1:(num*target_num)
    if (mirror_dis(distance_sort(2,k),distance_sort(3,k))~=2) && (mirror_count<3)
        mirror_dis(distance_sort(2,k),:)=2;
        mirror_dis(:,distance_sort(3,k))=2;
        mirror_dis(distance_sort(2,k),distance_sort(3,k))=1;
        mirror_count=mirror_count+1;
    end
end
for k=1:num
    for j=1:target_num
        if (mirror_dis(k,j)==1) || (mirror_dis(k,j)==0)
            goal_series(k)=j;
        end
        if mirror_dis(k,j)==0
            row=k;
            line=j;
        end
    end
end
if distance(row,line)>limited_init_distance
    for k=1:num%????????
        sum_dis(k)=(distance(k,line)+distance(row,goal_series(k)))-(distance(k,goal_series(k))+distance(row,line));
        if temp_sum_dis >= sum_dis(k)
            temp_sum_dis=sum_dis(k);
            second_row=k;
        end
    end
    goal_series(row)=goal_series(second_row);
    goal_series(second_row)=line;
end
for k=1:target_num
    temp_goal(k,:)=goal(goal_series(k),:);
end
goal=temp_goal;
end

author: heartxuxuxu


text:
function y = frunge ( x )

%*****************************************************************************80
%
%% FRUNGE sets the Runge data values.
%
%  Discussion:
%
%    Interpolation of the Runge function at evenly spaced points
%    in [-1,1] is a common test.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 January 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real FRUNGE, the value of the derivative of the Runge function.
%
  y = 1.0E+00 ./ ( 1.0E+00 + 25.0E+00 * x .* x );

  return
end

author: johannesgerer

repo: jburkardt-m