text:
%% Apparently it must be here!
% Written by Ali Akbar Eftekhari
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or
% without modification, are permitted provided that the following
% conditions are met:
%
%     *   Redistributions of source code must retain the above copyright notice,
%         this list of conditions and the following disclaimer.
%     *   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.

clc; clear; close all;
%% create a mesh and visualize it
% The first part of this demo shows you how to create a mesh and visualize
% it. The aim of the visualization is to introduce you to the mesh
% structure.
% Here, we are going to create and visualize a 1D uniform equidistant mesh
% and show the cell centers with a 'o' marker and the face of each cell
% with a '+' marker.
% We start by defining the length of the domain:
L = 1.0; % length of the domain
%%
% Then we define the number of cells in the domain:
%
Nx = 10; % number of cells in the domain
%%
% Now we call createMesh1D, one of the functions from a group of createMesh*
% functions:
%
m = createMesh1D(Nx, L);
%%
% This function generates a structure which keeps the domain and grid
% information, which will be used by almost every other function in the
% FVMtool.
%
% In general, 1D, 2D, 3D, 1D radial (axial symetry), and 2D cylindrical
% grids can be constructed. For more information, type
%%
%
%   help createMesh2D
%   help createMesh3D
%   help createMeshCylindrical1D
%   help createMeshCylindrical2D
%
%%
% Now, let's have a look at the structure of the created mesh:
%
disp(m);
%%
% You can easily get more information about this structure by typing
%%
%   help createMesh1D
%
%%
% As an example, we use the position of the cell centers to visualize the
% domain, the cells, and the interfaces between them:
%
figure(2);
plot(m.cellcenters.x, ones(size(m.cellcenters.x)), 'or', ...
     m.facecenters.x, ones(size(m.facecenters.x)), '-+b');
legend('cell centers', 'face centers');
title('Visualization of a 1D discretized domain');
%%
% The generated figure shows a 1D domain of length 1.0 [unit], discretized
% into 10 cells of the same size (dx = 1/10).
% I'm going to assume that you are familiar with the cell-centred finite
% volume method, discretization, and especially handling the boundary
% conditions using the ghost cells.
% You can also generate and visualize 2D and 3D meshes. To find some
% examples, type:
%
%   help createMesh2D;
%   help createMesh3D;
%
% You will end up with the following two simple examples:
%
Nx = 5;
Ny = 7;
Lx = 10;
Ly = 20;
m = createMesh2D(Nx, Ny, Lx, Ly);
[X, Y] = ndgrid(m.cellcenters.x, m.cellcenters.y);
[Xf,Yf]=ndgrid(m.facecenters.x, m.facecenters.y);
figure(3);
plot(X, Y, 'or', ...
     Xf, Yf, '-b', Xf', Yf', '-b');
%%
% that shows you a 2D grid, and
%
Nx = 2;
Lx = 1.0;
Ny = 3;
Ly = 2.0;
Nz = 4;
Lz = 3.0;
m = createMesh3D(Nx, Ny, Nz, Lx, Ly, Lz);
[X, Y, Z] = ndgrid(m.cellcenters.x, m.cellcenters.y, m.cellcenters.z);
[Xf, Yf, Zf] = ndgrid(m.facecenters.x, m.facecenters.y, m.facecenters.z);
figure(4);
plot3(X(:), Y(:), Z(:), 'or')
hold on;
plot3(Xf(:), Yf(:), Zf(:), '+b')
legend('cell centers', 'cell corners');
hold off;
%%
% that shows a 3D grid.
%% Boundary condition structure
% I had so many reasons to write this toy toolbox, of which, the most
% important one was to be able to implement different boundary conditions
% in the most convenient way! My final implementation makes the user able
% to define either a _periodic boundary condition_ or a _general boundary
% condition_ of the following form:
%%
%
% $$a (\nabla \phi .\mathbf{n}) + b \phi = c $$
%
%%
% In the above equation, $\phi$ is the unknown, and _a_, _b_, and _c_ are
% constants. In practice, this boundary condition equation will be discretized
% to the following system of algebraic equations:
%%
%
% $$M_{bc} \phi = {RHS}_{bc}$$
%
%%
% By adjusting the values of _a_ and _b_, onne can easily define
% one of the following well-known boundary conditions:
%%
%
% * Neumann (_a_ is nonzero; _b_ is 0)
% * Dirichlet (_a_ is zero; _b_ is nonzero)
% * Robin (_a_ and _b_ are both nonzero)
%
%%
% To clarify the above explanations, let us create a boundary condition
% structure for a 1D mesh.
%
Nx = 10; % number of cells in the domain
Lx = 1.0; % length of the domain
m = createMesh1D(Nx, Lx); % createMesh and createMesh are identical
BC = createBC(m); % creates a boundary condition structure
disp(BC); % display the BC structure
%%
% The BC structure has two substructures, i.e., _left_ that denotes the boundary
% at the left side of the domain (at x=0), and _right_ that denotes the
% boundary at the end of the domain (at x=Lx). Each of these substructures
% have three fields, i.e., _a_, _b_, and _c_. The default values are _a_=1,
% _b_=0, and _c_=0:
%
disp(BC.left); % show the values of the coefficients for the left boundary
%%
% There is one other field, i.e., _periodic_, which has a zero value. If
% you change it to one, then a periodic boundary condition will be created
% and the _a_, _b_, and _c_ values will be ignored.
% clearly, you can change the above boundary condition by assigning new
% values to the |BC.left| fields. For instance, you can define a Dirichlet
%  boundary (i.e., fixed value) with a value of 2.5 by typing
%
BC.left.a = 0;
BC.left.b = 1;
BC.left.c = 2.5;
%%
% You can define a periodic boundary condition simply by writing:
%
BC.left.periodic = 1;
%%
% For boundary condition structures created for 2D and 3D grids, we will
% have _left_, _right_, _bottom_, _top_, _back_, and _front_ boundaries
% and thus substructures. Let me show them to you in action:
%
m = createMesh2D(3,4, 1.0, 2.0);
BC = createBC(m);
disp(BC);
disp(BC.top);
%%
% Yes, that's right. _a_, _b_, and _c_ are vectors. It means that you can
% have different boundary conditions for
% different cell faces at each boundary. For instance, I can have a Neumann
% boundary condition for the first cell and a Dirichlet boundary condition
% for the last cell at the top boundary:
BC.top.a(1) = 1; BC.top.b(1) = 0; BC.top.c(1) = 0; % zero value Neumann
BC.top.a(end) =0; BC.top.b(end)=1; BC.top.c(end) = 0; % zero value Dirichlet
disp('top.  a     b     c'); % some fancy display!
disp('   ---------------');
disp([BC.top.a BC.top.b BC.top.c]);
%%
% The same procedure can be followed for a 3D grid. However, _a_, _b_, and
% _c_ values are 2D matrices for a 3D grid. I will discuss it in more
% details when we reach the practical examples.
% *Important note:* If you need to assign a boundary condition to the
% entire boundary, use (:) in your assignment. For instance, to define a
% Dirichlet boundary for the right boundary, you may write
BC.right.a(:)=0; BC.right.b(:)=1; BC.right.c(:)=0;
%% Solve a diffusion equation
% As the first example, we solve a steady-state diffusion equation of the
% following formform
%
% $$\nabla.\left(-D\nabla c\right)=0,$$
%
%%
% where _D_ is the diffusivity and _c_ is the concentration. Let me assume
% that we have a 1D domain, with Dirichlet boundary conditions at both
% boundaries, i.e.,
% at x=0, c=1; and at x=L, c=0.
% First of all, we need to define our domain, discretize it, and define the
% boundaries at the borders.
clc; clear; % clear the screen and memory
L = 0.01; % a 1 cm domain
Nx = 10; % number of cells
m = createMesh3D(Nx,Nx,Nx,L, L,L); % create the mesh
BC = createBC(m); % construct the BC structure (Neumann by default)
%%
% Now as you may remeber, we have to switch from Neumann to Dirichlet
% boundary conditions
%
BC.left.a(:) = 0; BC.left.b(:) = 1; BC.left.c(:) = 1; % Left boundary to Dirichlet
BC.right.a(:) = 0; BC.right.b(:) = 1; BC.left.c(:) = 1; % right boundary to Dirichlet
%%
% The next sep is to define the diffusivity coefficient. In this FVTool,
% the physical properties of the domain are defined for each cell, with the
% function createCellVariable
D = createCellVariable(m, 1e-5); % assign a constant value of 1e-5 to diffusivity value on each cell
%%
% However, the transfer coefficients must be known on the face of each cell.
% For this reason, we have a few averaging schemes implemented in the
% Utilities folder. For a 1D domain, we can use a harmonic mean scheme:
D_face = harmonicMean(D); % average diffusivity value on the cell faces
%%
% Now, we can convert the PDE to a algebraic system of linear equations,
% i.e.,
%%
%
% $$\nabla.\left(-D\nabla c\right) \approx Mc = 0$$
%
%%
% where M is the matrix of coefficient that is going to be calculated using
% this toolbox. The matrix of coefficient, _M_ has two parts. The diffusion
% equation and the boundary conditions. They are calculated by:
M_diff = diffusionTerm(D_face); % matrix of coefficients for diffusion term
[M_bc, RHS_bc] = boundaryCondition(BC); % matrix of coefficient and RHS vector for the boundary condition
%%
% A vector of right hand side values are always obtained during the
% discretization of the boundary conditions.
% Now that the PDE is discretized, we can solve it by a Matlab linear solver.
c = solvePDE(m, M_diff+M_bc, RHS_bc);
%%
% finally, the resut can be visualized:
visualizeCells(c);
%%
% Just to get excited a little bit, only change the mesh definition command
% from createMesh1D(Nx,L) to createMesh2D(Nx,Nx,L,L), run the code and see
% what happens. For even more excitement, change it to
% createMesh3D(Nx,Nx,Nx,L,L,L).
% This is usually the way we develop new mathematical models for a physical
% phenomenon. Write the equation, solve it in 1D, compare it to the
% analytical solution, then solve it numerically in 2D and 3D for more
% realistic cases with heterogeneous transfer coefficients and other
% nonidealities (and perhaps compare it to some experimental data)
%% Solve a convection-diffuison equation and compare it to analytical solution
% Here, I'm going to add a convection term to what we solved in the
% previous example. This tutorial is adopted from the fipy
% convection-diffusion example you can find at this address:
%%
% <http://www.ctcms.nist.gov/fipy/examples/convection/index.html>
%%
% The differential equation reads
%%
%
% $$\nabla.\left(\mathbf{u} \phi -D\nabla \phi \right)=0$$
%
%%
% Here, $\mathbf{u}$ is a velocity vector (face variable) and $D$ is the
% diffusion coefficient (again a face variable). Please see the PDF
% document for an explanation of cell and face variables. We use Dirichlet
% (constant value) boundary conditions on the left and right boundaries.
% It is zero at the left boundary and one at the right boundary. The
% analytical solution of this differential equation reads
%%
%
% $$c = \frac{1-exp(ux/D)}{1-exp(uL/D)}$$
%
%%
% We start the code as always with some cleaning up:
clc; clear;
%%
% Then we define the domain and mesh size:
L = 1;  % domain length
Nx = 25; % number of cells
meshstruct = createMesh1D(Nx, L);
x = meshstruct.cellcenters.x; % extract the cell center positions
%%
% The next step is to define the boundary condition:
BC = createBC(meshstruct); % all Neumann boundary condition structure
BC.left.a = 0; BC.left.b=1; % switch the left boundary to Dirichlet
BC.left.c=0; % value = 0 at the left boundary
BC.right.a = 0; BC.right.b=1; % switch the right boundary to Dirichlet
BC.right.c=1; % value = 1 at the right boundary
%%
% Now we define the transfer coefficients:
D_val = 1.0; % diffusion coefficient value
D = createCellVariable(meshstruct, D_val); % assign dif. coef. to all the cells
Dave = harmonicMean(D); % convert a cell variable to face variable
u = -10; % velocity value
u_face = createFaceVariable(meshstruct, u); % assign velocity value to cell faces
%%
% Now we discretize the differential equation into a system of linear
% algebraic equations:
%%
%
% $$(M_{conv}-M_{diff}+M_{bc})\phi={RHS}_{bc}$$
%
%%
% or if we use an upwind discretization scheme, we will obtain:
%%
%
% $$(M_{conv,uw}-M_{diff}+M_{bc})\phi={RHS}_{bc}$$
%
%%
%
Mconv =  convectionTerm(u_face); % convection term, central, second order
Mconvupwind = convectionUpwindTerm(u_face); % convection term, upwind, first order
Mdiff = diffusionTerm(Dave); % diffusion term
[Mbc, RHSbc] = boundaryCondition(BC); % boundary condition discretization
M = Mconv-Mdiff+Mbc; % matrix of coefficient for central scheme
Mupwind = Mconvupwind-Mdiff+Mbc; % matrix of coefficient for upwind scheme
RHS = RHSbc; % right hand side vector
c = solvePDE(meshstruct, M, RHS); % solve for the central scheme
c_upwind = solvePDE(meshstruct, Mupwind, RHS); % solve for the upwind scheme
c_analytical = (1-exp(u*x/D_val))/(1-exp(u*L/D_val)); % analytical solution
figure(5);
plot(x, c.value(2:Nx+1), x, c_upwind.value(2:Nx+1), '--',...
    x, c_analytical, '.');
