text:
%%
% \documentclass[12pt]{article}
%
% \title{ODEbox: A Toolbox for Ordinary Differential Equations\\
% Getting Started with the odeSolveIVP Tool}
%
% \author{Matthew Harker and Paul O'Leary\\
% Institute for Automation\\
% University of Leoben\\
% A-8700 Leoben,
% Austria\\
% URL: automation.unileoben.ac.at\\
% \\
% Original: January 9, 2013\\
% $\copyright$ 2013\\
% \\
% Last Modified: \today}
%%
% This is the documentation for the \lstinline{odeSolveIVP.m} file. This is
% a tool for the solution of initial value problems. It offers a very
% simple interface to solving such problems. Two different examples are
% documented here to show the different possible calls to the tool. The
% examples which have been selected are from Kreyszig~\cite{Kreyszig2010}
% and Adams~\cite{Adams}, since this
% give the analytical results enabling a comparison with the numerical
% solutions.
%
%%
% \section{Tidy up the Workspace}
clear all;
close all;
setUpGraphics;
%
%%
% \section{IVP Example 1}
%
% The equation being considered here is:
% %
% \begin{equation}
%   \dddot{y} + 3 \, \ddot{y} + 3 \, \dot{y} + y = 30 \, \mathrm{e}^{-x},
% \end{equation}
% %
% with the initial conditions
% %
% \begin{equation}
%   y(0) = 3,
%   \hspace{2mm}
%   \dot{y}(0) = -3,
%   \hspace{2mm}
%   \text{and}
%   \hspace{2mm}
%   \ddot{y}(0) = -47.
% \end{equation}
% %
% The analytical solution to this equation is,
% %
% \begin{equation}
%   (3 - 25 \, x^2 + 5 \, x^3)  \, \mathrm{e}^{-x}
% \end{equation}
%
% this example is from~\cite[Chapter 2]{Kreyszig2010}.
%%
% \subsection{Define the paramates for $x$.}
%
% Define the minimum, maximum x values and the number of points desired for
% the solution.
%
xMin = 0 ;
xMax = 10 ;
noPts = 73 ;
%%
% Setup the paramater vector
%
params = [xMin, xMax, noPts ];
%%
% \subsection{Call the \lstinline{odeSolveIVP}}
%
% A vector cof constant coefficients are used for this example equation.
%
[y, x, vals] = odeSolveIVP( [1;3;3;1] , '30*exp(-x)', [ 3 ; -3; -47 ] , params, 13);
%%
% Defining a MATLAB inline function for the analytical solution.
%
f = inline('(3 - 25 * x.^2).* exp(-x) + 5 * (x.^3).*exp(-x)') ;
%%
% \subsection{Plot the Results}
%
% The inline analytical solution is computed during plotting.
%
fig1 = figure;
plot( x, f(x), 'k');
hold on;
plot( x, y, 'ro');
grid on;
xlabel('$$x$$');
ylabel('$$y(x)$$');
legend('Analytic','Numerical','Location','SouthEast');
%
%\caption{Comparison of the analytical and the numerical solution for the first example.}
%%
% \section{IVP Example 2}
%
% The equation being considered here is:
% %
% \begin{equation}
%   x^2 \ddot{y} - 3 \, x\, \dot{y} + 13 y = 0,
% \end{equation}
% %
% with the initial conditions
% %
% \begin{equation}
%   y(0) = 5,
%   \hspace{2mm}
%   \text{and}
%   \hspace{2mm}
%   \dot{y}(0) = 0,
% \end{equation}
% %
% The analytical solution to this equation is,
% %
% \begin{equation}
%   5 \, x^2 \, \cos(3 \log(x)) - (10/3) \, x^2 \, \sin(3 \log(x))
% \end{equation}
%
% this example if from~\cite{Adams}.
%%
% \subsection{Define the paramates for $x$.}
%
% Define the minimum, maximum x values and the number of points desired for
% the solution.
%
noPts = 73 ;
xMin = 1 ;
xMax = 5 ;
%%
% Generate the vector of x values for which the equation should be solved.
% We are now artifically generating a nonuniformly spaced set of pointe to
% compute the solution. This is done to demonstrate the possibility of
% working with arbitrary nodes.
%
z = linspace( 0, 1, noPts )';
z = z.^2;
%
x = xMin + z * (xMax - xMin);
%
%%
% \subsection{Call the \lstinline{odeSolveIVP}}
%
% This is an equation with variable coefficients, i.e., the coefficients
% are functions of $x$. the functions are defines as strings which can be
% computed to inline functions. the notation [] is used to indicate that
% its is a homogeneous equation, i.e., the forcing function is identically
% zero. This call is with x as a vector of points where the solution should
% be computed.
%
[y, x, vals] = odeSolveIVP( {'x.^2';'-3*x';'13'} , [], [ 5 ; 0 ] , x, 13);
%%
% The inline for the analytical solution is, this is only needed for
% comparison purposes.
%
f = inline('5*x.^2.*cos(3*log(x)) - (10/3)*x.^2.*sin(3*log(x))') ;
%%
% \subsection{Plot the Results}
fig1 = figure;
plot( x, f(x), 'k');
hold on;
plot( x, y, 'ro');
grid on;
xlabel('$$x$$');
ylabel('$$y(x)$$');
legend('Analytic','Numerical','Location','NorthWest');
%
%\caption{Comparison of the analytical and the numerical solution for the
% second example. Note the non-uniform spacing of the nodes.}
%
%% Define the Bibliography
%
% \bibliographystyle{plain}
% \bibliography{odebib}


text:
function fx = p05_fx ( x )