legend('central', 'upwind', 'analytical');
%%
% As you see here, we obtain a more accurate result by using a central
% difference discretization scheme for the convection term compared to the
% first order upwind.
%% solve a transient diffusion equation
% This tutorial is adapted from the fipy 1D diffusion example
%%
% <http://www.ctcms.nist.gov/fipy/examples/diffusion/index.html FiPy diffusion tutorial>
% The transient diffusion equation reads
%%
%
% $$\alpha\frac{\partial c}{\partial t}+\nabla.\left(-D\nabla c\right)=0,$$
%
% where $c$ is the independent variable (concentration, temperature, etc)
% , $D$ is the diffusion coefficient, and $\alpha$ is a constant.
%%
% Once again, clean up:
clc; clear;
%%
% Define the domain and create a mesh structure
L = 50;  % domain length
Nx = 20; % number of cells
m = createMesh1D(Nx, L);
x = m.cellcenters.x; % cell centers position
%%
% Create the boundary condition structure:
BC = createBC(m); % all Neumann boundary condition structure
%%
% Switch the left and right boundaries to Dirichlet:
BC.left.a = 0; BC.left.b=1; BC.left.c=1; % left boundary
BC.right.a = 0; BC.right.b=1; BC.right.c=0; % right boundary
%%
% Define the transfer coefficients:
D_val = 1;
D = createCellVariable(m, D_val);
Dave = harmonicMean(D); % convert it to face variables
% Define alfa, the coefficient of the transient term:
alfa_val = 1;
alfa = createCellVariable(m, alfa_val);
%%
% Define the initial values:
c_init = 0;
c_old = createCellVariable(m, c_init, BC); % initial values
c = c_old; % assign the old value of the cells to the current values
%%
% Now define the time step and the final time:
dt = 0.1; % time step
final_t = 100;
%%
% Here, we first define the matrices of coefficients that will not change
% as we progress in time, viz. diffusion term and boundary condition:
Mdiff = diffusionTerm(Dave);
[Mbc, RHSbc] = boundaryCondition(BC);
%%
% The transitionTerm function gives a matrix of coefficient and a RHS
% vector. The matrix of coefficient does not change in each time step, but
% the RHS does (see the PDF documents). Therefore, we need to call the
% function inside the time loop.
% Start the loop here:
for t=dt:dt:final_t
    [M_trans, RHS_trans] = transientTerm(c_old, dt, alfa);
    M = M_trans-Mdiff+Mbc;
    RHS = RHS_trans+RHSbc;
    c = solvePDE(m,M, RHS);
    c_analytical = 1-erf(x/(2*sqrt(D_val*t)));
    c_old = c;
end
%%
% Now visualize the final results
figure(6)
plot(x, c.value(2:Nx+1), 'o', x, c_analytical);
xlabel('Length [m]'); ylabel('c');
legend('Numerical', 'Analytical');
%%
% you can visualize the results from each time step by moving the plot line
% inside the for loop.
%% convection equations; different discretization schemes
% If I want to highlight one special feature of this FVTool, I will point
% a finger on its various discretization schemes for a linear convection
% term, which includes central difference (second order), upwind (first
% order), and TVD scheme with various flux limiters.
%%
% Here, we are going to compare the performance of each scheme for solving
% two PDE's.
% First, a simple linear transient convection equation with an strange initial
% condition and later, we solve the well-known Burger's equation.
clc; clear;
% define a 1D domain and mesh
W = 1;
Nx = 500;
mesh1 = createMesh1D(Nx, W);
x = mesh1.cellcenters.x;
% define the boundaries
BC = createBC(mesh1); % all Neumann
BC.left.periodic=1;
BC.right.periodic =1;
% Initial values
phi_old = createCellVariable(mesh1, 0.0, BC);
phi_old.value(20:120) = 1;
phi_old.value(180:400)= sin(x(180:400)*10*pi());
% initial guess for phi
phi = phi_old;
phiuw_old=phi_old;
% initial values for upwind scheme
phiuw = phi;
% keep the initial values for visualization
phiinit=phi_old;
% velocity field
u = 0.3;
uf = createFaceVariable(mesh1, u);
% diffusion field
D = 1e-2;
Df = createFaceVariable(mesh1, D);
% transient term coefficient
alfa = createCellVariable(mesh1,1.0);
% upwind convection term
Mconvuw = convectionUpwindTerm1D(uf);
% define the BC term
[Mbc, RHSbc] = boundaryCondition(BC);
% choose a flux limiter
FL = fluxLimiter('Superbee');
% solver
dt = 0.001; % time step
final_t = W/u;
t = 0;
while t<final_t
    t = t+dt;
    % inner loop for TVD scheme
    for j = 1:5
        [Mt, RHSt] = transientTerm(phi_old, dt, alfa);
        [Mconv, RHSconv] = convectionTvdTerm1D(uf, phi, FL);
        M = Mconv+Mt+Mbc;
        RHS = RHSt+RHSbc+RHSconv;
        phi = solvePDE(mesh1, M, RHS);
    end
    [Mtuw, RHStuw] = transientTerm(phiuw_old, dt, alfa);
    Muw = Mconvuw+Mtuw+Mbc;
    RHSuw = RHStuw+RHSbc;
    phiuw = solvePDE(mesh1, Muw, RHSuw);
    phiuw_old = phiuw;
    phi_old = phi;
end
figure(7);plot(x, phiinit.value(2:Nx+1), x, phi.value(2:Nx+1), '-o', x, ...
        phiuw.value(2:Nx+1));

% %% method of lines: using Matlab's ODE solvers for adaptive time stepping
%
%
%
% %% solving a nonlinear PDE
%
% %% solving a system of linear PDE's: sequential and coupled methods
%
% %% solving a system of nonlinear PDE's: sequential and coupled olutions
%
% %% Real life cases: water-flooding in the production of oil
%
% %% and finally your examples?
%
%
%

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function loc_2d_shapes(NumOfSamples)
% Create Localized 2-d data from uniform distribution.
% Data are bount to the rectangle between points A(0,6) and B(4,10),
% a circle with center at C(7.5,2.5) and radius R = 2,
% and an isosceles triangle defined by the intesection points
% of the lines y=6, y=2x-5 and y=-2x+25.
% Every local set has the same number of samples (NumOfSamples/4).

clc;
input_dims = 2;

% Initialize Sample counters.
k1 = 0; k2 = 1; k3 = 0; k4 = 0;

% Initialize Data
Data = [7.5; 2.5;];

while size(Data,2)<NumOfSamples

In1 = 10*rand(input_dims,1);

% Apply Mask to input Data
if (In1(1,1)<4)&&(In1(1,1)>0)&&(In1(2,1)<10)&&(In1(2,1)>6)&&(k1<=NumOfSamples/4)
    Data = [Data In1];
    k1 = k1 + 1;
elseif norm(In1-[7.5; 2.5;])<=2&&(k2<=NumOfSamples/4)
    Data = [Data In1];
    k2 = k2 + 1;
elseif (In1(2,1)>6)&&(In1(2,1)<2*In1(1,1)-5)&&(In1(2,1)<-2*In1(1,1)+25)&&(k3<=NumOfSamples/4)
    Data = [Data In1];
    k3 = k3 + 1;
elseif (abs(In1(2,1)-In1(1,1))<.1)&&(In1(2,1)<5)&&(In1(2,1)>.5)&&(In1(1,1)<5)&&(In1(1,1)>.5)&&(k4<=NumOfSamples/4)
    Data = [Data In1];
    k4 = k4 + 1;
end
end

a = randperm(NumOfSamples);
Data = Data(:,a);

save local_uniform_2d Data;
plot(Data(1,:),Data(2,:),'.','MarkerSize',1);
grid on;
clear;

author: Sable

repo: mcbench-benchmarks


text:
%converts frequency to bark
%>
%> @param fInHz: frequency
%> @param cModel: 'Schroeder','Terhardt', 'Zwicker', or 'Traunmuller'
%>
%> @retval bark value
% ======================================================================
function [bark] = ToolFreq2Bark(fInHz, cModel)

    if (nargin < 2)
        cModel = 'Schroeder';
    end

    % set function handle
    hPitchFunc = str2func (['aca' cModel '_I']);

    bark = hPitchFunc(fInHz);
end

function [bark] = acaSchroeder_I(f)
    bark    = 7 * asinh(f/650);
end

function [bark] = acaTerhardt_I(f)
    bark    = 13.3 * atan(0.75 * f/1000);
end

function [bark] = acaZwicker_I(f)
    bark    = 13 * atan(0.76 * f/1000) + 3.5 * atan(f/7500);
end

function [bark] = acaTraunmuller_I(f)
    bark    = 26.81/(1+1960./f) - 0.53;
end

author: alexanderlerch

repo: ACA-Code


text:
function y = cropbyval(x, val)
%CROPBYVAL Crop an array by
%
%   Y = CROPBYVAL(X, VAL) returns a subarray identical to X but where any
%   values VAL on the borders of X have been removed.
%
%   Y = CROPBYVAL(X) will determine VAL automatically by looking at the
%   borders of X and picking the most common value.  En error is given if
%   there is no value which is more common that the other values on the
%   border.  Note that automatically chosing VAL slowes down CROPBYVAL.
%
%   For example, if
%
%     X = [ 0  0  0  0
%           1  2  3  0
%           4  5  6  0 ]
%
%   then both Y = CROPBYVAL(X, 0) and Y = CROPBYVAL(X) will return
%
%     Y = [ 1  2  3
%           4  5  6 ]

%   Author:      Peter J. Acklam
%   Time-stamp:  2002-03-03 13:51:33 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   % check number of input arguments
   error(nargchk(1, 2, nargin));

   % if VAL is give, check it, otherwise compute default value
   if nargin == 2
      if any(size(val) ~= 1)
         error('VAL must be a scalar.');
      end
   else
      vals = [];
      sx = size(x);
      dx = ndims(x);
      c = cell(dx, 1);
      for i = 1:dx
         c{i} = 2:sx(i)-1;
      end
      t(prod(sx)) = logical(uint8(1));          % use uint8 to save memory
      t = reshape(t, sx);
      t(:) = t(end);
      t(c{:}) = 0;
      [lens, vals] = rlencode(sort(x(t(:))));
      len = max(lens);
      k = find(len == lens);
      if length(k) > 1
         error('No unique most common value on boundary.');
      end
      val = vals(k);
   end

   % get linear index values of elements different from VAL
   if isnan(val)
      i = find(~isnan(x));
   elseif val == 0
      i = find(x);
   else
      i = find(x ~= val);
   end

   % quick exit if output will be empty
   if isempty(i)
      y = [];
      return
   end

   % The following is based on code by
   % Doug Schwarz <douglas.schwarz@kodak.com>

   sx     = size(x);            % size vector of x
   dx     = length(sx);         % number of dimensions of x
   c      = cell(dx, 1);        % initialize list of subscript
   [c{:}] = ind2sub(sx, i);     % compute the subscripts
   for j = 1:dx
      c{j} = min(c{j}) : max(c{j});
   end
   y = x(c{:});

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
%% cunique
% Below is a demonstration of the features of the |cunique| function

%%
clear; close all; clc;

%% Syntax
% |[A_uni,ind1,ind2,Ac]=cunique(A);|

%% Description
% The |imx| function provides a figure window based GUI for 3D image
% segmentation

%% Examples

%%
% Plot settings
fontSize=20;

%% Example 1: Getting unique entries and occurance counts for 1xN arrays

n=15;
A=round(25*rand(1,n)); %Rounded random set in range 0-25
A(1)=A(end); %Force at least one double occurance for this example
A

%Get unique set and counts
[A_uni,ind1,ind2,Ac]=cunique(A)

%% Example 2: Getting unique entries and occurance counts for NxM arrays

n=5;
m=6;
A=round(25*rand(n,m)); %Rounded random set in range 0-25
A(1)=A(end); %Force at least one double occurance for this example
A

%Get unique set and counts
[A_uni,ind1,ind2,Ac]=cunique(A)

%%
% Visualizing input array and occurange counts

cFigure;
subplot(1,2,1);
title('The input array')
hold on;
imagesc(A);
image_numeric(A,[],0,fontSize);
axis tight; axis equal;
colormap(gca,gjet(max(A(:))));
icolorbar;

subplot(1,2,2);
title('The occurance counts')
hold on;
imagesc(Ac);
image_numeric(Ac,[],0,fontSize);
axis tight; axis equal;
colormap(gca,gjet(max(Ac(:))));
icolorbar;
drawnow;

%% Example 3: Getting unique entries and occurance counts for NxMx... arrays

n=3;
m=4;
l=2;

A=round(25*rand(n,m,l)); %Rounded random set in range 0-25
A(1)=A(end); %Force at least one double occurance for this example
A

%Get unique set and counts
[A_uni,ind1,ind2,Ac]=cunique(A)

%% Example 4: Using 'rows' option

n=5;
m=3;

A=round(25*rand(n,m)); %Rounded random set in range 0-25
A(1,:)=A(end,:); %Force at least one double row for this example
A

%Get unique set and counts
[A_uni,ind1,ind2,Ac]=cunique(A,'rows')

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


text:
function [imageSeg,Rot]= getPlanSeg(XYZworldframeTest,SpaceTest,imgSize)
% input XYZworldframeTest is alinged in Z
%load([ '/n/fs/modelnet/NYUdataSet/NYUdatafeatureNew/' 'feature_' num2str(imageNum) '.mat'],'XYZworldframeTest','rgbTest','SpaceTest','imgDepth');
onPlaneThreshold =0.055;
sizethr =500;
normalAgreeThreshold =0.8;
removeNaN = find(~sum(isnan(XYZworldframeTest),2));
XYZworldframeTest = XYZworldframeTest(removeNaN,:);
[zhist,zc]=hist(XYZworldframeTest(:,3),round(20*(SpaceTest.Rz(2)-SpaceTest.Rz(1))));
[~,zpks] =findpeaks([0,zhist,0],'minpeakdistance',3);
zpks = zpks-1;