%*****************************************************************************80
%
%% P05_FX evaluates ( x + 3 ) * ( x - 1 )^2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(*), the point at which F is to be evaluated.
%
%    Output, real FX(*), the value of the function at X.
%
  fx = ( x + 3.0 ) .* ( x - 1.0 ) .* ( x - 1.0 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [e E] = xyzankurError(X1,X2,type)

% XYZANKURERROR Computes the error between two poses in xyz format
% FORMAT
% DESC Computes the error between two poses in xyz format
% data.
% ARG X1 : first set of poses
% ARG X2 : second set of poses
% RETURN E : mean error for the sequence
%
% SEEALSO : xyzankurDraw, xyzankurModify
%
% COPYRIGHT : Carl Henrik Ek, Andreas Damianou  and Neil Lawrence, 2011
%

% MOCAP

if(nargin<2)
    error('Too few arguments');
end
if(size(X1)~=size(X2))
    error('Dimensions mismatch');
end

E = zeros(size(X1,1),1);
for(i = 1:1:size(X1,1))
    E(i) = mean(sqrt(sum((xyzankur2joint(X1(i,:)) - xyzankur2joint(X2(i,:))).^2,2)));
end
e = mean(E);

return

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% Test file for @chebfun/assignColumns.m.

function pass = test_assignColumns(pref)

if ( nargin < 1 )
    pref = chebfunpref();
end

funlist = {@chebfun, @quasimatrix};

for k = 1:2
    myfun = funlist{k};

    % Generate a few random points to use as test values.
    seedRNG(6178);
    x = 2 * rand(10, 1) - 1;

    % Check a few examples.
    f = myfun(@(x) [sin(x) cos(x) exp(x)], [-1 0 1], pref);
    g = myfun(@(x) x, [-1 1], pref);

    h = assignColumns(f, 1, g);
    h_exact = @(x) [x cos(x) exp(x)];
    err = feval(h, x) - h_exact(x);
    pass(k,1) = norm(err(:), inf) < 10*vscale(h)*eps;

    h = assignColumns(f, 3, g);
    h_exact = @(x) [sin(x) cos(x) x];
    err = feval(h, x) - h_exact(x);
    pass(k,2) = norm(err(:), inf) < 10*vscale(h)*eps;

    g = myfun(@(x) [x x.^2], [-1 0 1], pref);

    h = assignColumns(f, [3 1], g);
    h_exact = @(x) [x.^2 cos(x) x];
    err = feval(h, x) - h_exact(x);
    pass(k,3) = norm(err(:), inf) < 10*vscale(h)*eps;

    h = assignColumns(f, [2 2], g);
    h_exact = @(x) [sin(x) x.^2 exp(x)];
    err = feval(h, x) - h_exact(x);
    pass(k,4) = norm(err(:), inf) < 10*vscale(h)*eps;

    g = myfun(@(x) [x x.^2 x.^3], [-1 -0.5 0.5 1], pref);

    h = assignColumns(f, [1 2 3], g);
    h_exact = @(x) [x x.^2 x.^3];
    err = feval(h, x) - h_exact(x);
    pass(k,5) = norm(err(:), inf) < 10*vscale(h)*eps;

    hc = assignColumns(f, ':', g);
    pass(k,6) = isequal(hc, h);

    h  = assignColumns(f, [2 1], [-0.5 0.5]);
    h_exact = @(x) [(0.5 + 0*x) (-0.5 + 0*x) exp(x)];
    err = feval(h, x) - h_exact(x);
    pass(k,7) = norm(err(:), inf) < 10*vscale(h)*eps;

    h = assignColumns(f.', [2 1], [-0.5 ; 0.5]);
    h_exact = @(x) [(0.5 + 0*x) (-0.5 + 0*x) exp(x)].';
    err = feval(h, x) - h_exact(x);
    pass(k,8) = norm(err(:), inf) < 10*vscale(h)*eps;

    % Check error conditions.
    try
        h = assignColumns(f, [1 2 3], g.');
        pass(k,9) = false;
    catch ME
        pass(k,9) = strcmp(ME.identifier, ...
            'CHEBFUN:CHEBFUN:assignColumns:numCols');
    end

    try
        h = assignColumns(f, [1 2], g);
        pass(k,10) = false;
    catch ME
        pass(k,10) = strcmp(ME.identifier, ...
            'CHEBFUN:CHEBFUN:assignColumns:numCols');
    end

    try
        g2 = myfun(@(x) [x x.^2 x.^3], [0 1], pref);
        h = assignColumns(f, [1 2 3], g2);
        pass(k,11) = false;
    catch ME
        pass(k,11) = strcmp(ME.identifier, ...
            'CHEBFUN:CHEBFUN:assignColumns:domain');
    end

    % Test assign outside of dimension of f:
    h = assignColumns(f, [1 2 4], g);
    h_exact = @(x) [x x.^2 exp(x) x.^3];
    err = feval(h, x) - h_exact(x);
    pass(k,12) = norm(err(:), inf) < 10*vscale(h)*eps;

    %% Test for function defined on unbounded domain:

    % Functions on [-inf b]:

    % Set the domain:
    dom = [-Inf -3*pi];
    domCheck = [-1e6 -3*pi];

    % Generate a few random points to use as test values:
    x = diff(domCheck) * rand(100, 1) + domCheck(1);

    % Array-valued function:
    op = @(x) [exp(x) x.*exp(x) (1-exp(x))./x];
    opg = @(x) exp(-x.^2);
    oph = @(x) [exp(x) exp(-x.^2) (1-exp(x))./x];

    f = myfun(op, dom);
    g = myfun(opg, dom);
    h = assignColumns(f, 2, g);
    err = feval(h, x) - oph(x);
    pass(k,13) = norm(err, inf) < 1e2*vscale(h)*eps;

    %% Test removing columns (#889)

    x = chebfun('x');
    xx = [x 2*x];
    if ( k == 2 )
        xx = cheb2quasi(xx);
    end
    xx(:,2) = [];
    pass(k,14) = numColumns(xx) == 1 && norm(x - xx) < eps;

    x = chebfun('x');
    xx = [x 2*abs(x)];
    if ( k == 1 )
        xx = quasi2cheb(xx);
    end
    xx(:,2) = [];
    pass(k,15) = numColumns(xx) == 1 && norm(x - xx) < eps;

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function oADEVArray = timeADevCalculation(fracFreq, sPeriod, tau)
%function for doing ovelapping Allan calculations array tau values

 phaseError = calculatePhaseError((1/sPeriod),fracFreq); %turn freq readings to frac freq values

 tCount = numel(tau); %get the size of the tau array
 oADEVArray = zeros(1,tCount); %allocate array for Allan Dev values

 %add a wait bar so user knows status
j = 1/tCount;
h = waitbar(0,'Performing ADEV calculations...','CreateCancelBtn','setappdata(gcbf,''canceling'',1)');
setappdata(h,'canceling',0)
%loop through each tau value and calculate ADEV
 for i = 1:tCount
     oADEVArray(i) = calculateADEV(tau(i),sPeriod,phaseError);
     % Check for Cancel button press
    if getappdata(h,'canceling')
        oADEVArray = [];
        oADEVArray = -42881;
        delete(h);
        return;
    end
     waitbar((i*j),h,'Performing ADEV calculations...');
 end

 waitbar(1.0,h,'ADEV calculations done');
 %plot(tau,allanArray);
 %semilogx(tau,allanArray);
 %loglog(tau,hadaArray);
 delete(h);

author: Sable

repo: mcbench-benchmarks


text:
function [D, dD1, dD2] = dist_earth(coord1, coord2)
% [D, dD] = dist_earth(coord1, coord2)
% coord1 is 2 x N
% coord2 is 2 x M
% Returns N x M matrix D of mutual distances.

if nargin < 2 || isempty(coord2)
  coord2 = coord1;
end

R = 6371.01; % Spherical Earth radius approximation

n1 = cols(coord1);
n2 = cols(coord2);

D = zeros(n1, n2);

coord1 = pi/180 * coord1;
coord2 = pi/180 * coord2;
if nargout >= 2
  dD1 = zeros([n1, n2, 2]);
end
if nargout >= 3
  error('Hmm.. maybe you shouldn''t use third output? :)');
  dD2 = zeros([n1, n2, 2]);
end

for i=1:n1
  lat1 = coord1(2,i);
  lat2 = coord2(2,:);
  dlon = (coord2(1,:) - coord1(1,i));
  f = sin(lat2).*sin(lat1) + cos(lat2).*cos(lat1).*cos(dlon);
  f(f>=1-eps) = 1-eps; % correction because of numerical errors, f should be [-1,1]
  f(f<=eps-1) = eps-1; % correction because of numerical errors, f should be [-1,1]
  D(i,:) = R * acos(f);
  if ~isreal(D)
    coord1
    f
%    D
    error('Oohps! Complex distance..');
  end
  %dD_df(isinf(dD_df)) = 0;%-R ./ sqrt(1-f.^2);

  if nargout >= 2
    dD_df = -R ./ sqrt(1-f.^2);
    dD1(i,:,1) = pi/180 * dD_df .* cos(lat1) .* cos(lat2) .* sin(dlon);
    dD1(i,:,2) = pi/180 * dD_df .* (cos(lat1).*sin(lat2) - sin(lat1).* ...
                                    cos(lat2).*cos(dlon));
  end
% $$$   if any(isnan(dD1))
% $$$     dD1
% $$$     error('WTF?');
% $$$   end
% $$$   if nargout >= 3
% $$$     dD2(i,:,1) = -dD1(i,:,1);
% $$$     dD2(i,:,2) = dD_df .* (sin(lat1).*cos(lat2) - cos(lat1).*sin(lat2).*cos(dlon));
% $$$   end
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [C,vol] = centroid(V,F,varargin)
  % CENTROID Compute the centroid of a closed polyhedron boudned by (V,F)
  %
  % C = centroid(V,F)
  % [C,vol] = centroid(V,F,'ParameterName',ParameterValue, ...)
  %
  % Inputs:
  %   V  #V by 3 list of mesh vertex positions
  %   F  #F by 3 list of triangle mesh indices
  %   Optional:
  %     'Robust' followed by whether to use more robust but costlier method for
  %       nearly closed input. {false}
  % Outputs:
  %   C  3-vector of centroid location
  %   vol  total volume of polyhedron
  %

  % default values
  robust = false;
  % Map of parameter names to variable names
  params_to_variables = containers.Map( {'Robust'}, {'robust'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  if robust
    assert(size(V,2) == 3,'Only 3d supported');
    [TV,TT] = cdt(V,F);
    BC = barycenter(TV,TT);
    w = winding_number(V,F,barycenter(TV,TT))/(4*pi);
    TT = TT(abs(w)>0.5,:);
    BC = BC(abs(w)>0.5,:);
    vol = volume(TV,TT);
    C = sum(bsxfun(@times,vol,BC))/sum(vol);
    vol = sum(vol);
  else
    % "Calculating the volume and centroid of a polyhedron in 3d" [Nuernberg 2013]
    % http://www2.imperial.ac.uk/~rn/centroid.pdf
    switch size(V,2)
    case 2
      % https://en.wikipedia.org/wiki/Centroid#Centroid_of_polygon
      % Rename corners
      A = V(F(:,1),:);
      B = V(F(:,2),:);
      D = A(:,1).*B(:,2) - B(:,1).*A(:,2);
      vol = 0.5*sum(D);
      C = (1/(6*vol)).*sum((A+B).*D);
    case 3
      % Rename corners
      A = V(F(:,1),:);
      B = V(F(:,2),:);
      C = V(F(:,3),:);
      % Needs to be **unnormalized** normals
      N = cross(B-A,C-A,2);
      % total volume via divergence theorem: ∫ 1
      vol = sum(sum(A.*N))/6;
      % centroid via divergence theorem and midpoint quadrature: ∫ x
      C = 1/(2*vol)*(1/24* sum(N.*((A+B).^2 + (B+C).^2 + (C+A).^2)));
    end
  end

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%computes the lag of the average magnitude difference function
%> called by ::ComputePitch
%>
%> @param x: audio signal
%> @param iBlockLength: block length in samples
%> @param iHopLength: hop length in samples
%> @param f_s: sample rate of audio data
%>
%> @retval f_0 amdf lag (in Hz)
%> @retval t time stamp of f_0 estimate (in s)
% ======================================================================
function [f_0, t] = PitchTimeAmdf(x, iBlockLength, iHopLength, f_s)

    % blocking
    [x_b, t] = ToolBlockAudio(x, iBlockLength, iHopLength, f_s);
    iNumOfBlocks = size(x_b, 1);

    % allocate memory
    f_0 = zeros(1, iNumOfBlocks);

    % initialization
    f_max = 2000;
    f_min = 50;
    eta_min = round(f_s/f_max);
    eta_max = round(f_s/f_min);

    for n = 1:iNumOfBlocks

        if (sum(abs(x_b(n, :))) == 0)
            f_0(n) = 0;
            continue;
        end

        % calculate the amdf_I minimum
        afAMDF = amdf_I(x_b(n, :), eta_max);
        [fDummy, T_0(n)] = min(afAMDF(1+eta_min:end));

        % convert to Hz
        f_0(n) = f_s ./ (T_0(n) + eta_min);
    end
end

function [AMDF] = amdf_I(x, eta_max)
    K = length(x);

    AMDF = ones(1, K);

    for eta=0:min(K-1,eta_max-1)
        AMDF(eta+1) = sum(abs(x(1:K-1-eta) - x(eta+2:end))) / K;
    end
end

author: alexanderlerch

repo: ACA-Code


text:
function [ net ] = replace_last_layer( net, lr_old_layers, lr_new_layer, nClasses, dropOutRatio )

  [net_info] = vl_simplenn_display(net);
  net_dimensions = net_info.dataSize(3,:);

  scal = 1 ;
  net.layers{end-1} = struct('type', 'conv', ...
                             'filters', 0.01/scal * randn(1,1,net_dimensions(end-2),nClasses,'single'), ...
                             'biases', zeros(1, nClasses, 'single'), ...
                             'stride', 1, ...
                             'pad', 0, ...
                             'filtersLearningRate',  lr_new_layer(1) , ...
                             'biasesLearningRate',  lr_new_layer(2) , ...
                             'filtersWeightDecay', 1, ...
                             'biasesWeightDecay', 0) ;

  net.layers{end} = struct('type', 'softmaxloss') ;

  if dropOutRatio > 0
    % inject dropout for the two FC layers:
    dropout = struct('type', 'dropout', ...
                               'rate', dropOutRatio) ;
    nL = numel(net.layers);
    net.layers = {net.layers{1:nL-4} dropout net.layers{nL-3:nL-2} dropout net.layers{nL-1:end}};
  end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% % multichannel superdirective beamformer
% % Ziteng Wang @ 201811
% % reference: Berkun R, Cohen I, Benesty J. A tunable beamformer for robust superdirective beamforming[C]//Acoustic Signal Enhancement (IWAENC), 2016 IEEE International Workshop on. IEEE, 2016: 1-5.
%

clear all
close all
warning('off')
addpath('..\STFT\')
addpath('..\Simulation\')
addpath('..\Simulation\RIR-Generator\')

%% simulation start
flatStart = 1;
prefix = '';   % for saving file

speechDir = '..\Simulation\Data\';
speechFile = 'fajw0_sa1.wav';
saveDir = 'GeneratedData\';
if ~exist(saveDir)
    mkdir(saveDir)
end
speech = audioread([speechDir speechFile]);

% configuration
cfg = [];
cfg.fs = 16000;                     % sampling rate
cfg.room = [6 5 3];                 % room dimension (m)
cfg.T60 = 0.3;                      % reverberation time (s)

cfg.Nch = 6;
cfg.micCenter = [2 3 1.5];          % array center (m)
cfg.micCoordinate =[0.0425,0.0,0.0;
        0.02125,0.03680608,0.0;
        -0.02125,0.03680608,0.0;
        -0.0425,0.0,0.0;
        -0.02125,-0.03680608,0.0;
        0.02125,-0.03680608,0.0;];  % microphone array coordinates

cfg.az = 180;
cfg.el = 0;
cfg.dist = 2;

cfg.SNR = 10;
if cfg.SNR ~= inf               % inf means no noise
    cfg.noiseType = 'diffuse';    % choice {'white' 'diffuse' 'recorded'}
    if strcmp(cfg.noiseType, 'recorded')
        % check first the noise is longer than speech!
        cfg.noiseFile = '';
    end
end

cfg.SIR = inf;
if cfg.SIR ~= inf
    % check the interference is longer than speech!
    cfg.interfFile = '';
    cfg.azITF = 180;
    cfg.elITF = 0;
    cfg.distITF = 2;
end

% setup room and collect data
setup_room
setup_noise

%% processing start
% the following parts are specific to the algorithm
Nfft = 1024;

Y = stft_multi_2(y, Nfft);
[Nframe, Nbin, Nch] = size(Y);

%% superdirective beamformer

%%% diagloading value !
diagLoad = 0.01;

% diffuse noise coherence
SDcov = ones(Nch, Nch, Nbin);

% in matlab: sinc(x) = sin(pi*x) / (pi*x)
for bin = 1:Nbin
    for chi=1:Nch-1
        for chj=chi+1:Nch
            SDcov(chi,chj,bin) = sinc(2*(bin-1)/Nfft*fs*norm(micPose(chi,:)-micPose(chj,:))/c);
            SDcov(chj,chi,bin) = SDcov(chi,chj,bin);
        end
    end
end

% free field steering vector
steerVec = zeros(Nch, Nbin);
for bin = 1:Nbin
    steerVec(:,bin)= exp(-2*1i*pi*(bin-1)/Nfft*fs*TDOA);
end

% filter coefficients:
% h = PhiN^-1 * d / d^H * PhiN^-1 * d
PhiNinv = zeros(Nch, Nch, Nbin);
for bin = 1:Nbin
    PhiNinv(:,:,bin) = inv(SDcov(:,:,bin) + diagLoad*eye(Nch));
end
tmp = squeeze(sum(bsxfun(@times, PhiNinv, permute(steerVec, [3,1,2])), 2));
hSuperDirective = bsxfun(@rdivide, tmp, sum(bsxfun(@times, conj(steerVec), tmp),1));

% apply the filter
Xest = sum(bsxfun(@times, conj(permute(hSuperDirective, [3,2,1])), Y), 3);
xest = istft_multi_2(Xest, length(speech));

%% record and plot
audiowrite([saveDir prefix 'SuperDirective' postfix '.wav'], xest, fs);

% plot the beampattern
dist = 2;
el = 0;
for az=1:360
    sourcePose = dist * [cos(el/180*pi)*cos(az/180*pi) cos(el/180*pi)*sin(az/180*pi) sin(el/180*pi)] + micCenter;
    TDOA = sqrt(sum((bsxfun(@minus, sourcePose, micPose)).^2, 2))/c;
    for bin=1:Nbin
        azVec = exp(-2*1i*pi*(bin-1)/Nfft*fs*TDOA);
        beamPattern(az,bin) = abs(sum(conj(hSuperDirective(:,bin)) .* azVec));
        beamPattern(az,bin) = max(20*log10(beamPattern(az,bin)), -30); % in dB
    end
end
figure;imagesc(beamPattern); axis xy; colorbar
title(['BeamPattern of SD steered towards azimuth ' num2str(cfg.az)]);
xlabel('frequency index'); ylabel('azimuth angle');

% white noise gain
for bin=1:Nbin
    whiteNoiseGain(bin) = 10*log10(1/(hSuperDirective(:,bin)'*hSuperDirective(:,bin)));
end
figure; plot(whiteNoiseGain)
title(['WNG of SD with diagnal loading (' num2str(diagLoad) ')']);
xlabel('frequency index'); ylabel('WNG in dB');

author: ZitengWang

repo: MASP

save_path: github-repos/MATLAB/ZitengWang-MASP


text:
function [class, Simil] = classifier(data, ideals, y)
%   INPUT:
%
%   data    =    datamatrix
%   ideals  =    idealvectors
%   y(1)    =    p-value in similarity measure
%   y(2)    =    alpha value for OWA weights
%   y(3)    =    used owa aggregator
%
%
%   OUTPUT:
%
%   class   =    column vector of the classes in which the samples are classified
%   Simil   =    similarity values for each class

[nc, v_dim] = size(ideals);
d_dim = size(data,1);
Simil = zeros(d_dim, nc);

if nargin==2
    y = [1, 1, 1];
end

for j = 1 : nc

    Ideal = repmat(ideals(j,:),d_dim,1);

    if y(3) == 1 %Using OWA with Linguistic quantifier 1
        w=owaw1(v_dim,y(2));
        tmpmatrix=(1-abs(data.^y(1)-Ideal.^y(1))).^(1/y(1));
        Simil(:,j)=owamatrix(tmpmatrix,w);
    elseif    y(3) == 2 %Using OWA with Linguistic quantifier 2
        w=owaw2(v_dim,y(2));
        tmpmatrix=(1-abs(data.^y(1)-Ideal.^y(1))).^(1/y(1));
        Simil(:,j)=owamatrix(tmpmatrix,w);
    elseif    y(3) == 3 %Using OWA with Linguistic quantifier 3
        w=owaw3(v_dim,y(2));
        tmpmatrix=(1-abs(data.^y(1)-Ideal.^y(1))).^(1/y(1));
        Simil(:,j)=owamatrix(tmpmatrix,w);
    elseif    y(3) == 4 %Using OWA with Linguistic quantifier 4
        w=owaw4(v_dim,y(2));
        tmpmatrix=(1-abs(data.^y(1)-Ideal.^y(1))).^(1/y(1));
        Simil(:,j)=owamatrix(tmpmatrix,w);
    elseif  y(3) == 5   %O'Hagan's method
        tmpmatrix=(1-abs(data.^y(1)-Ideal.^y(1))).^(1/y(1));
        weights=Ohaganw(v_dim,y(2));
        index= 1;
        w=weights(index,:);
        Simil(:,j)=owamatrix(tmpmatrix,w);
    end

end

[simil_val, class] = max(Simil');
class=class';


text:
function vw = roiSelectDisplay(vw,promptStr)
% Select a subset of the ROIS
%
%   vw = roiSelectDisplay(vw,promptStr)
%
% The selected ROIs are assigned to the variable vw.ui.showROIs
% They are positive if we show the whole ROI and negative if we show just
% the ROI outlines.
% The drawing mode is inherited from the current sign of ui.showROIs.  If
% that is zero, the mode is positive.
%
% Example:
%     tmp = roiSelectDisplay(vw,'Select subset of ROIs')
%

if notDefined('vw'),        vw = getCurView;            end
if notDefined('promptStr'), promptStr = 'Select ROIs '; end

R = viewGet(vw,'rois');

if isempty(R), disp('No ROIs.'); return; end

% There should be an roiGet/Set routine to get the names
% -ras: maybe; but probably not for this poor, molested codebase.
% there's a nice functional roiGet/Set in mrVista2.
roiNames = {R.name};
[sList ok] = listdlg('PromptString', promptStr, ...
    'SelectionMode', 'multiple', ...
    'ListString', roiNames);

if ok
    % Set the display mode
    % (ras 01/07: now the rendering method is handled separately)
    vw.ui.showROIs = sList;
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [tt,ind_left]=tt_crossl(d,n,fun,ind_right)
%One left-to-right sweep of the TT-cross method.
%   [TT,IND_LEFT]=TT_CROSSL(D,N,FUN,IND_RIGHT) Computes one QR-maxvol sweep
%   of the TT-cross method. The input is the pair (D,N) that determines the
%   size of the tensor, FUN is the function handle to evaluate a particular
%   element of the tensor, i.e., VAL=FUN(IND). To pass parameters, please
%   use anonymous function handles in MATLAB
%
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

if ( numel(n) == 1 )
  n=n*ones(d,1);
end
sz=n;
tt=cell(d,1);
ind_left=cell(d,1); %Computed ind_left
ind_r=ind_right{1};
r1=size(ind_r,2);
ncur=sz(1);
mat=zeros(ncur,r1);
for i=1:ncur
   for s=1:r1
      ind_f=[i,ind_r(:,s)'];
      val=fun(ind_f);
      mat(i,s)=val;
   end
end
[qs,~]=qr(mat,0);
ind=maxvol2(qs);
ind_left{1}=ind;
mat=qs/qs(ind,:);
tt{1}=mat;

for k=2:d-1
     ncur=sz(k);
     ind_l=ind_left{k-1};
     ind_r=ind_right{k};
     r2=size(ind_l,2);
     r3=size(ind_r,2);
     core=zeros(ncur,r2,r3);
     for i=1:ncur
         for s2=1:r2
             for s3=1:r3
                 ind_f=[ind_l(:,s2)',i,ind_r(:,s3)'];
                 val=fun(ind_f);
                 core(i,s2,s3)=val;
             end
         end
     end
     core=permute(core,[2,1,3]); core=reshape(core,[r2*ncur,r3]);
     rnew=min(r2*ncur,r3);
     [qs,rs]=qr(core,0);
     ind=maxvol2(qs);
     ind_old=ind_left{k-1};
     ind_new=zeros(k,rnew);
     ncur=sz(k);
     for s=1:rnew
        f_in=ind(s);
        w1=tt_ind2sub([r2,ncur],f_in);
        rs=w1(1); js=w1(2);
        ind_new(:,s)=[ind_old(:,rs)',js];
     end
     ind_left{k}=ind_new;
     core=qs/qs(ind,:);
     core=reshape(core,[r2,ncur,rnew]); core=permute(core,[2,1,3]);

    tt{k}=core;
end
ncur=sz(d);
ind_l=ind_left{d-1};
r=size(ind_l,2);
mat=zeros(ncur,r);
for j=1:ncur
    for s=1:r
      ind_f=[ind_l(:,s)',j];
      val=fun(ind_f);
      mat(j,s)=val;
    end
end
tt{d}=mat;
tt=tt_tensor(tt);
return
end

author: oseledets

repo: TT-Toolbox


text:
function vl_setupnn()
%VL_SETUPNN Setup the MatConvNet toolbox.
%   VL_SETUPNN() function adds the MatConvNet toolbox to MATLAB path.

% Copyright (C) 2014-15 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

root = vl_rootnn() ;
addpath(fullfile(root, 'matlab')) ;
addpath(fullfile(root, 'matlab', 'mex')) ;
addpath(fullfile(root, 'matlab', 'simplenn')) ;
addpath(fullfile(root, 'matlab', 'xtest')) ;
addpath(fullfile(root, 'examples')) ;

if ~exist('gather')
  warning('The MATLAB Parallel Toolbox does not seem to be installed. Activating compatibility functions.') ;
  addpath(fullfile(root, 'matlab', 'compatibility', 'parallel')) ;
end

if numel(dir(fullfile(root, 'matlab', 'mex', 'vl_nnconv.mex*'))) == 0
  warning('MatConvNet is not compiled. Consider running `vl_compilenn`.');
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function [plotMins,plotMaxs] = prtPlotUtilGetPlotLims(PrtClassifier,PrtDataSet)
% Internal function,
% xxx Need Help xxx

% Here we get the limits for plotting:
if nargin == 1 || isempty(PrtDataSet)
    testY = [];
    testX = [];
else
    testX = getObservations(PrtDataSet);
end

trainX = getObservations(PrtClassifier.PrtDataSet);
plotMins = min(cat(1,trainX,testX));
plotMaxs = max(cat(1,trainX,testX));

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [ energy x ] = neb_energies( directory )
%NEB_ENERGIES Extract energies and hyperdistance from NEB calculation.
%   [energy,x] = NEB_ENERGIES(directory) extracts the energies and
%   hyperdistance of each image from an NEB calculation. The optional
%   argument directory can be used to specify the directory containing the
%   subdirectories 00, 01, etc. If directory is specified, the current
%   directory is used. If no OUTCAR file is present in subdirectories 00
%   and N+1, then first and last value of energy of will be NaN.
%
%   See also NEB_SPLINE, PLOT_NEB_SPLINE, HYPERDISTANCE.

% shouldn't change directory

    startDir = pwd();
    if nargin == 1
        cd(directory);
    end

    nimg = num_images();

    energy = zeros(1,nimg+2);
    x = zeros(1,nimg+2);

    % read starting point
    if exist('00/OUTCAR','file')
        %buffer = import_oszicar('00/OSZICAR');
        energy(1) = import_outcar('00/OUTCAR','energy');
    else
        energy(1) = nan;
    end

    geo1 = import_poscar('00/POSCAR');
    x(1) = 0;

    % read images
    for i = 2:nimg+1
       %buffer = import_oszicar([sprintf('%02d',i-1) '/OSZICAR']);
       energy(i) = import_outcar([sprintf('%02d',i-1) '/OUTCAR'],'energy');
       geo2 = import_poscar([sprintf('%02d',i-1) '/CONTCAR']);
       x(i) = x(i-1) + hyperdistance(geo1,geo2);
       geo1=geo2;
    end

    % read ending point
    if exist([sprintf('%02d',nimg+1) '/OUTCAR'],'file')
        %buffer = import_oszicar([sprintf('%02d',nimg+1) '/OSZICAR']);
        energy(nimg+2) = import_outcar([sprintf('%02d',nimg+1) '/OUTCAR'],'energy');
    else
        energy(nimg+2) = nan;
    end

    geo2 = import_poscar([sprintf('%02d',nimg+1) '/POSCAR']);
    x(nimg+2) = x(nimg+1) + hyperdistance(geo1,geo2);

    cd(startDir);

end


text:
function [f,Ph]=haralick_n(qs,nL,q, nanFlag)
% Haralick textures measurements

Ph = coocurrance_alldir_mod(q);
% last row and column corresponds to outside ROI!
if nanFlag
    Ph=Ph(1:end-1,1:end-1);
    nL=nL-1;
end

%compute features
R=sum(Ph(:));
Ph=Ph/R;
% Energy (1)
f(1)=sum(sum(Ph.^2));
% Contrast (2)
f(2)=0.0;
for n=0:nL-1
   temp=0;
   for i=1:nL
      for j=1:nL
         if (abs(i-j) == n)
            temp=temp+Ph(i,j);
         end
      end
   end
   f(2)=f(2)+n^2*temp;
end
% Correlation
%using symmetry Ph'=Ph!
Px=sum(Ph);
Py=sum(Ph');
vec=[1:nL];
ux=sum(Px .*vec);
uy=sum(Py .*vec);

varx=sum(Px .* vec.^2)-ux^2;
sigx=sqrt(varx);
vary=sum(Py .* vec.^2)-uy^2;
sigy=sqrt(vary);
u=vec*Ph(i,j)*vec';
f(3)=(u-ux*uy)/(sigx*sigy);

%Entropy (3)
res = 1e-10; % realmin
f(4)=-sum(sum(Ph.*log2(Ph+res))); % log????

% variance
f(5)=0;
for i=1:nL
   for j=1:nL
      f(5)=f(5)+(i-u)^2*Ph(i,j);
   end
end
% P(x+y)
f(6)=0;  % sum of entropy...
for k=2:2*nL
   temp=0.0;
   for i=1:nL
      for j=1:nL
         if ((i+j) == k)
            temp=temp+Ph(i,j);
         end
      end
   end
   Pxpy(k)=temp;
   f(6)=f(6)-temp*log2(temp+res);
end

% inverse different moment..
f(7)=0;
f(8)=0;  %Homogeneity (4)
for i=1:nL
   for j=1:nL
      temp1=1/(1+(i-j)^2)*Ph(i,j);
      temp2=1/(1+abs(i-j))*Ph(i,j);
      f(7)=f(7)+temp1;
      f(8)=f(8)+temp2;
   end
end

return

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function indxV = getAssociatedDose(uid)
%"getAssociatedDose"
%   Returns a dose index for corresponding dose uid passed.
%
%	DK 12/07/2006
% copyright (c) 2001-2008, Washington University in St. Louis.
% Permission is granted to use or modify only for non-commercial,
% non-treatment-decision applications, and further only if this header is
% not removed from any file. No warranty is expressed or implied for any
% use whatever: use at your own risk.  Users can request use of CERR for
% institutional review board-approved protocols.  Commercial users can
% request a license.  Contact Joe Deasy for more information
% (radonc.wustl.edu@jdeasy, reversed).
%
%Usage:
%   indxV = getAssociatedDose(uid)

global planC
indexS = planC{end};

% preallocate memory
if length(planC{indexS.dose})< 1
    indxV = [];
    return
end

[doseUID{1:length(planC{indexS.dose})}] = deal(planC{indexS.dose}.doseUID);

if length(uid)==0
    indxV = [];
else
    [jnk,indxV] = ismember(uid,doseUID);
    if indxV == 0;
        indxV = [];
    end
end


text:
function a = pascal1 ( n )

%*****************************************************************************80
%
%% PASCAL1 returns the PASCAL1 matrix.
%
%  Formula:
%
%    if ( J = 1 )
%      A(I,J) = 1
%    elseif ( I = 0 )
%      A(1,J) = 0
%    else
%      A(I,J) = A(I-1,J) + A(I-1,J-1)
%
%  Example:
%
%    N = 5
%
%    1  0  0  0  0
%    1  1  0  0  0
%    1  2  1  0  0
%    1  3  3  1  0
%    1  4  6  4  1
%
%  Properties:
%
%    A is a "chunk" of the Pascal binomial combinatorial triangle.
%
%    A is generally not symmetric: A' /= A.
%
%    A is nonsingular.
%
%    A is lower triangular.
%
%    A is integral, therefore det ( A ) is integral, and
%    det ( A ) * inverse ( A ) is integral.
%
%    det ( A ) = 1.
%
%    A is unimodular.
%
%    LAMBDA(1:N) = 1.
%
%    (0, 0, ..., 0, 1) is an eigenvector.
%
%    The inverse of A is the same as A, except that entries in "odd"
%    positions have changed sign:
%
%      B(I,J) = (-1)^(I+J) * A(I,J)
%
%    The product A*A' is a Pascal matrix
%    of the sort created by subroutine PASCAL2.
%
%    Let the matrix C have the same entries as A, except that
%    the even columns are negated.  Then Inverse(C) = C, and
%    C' * C = the Pascal matrix created by PASCAL2.
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Gregory, David Karney,
%    Example 3.15,
%    A Collection of Matrices for Testing Computational Algorithms,
%    Wiley, 1969, page 43,
%    LC: QA263.G68.
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n

      if ( j == 1 )
        a(i,j) = 1.0;
      elseif ( i == 1 )
        a(i,j) = 0.0;
      else
        a(i,j) = a(i-1,j-1) + a(i-1,j);
      end

    end
  end

  return
end


text:
function [ mSynCat, vMain] = calc_SynCat(nNbkg,vSynCat,nSynCat_,nSynMode, vAfter,vEtas,mCatalog,bPSQ,vPSQ,mPSQ)
    % calculates a synthetic catalog with a backgroundrate based number of earthquakes in the catalog of the period before starting date.
    %
    % Example: [mSynCat, vMain] = calc_SynCat(2500,2.5,2.5,8,100,'January 1,1980','December 31,1990',6.5,1,0,a);
    %
    % This function (calc_SynCat) calculates a synthetic catalog with a
    % backgroundrate based number of earthquakes in the catalog of the period
    % before starting date. This function is doing the following:
    % 1) creates a synthetic catalog with background rate based on input parameters.
    % 2) calculate omori-type aftershock and add this to the background catalog,
    % that was created by the properties of the inputbackground catalog. Therefore
    % the program package ETESProject by karen Felzer was used.
    %
    % Author: van Stiphout, Thomas
    % Email: vanstiphout@sed.ethz.ch
    % Created: 2. Mar 2007
    % Changed: -
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %  Variables
    % Input:
    % nNbkg         Number of events in background catalog during assigned period
    % vSynCat       Vektor with contraints on synthetic catalog
    %                   1-4: fMinLon fMaxLon fMinLat fMaxLat in degree
    %                   5-6: fMinDepth fMaxDepth  in km
    %                   7-8: fMinTime fMaxTime (duration of catalog)
    %                   9:   fRate Reduction of earthquakes for entire catalog
    %                          (0.75 = 25% reduction) in period (fT-fTw) - fT
    %                   10:  fBValue           b-value of the catalog
    %                   11:  fMc Magnitude of Completness (minimum Magnitude in Catalog)
    %                   12:  fIncMagnitude increment steps
    % nSynCat_      1: only backgorund, 2: background + ETAS
    % nSynMode      Type of synthetic catalog 0:homogeneous distr; 1:based
    %                on real catalog 2:hypocenter based on real catalog, magnitude and focal
    %                time is randomly computed
    % vAfter        Definition/Constraints for Aftershocks
    %                   1: Mmin: Minimum earthquake magnitude to report and simulate
    %                   2: MminR: Minimum earthquake magnitude for reporting
    %                       in the output catalog
    %                   3: maxm: Maximum earthquake magnitude
    %                   4: DMax: Maximum distance for aftershocks
    %                   5: magb: Cutoff magnitude for the point source vs. plane source
    %                       representation of the mainshock
    % vEtas         Parameters for ETAS:[CX, P, A] standard=[0.095, 1.34, 0.008]
    % mCatalog      Declustered catalog needed only for nSynMode=1
    % bPSQ          with (1) or without (0) PSQ
    % vPSQ          values for PSQ (lon, lat, N, Tw, T)
    %
    % Output:
    % SynCat Simulated ETES catalog [yr lon lat month day mag depth hrs min sec]
    % vMain            Vector with assignments if main (true) or aftershock (false)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % initialize:
    clear FaultParam;

    % background:
    fBValue=vSynCat(10);
    fMc=vSynCat(11);
    fInc=vSynCat(12);
    fMinLat=vSynCat(3);
    fMaxLat=vSynCat(4);
    fMinLon=vSynCat(1);
    fMaxLon=vSynCat(2);
    fMinDepth=vSynCat(5);
    fMaxDepth=vSynCat(6);
    nRchange=vSynCat(9);
    [fYr1, nMn1, nDay1, nHr1, nMin1, nSec1]=decyear2mat(vSynCat(7));
    [fYr2, nMM2, nDD2, nHH2, nMN2, nSS2]=decyear2mat(vSynCat(8));

    % vAfter
    Mmin=vAfter(1);
    MminR=vAfter(2);
    maxm=vAfter(3);
    DMax=vAfter(4);
    magb=vAfter(5);
    FaultParam=[0 0 0 0 0 0 0];

    % vEtas
    CX=vEtas(1);
    P=vEtas(2);
    A=vEtas(3);

    switch nSynCat_
        case 1
            % create background rate
            mSynCat=syn_catalog(nNbkg, fBValue, fMc, fInc, fMinLat,fMaxLat,...
                fMinLon,fMaxLon,fMinDepth,fMaxDepth, fYr1,fYr2,nSynMode,...
                mCatalog);
            if bPSQ
                mSynCat=sr_makePSQ(mSynCat,vPSQ(2),vPSQ(3),vPSQ(4),...
                    vPSQ(1),mPSQ);
            end
            %       mSynCat(:,1)=mSynCat(randperm(size(mSynCat,1))',1);
            vMain=ones(nNbkg,1);
        case 2   % create background rate + ETAS
            nNbkgA=floor(nNbkg/(1+nRchange)*1);
            nNbkgB=ceil(nNbkg/(1+nRchange)*nRchange);
            if nRchange == 1
                mCat=syn_catalog(nNbkg, fBValue, fMc, fInc, fMinLat,fMaxLat,...
                    fMinLon,fMaxLon,fMinDepth,fMaxDepth, fYr1,fYr2,nSynMode,...
                    mCatalog);
            else
                mCatA=syn_catalog(nNbkgA, fBValue, fMc, fInc, fMinLat,fMaxLat,...
                    fMinLon,fMaxLon,fMinDepth,fMaxDepth,fYr1,(fYr1+fYr2)/2,...
                    nSynMode,mCatalog);
                mCatB=syn_catalog(nNbkgB, fBValue, fMc, fInc, fMinLat,fMaxLat,...
                    fMinLon,fMaxLon,fMinDepth,fMaxDepth,(fYr1+fYr2)/2,fYr2,...
                    nSynMode,mCatalog);
                mCat=[mCatA;mCatB];
            end
            %         mCat(:,1)=mCat(randperm(size(mCat,1))',1);
            [Y,Idx]=sort(mCat(:,3),1);
            mCat=mCat(Idx,:);
            if bPSQ
                mCat=sr_makePSQ(mCat,vPSQ(2),vPSQ(3),vPSQ(4),...
                    vPSQ(1),mPSQ);
            end
            mCat1=[fix(mCat(:,3)) mCat(:,4) mCat(:,5) mCat(:,8) mCat(:,9) rand(size(mCat,1),1)*60 mCat(:,2) mCat(:,1) mCat(:,7) mCat(:,6)];
            %         mCat1(:,1)=mCat1(:,1)-(max(mCat1(:,1))-min(mCat1(:,1)));
            %  [N,Mmax,catalog,IDAll,HistoryAll] = ETESProject2(cat1,RateGrid2,0.095,1.34,0.008,Mmin,MminR,maxm,DMax,sYr1,sYr2,6.5,FaultParam)
            [N,Mmax,catalog,IDAll,HistoryAll] = ETESProject3(mCat1,vEtas,vAfter,fYr1,fYr2,FaultParam);

            mSynCat(:,1) = catalog(:,8); % lon
            mSynCat(:,2) = catalog(:,7);  % lat
            mSynCat(:,3) = catalog(:,1);  % years
            mSynCat(:,4) = catalog(:,2);  % month
            mSynCat(:,5) = catalog(:,3);  % day
            mSynCat(:,6) = catalog(:,10); % mag
            mSynCat(:,7) = catalog(:,9);  % depth
            mSynCat(:,8) = catalog(:,4);  % hours
            mSynCat(:,9) = catalog(:,5);  % minutes
            mSynCat(:,10) = catalog(:,6);  % seconds

            mSynCat(:,3) = decyear([mSynCat(:,3) mSynCat(:,4) mSynCat(:,5) ...
                mSynCat(:,8) mSynCat(:,9) mSynCat(:,10)]); % yrs decimal
            size(mSynCat)
            vMain=(catalog(:,12)==0);
    end
    % disp('cat fertig')
end

author: CelsoReyes


text:
function a = stripe ( n )

%*****************************************************************************80
%
%% STRIPE returns the STRIPE matrix.
%
%  Example:
%
%    N = 7
%
%    5  2  1  1  .  .  .
%    2  6  3  1  1  .  .
%    1  3  6  3  1  1  .
%    1  1  3  6  3  1  1
%    .  1  1  3  6  3  1
%    .  .  1  1  3  6  2
%    .  .  .  1  1  2  5
%
%  Properties:
%
%    A is integral, therefore det ( A ) is integral, and
%    det ( A ) * inverse ( A ) is integral.
%
%    A is symmetric: A' = A.
%
%    A is banded, with bandwidth 7.
%
%    A is centrosymmetric: A(I,J) = A(N+1-I,N+1-J).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joan Westlake,
%    A Handbook of Numerical Matrix Inversion and Solution of
%    Linear Equations,
%    John Wiley, 1968,
%    ISBN13: 978-0471936756,
%    LC: QA263.W47.
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n

           if ( j == i - 3 )
        a(i,j) = 1.0;
      elseif ( j == i - 2 )
        a(i,j) = 1.0;
      elseif ( j == i - 1 )
        if ( j == 1 | j == n - 1 )
          a(i,j) = 2.0;
        else
          a(i,j) = 3.0;
        end
      elseif ( j == i )
        if ( i == 1 | i == n )
          a(i,j) = 5.0;
        else
          a(i,j) = 6.0;
        end
      elseif ( j == i + 1 )
        if ( j == 2 | j == n )
          a(i,j) = 2.0;
        else
          a(i,j) = 3.0;
        end
      elseif ( j == i + 2 )
        a(i,j) = 1.0;
      elseif ( j == i + 3 )
        a(i,j) = 1.0;
      else
        a(i,j) = 0.0;
      end

    end
  end

  return
end


text:
function [y, t] = Raw_to_DeltaConc_Example(filename, hbType, channels)
%
%    Example 1:   Plot channel 1 delta HbO concentration, derived from
%                 acquisition file test.snirf
%
%       cd '<root folder>\Homer3\UnitTests\Example_RunFuncOffline'
%       Raw_to_DeltaConc_Example('./test.snirf');
%
%    Example 2:   Plot channels 1 and 2 delta HbR concentration, derived
%                 from acquisition file test.snirf
%
%       cd '<root folder>\Homer3\UnitTests\Example_RunFuncOffline'
%       Raw_to_DeltaConc_Example('./test.snirf', 2, [1,2]);
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parse args
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Arg 1
[pname, fname] = fileparts(filename);

% Arg 2:  Hb type - {1 -> HbO,  2 -> HbR,  3 -> HbT}
if ~exist('hbType','var') || isempty(hbType)
    hbType   = 1;
end

% Arg 3:
if ~exist('channels','var') || isempty(channels)
    channels = 1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load acquisition file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
acquired = SnirfClass([pname, '/', fname, '.snirf']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Call sequence of user functions to get from acquired data
% to delta concentration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dod = hmrR_Intensity2OD(acquired.data);
dod = hmrR_BandpassFilt(dod, 0.01, 0.50);
dc  = hmrR_OD2Conc_new(dod, acquired.probe, [1,1]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Extract plot data from output of last function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t = dc.time;
y = dc.GetDataTimeSeries('reshape');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PlotData(t, y, hbType, channels);

author: BUNPC


text:
 function inv1 = de_ftab_inv1(fit, s1, varargin)
%function inv1 = de_ftab_inv1(fit, s1, [options])
%|
%| Build object that does 1D inverse of BH function for 1st material component,
%| (usually water), for conventional "water only" beam-hardening correction.
%|
%| in
%|	fit	strum	initialized by de_ftab_fit()
%|	s1	[N]	s_1 sample values for building table (ftab.sls.sl{1})
%|
%| option
%|	'll'	1	which component (default: 1)
%|	'type'	char	'interp' cubic interpolation (default)
%|			'exp'	log(exp()) form - under development
%|	'show'	0|1	plot it?
%|
%| out
%|	inv1	strum
%|	methods:
%|		inv1.fun(hf1)	map log values into corrected values
%|			(corresponding to line-integral of material density)
%|		inv1.plot(fit)	show fit
%|
%| Copyright 2008-09-28, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if nargin == 1 && streq(fit, 'test'), de_ftab_inv1_test, return, end

arg.show = false;
arg.type = 'interp'; % cubic interpolation
arg.ll = 1; % which material, one of {1, ..., LL}
arg = vararg_pair(arg, varargin);

% although we need only the 1D array of s1 values,
% fit.fmfun requires a [n1 1 LL] array.  here the other components are zero.
LL = fit.LL; % usually 1 anyway
sll = s1; sll(1,1,LL) = 0; % [n1 1 L]

% evaluate nonlinear BH function for each of the M spectra
f1 = fit.fmfun(sll); % [n1 1 M]
%f1 = squeeze(f1(:,1,:)); % [n1 M] % changed 2008-12-19
f1 = squeeze(f1); % [n1 M]

st.s1 = s1;
st.f1 = f1;

switch arg.type
case 'interp'
	inv1_fun = @de_ftab_inv1_interp;
case 'exp'
	st = de_ftab_inv1_exp_init(st, fit, arg);
	inv1_fun = @de_ftab_inv1_exp_eval;
otherwise
	fail('unknown type %s', arg.type)
end

meth = {'fun', inv1_fun, '(fh1 [N M]) -> [N M]';
	'plot', @de_ftab_inv1_plot, '(fit)'};

inv1 = strum(st, meth);

if arg.show && im
	inv1.plot(fit);
end

% de_ftab_inv1_exp_init()
function st = de_ftab_inv1_exp_init(st, fit, arg);

s1 = st.s1;
f1 = st.f1;

arg.fit0 = true; % todo
arg.fit0 = false;
arg.chat = true;
arg.thresh = 1e-16;
arg.nexp = 91;

y = exp(s1);
wt = ones(size(y)); % todo
Wh = spdiag(wt, 'nowarn');

MM = size(f1,2);
for mm=1:MM
	f_coef = fit.coef{mm};
	f_exp = fit.mac{mm}(:,1); % for ll=1 only
	emin = 1 / max(f_exp);
	emax = 1 / min(f_exp);
	inv_exp = linspace(emin, emax, arg.nexp)'; % [Ne 1] candidate exps
	Ab = exp(f1(:,mm) * inv_exp'); % [*Ns Ne] "over-complete" basis

	% initial coefficients for each candidate exponent
	if arg.fit0
		fder0 = fit.mac_eff(mm,1);
		inv_der0 = 1 / fder0;
		warg = {'inprodv', inv_exp' / inv_der0};
	else
		warg = {};
	end
	x = wls_simplex(Ab, y(:), Wh, [], warg{:}); % [Ne 1]

	if 1 % reduce to fewer essential terms
		ie = x > arg.thresh; % find key energies
		if arg.chat
			printm('%d coefficients of %d for m=%d', ...
				sum(ie), length(ie), mm)
		end
		inv_exp = inv_exp(ie);

		% final coefficients at key exponents
		if arg.fit0
			fder0 = fit.mac_eff(mm,1);
			inv_der0 = 1 / fder0;
			warg = {'inprodv', inv_exp' / inv_der0}; % shorter now
		else
			warg = {};
		end
		A = Ab(:,ie); % [*Ns Ne] final basis
		x = wls_simplex(A, y(:), Wh, [], warg{:}); % [Ne 1]
	end

	st.exp.coef{mm} = x; % [Ne 1]
	st.exp.exp{mm} = inv_exp;
end

% de_ftab_inv1_exp_eval()
function s1hat = de_ftab_inv1_exp_eval(st, fh1)

MM = size(st.f1,2);
dim = size(fh1);
fh1 = reshapee(fh1, [], MM); % [N M]
s1hat = zeros(size(fh1), 'single');
for mm=1:MM
	tmp = exp(fh1(:,mm) * st.exp.exp{mm}'); % [*N Ne]
	s1hat(:,mm) = log(tmp * st.exp.coef{mm}); % [*N 1]
end
s1hat = reshape(s1hat, dim);

warn 'untested'

% de_ftab_inv1_interp()
% inverse based on interpolation
function s1hat = de_ftab_inv1_interp(st, fh1)

s1 = st.s1;
f1 = st.f1;

MM = size(f1,2);
dim = size(fh1);
fh1 = reshapee(fh1, [], MM);
s1hat = zeros(size(fh1), 'single');
for mm=1:MM
	s1hat(:,mm) = interp1(f1(:,mm), s1, fh1(:,mm), 'pchip', 'extrap');
end
s1hat = reshape(s1hat, dim);

% de_ftab_inv1_plot()
function out = de_ftab_inv1_plot(st, fit)

if nargin ~= 2, fail 'need "fit" argument', end

s1 = st.s1;

s4 = max(s1);
s4 = linspace(-0.0*s4, 1.0*s4, 4*length(s1)+1)'; % [n4] fine sampling
LL = fit.LL;
MM = fit.MM;
sll = s4; sll(1,1,LL) = 0; % [n4 1 L]

f4 = fit.fmfun(sll); % [n4 1 M]
%f4 = squeeze(f4(:,1,:)); % [n4 M]
f4 = squeeze(f4); % [n4 M]

s4i = st.fun(f4); % [n4 M] one correction for each spectrum
err = s4i - repmat(s4, [1 MM]); % [n4 M]

clf, pl = @(n,m) subplot(200 + MM*10 + (n-1)*MM + m);
for mm=1:MM
	pl(1,mm)
	plot(s4, f4(:,mm), 'c-', s4i(:,mm), f4(:,mm), 'y:', ...
		s4, s4 * fit.mac_eff(mm), '--') % monenergetic line
	ir_legend({'true (fine)', 'fit', 'mono'})
	m = num2str(mm);
	axis tight, xlabel 's1', ylabel(['f' m]), title(['inv1 fit m=' m])

	pl(2,mm)
	plot(s4, err(:,mm), '.-')
	axis tight, xlabel 's1', ylabel(['err' m]), title 'inv1 error'
end
prompt

if nargout, out = []; end

% de_ftab_inv1_test()
function de_ftab_inv1_test
stype = 'ps1';
mtype = {'water', 'bone'};
xrs = xray_read_spectra(stype);
mas = xray_read_mac(mtype);
sls = de_ftab_sls;
fm = de_ftab_fm(sls.sll, mas.mac(xrs.en), xrs.Ide);
sl = sls.sl;
fit = de_ftab_fit(sl, fm, 'type', 'exp', 'mtype', mtype, 'show', 0);
itypes = {'exp', 'interp'};
for ii=1:length(itypes)
	itype = itypes{ii};
	inv1 = de_ftab_inv1(fit, sl{1}, 'type', itype);
%	inv1.fun([0 0])
	if im, inv1.plot(fit); end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function lmk = newLmk(Lmk)

% NEWLMK  Get a new landmark.
%   NEWLMK(Lmk) returns the index in structure array Lmk() of the first
%   non-used landmark.
%
%   See also NEWRANGE, INITNEWLMKS.

%   Copyright 2009-2009 Joan sola @ LAAS-CNRS.

lmk = find(~[Lmk.used],1,'first');

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
 function [sinit, bodycoil] = ir_mri_sensemap_init(ykj, varargin)
%function [sinit, bodycoil] = ir_mri_sensemap_init(ykj, varargin)
%|
%| Form initial sensitivity map estimates from surface coils and bodycoil.
%|
%| in
%| ykj [(N) ncoil] surface coil data (2D or 3D)
%|
%| option
%| 'type' initialization method (def: '' = 'ratio')
%|		options: 'median', 'order1', 'order2', ...
%|		if numeric, then simply return as "sinit"
%| 'bodycoil' [(N)] body coil data (def: SSoS of ykj with phase of 1st coil)
%| 'thresh' fraction of max |bodycoil| that is "good" (def: 0.05)
%| 'mask' [(N)] optional mask to pick good values instead of thresholding
%| 'chat' 0|1 verbosity
%|
%| out
%| sinit [(N) ncoil] initial sensitivity map estimates
%| bodycoil [(N)] computed SSoS bodycoil image (if not provided)
%|
%| used in ir_mri_sensemap_admm - see that m-file for details.
%|
%| Michael Allison
%| 2015-08-10 Jeff Fessler revise args, add tests, support 3D

if nargin < 1, ir_usage(), end
if streq(ykj, 'test', 4), ir_mri_sensemap_init_test(ykj), return, end

arg.type = '';
arg.thresh = 0.05;
arg.mask = [];
arg.bodycoil = [];
arg.chat = false;
arg = vararg_pair(arg, varargin);

bodycoil = arg.bodycoil;
if isempty(bodycoil)
	if ndims(ykj) < 3, fail('ykj must be 3D or 4D to infer bodycoil'), end
	ndim = ndims(ykj) - 1; % 2 or 3
	bodycoil = sqrt(sum(abs(ykj).^2, ndim+1)); % sum-of-squares
        tmp = stackpick(ykj, 1); % ykj(:,...,:,1)
        tmp = angle(tmp); % crude estimate of phase of image f_j
        bodycoil = bodycoil .* exp(1i * tmp);
else
	ndim = ndims(bodycoil);
	dim = size(bodycoil);
end

if ~isempty(arg.type) && isnumeric(arg.type) % trivial case
	sinit = arg.type;
return
end

dim = size(bodycoil);
ncoil = numel(ykj) / prod(dim);
if round(ncoil) ~= ncoil, fail('size'), end

if ~isempty(arg.mask)
	jf_equal(size(arg.mask), dim)
	if arg.chat
		warn('using mask instead of threshold to determine good pixels')
	end
	good = mask;
end

if arg.chat
	tell = @(s) printm(s);
else
	tell = @(s) nop(s);
end

sinit = zeros([prod(dim) ncoil], 'single');
for ic = 1:ncoil
	% cannot use stackpick because of case of just one coil
	if ndim == 2
		zj = ykj(:,:,ic);
	elseif ndim == 3
		zj = ykj(:,:,:,ic);
	end

	tmp = div0(zj, bodycoil); % usual ratio

	% determine "good" pixels
	if isempty(arg.mask)
		good = abs(bodycoil) > arg.thresh * max(abs(bodycoil(:)));
	end

	switch arg.type
	case {'', 'ratio'}
		% set all uncertain map values to median of good ones
		tell('Using zero background.')
		tmp(~good) = 0;

	case 'median'
		% set all uncertain map values to median of good ones
		tell('Using median background.')
		tmp1 = median(abs(tmp(good)));
		tmp2 = median(angle(tmp(good))); % dubious
		tmp(~good) = tmp1 .* exp(1i*tmp2);

	case 'avg'
		% set all uncertain map values to mean of good ones
		tell('Using mean background.')
		tmp1 = mean(abs(tmp(good)));
		tmp2 = mean(angle(tmp(good))); % dubious
		tmp(~good) = tmp1 .* exp(1i*tmp2);

	case 'zeros'
		tmp = zeros(dim);

	case 'order1'
		tell('Using 1st-order fit for background.')
		expo = [0 0; 1 0; 0 1];
		tmp = ortho_init(dim, expo, zj, bodycoil, good);

	case 'order2'
		tell('Using 2nd order fit for background.')
		expo = [0 0; 1 0; 0 1; 1 1; 2 0; 0 2];
		tmp = ortho_init(dim, expo, zj, bodycoil, good);

	case 'order3'
		tell('Using 3rd order fit for background.')
		expo = [0 0; 1 0; 0 1; 1 1; 2 0; 2 1; 0 2; 1 2; ...
			3 0; 0 3];
		tmp = ortho_init(dim, expo, zj, bodycoil, good);

	case 'order4'
		tell('Using 4th order fit for background.')
		expo = [0 0; 1 0; 0 1; 1 1; 2 0; 2 1; 0 2; 1 2; ...
			3 0; 0 3; 2 2; 3 1; 1 3; 4 0; 0 4];
		tmp = ortho_init(dim, expo, zj, bodycoil, good);

	otherwise
		fail('Unknown initialization type "%s"', init)
	end

	sinit(:,ic) = single(tmp(:));
end % ic
sinit = reshape(sinit, [dim ncoil]);

end %

function nop(varargin)
end % nop

% ortho_init()
% initialize using orthogonal polynomial basis fit
function init = ortho_init(dim, expo, zj, bodycoil, good)
switch numel(dim)
case 2
	%xx = ndgrid_jf('mat', 0:nx-1, 0:ny-1);
	A = cheby2d(dim(1), dim(2), expo);
case 3
	nx = dim(1);
	ny = dim(2);
	ny = dim(2);
	A = cheby3d(dim(1), dim(2), dim(3), expo);
otherwise
	fail('#dim = %d not done', numel(dim))
end

A = reshape(A, [prod(dim), size(expo,1)]);
maskBodyI = good .* bodycoil;
Aw = repmat(maskBodyI(:), [1 size(expo,1)]) .* A;
theta = pinv(Aw) * zj(:);
clear Aw
init = A * theta; % not Aw
init = reshape(init, size(bodycoil));

end % ortho_init()

% cheby2d()
function polys = cheby2d(nx,ny,expo)
nbases = length(expo);
x = linspace(-1,1,nx);
y = linspace(-1,1,ny);
polys = zeros(nx,ny,nbases);
for n = 1:nbases
	tmp = polyval(ir_chebyshev_poly(expo(n,1)),x);
	tmp = tmp';
	X = repmat(tmp,[1 ny]);
	tmp = polyval(ir_chebyshev_poly(expo(n,2)),y);
	Y = repmat(tmp,[nx 1]);
	polys(:,:,n) = X .* Y;
end

end % cheby2d()

% cheby3d()
function polys = cheby3d(nx,ny,expo)
nbases = size(expo,1);
x = linspace(-1,1,nx);
y = linspace(-1,1,ny);
z = linspace(-1,1,nz);
polys = zeros(nx, ny, nz, nbases, 'single');
for n = 1:nbases
	tmp = polyval(ir_chebyshev_poly(expo(n,1)), x);
	X = repmat(tmp', [1 ny]);
	tmp = polyval(ir_chebyshev_poly(expo(n,2)), y);
	Y = repmat(tmp, [nx 1]);
	tmp = polyval(ir_chebyshev_poly(expo(n,3)), z);
	Z = repmat(tmp, [nx 1]);
	polys(:,:,:,n) = X .* Y .* Z;
end

end % cheby3d()

% ir_mri_sensemap_init_test2()
% 2D
function ir_mri_sensemap_init_test2
nx = 28; ny = 32;
ig = image_geom('nx', nx, 'ny', ny, 'fov', 22);
tmp = ellipse_im(ig, 'shepplogan-emis', 'oversample', 2);
%tmp = ones(nx,ny);
xtrue = tmp .* exp(1i * (tmp-2));
im(abs(xtrue))
ncoil = 4;
strue = ir_mri_sensemap_sim('nx', ig.nx, 'ny', ig.ny, 'dx', ig.dx, ...
	'ncoil', ncoil);
scale = 1 / sqrt(sum(abs(strue(end/2,end/2,:).^2)));
strue = strue * scale; % normalize so SSoS=1 near center
rng(0)
ykj = strue .* repmat(xtrue, [1 1 ncoil]) ...
	+ 2^-8 * (randn(size(strue)) + 1i * randn(size(strue)));
im plc 5 2
clim = minmax(abs(strue))';
im(1, 'row', 1, abs(strue), clim, 'true'), cbar
im(2, 'row', 1, abs(ykj)), cbar
sinit1 = ir_mri_sensemap_init(ykj);
im(3, 'row', 1, abs(sinit1), clim), cbar
slist = {'ratio', 'median', 'avg', 'order1', 'order2', 'order3', 'order4', 'zeros'};
for is=1:numel(slist)
	stype = slist{is};
	sinit2 = ir_mri_sensemap_init(ykj, 'type', stype);
	im(2+is, 'row', 1, abs(sinit2), clim, stype), cbar
	drawnow
end

end % ir_mri_sensemap_init_test2()

% ir_mri_sensemap_init_test3()
% 3D
function ir_mri_sensemap_init_test3
nx = 28; ny = 32; nz=6;
ig = image_geom('nx', nx, 'ny', ny, 'nz', nz, 'fov', 22, 'zfov', 10);
tmp = ellipsoid_im(ig, 'shepp-logan-e3d', 'oversample', 2);
%tmp = ig.ones;
xtrue = tmp .* exp(1i * tmp); % arbitrary phase
%im(abs(xtrue))
nring = 3;
ncoil = 4 * nring;
strue = ir_mri_sensemap_sim('nx', ig.nx, 'ny', ig.ny, 'nz', ig.nz, ...
	'dx', ig.dx, 'nring', nring, 'ncoil', ncoil);
%im('row', 3, abs(strue))
scale = 1 / sqrt(sum(abs(strue(end/2,end/2,end/2,:).^2)));
strue = strue * scale; % normalize so SSoS=1 near center
rng(0)
ykj = strue .* repmat(xtrue, [1 1 1 ncoil]) ...
	+ 1 * 2^-12 * (randn(size(strue)) + 1i * randn(size(strue)));
im plc 4 1
clim = minmax(abs(strue))';
im(1, 'row', nring, abs(strue), clim, 'true'), cbar
titlef('True smaps for 3D [%d,%d,%d] with %d coils (%d rings of %d coils)', ...
	nx, ny, nz, ncoil, nring, ncoil/nring)
im(2, 'row', nring, abs(ykj), 'surface coils'), cbar
sinit1 = ir_mri_sensemap_init(ykj);
im(3, 'row', nring, abs(sinit1), clim, 'ratio'), cbar
slist = {'median'}; % 'avg', 'order1', 'order2', 'order3', 'order4', 'zeros'};
for is=1:numel(slist)
	stype = slist{is};
	sinit2 = ir_mri_sensemap_init(ykj, 'type', stype);
	im(3+is, 'row', nring, abs(sinit2), clim, stype), cbar
	drawnow
end
end % ir_mri_sensemap_init_test3()

% ir_mri_sensemap_init_test()
function ir_mri_sensemap_init_test(arg)
switch arg
case 'test2'
	ir_mri_sensemap_init_test2
case 'test3'
	ir_mri_sensemap_init_test3
case 'test'
	ir_mri_sensemap_init_test2
	if im, prompt, end
	ir_mri_sensemap_init_test3
otherwise
	fail('unknown test "%s"', arg)
end
end % ir_mri_sensemap_init_test()


text:
function gpsdata_all=loadAllGPSData(gpsfile, gpsSE, coordinateStyle)
% load rtklib output of gps solutions, input gpsfile name, gps start and
% end time of gps TOW, gpspostype are as described in readgpsheader
% output: Nx9 gps entries in the session gpsSE, each entry is
% GPS TOW, ECEF XYZ, Q and no of satels and sdx sdy sdz

% if coordinateStyle == lla, the first four columns will be
% GPS TOW, lat lon ellipsoid height.
if nargin < 3
    coordinateStyle = 'ecef';
end

gpsdata_all=zeros( ceil((gpsSE(2)- gpsSE(1)+100)*5), 12);
num_entries=0;
[fgps, gpsdata, gpspostype]=readgpsheader(gpsfile, gpsSE(1));
num_entries=num_entries+1;
gpsdata_all(num_entries, :) = gpsdata';
while(1)
    [fgps, gpsdata]=grabnextgpsdata(fgps, gpspostype);
    if(gpsdata(1)==inf || gpsdata(1)> gpsSE(2))
        break;
    else
        num_entries=num_entries+1;
        gpsdata_all(num_entries, :) = gpsdata';
    end
end
fclose(fgps);
gpsdata_all= gpsdata_all(1:num_entries, :);
if strcmp(coordinateStyle, 'lla')
    gpsdata_all(:, 2:4) = ecef2lla(gpsdata_all(:, 2:4));
    gpsdata_all(:, 2:3) = gpsdata_all(:, 2:3) * pi / 180;
end
end

repo: ekfmonoslam


text:
function plot_strike_dip(catalog,ax)
    % plot strike and dip for a catalog using DipDirection and Dip
    % specify the axes

    dipdir=catalog.DipDirection;
    dip=catalog.Dip;
    x0=catalog.Longitude;
    y0=catalog.Latitude;

    strikescale=0.5;
    dipscale=0.3;
    linecolor = [1 0 0];
    limit_to_axes = true;

    % figure out some scaling stuff
    axun = get(ax,'Units');
    set(ax,'Units','pixels')
    p=fix(get(ax,'Position'));
    set(ax,'Units',axun);

    % TODO: do something with p and the xlim to come up with an appropriately sized symbol

    delete(findobj(gcf,'Tag','diptest'))

    strike = wrapTo360(dipdir-90);
    scaling = ax.DataAspectRatio;
    labels = cellstr( num2str(dip) );
    labels="  "+labels;

    is_horiz = dip==0;
    is_vert = dip==90;

    dx = cosd(strike) .* strikescale ./scaling(2);
    dy = sind(strike) .* strikescale ./ scaling(1);
    xx = ([x0,x0,x0] + [-dx, dx, nan(size(x0))])';
    yy = ([y0,y0,y0] + [-dy, dy, nan(size(y0))])';
    xx=xx(:);
    yy=yy(:);

    hold on;
    plot(ax,xx,yy,'r','linewidth',1,'Tag','diptest','DisplayName','Strikes')

    dipx=cosd(dipdir(:))* dipscale./scaling(2);
    dipy=sind(dipdir(:))* dipscale./scaling(1);
    xdx = ([x0,x0,x0] + [zeros(size(x0)) - is_vert .* dipx, dipx, nan(size(x0))])';
    ydy = ([y0,y0,y0] + [zeros(size(y0)) - is_vert .* dipy, dipy, nan(size(y0))])';

    xdx=xdx(:);
    ydy=ydy(:);
    plot(ax,xdx,ydy,'Color',linecolor,'linewidth',1.5,'Tag','diptest','DisplayName','Dips')
    if any(is_horiz)
        plot(ax,x0(is_horiz),y0(is_horiz),'o','Color',linecolor,'linewidth',1.5,'Tag','diptest','DisplayName','HorizDips')
    end
    if limit_to_axes
        rangeidx=in_range(x0,xlim) & in_range(y0,ylim);
        text(ax,x0(rangeidx), y0(rangeidx), labels(rangeidx), 'VerticalAlignment','top',...
            'HorizontalAlignment','left',...
            'Tag','diptest',...
            'Color',linecolor .* 0.66);

    else
        text(ax,x0, y0, labels, 'VerticalAlignment','top',...
            'HorizontalAlignment','left',...
            'Tag','diptest',...
            'Color',linecolor .* 0.66);
    end
end


text:
% system process model in local s0 frame with dcm formulation, Cs0e is the
% transformation between the sensor local s0 frame and e-frame
% written by Huai
% rqs02e, r s0 in e, q s0 2e, rvqs0, rs in s0, vs in s0, qs0 2s
% acc, m/s^2, acc by imu in s frame, gyro, rad/s, angular rate by imu
% w.r.t i frame coordinated in s frame, dt, time interval for covariance update
function [STM Qd]=sys_local_dcm_v000(rqs02e, rvqs0, acc, gyro, dt, imutype, modelNo)
% the covariance corresponds to states, rs0 in e, q s02e,
% rs in s0, v s in s0, q s02s, ba, bg, sa, sg, qs2c, Ts in c,
% for each group frame, qs02si and Tsi in s0, for each
% point, its inverse depth, here we are only interested up to sa and sg
% for modelNo=3 and imutype=5, MEMS
% acc bias drift, random walk
% gyro bias drift, random walk
% acc scale factor, random walk
% gyro scale factor, random walk

% for modelNo=1 and imutype=5 the acc and gyro turn on bias are removed
R=6317000; % earth average radius
%system disturbance coefs
N=zeros(15,6);
N(13:15,4:6)=eye(3); %attitude
Cs02s=quat2dcm_v000(rvqs0(7:10));
N(10:12,1:3)=Cs02s'; %velocity

%system matrix
A=zeros(15);
% 0s for r s0 in e and q s0 to e
A(7:9, 10:12)=eye(3); % rs in s0
%Velocity
xyz_imu=rqs02e(1:3)+quatrot_v000(rqs02e(4:7),rvqs0(1:3),0);
unit=xyz_imu/norm(xyz_imu,2);
Llh=ecef2geo_v000(xyz_imu,0);
[Rn, Re, g, sL, cL, WIE_E]=geoparam_v000(Llh);
geCoeff=-g*R^2/norm(xyz_imu,2)^3*(eye(3)-3*(unit*unit'))-skew([0;0;WIE_E])*skew([0;0;WIE_E]);
Cs02e=quat2dcm_v000(rqs02e(4:7));
A(10:12,1:3)=Cs02e'*geCoeff;
A(10:12,7:9)=Cs02e'*geCoeff*Cs02e;

Cne=pos2Cne_v000(Llh(1), Llh(2));
ge=Cne(:,3)*g-[xyz_imu(1:2);0]*WIE_E^2; % gravitation in e frame
A(10:12,4:6)=Cs02e'*(-skew(ge)+geCoeff*skew(quatrot_v000(rqs02e(4:7),rvqs0(1:3),0))+...
    skew(skew([0;0;WIE_E])*skew([0;0;WIE_E])*xyz_imu)-2*skew(quatrot_v000(rqs02e(4:7),rvqs0(4:6),0))*skew([0;0;WIE_E]));

A(10:12, 10:12)=-2*skew(quatrot_v000(rqs02e(4:7),[0;0;WIE_E], 1));
A(10:12, 13:15)=-Cs02s'*skew(acc);
%Attitude
A(13:15,4:6)=Cs02s*Cs02e'*skew([0;0;WIE_E]);
A(13:15,13:15)=skew(-gyro);
Anav=A;
Nnav=N;
% X(k+1) = ffun[X(k),U(k),V(k)]
% X(k+1) = Ak*X(k)+Gk*Vk

%%%%Imu error model parameters
[Aimu_d, Qimu_d, Cimu, Rimu]=imu_err_model_v001(acc, gyro, dt, imutype, modelNo);

%%%%Combine and discretize nav and imu models
% this discretization can also be accomplished by Loan's matrix exponential
% method, see sys_metric_phipsi_v000.m
navStates=15;
Anav_d=eye(navStates)+dt*Anav;  %Use 1st order taylor series to discretize Anav
Qnav=Nnav*Rimu*Nnav';
Qnav_d=dt/2*(Anav_d*Qnav+Qnav*Anav_d');      %Use trapezoidal rule to discretize Rimu

STM=zeros(navStates+size(Aimu_d,1));
STM(1:navStates,1:navStates)=Anav_d;
STM(1:navStates,1+navStates:end)=Nnav*Cimu*dt;
STM(1+navStates:end,1+navStates:end)=Aimu_d;

Qd=zeros(navStates+size(Aimu_d,1));
Qd(1:navStates,1:navStates)=Qnav_d;
Qd(1+navStates:end,1+navStates:end)=Qimu_d;
Qd(1:navStates,1+navStates:end)=Nnav*Cimu*Qimu_d*dt/2;
Qd(1+navStates:end,1:navStates)=Qd(1:navStates,1+navStates:end)';
end


text:
function prob = precalcgpstruct(prob)
prob.b = full(prob.b);
try
    % In 2013a they changed the logic in the output. Tell MATLAB to use old
    % style
    [aa,bb,cc] = unique(prob.map,'legacy');
catch
    % If we have an old version, MATLAB would barf at the legacy flag
    [aa,bb,cc] = unique(prob.map);
end

if aa(1)==1
    bb=[0 bb(1:end)'];
end
indsi = [];
indsj = [];
vals = [];
for i = 1:max(prob.map)
    ind = [(bb(i)+1):(bb(i+1))];
    indsj = [indsj ind];
end
indsi = prob.map(prob.map>0);
vals = prob.b(prob.map>0);
prob.B = sparse(indsi,indsj,vals,max(prob.map),size(prob.A,1));

ind = find(prob.map==0);
prob.Afun = prob.A(ind,:);
prob.bfun = prob.b(ind);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [plots,rssi1,rssi2]=ad9361_LTE_2ch(ip,LTEmode,channel)

clc;

%% Pick up LTE parameters according to LTE Mode

if strcmp(LTEmode,'LTE1.4') == 1
    configuration = 'R.4';
    samplingrate = 1.92e6;
    bandwidth = 1.08e6;
    fir_data_file = 'LTE1p4_MHz.ftr';
elseif strcmp(LTEmode,'LTE3') == 1
    configuration = 'R.5';
    samplingrate = 3.84e6;
    bandwidth = 2.7e6;
    fir_data_file = 'LTE3_MHz.ftr';
elseif strcmp(LTEmode,'LTE5') == 1
    configuration = 'R.6';
    samplingrate = 7.68e6;
    bandwidth = 4.5e6;
    fir_data_file = 'LTE5_MHz.ftr';
elseif strcmp(LTEmode,'LTE10') == 1
    configuration = 'R.7';
    samplingrate = 15.36e6;
    bandwidth = 9e6;
    fir_data_file = 'LTE10_MHz.ftr';
else
    error('Please input LTE1.4, LTE3, LTE5 or LTE10.');
end

%% Generate LTE Signal using LTE System Toolbox

% Check for LST presence
if isempty(ver('lte'))
    error('ad9361_LTE:NoLST','Please install LTE System Toolbox to run this example.');
end

% Generate the LTE signal
txsim.RC = configuration; % Base RMC configuration
txsim.NCellID = 17;       % Cell identity
txsim.NFrame = 700;       % Initial frame number
txsim.TotFrames = 1;      % Number of frames to generate
txsim.RunTime = 20;       % Time period to loop waveform in seconds
txsim.DesiredCenterFrequency = 2.45e9; % Center frequency in Hz

% Generate RMC configuration and customize parameters
rmc = lteRMCDL(txsim.RC);
rmc.NCellID = txsim.NCellID;
rmc.NFrame = txsim.NFrame;
rmc.TotSubframes = txsim.TotFrames*10; % 10 subframes per frame
rmc.OCNG = 'On'; % Add noise to unallocated PDSCH resource elements

% Generate RMC waveform
trData = [1;0;0;1]; % Transport data
[eNodeBOutput,~,rmc] = lteRMCDLTool(rmc,trData);
txsim.SamplingRate = rmc.SamplingRate;

% Scale the signal for better power output and cast to int16. This is the
% native format for the SDR hardware. Since we are transmitting the same
% signal in a loop, we can do the cast once to save processing time.
powerScaleFactor = 0.7;
eNodeBOutput = eNodeBOutput.*(1/max(abs(eNodeBOutput))*powerScaleFactor);
eNodeBOutput = int16(eNodeBOutput*2^15);

%% Transmit and Receive using MATLAB libiio

% System Object Configuration
s = iio_sys_obj_matlab; % MATLAB libiio Constructor
s.ip_address = ip;
s.dev_name = 'ad9361';
s.in_ch_no = 4;
s.out_ch_no = 4;
s.in_ch_size = length(eNodeBOutput);
s.out_ch_size = length(eNodeBOutput)*2;

s = s.setupImpl();

% Configure the FIR filter on AD9361
s.writeFirData(fir_data_file);

input = cell(1, s.in_ch_no + length(s.iio_dev_cfg.cfg_ch));
output = cell(1, s.out_ch_no + length(s.iio_dev_cfg.mon_ch));

% Set the attributes of AD9361
input{s.getInChannel('RX_LO_FREQ')} = 2.45e9;
input{s.getInChannel('RX_SAMPLING_FREQ')} = samplingrate;
input{s.getInChannel('RX_RF_BANDWIDTH')} = bandwidth;
input{s.getInChannel('RX1_GAIN_MODE')} = 'slow_attack';
input{s.getInChannel('RX1_GAIN')} = 0;
input{s.getInChannel('RX2_GAIN_MODE')} = 'slow_attack';
input{s.getInChannel('RX2_GAIN')} = 0;
input{s.getInChannel('TX_LO_FREQ')} = 2.45e9;
input{s.getInChannel('TX_SAMPLING_FREQ')} = samplingrate;
input{s.getInChannel('TX_RF_BANDWIDTH')} = bandwidth;

% Keep transmiting and receiving the LTE signal
fprintf('Starting transmission at Fs = %g MHz\n',txsim.SamplingRate/1e6);
for i = 1:5
    fprintf('Transmitting Data Block %i ...\n',i);
    input{1} = real(eNodeBOutput);
    input{2} = imag(eNodeBOutput);
    input{3} = real(eNodeBOutput);
    input{4} = imag(eNodeBOutput);
    output = stepImpl(s, input);
end
fprintf('Transmission finished\n');

% Read the RSSI attributes of both channels
rssi1 = output{s.getOutChannel('RX1_RSSI')};
rssi2 = output{s.getOutChannel('RX2_RSSI')};

s.releaseImpl();

%% Post Processing of Caputured Data

if channel == 1
    I = output{1};
    Q = output{2};
elseif channel == 2
    I = output{3};
    Q = output{4};
else
    error('Please select a channel: 1 or 2.');
end
Rx = I+1i*Q;

% Call LTE Reciever Function
[plots]=LTEReceiver(Rx,samplingrate,configuration);

author: analogdevicesinc


text:
function w = anmp (a)

% normalize angle into the range -pi <= a < +pi.

% ported from NOVAS 3.0

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

w = mod(a, 2.0 * pi);

if (abs(w) >= pi)

    w = w - sign(a) * (2.0 * pi);

end

author: goGPS-Project


text:
%% DEMO_febio_0077_cube_perfect_osmometer_01.m
% Below is a demonstration for:
%
% * Building geometry for a cube with hexahedral elements
% * Defining the boundary conditions
% * Coding the febio structure
% * Running the model
% * Importing and visualizing the displacement and stress results

%% Keywords
%
% * febio_spec version 3.0
% * febio, FEBio
% * perfect osmometer
% * hexahedral elements, hex8
% * cube, box, rectangular
% * static, solid
% * displacement logfile
% * stress logfile

%%

clear; close all; clc;

%% Plot settings
fontSize=20;
faceAlpha1=0.8;
markerSize=40;
markerSize2=35;
lineWidth=3;
cMap=viridis(250); %colormap

%% Control parameters

% Path names
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');

% Defining file names
febioFebFileNamePart='tempModel';
febioFebFileName=fullfile(savePath,[febioFebFileNamePart,'.feb']); %FEB file name
febioLogFileName=[febioFebFileNamePart,'.txt']; %FEBio log file name
febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacement
febioLogFileName_vol=[febioFebFileNamePart,'_vol_out.txt']; %Log file name for exporting stress
febioLogFileName_stress_prin=[febioFebFileNamePart,'_stress_prin_out.txt']; %Log file name for exporting principal stress

%Specifying dimensions and number of elements
cubeSize=1;
sampleWidth=cubeSize; %Width
sampleThickness=cubeSize; %Thickness
sampleHeight=cubeSize; %Height
pointSpacings=0.3*ones(1,3); %Desired point spacing between nodes
numElementsWidth=round(sampleWidth/pointSpacings(1)); %Number of elemens in dir 1
numElementsThickness=round(sampleThickness/pointSpacings(2)); %Number of elemens in dir 2
numElementsHeight=round(sampleHeight/pointSpacings(3)); %Number of elemens in dir 3

%Material parameter set
E_youngs=1;
v_pois=0.3;

iosm=300;
bosm_ini=300;
bosm_diff_amp=200;

% FEA control settings
numTimeSteps=50; %Number of time steps desired
max_refs=25; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=12; %Optimum number of iterations
max_retries=5; %Maximum number of retires
dtmin=(1/numTimeSteps)/100; %Minimum time step size
dtmax=1/numTimeSteps; %Maximum time step size
runMode='internal';

%% Creating model geometry and mesh
% A box is created with tri-linear hexahedral (hex8) elements using the
% |hexMeshBox| function. The function offers the boundary faces with
% seperate labels for the top, bottom, left, right, front, and back sides.
% As such these can be used to define boundary conditions on the exterior.

% Create a box with hexahedral elements
cubeDimensions=[sampleWidth sampleThickness sampleHeight]; %Dimensions
cubeElementNumbers=[numElementsWidth numElementsThickness numElementsHeight]; %Number of elements
outputStructType=2; %A structure compatible with mesh view
[meshStruct]=hexMeshBox(cubeDimensions,cubeElementNumbers,outputStructType);

%Access elements, nodes, and faces from the structure
E=meshStruct.elements; %The elements
V=meshStruct.nodes; %The nodes (vertices)
Fb=meshStruct.facesBoundary; %The boundary faces
Cb=meshStruct.boundaryMarker; %The "colors" or labels for the boundary faces
elementMaterialIndices=ones(size(E,1),1); %Element material indices

%%
% Plotting model boundary surfaces and a cut view

hFig=cFigure;

subplot(1,2,1); hold on;
title('Model boundary surfaces and labels','FontSize',fontSize);
gpatch(Fb,V,Cb,'k',faceAlpha1);
colormap(gjet(6)); icolorbar;
axisGeom(gca,fontSize);

hs=subplot(1,2,2); hold on;
title('Cut view of solid mesh','FontSize',fontSize);
optionStruct.hFig=[hFig hs];
meshView(meshStruct,optionStruct);
axisGeom(gca,fontSize);

drawnow;

%% Defining the boundary conditions
% The visualization of the model boundary shows colors for each side of the
% cube. These labels can be used to define boundary conditions.

%Define supported node sets
logicFace=Cb==1; %Logic for current face set
Fr=Fb(logicFace,:); %The current face set
bcSupportList_X=unique(Fr(:)); %Node set part of selected face

logicFace=Cb==3; %Logic for current face set
Fr=Fb(logicFace,:); %The current face set
bcSupportList_Y=unique(Fr(:)); %Node set part of selected face

logicFace=Cb==5; %Logic for current face set
Fr=Fb(logicFace,:); %The current face set
bcSupportList_Z=unique(Fr(:)); %Node set part of selected face

%Prescribed displacement nodes
logicPrescribe=Cb==6; %Logic for current face set
Fr=Fb(logicPrescribe,:); %The current face set
bcPrescribeList=unique(Fr(:)); %Node set part of selected face

%%
% Visualizing boundary conditions. Markers plotted on the semi-transparent
% model denote the nodes in the various boundary condition lists.

hf=cFigure;
title('Boundary conditions','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;

gpatch(Fb,V,'kw','k',0.5);

hl(1)=plotV(V(bcSupportList_X,:),'r.','MarkerSize',markerSize);
hl(2)=plotV(V(bcSupportList_Y,:),'g.','MarkerSize',markerSize);
hl(3)=plotV(V(bcSupportList_Z,:),'b.','MarkerSize',markerSize);
hl(4)=plotV(V(bcPrescribeList,:),'k.','MarkerSize',markerSize);

legend(hl,{'BC x support','BC y support','BC z support','BC z prescribe'});

axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Defining the FEBio input structure
% See also |febioStructTemplate| and |febioStruct2xml| and the FEBio user
% manual.

%Get a template with default settings
[febio_spec]=febioStructTemplate;

%Set globals/constants
febio_spec.Globals.Constants.R=8.314e-6;
febio_spec.Globals.Constants.T=310;

%febio_spec version
febio_spec.ATTR.version='3.0';

%Module section
febio_spec.Module.ATTR.type='biphasic';

febio_spec.Output.plotfile.var{end+1}.ATTR.type='fluid pressure';
febio_spec.Output.plotfile.var{end+1}.ATTR.type='effective fluid pressure';

%Control section
febio_spec.Control.analysis='STATIC';
febio_spec.Control.time_steps=numTimeSteps;
febio_spec.Control.step_size=1/numTimeSteps;
febio_spec.Control.solver.max_refs=max_refs;
febio_spec.Control.solver.max_ups=max_ups;
febio_spec.Control.time_stepper.dtmin=dtmin;
febio_spec.Control.time_stepper.max_retries=max_retries;
febio_spec.Control.time_stepper.opt_iter=opt_iter;

%Use must-points by specifying
febio_spec.Control.time_stepper=rmfield(febio_spec.Control.time_stepper,'dtmax'); %Remove existing template dtmax definition
febio_spec.Control.time_stepper.dtmax.ATTR.lc=1; %Set load curve id for dtmax
febio_spec.Control.time_stepper.dtmax.VAL=1; %Set value

%Material section
materialName1='Material1';
febio_spec.Material.material{1}.ATTR.name=materialName1;
febio_spec.Material.material{1}.ATTR.type='solid mixture';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.mat_axis.ATTR.type='vector';
febio_spec.Material.material{1}.mat_axis.a=[1 0 0];
febio_spec.Material.material{1}.mat_axis.d=[0 1 0];

febio_spec.Material.material{1}.solid{1}.ATTR.type='perfect osmometer';
febio_spec.Material.material{1}.solid{1}.phiw0=0.8;
febio_spec.Material.material{1}.solid{1}.iosm=iosm;
febio_spec.Material.material{1}.solid{1}.bosm.ATTR.lc=2;
febio_spec.Material.material{1}.solid{1}.bosm.VAL=1;

febio_spec.Material.material{1}.solid{2}.ATTR.type='natural neo-Hookean';
febio_spec.Material.material{1}.solid{2}.E=E_youngs;
febio_spec.Material.material{1}.solid{2}.v=v_pois;

% Mesh section
% -> Nodes
febio_spec.Mesh.Nodes{1}.ATTR.name='Object1'; %The node set name
febio_spec.Mesh.Nodes{1}.node.ATTR.id=(1:size(V,1))'; %The node id's
febio_spec.Mesh.Nodes{1}.node.VAL=V; %The nodel coordinates

% -> Elements
partName1='Part1';
febio_spec.Mesh.Elements{1}.ATTR.name=partName1; %Name of this part
febio_spec.Mesh.Elements{1}.ATTR.type='hex8'; %Element type
febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(E,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E; %The element matrix

% -> NodeSets
nodeSetName1='bcSupportList_X';
nodeSetName2='bcSupportList_Y';
nodeSetName3='bcSupportList_Z';
nodeSetName4='bcPrescribeList';

febio_spec.Mesh.NodeSet{1}.ATTR.name=nodeSetName1;
febio_spec.Mesh.NodeSet{1}.node.ATTR.id=bcSupportList_X(:);

febio_spec.Mesh.NodeSet{2}.ATTR.name=nodeSetName2;
febio_spec.Mesh.NodeSet{2}.node.ATTR.id=bcSupportList_Y(:);

febio_spec.Mesh.NodeSet{3}.ATTR.name=nodeSetName3;
febio_spec.Mesh.NodeSet{3}.node.ATTR.id=bcSupportList_Z(:);

%MeshDomains section
febio_spec.MeshDomains.SolidDomain.ATTR.name=partName1;
febio_spec.MeshDomains.SolidDomain.ATTR.mat=materialName1;

%Boundary condition section
% -> Fix boundary conditions
febio_spec.Boundary.bc{1}.ATTR.type='fix';
febio_spec.Boundary.bc{1}.ATTR.node_set=nodeSetName1;
febio_spec.Boundary.bc{1}.dofs='x';

febio_spec.Boundary.bc{2}.ATTR.type='fix';
febio_spec.Boundary.bc{2}.ATTR.node_set=nodeSetName2;
febio_spec.Boundary.bc{2}.dofs='y';

febio_spec.Boundary.bc{3}.ATTR.type='fix';
febio_spec.Boundary.bc{3}.ATTR.node_set=nodeSetName3;
febio_spec.Boundary.bc{3}.dofs='z';

%LoadData section
% -> load_controller
febio_spec.LoadData.load_controller{1}.ATTR.id=1;
febio_spec.LoadData.load_controller{1}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{1}.interpolate='STEP';
febio_spec.LoadData.load_controller{1}.points.point.VAL=[0 dtmax; 0.25 dtmax; 0.5 dtmax; 0.75 dtmax; 1 dtmax;]; %Defines must-points

febio_spec.LoadData.load_controller{2}.ATTR.id=2;
febio_spec.LoadData.load_controller{2}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{2}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{2}.points.point.VAL=[0 bosm_ini; 0.25 bosm_ini+bosm_diff_amp; 0.5 bosm_ini; 0.75 bosm_ini-bosm_diff_amp; 1 bosm_ini;];

%Output section
% -> log file
febio_spec.Output.logfile.ATTR.file=febioLogFileName;
febio_spec.Output.logfile.node_data{1}.ATTR.file=febioLogFileName_disp;
febio_spec.Output.logfile.node_data{1}.ATTR.data='ux;uy;uz';
febio_spec.Output.logfile.node_data{1}.ATTR.delim=',';

febio_spec.Output.logfile.element_data{1}.ATTR.file=febioLogFileName_vol;
febio_spec.Output.logfile.element_data{1}.ATTR.data='J';
febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';

febio_spec.Output.logfile.element_data{2}.ATTR.file=febioLogFileName_stress_prin;
febio_spec.Output.logfile.element_data{2}.ATTR.data='s1;s2;s3';
febio_spec.Output.logfile.element_data{2}.ATTR.delim=',';

%% Quick viewing of the FEBio input file structure
% The |febView| function can be used to view the xml structure in a MATLAB
% figure window.

%%
% |febView(febio_spec); %Viewing the febio file|

%% Exporting the FEBio input file
% Exporting the febio_spec structure to an FEBio input file is done using
% the |febioStruct2xml| function.

febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode
%system(['gedit ',febioFebFileName,' &']);

%% Running the FEBio analysis
% To run the analysis defined by the created FEBio input file the
% |runMonitorFEBio| function is used. The input for this function is a
% structure defining job settings e.g. the FEBio input file name. The
% optional output runFlag informs the user if the analysis was run
% succesfully.

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.runMode=runMode;
febioAnalysis.maxLogCheckTime=10; %Max log file checking time

[runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!

%% Import FEBio results

if runFlag==1 %i.e. a succesful run

    %%

    % Importing nodal displacements from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp),1,1);

    %Access data
    N_disp_mat=dataStruct.data; %Displacement
    timeVec=dataStruct.time; %Time

    %Create deformed coordinate set
    V_DEF=N_disp_mat+repmat(V,[1 1 size(N_disp_mat,3)]);

    %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations

    DN_magnitude=sqrt(sum(N_disp_mat(:,:,end).^2,2)); %Current displacement magnitude

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    title('Displacement magnitude [mm]','Interpreter','Latex')
    hp=gpatch(Fb,V_DEF(:,:,end),DN_magnitude,'k',1,2); %Add graphics object to animate
    hp.Marker='.';
    hp.MarkerSize=markerSize2;
    hp.FaceColor='interp';
    gpatch(Fb,V,0.5*ones(1,3),'none',0.25); %A static graphics object

    axisGeom(gca,fontSize);
    colormap(cMap); colorbar;
    caxis([0 max(DN_magnitude)]); caxis manual;
    axis(axisLim(V_DEF)); %Set axis limits statically
    view(140,30);
    camlight headlight;

    % Set up animation features
    animStruct.Time=timeVec; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments
        DN_magnitude=sqrt(sum(N_disp_mat(:,:,qt).^2,2)); %Current displacement magnitude

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp hp]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData'}; %Properties of objects to animate
        animStruct.Set{qt}={V_DEF(:,:,qt),DN_magnitude}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    drawnow;

    %%
    % Importing element stress from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_vol),1,1);

    %Access data
    E_J_mat=dataStruct.data;

    %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations

    [CV]=faceToVertexMeasure(E,V,E_J_mat(:,:,end));

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure  /usr/local/MATLAB/R2020a/bin/glnxa64/jcef_helper: symbol lookup error: /lib/x86_64-linux-gnu/libpango-1.0.so.0: undefined symbol: g_ptr_array_copy

    gtitle([febioFebFileNamePart,': Press play to animate']);
    title('$J$ [.]','Interpreter','Latex')
    hp=gpatch(Fb,V_DEF(:,:,end),CV,'k',1,2); %Add graphics object to animate
    hp.Marker='.';
    hp.MarkerSize=markerSize2;
    hp.FaceColor='interp';
    gpatch(Fb,V,0.5*ones(1,3),'none',0.25); %A static graphics object

    axisGeom(gca,fontSize);
    colormap(cMap); colorbar;
    caxis([min(E_J_mat(:)) max(E_J_mat(:))]);
    axis(axisLim(V_DEF)); %Set axis limits statically
    view(140,30);
    camlight headlight;

    % Set up animation features
    animStruct.Time=timeVec; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments

        [CV]=faceToVertexMeasure(E,V,E_J_mat(:,:,qt));

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp hp]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData'}; %Properties of objects to animate
        animStruct.Set{qt}={V_DEF(:,:,qt),CV}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    drawnow;

    %%
%     % Importing element stresses from a log file
%     dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_stress_prin),1,1);
%
%     %Access data
%     E_stress_prin_mat=dataStruct.data;
%
%     %Compute pressure
%     P = squeeze(sum(sum(E_stress_prin_mat,1),2));%squeeze(-1/3*mean(sum(E_stress_prin_mat,2),1));
%
%     %%
%     % Visualize pressure-stretch curve
%
%     cFigure; hold on;
%     title('Pressure-time curve','FontSize',fontSize);
%     xlabel('$t$ [s]','FontSize',fontSize,'Interpreter','Latex');
%     ylabel('$p$ [MPa]','FontSize',fontSize,'Interpreter','Latex');
%
%     plot(timeVec(:),P(:),'r-','lineWidth',lineWidth);
%
%     view(2); axis tight;  grid on; axis square; box on;
%     set(gca,'FontSize',fontSize);
%     drawnow;

end

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*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