range =[max(min(XYZworldframeTest(:,3)),zc(max(1,zpks-1))'),...
        min(zc(min(size(zc,2),zpks+1)),max(XYZworldframeTest(:,3)))';
        zc(end-2),max(XYZworldframeTest(:,3));
        min(XYZworldframeTest(:,3)),zc(3)];
removepts = inrange(XYZworldframeTest(:,3),range);
XYZworldframeTestRemove=XYZworldframeTest;
XYZworldframeTestRemove(removepts>0,:) =[];

% hough lines
[H,T,R] = votelines(XYZworldframeTestRemove(:,1),XYZworldframeTestRemove(:,2));
% blur it
allangleVote = max(H);
allangleVote = [allangleVote,allangleVote];
shift= 0;
clear rightAngle
for i =1:length(shift)
    rightAngle(i,:) = allangleVote(91+shift(i):180+shift(i))+allangleVote(1:90);
end
[~,AnglepickedLin]= max(rightAngle(:));
[s,Anglepicked]=ind2sub(size(rightAngle),AnglepickedLin);
[~,Rpicked1] = max(H(:,Anglepicked));
[~,Rpicked2] = max(H(:,Anglepicked+90+shift(s)));
P = [R(Rpicked1),T(Anglepicked);R(Rpicked2),T(Anglepicked)+90+shift(s)];

rotation = min(P(:,2));
Rot = getRotationMatrix('z',-1*rotation/180*pi);
XYZworldframeTestNew = [[Rot(1:2,1:2)*XYZworldframeTest(:,[1,2])']', XYZworldframeTest(:,3)];
normals = points2normals(XYZworldframeTestNew);

% find  corner
[xhist,xc]=hist(XYZworldframeTestNew(:,1),round((SpaceTest.Rx(2)-SpaceTest.Rx(1))/0.05));
[~,xpks] =findpeaks([0,xhist,0],'MINPEAKHEIGHT',50);
xpks = xpks-1;
[yhist,yc]=hist(XYZworldframeTestNew(:,2),round((SpaceTest.Ry(2)-SpaceTest.Ry(1))/0.05));
[~,ypks] =findpeaks([0,yhist,0],'MINPEAKHEIGHT',50);
ypks = ypks-1;
[zhist,zc]=hist(XYZworldframeTestNew(:,3),round((SpaceTest.Rz(2)-SpaceTest.Rz(1))/0.05));
[~,zpks] =findpeaks([0,zhist,0],'MINPEAKHEIGHT',50);
zpks = zpks-1;

%project points onto this plan and caculate conected component

gid =1;
imageSeg = zeros(imgSize);
for p =1:size(xpks,2)
    inLine =find(abs(XYZworldframeTestNew(:,1)-xc(xpks(p)))<onPlaneThreshold&abs(normals(1,:)')>normalAgreeThreshold);
    conthisSegMask = zeros(imgSize);
    conthisSegMask(removeNaN(inLine)) =1;
    label = bwlabel(conthisSegMask,8);
    unique_label = unique(label);
    unique_label(unique_label==0)=[];
    for i =1:length(unique_label),
         if sum(label(:)==unique_label(i))>sizethr,
            imageSeg(label==unique_label(i))=gid;
            gid = gid+1;
         end
    end
end

for p =1:size(ypks,2)
    inLine =find(abs(XYZworldframeTestNew(:,2)-yc(ypks(p)))<onPlaneThreshold&abs(normals(2,:)')>normalAgreeThreshold);
    conthisSegMask = zeros(imgSize);
    conthisSegMask(removeNaN(inLine)) =1;
    label = bwlabel(conthisSegMask,8);
    unique_label = unique(label);
    unique_label(unique_label==0)=[];
    for i =1:length(unique_label),
         if sum(label(:)==unique_label(i))>sizethr,
            imageSeg(label==unique_label(i))=gid;
            gid = gid+1;
         end
    end
end

for p =1:size(zpks,2)
    inLine =find(abs(XYZworldframeTestNew(:,3)-zc(zpks(p)))<onPlaneThreshold&abs(normals(3,:)')>normalAgreeThreshold);
    conthisSegMask = zeros(imgSize);
    conthisSegMask(removeNaN(inLine)) =1;
    label = bwlabel(conthisSegMask,8);
    unique_label = unique(label);
    unique_label(unique_label==0)=[];
    for i =1:length(unique_label),
         if sum(label(:)==unique_label(i))>sizethr,
            imageSeg(label==unique_label(i))=gid;
            gid = gid+1;
         end
    end
end
%%

% rotate back
%{
corner_r = [xc(xpks(1)) yc(ypks(1));xc(xpks(end)) yc(ypks(end))];
corner_r = get4corner(corner_r);
corner = [[Rot(1:2,1:2)'*corner_r(:,[1:2])']'];
D1 = cos(P(1,2)*pi/180)*corner(1,1)+sin(P(1,2)*pi/180)*corner(1,2);
D2 = cos(P(2,2)*pi/180)*corner(1,1)+sin(P(2,2)*pi/180)*corner(1,2);
D3 = cos(P(1,2)*pi/180)*corner(3,1)+sin(P(1,2)*pi/180)*corner(3,2);
D4 = cos(P(2,2)*pi/180)*corner(3,1)+sin(P(2,2)*pi/180)*corner(3,2);
P_new = [D1,P(1,2);D2,P(2,2);D3,P(1,2);D4,P(2,2)];
f = figure,
vis_point_cloud(XYZworldframeTest,rgbTest,10,5000);hold on;
for j =1:3
    plot3(corner([j,j+1],1),corner([j,j+1],2),[max(XYZworldframeTest(:,3));max(XYZworldframeTest(:,3))],'-xr','LineWidth',10)
end
plot3(corner([1,4],1),corner([1,4],2),[max(XYZworldframeTest(:,3));max(XYZworldframeTest(:,3))],'-xr','LineWidth',10)

for i =1:gid
    hold on;
    plot3(XYZworldframeTestOrg(imageSeg(:)==i,1),XYZworldframeTestOrg(imageSeg(:)==i,2),XYZworldframeTestOrg(imageSeg(:)==i,3),'+','Color',rand([1,3]));
end

axis equal;
axis tight;
view(30,50)
saveas(f,['./result/' num2str(imageNum) '.fig']);
saveas(f,['./result/' num2str(imageNum) '.jpg']);

for gid =1:length(gropuind)
    imageSeg(removeNaN(gropuind{gid})) = gid;
end
%}

%figure(1),imagesc(imageSeg)
%Boudary = [D1 D2 D3 D4];
%{
if imageNum> 0,
    im = getImagesc(imageSeg);
    mkdir(segpath);
    imwrite(im,sprintf('%s/%04d.jpg',segpath,imageNum));
end
%}
end
function [hough_transform,T,R] = votelines(X,Y)
         thetaResolution = 1;
         rhoResolution = 0.1;
         T = [1:180];
         theta = T/180*pi;
         rho = X(:)*cos(theta)+ Y(:)*sin(theta);
         % quantize the rho
         rhoNorm = max(1,round((rho-min(rho(:)))/rhoResolution));
         R =[1:max(rhoNorm(:))]*rhoResolution+min(rho(:));
         hough_transform = zeros(max(rhoNorm(:)),size(T,2));
         TT = repmat(T,[size(rhoNorm,1),1]);
         hough_transform = accumarray([rhoNorm(:),TT(:)],1,size(hough_transform));
end

author: thusiyuan


text:
function [kden,N,K] = density_dir(CIJ)
%DENSITY        Density
%
%   kden = density_dir(CIJ);
%   [kden,N,K] = density_dir(CIJ);
%
%   Density is the fraction of present connections to possible connections.
%
%   Input:      CIJ,    directed (weighted/binary) connection matrix
%
%   Output:     kden,   density
%               N,      number of vertices
%               K,      number of edges
%
%   Notes:  Assumes CIJ is directed and has no self-connections.
%           Weight information is discarded.
%
%
%   Olaf Sporns, Indiana University, 2002/2007/2008

N = size(CIJ,1);
K = nnz(CIJ);
kden = K/(N^2-N);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function boolean = CalculateTimeForTakingOutBillet(T_new,reference_Temperature,x_intervals,y_intervals)
boolean = 0;
count = 0;

for x_index = 2:1:(x_intervals-1)
    for y_index = 2:1:(y_intervals-1)
        if(abs(T_new(x_index,y_index))-reference_Temperature < 01 )
        count = 1;
        end
    end
    if(count >= ((x_intervals*y_intervals)/2))
        boolean = 1;
        break;
    end
end

end


text:
function [ a, done, change ] = vec_gray_next ( n, base, a, done )

%*****************************************************************************80
%
%% VEC_GRAY_NEXT computes the elements of a product space.
%
%  Discussion:
%
%    The elements are produced one at a time.
%
%    This routine handles the case where the number of degrees of freedom may
%    differ from one component to the next.
%
%    A method similar to the Gray code is used, so that successive
%    elements returned by this routine differ by only a single element.
%
%    The routine uses internal static memory.
%
%  Example:
%
%    N = 2, BASE = ( 2, 3 ), DONE = TRUE
%
%     A    DONE  CHANGE
%    ---  -----  ------
%    0 0  FALSE    1
%    0 1  FALSE    2
%    0 2  FALSE    2
%    1 2  FALSE    1
%    1 1  FALSE    2
%    1 0  FALSE    2
%    1 0   TRUE   -1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dennis Stanton, Dennis White,
%    Constructive Combinatorics,
%    Springer, 1986,
%    ISBN: 0387963472,
%    LC: QA164.S79.
%
%  Parameters:
%
%    Input, integer N, the number of components.
%
%    Input, integer BASE(N), contains the number of degrees of
%    freedom of each component.  The output values of A will
%    satisfy 0 <= A(I) < BASE(I).
%
%    Input, integer A(N).  On the first call, the input value
%    of A doesn't matter.  Thereafter, it should be the same as
%    its output value from the previous call.
%
%    Input, logical DONE.  On the first call, the user must
%    set DONE to TRUE.  Thereafter, DONE should be set to the output
%    value of DONE on the previous call.
%
%    Output, integer A(N), the next element of the space.
%
%    Output, logical DONE.  If DONE is FALSE, the program has computed
%    another entry, which is contained in A.  If DONE is TRUE,
%    then there are no more entries.
%
%    Output, integer CHANGE, is set to the index of the element whose
%    value was changed.  On return from the first call, CHANGE
%    is 1, even though all the elements have been "changed".  On
%    return with DONE equal to TRUE, CHANGE is -1.
%
  persistent active;
  persistent dir;
%
%  The user is calling for the first time.
%
  if ( done )

    done = 0;
    a(1:n) = 0;

    dir(1:n) = 1;
    active(1:n) = 1;

    for i = 1 : n

      if ( base(i) < 1 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'VEC_GRAY_NEXT - Warning!\n' );
        fprintf ( 1, '  For index I = %d\n', i );
        fprintf ( 1, '  the nonpositive value of BASE(I) = %d\n',  base(i) );
        fprintf ( 1, '  which was reset to 1!\n' );
        base(i) = 1;
        active(i) = 0;
      elseif ( base(i) == 1 )
        active(i) = 0;
      end

    end

    change = 1;

    return

  end
%
%  Find the maximum active index.
%
  change = -1;

  for i = 1 : n
    if ( active(i) ~= 0 )
      change = i;
    end
  end
%
%  If there are NO active indices, we have generated all vectors.
%
  if ( change == -1 )
    done = 1;
    return
  end
%
%  Increment the element with maximum active index.
%
  a(change) = a(change) + dir(change);
%
%  If we attained a minimum or maximum value, reverse the direction
%  vector, and deactivate the index.
%
  if ( a(change) == 0 | a(change) == base(change) - 1 )
    dir(change) = -dir(change);
    active(change) = 0;
  end
%
%  Activate all subsequent indices.
%
  for i = change + 1 : n
    if ( 1 < base(i) )
      active(i) = 1;
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function I_s = segImageRegion(I,S)

I_s = I;
I_s1 = I(:,:,1);
I_s2 = I(:,:,2);
I_s3 = I(:,:,3);

for i=1:max(S(:))

    mask = S==i;

%     I_s(mask) = mean(I(mask));
    I_s1(mask) = mean(I_s1(mask));
    I_s2(mask) = mean(I_s2(mask));
    I_s3(mask) = mean(I_s3(mask));

    I_s(:,:,1) = I_s1;
    I_s(:,:,2) = I_s2;
    I_s(:,:,3) = I_s3;
end

end

author: andrewssobral

repo: lrslibrary


text:
function l = invgammaPriorLogProb(prior, x)

% INVGAMMAPRIORLOGPROB Log probability of inverse gamma prior.

% PRIOR

% PRIOR

% Compute log prior
D = length(x);
l = D*prior.a*log(prior.b)-D*gammaln(prior.a)-(prior.a+1)*sum(log(x))-prior.b*sum(1./x);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef Dynamics < handle

% Copyright (C) 2013, Georgia Tech Research Corporation
% see the LICENSE file included with this software

    properties

    end

    methods
        function pose_t_1 = apply_dynamics(obj, pose_t, dt)

	end
    end

end

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
function [lf] = eeg_halfspace_medium_leadfield(rd, elc, vol)

% HALFSPACE_MEDIUM_LEADFIELD calculate the halfspace medium leadfield
% on positions pnt for a dipole at position rd and conductivity cond
% The halfspace solution requires a plane dividing a conductive zone of
% conductivity cond, from a non coductive zone (cond = 0)
%
% [lf] = halfspace_medium_leadfield(rd, elc, cond)

% Copyright (C) 2011, Cristiano Micheli and Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: eeg_halfspace_medium_leadfield.m 2781 2011-02-03 10:48:53Z roboos $

siz = size(rd);
if any(siz==1)
  % positions are specified as a single vector
  Ndipoles = prod(siz)/3;
  rd = rd(:)'; % ensure that it is a row vector
elseif siz(2)==3
  % positions are specified as a Nx3 matrix -> reformat to a single vector
  Ndipoles = siz(1);
  rd = rd';
  rd = rd(:)'; % ensure that it is a row vector
else
  error('incorrect specification of dipole locations');
end

Nelc     = size(elc,1);
lf       = zeros(Nelc,3*Ndipoles);

for i=1:Ndipoles
  % this is the position of dipole "i"
  dip1 = rd((1:3) + 3*(i-1));

  % distances electrodes - dipole
  r1 = elc - ones(Nelc,1) * dip1;

  % Method of mirror dipoles:
  % Defines the position of mirror dipoles being symmetric to the plane
  dip2 = get_mirror_pos(dip1,vol);

  % distances electrodes - mirror dipole
  r2 = elc - ones(Nelc,1) * dip2;

  % denominator
  R1 =  (4*pi*vol.cond) * (sum(r1' .^2 ) .^ 1.5)';
  % denominator, mirror term
  R2 = -(4*pi*vol.cond) * (sum(r2' .^2 ) .^ 1.5)';

  % condition of dipoles falling in the non conductive halfspace
  invacuum = acos(dot(vol.ori,(dip1-vol.pnt)./norm(dip1-vol.pnt))) < pi/2;

  if invacuum
    warning('dipole lies on the vacuum side of the plane');
    lf(:,(1:3) + 3*(i-1)) = NaN(Nelc,3);
  elseif any(R1)==0
    warning('dipole coincides with one of the electrodes');
    lf(:,(1:3) + 3*(i-1)) = NaN(Nelc,3);
  else
    lf(:,(1:3) + 3*(i-1)) = (r1 ./ [R1 R1 R1]) + (r2 ./ [R2 R2 R2]);
  end
end

author: PatternRecognition

repo: OpenBMI


text:
function sMat = in_bst_matrix(MatFile, varargin)
% IN_BST_MATRIX: Read a "matrix" file in Brainstorm format.
%
% USAGE:  sMat = in_bst_matrix(MatFile, FieldsList) : Read the specified fields
%         sMat = in_bst_matrix(MatFile)             : Read all the fields
%
% INPUT:
%    - MatFile    : Absolute or relative path to the file to read
%    - FieldsList : List of fields to read from the file
% OUTPUT:
%    - sMat : Brainstorm matrix file structure

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

%% ===== PARSE INPUTS =====
% Get protocol folders
ProtocolInfo = bst_get('ProtocolInfo');
% Filename: Relative / absolute
if ~file_exist(MatFile)
    MatFile = bst_fullfile(ProtocolInfo.STUDIES, MatFile);
    if ~file_exist(MatFile)
        error('File not found.');
    end
end

% Specific fields
if (nargin < 2)
    % Read all fields
    sMat = load(MatFile);
    % Get all the fields
    FieldsToRead = fieldnames(sMat);
    FieldsToRead{end+1} = 'ZScore';
else
    % Get fields to read
    FieldsToRead = varargin;
    % Add other necessary fields
    if ismember('Value', FieldsToRead)
        FieldsToRead{end+1} = 'ZScore';
    end
    % When reading Leff, make sure nAvg is read as well
    if ismember('Leff', FieldsToRead) && ~ismember('nAvg', FieldsToRead)
        FieldsToRead{end+1} = 'nAvg';
    end
    % Read each field only once
    FieldsToRead = unique(FieldsToRead);
    % Read specified files only
    warning off MATLAB:load:variableNotFound
    sMat = load(MatFile, FieldsToRead{:});
    warning on MATLAB:load:variableNotFound
end

%% ===== FILL OTHER MISSING FIELDS =====
for i = 1:length(FieldsToRead)
    if ~isfield(sMat, FieldsToRead{i}) || isempty(sMat.(FieldsToRead{i}))
        switch(FieldsToRead{i})
            case 'nAvg'
                sMat.(FieldsToRead{i}) = 1;
            case 'Leff'
                if isfield(sMat, 'nAvg') && ~isempty(sMat.nAvg)
                    sMat.Leff = sMat.nAvg;
                else
                    sMat.Leff = 1;
                end
            otherwise
                sMat.(FieldsToRead{i}) = [];
        end
    end
end

%% ===== APPLY DYNAMIC ZSCORE =====
% DEPRECATED
% Check for structure integrity
if ismember('ZScore', FieldsToRead) && ~isempty(sMat.ZScore) && (~isfield(sMat.ZScore, 'mean') || ~isfield(sMat.ZScore, 'std') || ~isfield(sMat.ZScore, 'abs') || ~isfield(sMat.ZScore, 'baseline') || isempty(sMat.ZScore.abs))
    sMat.ZScore = [];
end
% Apply to file values
if ismember('ZScore', FieldsToRead) && ismember('Value', FieldsToRead) && ~isempty(sMat.ZScore) && ~isempty(sMat.Value)
    sMat.Value = process_zscore_dynamic('Compute', sMat.Value, sMat.ZScore);
    sMat = rmfield(sMat, 'ZScore');
end

% ===== FIX TRANSPOSED TIME VECTOR =====
if isfield(sMat, 'Time') && (size(sMat.Time,1) > 1)
    sMat.Time = sMat.Time';
end

author: brainstorm-tools


text:
% This file is part of TREEQSM.
%
% TREEQSM 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.
%
% TREEQSM 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 TREEQSM.  If not, see <http://www.gnu.org/licenses/>.

function Pass = filtering(P,inputs)

% ---------------------------------------------------------------------
% FILTERING.M       Filters noise from point clouds.
%
% Version 3.0.0
% Latest update     3 May 2022
%
% Copyright (C) 2013-2022 Pasi Raumonen
% ---------------------------------------------------------------------

% Filters the point cloud as follows:
%
% 1) the possible NaNs are removed.
%
% 2) (optional, done if filter.k > 0) Statistical kth-nearest neighbor
% distance outlier filtering based on user defined "k" (filter.k) and
% multiplier for standard deviation (filter.nsigma): Determines the
% kth-nearest neighbor distance for all points and then removes the points
% whose distances are over average_distance + nsigma*std. Computes the
% statistics for each meter layer in vertical direction so that the
% average distances and SDs can change as the point density decreases.
%
% 3) (optional, done if filter.radius > 0) Statistical point density
% filtering based on user defined ball radius (filter.radius) and multiplier
% for standard deviation (filter.nsigma): Balls of radius "filter.radius"
% centered at each point are defined for all points and the number of
% points included ("point density") are computed and then removes the points
% whose density is smaller than average_density - nsigma*std. Computes the
% statistics for each meter layer in vertical direction so that the
% average densities and SDs can change as the point density decreases.
%
% 4) (optional, done if filter.ncomp > 0) Small component filtering based
% on user defined cover (filter.PatchDiam1, filter.BallRad1) and threshold
% (filter.ncomp): Covers the point cloud and determines the connected
% components of the cover and removes the points from the small components
% that have less than filter.ncomp cover sets.
%
% 5) (optional, done if filter.EdgeLength > 0) cubical downsampling of the
% point cloud based on user defined cube size (filter.EdgeLength):
% selects randomly one point from each cube
%
% Does the filtering in the above order and thus always applies the next
% fitering to the point cloud already filtered by the previous methods.
% Statistical kth-nearest neighbor distance outlier filtering and the
% statistical point density filtering are meant to be exlusive to each
% other.
%
% Inputs:
% P         Point cloud
% inputs    Inputs structure with the following subfields:
%   filter.EdgeLength   Edge length of the cubes in the cubical downsampling
%   filter.k            k of knn method
%   filter.radius       Radius of the balls in the density filtering
%   filter.nsigma       Multiplier for standard deviation, determines how
%                         far from the mean the threshold is in terms of SD.
%                         Used in both the knn and the density filtering
%   filter.ncomp        Threshold number of components in the small
%                         component filtering
%   filter.PatchDiam1   Defines the patch/cover set size for the component
%                         filtering
%   filter.BallRad1     Defines the neighbors for the component filtering
%   filter.plot         If true, plots the filtered point cloud
% Outputs:
% Pass      Logical vector indicating points passing the filtering
% ---------------------------------------------------------------------

% Changes from version 2.0.0 to 3.0.0, 3 May 2022:
% Major changes and additions.
% 1) Added two new filtering options: statistical kth-nearest neighbor
%    distance outlier filtering and cubical downsampling.
% 2) Changed the old point density filtering, which was based on given
%    threshold, into statistical point density filtering, where the
%    threshold is based on user defined statistical measure
% 3) All the input parameters are given by "inputs"-structure that can be
%    defined by "create_input" script
% 4) Streamlined the coding and what is displayed

%% Initial data processing
% Only double precision data
if ~isa(P,'double')
  P = double(P);
end
% Only x,y,z-data
if size(P,2) > 3
  P = P(:,1:3);
end
np = size(P,1);
np0 = np;
ind = (1:1:np)';
Pass = false(np,1);

disp('----------------------')
disp(' Filtering...')
disp(['  Points before filtering:  ',num2str(np)])

%% Remove possible NaNs
F = ~any(isnan(P),2);
if nnz(F) < np
  disp(['  Points with NaN removed:  ',num2str(np-nnz(Pass))])
  ind = ind(F);
end

%% Statistical kth-nearest neighbor distance outlier filtering
if inputs.filter.k > 0
  % Compute the knn distances
  Q = P(ind,:);
  np = size(Q,1);
  [~, kNNdist] = knnsearch(Q,Q,'dist','euclidean','k',inputs.filter.k);
  kNNdist = kNNdist(:,end);

  % Change the threshold kNNdistance according the average and standard
  % deviation for every vertical layer of 1 meter in height
  hmin = min(Q(:,3));
  hmax = max(Q(:,3));
  H = ceil(hmax-hmin);
  F = false(np,1);
  ind = (1:1:np)';
  for i = 1:H
    I = Q(:,3) < hmin+i & Q(:,3) >= hmin+i-1;
    points = ind(I);
    d = kNNdist(points);
    J = d < mean(d)+inputs.filter.nsigma*std(d);
    points = points(J);
    F(points) = 1;
  end
  ind = ind(F);
  disp(['  Points removed as statistical outliers:  ',num2str(np-length(ind))])
end

%% Statistical point density filtering
if inputs.filter.radius > 0
  Q = P(ind,:);
  np = size(Q,1);

  % Partition the point cloud into cubes
  [partition,CC] = cubical_partition(Q,inputs.filter.radius);

  % Determine the number of points inside a ball for each point
  NumOfPoints = zeros(np,1);
  r1 = inputs.filter.radius^2;
  for i = 1:np
    if NumOfPoints(i) == 0
      points = partition(CC(i,1)-1:CC(i,1)+1,CC(i,2)-1:CC(i,2)+1,CC(i,3)-1:CC(i,3)+1);
      points = vertcat(points{:,:});
      cube = Q(points,:);
      p = partition{CC(i,1),CC(i,2),CC(i,3)};
      for j = 1:length(p)
        dist = (Q(p(j),1)-cube(:,1)).^2+(Q(p(j),2)-cube(:,2)).^2+(Q(p(j),3)-cube(:,3)).^2;
        J = dist < r1;
        NumOfPoints(p(j)) = nnz(J);
      end
    end
  end

  % Change the threshold point density according the average and standard
  % deviation for every vertical layer of 1 meter in height
  hmin = min(Q(:,3));
  hmax = max(Q(:,3));
  H = ceil(hmax-hmin);
  F = false(np,1);
  ind = (1:1:np)';
  for i = 1:H
    I = Q(:,3) < hmin+i & Q(:,3) >= hmin+i-1;
    points = ind(I);
    N = NumOfPoints(points);
    J = N > mean(N)-inputs.filter.nsigma*std(N);
    points = points(J);
    F(points) = 1;
  end
  ind = ind(F);
  disp(['  Points removed as statistical outliers:  ',num2str(np-length(ind))])
end

%% Small component filtering
if inputs.filter.ncomp > 0
  % Cover the point cloud with patches
  input.BallRad1 = inputs.filter.BallRad1;
  input.PatchDiam1 = inputs.filter.PatchDiam1;
  input.nmin1 = 0;
  Q = P(ind,:);
  np = size(Q,1);
  cover = cover_sets(Q,input);

  % Determine the separate components
  Components = connected_components(cover.neighbor,0,inputs.filter.ncomp);

  % The filtering
  B = vertcat(Components{:}); % patches in the components
  points = vertcat(cover.ball{B}); % points in the components
  F = false(np,1);
  F(points) = true;
  ind = ind(F);
  disp(['  Points with small components removed:  ',num2str(np-length(ind))])
end

%% Cubical downsampling
if inputs.filter.EdgeLength > 0
  Q = P(ind,:);
  np = size(Q,1);
  F = cubical_downsampling(Q,inputs.filter.EdgeLength);
  ind = ind(F);
  disp(['  Points removed with downsampling:  ',num2str(np-length(ind))])
end

%% Define the output and display summary results
Pass(ind) = true;
np = nnz(Pass);
disp(['  Points removed in total: ',num2str(np0-np)])
disp(['  Points removed in total (%): ',num2str(round((1-np/np0)*1000)/10)])
disp(['  Points left: ',num2str(np)])

%% Plot the filtered and unfiltered point clouds
if inputs.filter.plot
  plot_comparison(P(Pass,:),P(~Pass,:),1,1,1)
  plot_point_cloud(P(Pass,:),2,1)
end

author: InverseTampere

repo: TreeQSM


text:
function  [cent, varargout]=FastPeakFind(d, threshold, filt ,edg, fid)
% Analyze noisy 2D images and find peaks to 1 pixel accuracy.
% The code is desinged to be as fast as possible, so I kept it pretty basic.
% The code assumes that the peaks are relatively sparse, test whether there
% is too much pile up and set threshold or user defined filter accordingly.
%
% Inputs:
%   d           The 2D data raw image - assumes a Double\Single-precision
%               floating-point, uint8 or unit16 array. Please note that the code
%               casts the raw image to uint16 if needed.  If the image dynamic range is
%               between 0 and 1, I multiplied to fit uint16. This might not be optimal for
%               generic use, so modify according to your needs.
%   threshold   A number between 0 and max(raw_image(:)) to remove  background
%   filt        A filter matrix used to smooth the image. The filter size
%               should correspond the characteristic size of the peaks
%   edg         A number>1 for skipping the first few and the last few 'edge' pixels
%   fid         In case the user would like to save the peak positions to
%               a file, the code assumes a "fid = fopen([filename], 'w+');" line
%               in the script that uses this function.
%
% Optional Outputs:
%   cent        a 1xN vector of coordinates of peaks (x1,y1,x2,y2,...
%   [cent cm]   in addition to cent, cm is a binary matrix  of size(d)
%               with 1's for peak positions.
%
%   Example:
%
%   p=FastPeakFind(image);
%   imagesc(image); hold on
%   plot(p(2:2:end),p(1:2:end),'r+')
%
%
%   Natan (nate2718281828@gmail.com)
%   Ver 1.61 , Date: June 5th 2013
%
%
%% defaults
if (nargin < 1)
    d=uint16(conv2(reshape(single( 2^14*(rand(1,1024*1024)>0.99995) ),[1024 1024]) ,fspecial('gaussian', 15,3),'same')+2^8*rand(1024));
    imagesc(d);
end

if ndims(d)>2 %I added this in case one uses imread (JPG\PNG\...).
    d=uint16(rgb2gray(d));
end

if isfloat(d) %For the case the input image is double, casting to uint16 keeps enough dynamic range while speeds up the code.
    if max(d(:))<=1
        d =  uint16( d.*2^16./(max(d(:))));
    else
        d = uint16(d);
    end
end

if (nargin < 2)
    threshold = (max([min(max(d,[],1))  min(max(d,[],2))])) ;
end

if (nargin < 3)
    filt = (fspecial('gaussian', 7,1)); %if needed modify the filter according to the expected peaks sizes
end

if (nargin < 4)
    edg =3;
end

if (nargin < 5)
    savefileflag = false;
else
    savefileflag = true;
end

%% Analyze image
if any(d(:))  ; %for the case of non zero raw image

    d = medfilt2(d,[3,3]);

    % apply threshold
    if isa(d,'uint8')
    d=d.*uint8(d>threshold);
    else
    d=d.*uint16(d>threshold);
    end

    if any(d(:))   ; %for the case of the image is still non zero

        % smooth image
        d=conv2(single(d),filt,'same') ;

        % Apply again threshold (and change if needed according to SNR)
        d=d.*(d>0.9*threshold);

        % peak find - using the local maxima approach - 1 pixle resolution
        % d will be noisy on the edges, since no hits are expected there anyway we'll skip 'edge' pixels.
        sd=size(d);
        [x y]=find(d(edg:sd(1)-edg,edg:sd(2)-edg));

        % initialize outputs
          cent=[];%
          cent_map=zeros(sd);

        x=x+edg-1;
        y=y+edg-1;
        for j=1:length(y)
            if (d(x(j),y(j))>=d(x(j)-1,y(j)-1 )) &&...
                    (d(x(j),y(j))>d(x(j)-1,y(j))) &&...
                    (d(x(j),y(j))>=d(x(j)-1,y(j)+1)) &&...
                    (d(x(j),y(j))>d(x(j),y(j)-1)) && ...
                    (d(x(j),y(j))>d(x(j),y(j)+1)) && ...
                    (d(x(j),y(j))>=d(x(j)+1,y(j)-1)) && ...
                    (d(x(j),y(j))>d(x(j)+1,y(j))) && ...
                    (d(x(j),y(j))>=d(x(j)+1,y(j)+1));

                  %All these alternatives were slower...
                %if all(reshape( d(x(j),y(j))>=d(x(j)-1:x(j)+1,y(j)-1:y(j)+1),9,1))
                %if  d(x(j),y(j)) == max(max(d((x(j)-1):(x(j)+1),(y(j)-1):(y(j)+1))))
                %if  d(x(j),y(j))  == max(reshape(d(x(j),y(j))  >=  d(x(j)-1:x(j)+1,y(j)-1:y(j)+1),9,1))

                        cent = [cent ; x(j) ; y(j)];
                        cent_map(x(j),y(j))=cent_map(x(j),y(j))+1; % if a binary matrix output is desired

            end
        end

        if savefileflag
            % previous version used dlmwrite, which can be slower than  fprinf
            %             dlmwrite([filename '.txt'],[cent],   '-append', ...
            %                 'roffset', 0,   'delimiter', '\t', 'newline', 'unix');+

            fprintf(fid, '%f ', cent(:));
            fprintf(fid, '\n');

        end

    else % in case image after threshold is all zeros
        cent=[];
        cent_map=zeros(size(d));
        if nargout>1 ;  varargout{1}=cent_map; end
        return
    end

else % in case raw image is all zeros (dead event)
    cent=[];
    cent_map=zeros(size(d));
    if nargout>1 ;  varargout{1}=cent_map; end
    return
end

%demo mode - no input to the function
if (nargin < 1); colormap(bone);hold on; plot(cent(2:2:end),cent(1:2:end),'rs');hold off; end

% return binary mask of centroid positions if asked for
if nargout>1 ;  varargout{1}=cent_map; end


text:
function varargout = process_baseline_norm2( varargin )
% PROCESS_BASELINE_NORM2: Normalization with respect to a baseline (A=baseline).

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Baseline normalization (A=baseline)';
    sProcess.FileTag     = @GetFileTag;
    sProcess.Category    = 'Filter2';
    sProcess.SubGroup    = 'Standardize';
    sProcess.Index       = 205;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/SourceEstimation#Z-score';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'results', 'timefreq', 'matrix'};
    sProcess.OutputTypes = {'data', 'results', 'timefreq', 'matrix'};
    sProcess.nInputs     = 2;
    sProcess.nMinFiles   = 1;
    % Default values for some options
    sProcess.isSourceAbsolute = 0;
    sProcess.processDim       = 1;    % Process channel by channel
    sProcess.isPaired         = 1;

    % === Process description
    sProcess.options.label1.Comment = ['This process normalizes each signal and frequency bin (FilesB)
' ...
                                       'with respect to a baseline (FilesA). In the formulas below:
'...
                                       '  x = data to normalize (FilesB)
' ...
                                       '  μ = mean over the baseline (FilesA) [mean(x(iBaseline))]
' ...
                                       '  σ = standard deviation over the baseline (FilesA) [std(x(iBaseline))]

' ...
                                       'Data selection:'];
    sProcess.options.label1.Type = 'label';
    % Common options
    sProcess = process_baseline_norm('DefineOptions', sProcess);
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = process_baseline_norm('FormatComment', sProcess);
end

%% ===== GET FILE TAG =====
function fileTag = GetFileTag(sProcess)
    fileTag = sProcess.options.method.Value;
end

%% ===== RUN =====
function sInputB = Run(sProcess, sInputA, sInputB) %#ok<DEFNU>
    % Call the corresponding process1
    sInputB = process_baseline_norm('Run', sProcess, sInputA, sInputB);
end


text:
filename='Bridge_triangle_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'PROJECTED GRADIENT'; incrementFactor = 1;
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.2;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [y,fs,wrd,phn,ffx]=v_readsph(filename,mode,nmax,nskip)
%V_READSPH  Read a SPHERE/TIMIT format sound file [Y,FS,WRD,PHN,FFX]=(FILENAME,MODE,NMAX,NSKIP)
%
% Input Parameters:
%
%	FILENAME gives the name of the file (with optional .SPH extension) or alternatively
%                 can be the FFX output from a previous call to READSPH having the 'f' mode option
%	MODE		specifies the following (*=default):
%
%    Scaling: 's'    Auto scale to make data peak = +-1 (use with caution if reading in chunks)
%             'r'    Raw unscaled data (integer values)
%             'p' *	Scaled to make +-1 equal full scale
%             'o'    Scale to bin centre rather than bin edge (e.g. 127 rather than 127.5 for 8 bit values)
%                     (can be combined with n+p,r,s modes)
%             'n'    Scale to negative peak rather than positive peak (e.g. 128.5 rather than 127.5 for 8 bit values)
%                     (can be combined with o+p,r,s modes)
%   Format    'l'    Little endian data (Intel,DEC) (overrides indication in file)
%             'b'    Big endian data (non Intel/DEC) (overrides indication in file)
%
%   File I/O: 'f'    Do not close file on exit
%             'd'    Look in data directory: v_voicebox('dir_data')
%             'w'    Also read the annotation file *.wrd if present (as in TIMIT)
%             't'    Also read the phonetic transcription file *.phn if present (as in TIMIT)
%                    Eac line of the annotation and transcription files is of the form: m n token
%                    where m and n are start end end times in samples and token is a word or phoneme test descriptor
%                    The corresponding cell arrays WRD and PHN contain two elements per row: {[m n]/fs 'token'}
%                    These outputs are only present if the corresponding 'w' and 't' options are selected
%
%	NMAX     maximum number of samples to read (or -1 for unlimited [default])
%	NSKIP    number of samples to skip from start of file
%               (or -1 to continue from previous read when FFX is given instead of FILENAME [default])
%ff
% Output Parameters:
%
%	Y          data matrix of dimension (samples,channels)
%	FS         sample frequency in Hz
%	WRD{*,2}   cell array with word annotations: WRD{*,:)={[t_start t_end],'text'} where times are in seconds
%              with the first sample at t=0 [only present if 'w' option is selected]
%	PHN{*,2}   cell array with phoneme annotations: PHN{*,:)={[t_start	t_end],'phoneme'} where times are in seconds
%              with the first sample at t=0 [only present if 't' option is selected]
%	FFX        Cell array containing
%
%     {1}     filename
%     {2}     header information
%        {1}  first header field name
%        {2}  first header field value
%     {3}     format string (e.g. NIST_1A)
%     {4}(1)  file id
%        (2)  current position in file
%        (3)  dataoff	byte offset in file to start of data
%        (4)  order  byte order (l or b)
%        (5)  nsamp	number of samples
%        (6)  number of channels
%        (7)  nbytes	bytes per data value
%        (8)  bits	number of bits of precision
%        (9)  fs	sample frequency
%		 (10) min value
%        (11) max value
%        (12) coding: 0=PCM,1=uLAW + 0=no compression,10=shorten,20=wavpack,30=shortpack
%        (13) file not yet decompressed
%     {5}     temporary filename
%
%   If no output parameters are specified, header information will be printed.
%   To decode shorten-encoded files, the program shorten.exe must be in the same directory as this m-file
%
%  Usage Examples:
%
% (a) Draw an annotated spectrogram of a TIMIT file
%           filename='....TIMIT/TEST/DR1/FAKS0/SA1.WAV';
%           [s,fs,wrd,phn]=v_readsph(filename,'wt');
%           v_spgrambw(s,fs,'Jwcpta',[],[],[],[],wrd);

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

persistent BYTEORDER
codes={'sample_count'; 'channel_count';  'sample_n_bytes';'sample_sig_bits'; 'sample_rate'; 'sample_min'; 'sample_max'};
codings={'pcm'; 'ulaw'};
compressions={',embedded-shorten-';',embedded-wavpack-'; ',embedded-shortpack-'};
if isempty(BYTEORDER), BYTEORDER='l'; end
if nargin<1, error('Usage: [y,fs,hdr,fidx]=READSPH(filename,mode,nmax,nskip)'); end
if nargin<2, mode='p';
else mode = [mode(:).' 'p'];
end
k=find((mode>='p') & (mode<='s'));
mno=all(mode~='o');                      % scale to input limits not output limits
sc=mode(k(1));
if any(mode=='l'), BYTEORDER='l';
elseif any(mode=='b'), BYTEORDER='b';
end
if nargout
    ffx=cell(5,1);
    if ischar(filename)
        if any(mode=='d')
            filename=fullfile(v_voicebox('dir_data'),filename);
        end
        fid=fopen(filename,'rb',BYTEORDER);
        if fid == -1
            fn=[filename,'.sph'];
            fid=fopen(fn,'rb',BYTEORDER);
            if fid ~= -1, filename=fn; end
        end
        if fid == -1
            error('Can''t open %s for input',filename);
        end
        ffx{1}=filename;
    else
        if iscell(filename)
            ffx=filename;
        else
            fid=filename;
        end
    end

    if isempty(ffx{4})
        fseek(fid,0,-1);
        str=char(fread(fid,16)');
        if str(8) ~= 10 || str(16) ~= 10, fclose(fid); error('File does not begin with a SPHERE header'); end
        ffx{3}=str(1:7);
        hlen=str2double(str(9:15));
        hdr={};
        while 1
            str=fgetl(fid);
            if str(1) ~= ';'
                [tok,str]=strtok(str);
                if strcmp(tok,'end_head'), break; end
                hdr(end+1,1)={tok};
                [tok,str]=strtok(str);
                if tok(1) ~= '-', error('Missing ''-'' in SPHERE header'); end
                if tok(2)=='s'
                    hdr(end,2)={str(2:str2num(tok(3:end))+1)};
                elseif tok(2)=='i'
                    hdr(end,2)={sscanf(str,'%d',1)};
                else
                    hdr(end,2)={sscanf(str,'%f',1)};
                end
            end
        end
        i=find(strcmp(hdr(:,1),'sample_byte_format'));
        if ~isempty(i)
            bord=char('b'+('l'-'b')*(hdr{i,2}(1)=='0'));
            if bord ~= BYTEORDER && all(mode~='b') && all(mode ~='l')
                BYTEORDER=bord;
                fclose(fid);
                fid=fopen(filename,'rb',BYTEORDER);
            end
        end
        i=find(strcmp(hdr(:,1),'sample_coding'));
        icode=0;                % initialize to PCM coding
        if ~isempty(i)
            icode=-1;                   % unknown code
            scode=hdr{i,2};
            nscode=length(scode);
            for j=1:length(codings)
                lenj=length(codings{j});
                if strcmp(scode(1:min(nscode,lenj)),codings{j})
                    if nscode>lenj
                        for k=1:length(compressions)
                            lenk=length(compressions{k});
                            if strcmp(scode(lenj+1:min(lenj+lenk,nscode)),compressions{k})
                                icode=10*k+j-1;
                                break;
                            end
                        end
                    else
                        icode=j-1;
                    end
                    break;
                end
            end
        end

        info=[fid; 0; hlen; double(BYTEORDER); 0; 1; 2; 16; 1 ; 1; -1; icode];
        for j=1:7
            i=find(strcmp(hdr(:,1),codes{j}));
            if ~isempty(i)
                info(j+4)=hdr{i,2};
            end
        end
        if ~info(5)
            fseek(fid,0,1);
            info(5)=floor((ftell(fid)-info(3))/(info(6)*info(7)));
        end
        ffx{2}=hdr;
        ffx{4}=info;
    end
    info=ffx{4};
    if nargin<4, nskip=info(2);
    elseif nskip<0, nskip=info(2);
    end

    ksamples=info(5)-nskip;
    if nargin>2
        if nmax>=0
            ksamples=min(nmax,ksamples);
        end
    end

    if ksamples>0
        fid=info(1);
        if icode>=10 && isempty(ffx{5})
            fclose(fid);
            dirt=v_voicebox('dir_temp');
            filetemp=fullfile(dirt,'shorten.wav');
            cmdtemp=fullfile(dirt,'shorten.bat');               % batch file needed to convert to short filenames
            % if ~exist(cmdtemp,'file')                   % write out the batch file if it doesn't exist
                cmdfid=fopen(cmdtemp,'wt');
                fprintf(cmdfid,'@"%s" -x -a %%1 "%%~s2" "%%~s3"\n',v_voicebox('shorten'));
                fclose(cmdfid);
            % end
            if exist(filetemp,'file')                          % need to explicitly delete old file since shorten makes read-only
                doscom=['del /f "' filetemp '"'];
                if dos(doscom) % run the program
                    error('Error running DOS command: %s',doscom);
                end
            end
            if floor(icode/10)==1               % shorten
                doscom=['"' cmdtemp '" ' num2str(info(3)) ' "' filename '" "' filetemp '"'];
                %                     fprintf(1,'Executing: %s\n',doscom);
                if dos(doscom) % run the program
                    error('Error running DOS command: %s',doscom);
                end
            else
                error('unknown compression format');
            end
            ffx{5}=filetemp;
            fid=fopen(filetemp,'r',BYTEORDER);
            if fid<0, error('Cannot open decompressed file %s',filetemp); end
            info(1)=fid;                            % update fid
        end
        info(2)=nskip+ksamples;
        pk=pow2(0.5,8*info(7))*(1+(mno/2-all(mode~='n'))/pow2(0.5,info(8)));  % use modes o and n to determine effective peak
        fseek(fid,info(3)+info(6)*info(7)*nskip,-1);
        nsamples=info(6)*ksamples;
        if info(7)<3
            if info(7)<2
                y=fread(fid,nsamples,'uchar');
                if mod(info(12),10)==1
                    y=v_pcmu2lin(y);
                    pk=2.005649;
                elseif mod(info(12),10)==2,
                    y=v_pcma2lin(y);
                    pk=2.005649;
                else
                    y=y-128;
                end
            else
                y=fread(fid,nsamples,'short');
            end
        else
            if info(7)<4
                y=fread(fid,3*nsamples,'uchar');
                y=reshape(y,3,nsamples);
                y=[1 256 65536]*y-pow2(fix(pow2(y(3,:),-7)),24);
            else
                y=fread(fid,nsamples,'long');
            end
        end
        if sc ~= 'r'
            if sc=='s'
                if info(10)>info(11)
                    info(10)=min(y);
                    info(11)=max(y);
                end
                sf=1/max(max(abs(info(10:11))),1);
            else sf=1/pk;
            end
            y=sf*y;
        end
        if info(6)>1, y = reshape(y,info(6),ksamples).'; end
    else
        y=[];
    end

    if mode~='f'
        fclose(fid);
        info(1)=-1;
        if ~isempty(ffx{5})
            doscom=['del /f ' ffx{5}];
            if dos(doscom) % run the program
                error('Error running DOS command: %s',doscom);
            end
            ffx{5}=[];
        end
    end
    ffx{4}=info;
    fs=info(9);
    wrd=ffx;        % copy ffx into the other arguments in case 'w' and/or 't' are not specified
    phn=ffx;
    if any(mode=='w')
        wrd=cell(0,0);
        fidw=fopen([filename(1:end-3) 'wrd'],'r');
        if fidw>0
            while 1
                tline = fgetl(fidw); % read an input line
                if ~ischar(tline)
                    break
                end
                [wtim, ntim, ee, nix] = sscanf(tline,'%d%d',2);
                if ntim==2
                    wrd{end+1,1}=wtim(:)'/fs;
                    wrd{end,2}=strtrim(tline(nix:end));
                end
            end
            fclose(fidw);
        end
    end
    if any(mode=='t')
        ph=cell(0,0);
        fidw=fopen([filename(1:end-3) 'phn'],'r');
        if fidw>0
            while 1
                tline = fgetl(fidw); % read an input line
                if ~ischar(tline)
                    break
                end
                [wtim, ntim, ee, nix] = sscanf(tline,'%d%d',2);
                if ntim==2
                    ph{end+1,1}=wtim(:)'/fs;
                    ph{end,2}=strtrim(tline(nix:end));
                end
            end
            fclose(fidw);
        end
        if any(mode=='w')
            phn=ph;             % copy into 4th argument
        else
            wrd=ph;             % copy into 3rd argument
        end
    end
else
    [y1,fs,ffx]=v_readsph(filename,mode,0);
    info=ffx{4};
    icode=info(12); % could convert this into text
    if ~isempty(ffx{1}), fprintf(1,'Filename: %s\n',ffx{1}); end
    fprintf(1,'Sphere file type: %s, coding %d\n',ffx{3}, icode);
    fprintf(1,'Duration = %ss: %d channel * %d samples @ %sHz\n',v_sprintsi(info(5)/info(9)),info(6),info(5),v_sprintsi(info(9)));
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
classdef VARtypeSpecificFAVARsettings < bear.settings.favar.FAVARsettings

    properties
        onestep (1,1) logical = 1; % Bayesian estimation of factors and the model in an one-step estimation (1=yes, 0=no (two-step))
        % thining of Gibbs draws
        thin = 1; % (=1 default, no thinning)
        % priors on factor equation
        % Loadings L~N(0,L0*eye)
        L0 = 1; %BBE set-up
        % Covariance Sigma~IG(a,b)
        a0 = 3; %BBE set-up
        b0 = 0.001; %BBE set-up
    end

end

author: european-central-bank

repo: BEAR-toolbox


text:
function [month, day, year] = gdate (jdate)

% convert Julian date to Gregorian (calendar) date

% input

%  jdate = julian day

% output

%  month = calendar month [1 - 12]
%  day   = calendar day [1 - 31]
%  year  = calendar year [yyyy]

%  note: day may include fractional part

% Orbital Mechanics with MATLAB

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

jd = jdate;

z = fix(jd + .5);
fday = jd + .5 - z;

if (fday < 0)
   fday = fday + 1;
   z = z - 1;
end

if (z < 2299161)
   a = z;
else
   alpha = floor((z - 1867216.25) / 36524.25);
   a = z + 1 + alpha - floor(alpha / 4);
end

b = a + 1524;
c = fix((b - 122.1) / 365.25);
d = fix(365.25 * c);
e = fix((b - d) / 30.6001);
day = b - d - fix(30.6001 * e) + fday;

if (e < 14)
   month = e - 1;
else
   month = e - 13;
end

if (month > 2)
   year = c - 4716;
else
   year = c - 4715;
end

author: goGPS-Project


text:
function varargout = x
%X   A chebfun of the identity on [-1,1].
%   X = CHEB.X returns a chebfun object for the function @(x)x on [-1,1].
%
%   CHEB.X is shorthand for the expression X = CHEBFUN(@(X) X).

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

x = chebfun(@(x) x);

if ( nargout > 0 )

    % For the syntax x = cheb.x:
    varargout{1} = x;

    if ( nargout > 1 )
        error('CHEB:X:TooManyOutputs',...
            'Too many output arguments. CHEB.X only returns "x".')
    end

else

   % Put 'x' into the workspace:
   assignin('base', 'x', x)

end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [X0,isYout] = VBA_spm_Xadjust(SPMfile,VOIfile,varthresh)
% gets the effects of no interest, wrt to which the data has been adjusted
% [X0,isYout] = spm_Xadjust(SPMfile,VOIfile,varthresh)
% IN:
%   - SPMfile: name of the SPM file
%   - VOIfile: name of the VOI file
%   - varthresh: threshold for PCA on X0 (fraction of explained variance).
%   Default is 0.95.
% OUT:
%   - X0: confounds matrix
%   - isYout: vector of indices of scans that were effectively removed from
%   the GLM using scan-nulling regressors included in the original
%   confounds matrix.

try, varthresh; catch, varthresh = 0.95; end

load(SPMfile)
load(VOIfile)

if isequal(xY.Ic,0)
    X0 = [];
    isYout = [];
    return
end
Fc = SPM.xCon(xY.Ic);
ind = sum(abs(Fc.c),2)==0;
X0 = SPM.xX.X(:,ind);

% remove scan-nulling regressors
n0 = size(X0,2);
remove = [];
isYout = [];
for i=1:n0
    X0i = X0(:,i) - X0(1,i);
    i1 = find(X0i~=0);
    tmp = zeros(size(X0,1),1);
    tmp(i1) = 1;
    if isequal(X0i,tmp)
        remove = [remove;i];
        isYout = [isYout,i1];
    end
end
X0 = X0(:,setdiff(1:size(X0,2),remove));
X0 = bsxfun(@minus,X0,mean(X0,1));
X0 = bsxfun(@rdivide,X0,std(X0,[],1));
if varthresh < 1
    [u,s,v] = svd(X0);
    s2 = diag(s).^2;
    ev = cumsum(s2./sum(s2));
    it = find(ev>=varthresh, 1 );
    X0 = [u(:,1:it),ones(size(X0,1),1)./sqrt(size(X0,1))];
end

return

% the following reproduces the code in spm_regions.m:

% sX = SPM.xX.xKXs;
% X0 = sX.X*(eye(spm_sp('size',sX,2)) - spm_sp('xpx-',sX)*sf_H(Fc,sX));
%
% function H = sf_H(Fc,sX)
% if sf_ver(Fc) > 1,
%    hsqr = sf_Hsqr(Fc,sX);
%    H    = hsqr' * hsqr;
% else
%    H = Fc.c * pinv(Fc.X1o' * Fc.X1o) * Fc.c';
% end
%
% function hsqr = sf_Hsqr(Fc,sX)
% if sf_ver(Fc) > 1,
%    hsqr = spm_sp('ox',spm_sp('set',Fc.X1o.ukX1o))' * spm_sp('cukx',sX);
% else
%    hsqr = spm_sp('ox',spm_sp('set',Fc.X1o))'*spm_sp('x',sX);
% end
%
% function v = sf_ver(Fc)
% if isstruct(Fc.X0), v = 2; else v = 1; end

author: MBB-team

repo: VBA-toolbox


text:
function [qx,qP,qD,xhist] = spm_pf(M,y,U)
% Particle Filtering for dynamic models
% FORMAT [qx,qP,qD,xhist] = spm_pf(M,y)
% M - model specification structure
% y - output or data (N x T)
% U - exogenous input
%
% M(1).x                            % initial states
% M(1).f  = inline(f,'x','v','P')   % state equation
% M(1).g  = inline(g,'x','v','P')   % observer equation
% M(1).pE                           % parameters
% M(1).V                            % observation noise precision
%
% M(2).v                            % initial process noise
% M(2).V                            % process noise precision
%
% qx - conditional expectation of states
% qP - {1 x T} conditional covariance of states
% qD - full sample
%__________________________________________________________________________
% See notes at the end of this script for details and a demo.  This routine
% is based on:
%
% var der Merwe R, Doucet A, de Freitas N and Wan E (2000). The
% unscented particle filter.  Technical Report CUED/F-INFENG/TR 380
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_pf.m 1143 2008-02-07 19:33:33Z spm $

% check model specification
%--------------------------------------------------------------------------
M  = spm_DEM_M_set(M);
dt = M(1).E.dt;
if length(M) ~=2
    errordlg('spm_pf requires a two-level model')
    return
end

% INITIALISATION:
%==========================================================================
T    = length(y);                          % number of time points
n    = M(2).l;                             % number of innovations
N    = 200;                                % number of particles.

% precision of measurement noise
%--------------------------------------------------------------------------
R    = M(1).V;
for i = 1:length(M(1).Q)
    R = R + M(1).Q{i}*exp(M(1).h(i));
end
P  = M(1).pE;                              % parameters
Q  = M(2).V.^-.5;                          % root covariance of innovations
v  = kron(ones(1,N),M(2).v);               % innovations
x  = kron(ones(1,N),M(1).x);               % hidden states
v  = v + 128*Q*randn(size(v));

% inputs
%--------------------------------------------------------------------------
if nargin < 3
    U = sparse(n,T);
end

for t = 1:T

    % PREDICTION STEP: with the (8x) transition prior as proposal
    %----------------------------------------------------------------------
    for i = 1:N
        v(:,i)     = 8*Q*randn(n,1) + U(:,t);
        f          = M(1).f(x(:,i),v(:,i),P);
        dfdx       = spm_diff(M(1).f,x(:,i),v(:,i),P,1);
        xPred(:,i) = x(:,i) + spm_dx(dfdx,f,dt);
    end

    % EVALUATE IMPORTANCE WEIGHTS: and normalise
    %----------------------------------------------------------------------
    for i = 1:N
        yPred  = M(1).g(xPred(:,i),v(:,i),P);
        ePred  = yPred - y(:,t);
        w(i)   = ePred'*R*ePred;
    end
    w   = w - min(w);
    w   = exp(-w/2);
    w   = w/sum(w);

    % SELECTION STEP: multinomial resampling.
    %----------------------------------------------------------------------
    x   = xPred(:,multinomial(1:N,w));

    % report and record moments
    %----------------------------------------------------------------------
    qx(:,t)  = mean(x,2);
    qP{t}    = cov(x');
    qX(:,t)  = x(:);
    fprintf('PF: time-step = %i : %i\n',t,T);
end

% sample density
%==========================================================================
if nargout > 3
    xhist = linspace(min(qX(:)),max(qX(:)),32);
    for i = 1:T
        q = hist(qX(:,i),xhist);
        qD(:,i) = q(:);
    end
end

return

function I = multinomial(inIndex,q);
%==========================================================================
% PURPOSE : Performs the resampling stage of the SIR
%           in order(number of samples) steps.
% INPUTS  : - inIndex = Input particle indices.
%           - q       = Normalised importance ratios.
% OUTPUTS : - I = Resampled indices.
% AUTHORS : Arnaud Doucet and Nando de Freitas

% MULTINOMIAL SAMPLING:
% generate S ordered random variables uniformly distributed in [0,1]
% high speed Niclas Bergman Procedure
%--------------------------------------------------------------------------
q        = q(:);
S        = length(q);  % S = Number of particles.
N_babies = zeros(1,S);
cumDist  = cumsum(q');

u = fliplr(cumprod(rand(1,S).^(1./(S:-1:1))));
j = 1;
for i = 1:S
    while (u(1,i) > cumDist(1,j))
        j = j + 1;
    end
    N_babies(1,j) = N_babies(1,j) + 1;
end;

% COPY RESAMPLED TRAJECTORIES:
%--------------------------------------------------------------------------
index = 1;
for i = 1:S
    if (N_babies(1,i)>0)
        for j=index:index+N_babies(1,i)-1
            I(j) = inIndex(i);
        end;
    end;
    index = index + N_babies(1,i);
end

return
%==========================================================================

% notes and demo:
%==========================================================================
% The code below generates a nonlinear, non-Gaussian problem (S) comprising
% a model S.M and data S.Y (c.f. van der Merwe et al 2000))
%
% The model is   f(x) = dxdt
%                     = 1 + sin(0.04*pi*t) - log(2)*x + n
%                y    = g(x)
%                     = (x.^2)/5  : if t < 30
%                       -2 + x/2  : otherwise
% i.e. the output nonlinearity becomes linear after 30 time steps.  In this
% implementation time is modelled as an auxiliary state variable.  n is
% the process noise, which is modelled as a log-normal variate.  e is
% Gaussian observation noise.

% model specification
%--------------------------------------------------------------------------
f       = '[1; (1 + sin(P(2)*pi*x(1)) - P(1)*x(2) + exp(v))]';
g       = '(x(1) > 30)*(-2 + x(2)/2) + ~(x(1) > 30)*(x(2).^2)/5';
M(1).x  = [1; 1];                  % initial states
M(1).f  = inline(f,'x','v','P');   % state equation
M(1).g  = inline(g,'x','v','P');   % observer equation
M(1).pE = [log(2) 0.04];           % parameters
M(1).V  = exp(4);                  % observation noise precision

M(2).v  = 0;                       % initial process log(noise)
M(2).V  = 2.4;                     % process log(noise) precision

% generate data (output)
%--------------------------------------------------------------------------
T       = 60;                      % number of time points
S       = spm_DEM_generate(M,T);

% Particle filtering
%--------------------------------------------------------------------------
pf_x    = spm_pf(M,S.Y);

% plot results
%--------------------------------------------------------------------------
x       = S.pU.x{1};
plot([1:T],x(2,:),[1:T],pf_x(2,:))
legend({'true','PF'})

author: spm


text:
classdef MatrixVectorizedInverter_2x2 < MatrixVectorizedInverter_Interface

    methods (Access = public)

        function B = computeInverse(obj,A)
            d(1,1,:) = A(2,2,:);
            d(1,2,:) = -A(1,2,:);
            d(2,1,:) = -A(2,1,:);
            d(2,2,:) = A(1,1,:);

            det = obj.computeDeterminant(A);

            B = zeros(size(A));
            for i = 1:2
                for j = 1:2
                    B(i,j,:) = squeeze(d(i,j,:))./det;
                end
            end
        end

        function detA = computeDeterminant(~,A)
            detA = squeeze(A(1,1,:).*A(2,2,:)-A(1,2,:).*A(2,1,:));
        end

    end

end


text:
function y = cvx_value( x )
y = apply( @cvx_value, x );

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [varargout] = acq_makeDataFolder(varargin)
%ACQ_MAKEDATAFOLDER - Create a folder for saving EEG data
%
%Synopsis:
% folder= acq_makeDataFolder;
% folder= acq_makeDataFolder(OPT);
%
%Arguments:
% 'tmp' : creates temporal directory
% OPT: struct or property/value list of optinal properties:
%
%Returns:
% folder: name of the folder in which EEG signals will be saved
%Side effect:
% Set global variables BTB.Tp.Code, BTB.Tp.Dir
%
% XX-XX Benjamin Blankertz
% 06-12 Javier Pascual. Updated. Created new function getSubjectcode

global BTB

%BTB.Tp.Code = [];

%% Get the date
today_vec= clock;
today_str= sprintf('%02d_%02d_%02d', today_vec(1)-2000, today_vec(2:3));

if length(varargin)==1 && isequal(varargin{1},'tmp'),
  BTB.Tp.Dir= [BTB.RawDir 'Temp_' today_str '\'];
  if ~exist(BTB.Tp.Dir, 'dir'),
    mkdir_rec(BTB.Tp.Dir);
  end
  return;
end

if isempty(BTB.Acq.StartLetter),
  BTB.Acq.StartLetter= 'a';
end

props = {'MultipleFolders'	0   'BOOL'
         'Interactive'      1   'BOOL'};

% if nargin==0,
%   varargout{1} = opt_catProps(props, acq_getSubjectCode);
%   return;
% end;

opt= opt_proplistToStruct(varargin{:});
[opt,isdefault] = opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

%% Check whether a directory exists that is to be used
dd= dir([BTB.RawDir 'VP???_' today_str '*']);

FolderList{1}= strcat('Temp_', today_str);
Comment{1}= 'for testing purpose';
NewCode= acq_getSubjectCode;
FolderList{2}= strcat(NewCode, '_', today_str);
Comment{2}= 'to be generated as new folder';
for k= 1:length(dd),
  FolderList{2+k}= dd(k).name;
  de= dir([BTB.RawDir dd(k).name '\*']);
  deeg= dir([BTB.RawDir dd(k).name '\*.eeg']);
  if length(de)<=2,
    Comment{2+k}= 'empty folder';
  elseif isempty(deeg),
    Comment{2+k}= 'folder without *.eeg files';
  else
    Comment{2+k}= 'folder with *.eeg files';
  end
end

if opt.Interactive,
  for k= 1:length(FolderList),
    fprintf('(%d) %s  (%s)\n', k, FolderList{k}, Comment{k});
  end
  choice= 0;
  while choice<1 || choice>length(FolderList),
    msg= sprintf(' -> Input your choice (1-%d): ', length(FolderList));
    choice= input(msg);
    if isempty(choice),
      choice= 1;
    end
  end
else
  choice= 2;
end
BTB.Tp.Dir= fullfile(BTB.RawDir, FolderList{choice});
is= find(FolderList{choice}=='_', 1, 'first');
BTB.Tp.Code= FolderList{choice}(1:is-1);

if ~exist(BTB.Tp.Dir, 'dir'),
  mkdir(BTB.Tp.Dir);
end

fprintf('EEG data will be saved in <%s>.\n', BTB.Tp.Dir);

author: bbci


text:
function storeLabelListGT(varargin)
% storeLabelListGT(varargin)
%
% Create the labelListGT entry for each ground-truth segmentation file of the dataset.
%
% Copyright by Holger Caesar, 2015

% Parse input
p = inputParser;
addParameter(p, 'projectName', 'WeaklySupervisedLearning');
addParameter(p, 'dataset', SiftFlowDatasetMC());
addParameter(p, 'proposalNameGT', 'GroundTruth');
addParameter(p, 'subset', 'all');
parse(p, varargin{:});

projectName = p.Results.projectName;
dataset = p.Results.dataset;
proposalNameGT = p.Results.proposalNameGT;
subset = p.Results.subset;

% Create paths
global glFeaturesFolder;
segmentFolder = fullfile(glFeaturesFolder, projectName, dataset.name, 'segmentations', proposalNameGT);

% Get list of images in dataset
if strcmp(subset, 'all'),
    [imageList, imageCount] = dataset.getImageList();
elseif strcmp(subset, 'train'),
    [imageList, imageCount] = dataset.getImageListTrn();
elseif strcmp(subset, 'test'),
    [imageList, imageCount] = dataset.getImageListTst();
else
    error('Error: Unknown subset: %s', subset);
end;

for imageIdx = 1 : imageCount,
    printProgress('Processing image', imageIdx, imageCount, 50);

    % Get segmentation
    imageName = imageList{imageIdx};
    segmentPath = [segmentFolder, filesep, imageName, '.mat'];
    segmentStruct = load(segmentPath);
    if isfield(segmentStruct, 'labelListGT'),
        fprintf('Skipping existing file %s...\n', segmentPath);
        continue;
    end;
    propBlobs = segmentStruct.propBlobs;
    blobCount = numel(propBlobs);

    % Get labelMap
    labelMap = dataset.getImLabelMap(imageName);
    labelListGT = nan(blobCount, 1);

    for blobIdx = 1 : blobCount,
        % Get labels for that blob
        blob = propBlobs(blobIdx);
        blobInds = blobToImageInds(blob, size(labelMap));
        blobLabel = unique(labelMap(blobInds));
        assert(numel(blobLabel) == 1 && blobLabel > 0);
        labelListGT(blobIdx) = blobLabel;
    end;

    % Sanity check
    assert(~any(isnan(labelListGT)));

    % Append to file
    save(segmentPath, 'labelListGT', '-append');
end;

author: nightrome

repo: matconvnet-calvin


text:
function rs = fft(s)

%tstoolbox/@signal/fft
%   Syntax:
%     * f = fft(s)
%
%   Output arguments:
%     * f - n by 2 array, the first column contains the magnitudes, the
%       second one the phases.
%
%   Fourier transform of scalar signal s.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(1,1);

if ndim(s) > 1
	error('Only for scalar signals');
end

x = fftshift(fft(data(s)));
rs = signal(core([abs(x)/dlens(s,1) (abs(x) > 10 * sqrt(eps)) .* angle(x)]), s);
a = getaxis(rs, 1);
rs = setaxis(rs, 1, achse(unit(a)^(-1), -samplerate(a)/2, samplerate(a)/dlens(s,1)));
rs = setaxis(rs, 2, achse);
rs = setplothint(rs, 'subplotgraph');
rs = addhistory(rs, ['Fourier transform']);
rs = addcommandlines(rs, 's = fft(s');

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [bndry_pairs, edgelet_ids, edgelet_sp] = ...
        generate_sp_neighborhood(sp_seg)
% finds all pairs of superpixels which share a boundary.
%
% @authors:     Ahmad Humayun
% @contact:     ahumayun@cc.gatech.edu
% @affiliation: Georgia Institute of Technology
% @date:        Fall 2013 - Summer 2014

    % create offset images to find superpixel boundary pixels
    south_seg = [sp_seg(2:end,:); sp_seg(end,:)];
    east_seg = [sp_seg(:,2:end), sp_seg(:,end)];

    % find the actual boundary pixels
    south_bndry = south_seg ~= sp_seg;
    east_bndry = east_seg ~= sp_seg;

    % find the src and the dst superpixel for each boundary pixel
    south_sp_src = sp_seg(south_bndry);
    south_sp_dst = south_seg(south_bndry);
    east_sp_src = sp_seg(east_bndry);
    east_sp_dst = east_seg(east_bndry);

    south_bndry_pairs = find(south_bndry);
    east_bndry_pairs = find(east_bndry);
    bndry_pairs = [south_bndry_pairs, south_bndry_pairs+1;
                   east_bndry_pairs, ...
                   east_bndry_pairs+size(sp_seg,1)];

    src_dst_pairs = [south_sp_src south_sp_dst; ...
                     east_sp_src east_sp_dst];

    [src_dst_pairs, sorted_ind] = sort(src_dst_pairs, 2);

    sorted_ind = sub2ind(size(src_dst_pairs), ...
        repmat((1:size(sorted_ind,1))', 1, 2), sorted_ind);

    bndry_pairs = bndry_pairs(sorted_ind);

    [edgelet_sp, ~, edgelet_ids] = unique(src_dst_pairs, 'rows');
end

author: Cloud-CV

repo: object-proposals


text:
function estimate_rr(up)
%ESTIMATE_RR estimates RR from respiratory signals using each possible set
% of RR estimation options, as specified in PC's literature review.
%	            estimate_rr(up)
%
%	Inputs:
%		data            data files should be stored in the specified format
%       up              universal parameters structure
%
%	Outputs:
%       	for each subject, n:
%       n_rrEsts.m      - a file of RR estimates
%

fprintf('\n--- Estimating RRs ');
%% Extract list of resp signals from the first pt:
subj = up.paramSet.subj_list(1);
loadpath = [up.paths.data_save_folder, num2str(subj), up.paths.filenames.respSigs];
if exist(loadpath, 'file')
    filecontents = whos('-file', loadpath);
    respSigs = extractfield(filecontents, 'name');
else
    warning(['No respiratory signals found for Subject ', num2str(subj) '.'])
end

%% Cycle through each patient
for subj = up.paramSet.subj_list
    loaded_this_subj_respSigs = 0;
    %% Make window timings if necessary
    identify_subj_wins(subj, up);
    %% Cycle through each resp signal
    for respSig_no = 1:length(respSigs)
        for option_no = 1 : length(up.al.options.estimate_rr)
            % Skip if this processing has been done previously
            if strcmp(up.al.options.estimate_rr{option_no}, 'GCE')
                save_name = [ respSigs{respSig_no}(1:3) '_' up.al.options.estimate_rr{option_no} ];
            else
                save_name = [ respSigs{respSig_no} '_' up.al.options.estimate_rr{option_no} ];
            end
            savepath = [up.paths.data_save_folder, num2str(subj), up.paths.filenames.rrEsts, '.mat'];
            exist_log = check_exists(savepath, save_name);
            if exist_log
                continue
            end

            % load data if it hasn't yet been loaded
            if ~loaded_this_subj_respSigs
                % Signals
                load([up.paths.data_save_folder, num2str(subj), up.paths.filenames.respSigs]);
                % Window timings
                load([up.paths.data_save_folder, num2str(subj), up.paths.filenames.win_timings, '.mat']);
                loaded_this_subj_respSigs = 1;
            end
            % Identify the relevant respSig data
            eval(['rel_data = ' respSigs{respSig_no} ';']);
            % add current subject and resp sig for any methods that do not use a resp sig.
            rel_data.subj = subj;
            rel_data.respSig = respSigs{respSig_no};
            %% Calculate RR from this resp sig using each option for estimating RR
            if (length(rel_data.t) == 1 && isnan(rel_data.t)) || sum(isnan(rel_data.v))==length(rel_data.v)
                temp_rr.t = mean([wins.t_start(:)' ; wins.t_end(:)']); temp_rr.t = temp_rr.t(:);
                temp_rr.v = nan(length(temp_rr.t),1);
            else
                temp_rr = feval(up.al.options.estimate_rr{option_no}, rel_data, wins, up);
            end
            % store this series of rrs:
            if strcmp(up.al.options.estimate_rr{option_no}, 'GCE')
                eval([respSigs{respSig_no}(1:3) '_' up.al.options.estimate_rr{option_no} ' = temp_rr;']);
            else
                eval([respSigs{respSig_no} '_' up.al.options.estimate_rr{option_no} ' = temp_rr;']);
            end

            clear temp_rr
            %% Save RRs to file
            save_or_append_data
        end
    end
    clear ekg* ppg* wins            % clear resp sigs
end

end

author: peterhcharlton

repo: RRest