text:
function mv = mv_exportSelectivity(mv, saveFlag, threshold, conditions);
%
% mv = mv_exportSelectivity(mv, [saveFlag], [threshold], conditions);
%
% Export, from a MultiVoxel UI to a mrVista view,
% a  parameter map of selectivity indices for each voxel,
% this map is saved as Selectivity-ROI.mat
% a corresponding colormap is save SelectivityCmap.mat
% which can be loaded from the Color Map UI
%
% Category indepdendent selecitivities range from 0-1
% The color indicates both the selectivity and category preference
% Thus, selectivity values are mapped to: category+selectivity
% where category are integers ranging from 0:Ncats-1
%
%  A second map is saved as a corAnal
%  note that this will overwrite the default corAnal performed on these data
%  This map is used for thresholding the selectivity map
%  co field contains voxel reliability (note that cothres has values of 0-1
%  thus, all negative reliabilities are going to be thresholded by default)
%  ph contains category independent selectivities
%  The co and ph fields can used to threshold the selectivity maps
%
% When loading a saved selectivity map need to load 3 files:
% File -> load Parmater Map-> Selectvity-ROIname
% File-> load CorAnal
% Color Map -> Parameter Map Mode -> Load colormap from file -> SelectivityCamp
%
% ras, 09/2005
% kgs 10/2005

if ieNotDefined('mv'), mv = get(gcf,'UserData'); end
if ~exist('threshold','var')
    threshold=0.1;
end
%%%%%check that a view exists
mrGlobals;
switch mv.roi.viewType
    case 'Inplane', view = getSelectedInplane;
    case 'Volume', view = getSelectedVolume;
    case 'Gray', view = getSelectedGray;
    case 'Flat', view = getSelectedFlat;
end
if isempty(view), error('No mrVista view opened.'); return; end

% get relevant params about the view
mapdims = viewGet(view,'dataSize');
nScans = viewGet(view,'numScans');
scan = mv.params.scans(1);

%%%%%get amplitudes for selected conditions
whichConds = find(tc_selectedConds(mv));
whichConds = whichConds(whichConds>1)-1; % remove null
nConds = length(whichConds);

%%%%%compute selectivity index
[scaledSel sel]= mv_selectivity(mv, whichConds, threshold);

% plug in the values to the map volume:
mapvol = zeros(mapdims);
ind = roiIndices(view,mv.coords);
mapvol(ind) = scaledSel;

%%%%%get voxel reliabilities, map to co map
if ~isfield(mv,'wta'), mv=mv_reliability(mv,'plotFlag',0); end
covol = zeros(mapdims);
covol(ind) = sel;
% % hack: move to range 0:1 instead of -1:1 so cothres can apply to the
% % whole range and not only positive numbers
% covol = 0.5*ones(mapdims);
% covol(ind) = mv.wta.voxR/2+covol(ind);

%%%%%map the reliabilityinto the ph map
phvol = zeros(mapdims);
% phvol(ind) = mod(sel, 1);
phvol(ind) = normalize(mv.wta.voxR, 0, 2*pi);

%%%%%create a color map
M = view.ui.mapMode;
colorsPerCond = floor(128/nConds);
M.numColors = colorsPerCond*nConds;
colors=[];
for i=1:nConds
    % make each condition have a gradient of colors up to
    % the full color specified by the color order
    col = mv.trials.condColors{whichConds(i)+1}; % full color
    for j = 1:colorsPerCond
        w = 0.3+ 0.5*(j-1)/colorsPerCond; % weight of color
        colors(end+1,:)=(w*col+ ((colorsPerCond-j+1)/(colorsPerCond))*[.7 .7 .7]);
    end
end
M.cmap = [gray(M.numGrays); colors];
max(M.cmap)
M.clipMode = [0.01 nConds+0.02]; % set to manual clip mode

% save the color map
cmapPath = fullfile(dataDir(view),'SelectivityCmap.mat');
cmap = M.cmap(M.numGrays+1:end,:);
save(cmapPath,'cmap');
fprintf('Saved color map in %s\n',cmapPath);

%%%%%set in view
map = cell(1,nScans); co = cell(1,nScans); ph = cell(1,nScans);
mapName = sprintf('Selectivity_%s',mv.roi.name);
map{scan} = mapvol; co{scan} = covol; ph{scan} = phvol;
view.map = map; view.mapName = mapName;
view.co = co; view.ph = ph; view.amp=map;
view.ui.mapMode = M;
try
    view=setParameterMap(view,map,mapName);
end
if saveFlag==1,
    saveParameterMap(view,[],1);
    saveCorAnal(view,1)
end
try
    refreshScreen(view);
catch
    return
end

% evaluate this in the workspace, so the view
% itself is updated
assignin('base','map',map);
assignin('base','co',co);
assignin('base','ph',ph);
assignin('base','tmp',M);

evalin('base',sprintf('%s=setParameterMap(%s,map,''%s'');',...
    view.name,view.name,mapName));
evalin('base',sprintf('%s=setDisplayMode(%s,''map'');',...
    view.name,view.name));
evalin('base',sprintf('%s.map = map;',view.name));
evalin('base',sprintf('%s.ui.mapMode=tmp;',view.name));
evalin('base',sprintf('%s.co=co;',view.name));
evalin('base',sprintf('%s.ph=ph;',view.name));
evalin('base',sprintf('%s=refreshScreen(%s);',view.name,view.name));

disp('Exported voxel selectivity to mrVista view.')

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [U, R, V] = cod(A, tol)
%COD    Complete orthogonal decomposition.
%       [U, R, V] = COD(A, TOL) computes a decomposition A = U*T*V,
%       where U and V are unitary, T = [R 0; 0 0] has the same dimensions as
%       A, and R is upper triangular and nonsingular of dimension rank(A).
%       Rank decisions are made using TOL, which defaults to approximately
%       LENGTH(A)*NORM(A)*EPS.
%       By itself, COD(A, TOL) returns R.

%       Reference:
%       G. H. Golub and C. F. Van Loan, Matrix Computations, third
%       edition, Johns Hopkins University Press, Baltimore, Maryland,
%       1996; sec. 5.4.2.

[m, n] = size(A);

% QR decomposition.
[U, R, P] = qr(A);    % AP = UR
V = P';               % A = URV;
if nargin == 1, tol = max(m,n)*eps*abs(R(1,1)); end  % |R(1,1)| approx NORM(A).

% Determine r = effective rank.
r = sum(abs(diag(R)) > tol);
r = r(1);             % Fix for case where R is vector.
R = R(1:r,:);         % Throw away negligible rows (incl. all zero rows, m>n).

if r ~= n

   % Reduce nxr R' =  r  [L]  to lower triangular form: QR' = [Lbar].
   %                 n-r [M]                                  [0]

   [Q, R] = trap2tri(R');
   V = Q*V;
   R = R';

end

if nargout <= 1, U = R; end

author: opencobra

repo: cobratoolbox


text:
clear all
close all
path(path,'..\..\FUZZCLUST')
colors={'r.' 'gx' 'b+' 'ys' 'md' 'cv' 'k.' 'r*' 'g*' 'b*' 'y*' 'm*' 'c*' 'k*' };

%the data
data.X=nDexample(5,250,2,1);
%normalization
data=clust_normalize(data,'range');

%parameters
param.c=3;
param.m=2;
param.e=1e-3;
param.ro=ones(1,param.c);
param.val=1;

%Gustafson Kessel-clustering
result = GKclust(data,param);

%validation
result = validity(result,data,param);
plot(data.X(:,1),data.X(:,2),'b.',result.cluster.v(:,1),result.cluster.v(:,2),'ro');
hold on
plot(result.cluster.v(:,1),result.cluster.v(:,2),'ro');

%evaluation
new.X=data.X;
eval = clusteval(new,result,param);
result.validity

author: Sable

repo: mcbench-benchmarks


text:
function [model,delta] = getPredictors(stimList, HRF,nsess,varargin)
% function [model,delta] = getPredictorsOddEven(stimList, HRF, nsess,[downsample factor])
%
% Build predictors and delta functions, given a condition function or delta
% function and either a convolution matrix or vector.
%
% IMPORTANT: The intercept in this function is added.
%
% stimList: condition function or delta function
% HRF:      hemodynamic response function, or convolution matrix (columns
%           are HRF)
% [downsamp]    takes every nth element of the design matrix
%
%
%inputs:
% 1     a col. vector of stimulus conditions OR a delta function matrix
% 2	    an HRF vector sampled at the frequency of the stimulus vector, OR
%	    a convolution matrix H
% 3     Number of sessions (must be a divisor of the stimList length)
%
%outputs:
% 1     a n x 2 matrix of regressors (cols) for each condition
% 2     a n x k delta matrix with onsets
%
% Tor Wager, last modified 2/22/04 to center predictors
%            modified to optionally take H convolution matrix as input
%            in place of HRF vector.  See convmtx.m (Buracas,mseq toolbox)
%
% stimList can be condition function e.g., [1 3 2 4 3 2 1]' or
% delta matrix (n x k), n samples and k conditions, e.g., [1 0 0 0 1 0 1]'
%
% Resampling: the default N in matlab resample has built-in antialiasing,
% but may not be good for fmri designs!  The appropriate downsampling
% is expected to be res*TR (res is units of samples/s), but we use 0
% because the model will depend on the analysis method used, and this is
% the most veridical approach.  With N = 0, every ith sample is used, where
% i is the downsampling factor you input.  Popular choices are 16*TR (for
% onsets2delta.m), using the SPM default res of 16.
% Delta is NOT resampled.
%
% example: TR = 2, 16 samples per second in hi-res delta dhr
% [tmp,d] = downsample_delta(dhr,16*2); X=getPredictors(d,hrf);

model = [];

% build odd/even session list
div = size(stimList,1) ./ nsess;
if div ~= round(div), error('The stimList length must be divisible by nsess!'),end
tmp = repmat([ones(div,1); 2*ones(div,1)],ceil(nsess./2),1);
odde = tmp(1:size(stimList,1));

if min(size(stimList)) > 1 % delta matrix

    delta = stimList;
    odde = repmat(odde,1,size(delta,2));
    delta = [delta delta];
    odde = [odde==1 odde==2];
    delta = delta .* odde;

    if min(size(HRF)) == 1
        for i = 1:size(delta,2)
            model(:,i) = conv(delta(:,i), HRF);
        end
    end

else

    for i = 1:max(stimList(:,1)) % condition function

        delta(:,i) = (stimList == i) & odde==1;
        d2(:,i) = (stimList == i) & odde==2;

        if min(size(HRF)) == 1
            model(:,i) = conv(delta(:,i), HRF);
            model2(:,i) = conv(d2(:,i), HRF);
        end

    end

    delta = [delta d2];
    model = [model model2];
end

if min(size(HRF)) > 1 % convolution matrix
    model = HRF * delta;
end

model = model(1:size(stimList,1),:);      	% eliminate extra values

% downsample, if necessary
if length(varargin) > 1
    model = model(1:varargin{1}:end,:); % equivalent to resample(X,1,varargin{1},0)
end

model = model - repmat(mean(model),size(model,1),1);    % center predictors

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef MimMatNatDatabaseSeries < handle
    % MimMatNatDatabaseSeries.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (SetAccess = private)
        Name
        StudyName
        Modality
        Date
        Time
        SeriesUid
        StudyUid
        NumberOfImages
    end

    methods
        function obj = MimMatNatDatabaseSeries(series_id, single_image_metainfo, number_of_images)
            if nargin > 0
                obj.Name = single_image_metainfo.SeriesDescription;
                obj.StudyName = single_image_metainfo.StudyDescription;
                obj.Modality = single_image_metainfo.Modality;
                obj.Date = single_image_metainfo.Date;
                obj.Time = single_image_metainfo.Time;
                obj.SeriesUid = series_id;
                obj.StudyUid = single_image_metainfo.StudyUid;
                obj.NumberOfImages = number_of_images;
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function savefast(filename, varargin)
% savefast: fast saves of large arrays to .mat files
%
% Matlab's 'save' command can be very slow when saving large arrays,
% because by default Matlab attempts to use compression. This function
% provides a much faster alternative, at the cost of larger files.
%
% The syntax is identical to that of the Matlab save command.
%
% Example:
% >> ops = struct('algorithm', 'greedy');
% >> A = int32(randi(20, 1000, 1200, 40));
% >> B = randn(500, 1800, 60);
% >> tic; save /tmp/test ops A B; toc
% Elapsed time is 22.980294 seconds.
% >> tic; savefast /tmp/test ops A B; toc
% Elapsed time is 0.571098 seconds.

% Copyright 2013 by Timothy E. Holy

  % Extract the variable values
  vars = cell(size(varargin));
  for i = 1:numel(vars)
    vars{i} = evalin('caller', varargin{i});
  end

  % Separate numeric arrays from the rest
  isnum = cellfun(@(x) isa(x, 'numeric'), vars);

  % Append .mat if necessary
  [filepath, filebase, ext] = fileparts(filename);
  if isempty(ext)
    filename = fullfile(filepath, [filebase '.mat']);
  end

  create_dummy = false;
  if all(isnum)
    % Save a dummy variable, just to create the file
    dummy = 0; %#ok<NASGU>
    save(filename, '-v7.3', 'dummy');
    create_dummy = true;
  else
    s = struct;
    for i = 1:numel(isnum)
      if ~isnum(i)
        s.(varargin{i}) = vars{i};
      end
    end
    save(filename, '-v7.3', '-struct', 's');
  end

  % Delete the dummy, if necessary, just in case the user supplied a
  % variable called dummy
  if create_dummy
    fid = H5F.open(filename,'H5F_ACC_RDWR','H5P_DEFAULT');
    H5L.delete(fid,'dummy','H5P_DEFAULT');
    H5F.close(fid);
  end

  % Save all numeric variables
  for i = 1:numel(isnum)
    if ~isnum(i)
      continue
    end
    varname = ['/' varargin{i}];
    h5create(filename, varname, size(vars{i}), 'DataType', class(vars{i}));
    h5write(filename, varname, vars{i});
  end
end

author: rksltnl


text:
function [success] = modelAddOutsideDependencies(dependencyName,varargin)
	% Used to request certain outside dependencies from users.
	% Biafra Ahanonu
	% started: 2017.11.16 [16:50:28]
	% inputs
		%
	% outputs
		%

	[success] = ciapkg.io.modelAddOutsideDependencies(dependencyName,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function center = imGeodesicExtremities(img, varargin)
%IMGEODESICEXTREMITIES Compute geodesic extremities of a binary particle
%
%   RES = imGeodesicExtremities(IMG);
%   IMG is a binary image representing a connected particle. The result RES
%   is a binary image with foreground pixels corresponding to the geodesic
%   extremities of the particle. The geodesic extremities are the pixels
%   belonging to the particle whose geodesic propagation equal the geodesic
%   length of the particle.
%
%   RES = imGeodesicExtremities(IMG, WEIGHTS);
%   use different weights for the computation of distances. See
%   imChamferDistance for further details.
%
%   Example
%   imGeodesicExtremities
%
%   Note:
%   As the algorithm propagates geodesic distances from each foreground
%   pixel, the computation time may be expensive.
%
%   See also
%   imGeodesics, imChamferDistance, imGeodesicPropagation
%   imGeodesicRadius, imGeodesicCenter, imGeodesicDiameter
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-05-22,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

propag = imGeodesicPropagation(img, varargin{:});
maxVal = max(propag(img));
center = propag==maxVal;


text:
% Test file for @deltafun/zeroDeltaFun.m.

function pass = test_zeroDeltaFun(pref)

if (nargin < 1)
    pref = chebfunpref();
end
%%
dTol = pref.deltaPrefs.deltaTol;

d = deltafun.zeroDeltaFun();
pass(1) = iszero(d.funPart);
pass(2) = all( d.funPart.domain == [-1, 1] );
pass(3) = isempty(d.deltaMag);
pass(4) = isempty(d.deltaLoc);

d = deltafun.zeroDeltaFun([4, 5]);
pass(5) = iszero(d.funPart);
pass(6) = all( d.funPart.domain == [4, 5] );
pass(7) = isempty(d.deltaMag);
pass(8) = isempty(d.deltaLoc);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% FEM2D_FSU uses the fsuClusterMatlab command to run the FEM2D computation.
%
%  Discussion:
%
%    The arguments to fsuClusterMatlab mean:
%
%    [] allows us to specify an output directory;
%    [] allows us to specify switches to the MOAB queue handler, such as more time.
%    'm" indicates this is an "MPI-like" job;
%    'w' indicates that the MATLAB session should wait for completion;
%    4 is the number of processors requested;
%    @fem2d_fun indicates the function we want to execute;
%    { param } contains input arguments.
%
%    The function "fem2d_fun" must correspond to a MATLAB M-file "fem2d_fun.m"
%    and that file must be in MATLAB's path.
%
%    The easiest way to ensure that the file is in MATLAB's path is to
%    create a subdirectory called "matlab" immediately below your home directory,
%    and place the function file there (as well as any other files needed by the job).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 March 2010
%
%  Author:
%
%    John Burkardt
%
  param = p_data ( );

  fsuClusterMatlab ( [], [], 'm', 'w', 4, @fem2d_fun, { param } )

author: johannesgerer

repo: jburkardt-m


text:
% Function to estimate the glottal parameters: NAQ, QOQ, H1-H2, HRF and PSP
%
% Octave compatible
%
% Description
%  This function can be used to estimate a range of conventional glottal
%  source parameters often used in the literature. This includes: the
%  normalised amplitude quotient (NAQ), the quasi-open quotient (QOQ), the
%  difference in amplitude of the first two harmonics of the differentiated
%  glottal source spectrum (H1-H2), the harmonic richness factor (HRF) and
%  the parabolic spectral parameter (PSP)
%
% Inputs
%  gf       : [samples] [Nx1] Glottal flow estimation
%  gfd      : [samples] [Nx1] Glottal flow derivative estimation
%  fs       : [Hz]      [1x1] sampling frequency
%  GCI      : [s]       [Mx1] Glottal closure instants
%
% Outputs
%  NAQ      : [s,samples] [Mx2] Normalised amplitude quotient
%  QOQ      : [s,samples] [Mx2] Quasi-open quotient
%  H1H2     : [s,dB]      [Mx2] Difference in glottal harmonic amplitude
%  HRF      : [s,samples] [Mx2] Harmonic richness factor
%  PSP      : [s,samples] [Mx2] Parabolic spectral parameter
%
% Example
%  Please see the HOWTO_glottalsource.m example file.
%
% References
%  [1] Alku, P., B ackstrom, T., and Vilkman, E. Normalized amplitude quotient
%     for parameterization of the glottal flow. Journal of the Acoustical
%     Society of America, 112(2):701–710, 2002.
%  [2] Hacki, T. Klassifizierung von glottisdysfunktionen mit hilfe der
%     elektroglottographie. Folia Phoniatrica, pages 43–48, 1989.
%  [3] Alku, P., Strik, H., and Vilkman, E. Parabolic spectral parameter -
%     A new method for quantification of the glottal flow. Speech
%     Communication, 22(1):67–79, 1997.
%  [4] Hanson, H. M. Glottal characteristics of female speakers: Acoustic
%     correlates. Journal of the Acoustical Society of America,
%     10(1):466–481, 1997.
%  [5] Childers, D. G. and Lee, C. K. Voice quality factors: Analysis,
%     synthesis and perception. Journal of the Acoustical Society of
%     America, 90(5):2394–2410, 1991.
%
% Copyright (c) 2013 Trinity College Dublin
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  John Kane kanejo@tcd.ie

function [NAQ,QOQ,H1H2,HRF,PSP] = get_vq_params(gf,gfd,fs,GCI)

GCI = round(GCI*fs)+1;

%% Initial settings
F0min=20;
F0max=500;
NAQ=zeros(1,length(GCI));
QOQ=zeros(1,length(GCI));
H1H2=zeros(1,length(GCI));
HRF=zeros(1,length(GCI));
PSP=zeros(1,length(GCI));

glot_shift=round(0.5/1000*fs);
qoq_level=0.5; % Threshold for QOQ estimation
T0_num=3; % Number of local glottal pulses to be used for harmonic spectrum
min_harm_num=5;
HRF_freq_max=5000; % Maximum frequency used for harmonic measurement
PSP_fft_size=2048; % As per Alku et al (1997)

%% Do processing
for n=1:length(GCI)

    % Get glottal pulse compensated for zero-line drift
    if n==1
        start=1;
        stop=GCI(n);
        T0=GCI(n+1)-GCI(n);
    else start=GCI(n-1);
        stop=GCI(n);
        T0=GCI(n)-GCI(n-1);
    end
    F0=fs/T0;

    if isinf(F0)==0 && T0~=0 && F0 > F0min && F0<F0max

        % Compensate for zero-line drift in the glottal flow pulse
        gf_comb=[gf(start) gf(stop)];
        if start~=stop
            if length(gf_comb)>1
                line=interp1(linspace(1,stop-start+1,2),gf_comb, ...
                    1:stop-start+1);
            else line=gf_comb;
            end
        else line=0;
        end
        gf_seg=gf(start:stop);
        gf_seg_comp=gf_seg(:)-line(:);

        if stop+glot_shift <= length(gfd)
            stop2=stop+glot_shift;
        else
            stop2=stop;
        end
        gfd_seg=gfd(start:stop2);

        % Get NAQ and QOQ
        d_peak = max(abs(gfd_seg));
        [f_ac,max_idx] = max(gf_seg_comp);
        Amid=f_ac*qoq_level;
        [T1,T2] = findAmid_t(gf_seg_comp,Amid,max_idx);

        NAQ(n) = (f_ac/d_peak)/T0;
        QOQ(n)=(T2-T1)/(fs/F0);

         % Get frame positions for H1-H2 parameter
        if GCI(n)-round((T0*T0_num)/2) > 0
            f_start = GCI(n)-round((T0*T0_num)/2);
        else f_start=1;
        end
        if GCI(n)+round((T0*T0_num)/2) <= length(gfd)
            f_stop = GCI(n)+round((T0*T0_num)/2);
        else f_stop=length(gfd);
        end
        f_frame=gfd(f_start:f_stop);
        f_frame=f_frame.*(2^15);
        f_win=f_frame(:).*hamming(length(f_frame));
        f_spec = mag2db(abs(fft(f_win,fs)));

         % Get H1-H2/HRF measurements
        [h_idx,h_amp]=v_findpeaks(f_spec,[],F0/2);
        HRF_harm_num=floor(HRF_freq_max/F0);
        if length(h_idx) >= min_harm_num
            [~,f0_idx] = min(abs(bsxfun(@minus,h_idx,(1:HRF_harm_num)*F0)),[],1);
            H1H2(n)=h_amp(f0_idx(1))-h_amp(f0_idx(2));
            HRF(n) = sum(h_amp(f0_idx(2:end)))/h_amp(f0_idx(1));
        end

        % Get PSP
        if n>1
           start=GCI(n-1);
           stop=GCI(n);
           gf_seg=gf(start:stop);
           gf_seg=gf_seg-min(gf_seg); % Adjust minimum to zero
           gf_seg=gf_seg/max(gf_seg); % Scale to unity
           gf_seg(PSP_fft_size)=0; % Zero-pad
           X = mag2db(abs(fft(gf_seg))); % Spectrum
           X = X(linspace(1,fs,length(X))<=fs/2); % Use up to the nyquist
           a = psp_get_a(X);

           % Estimate measure of maximal spectral decay
           a_max_sig=1/round(fs/F0)*ones(1,round(fs/F0));
           a_max_sig(PSP_fft_size)=0; % Zero-pad
           X_a_max = mag2db(abs(fft(a_max_sig))); % Spectrum
           X_a_max = X_a_max(linspace(1,fs,length(X_a_max))<=fs/2); % Use up to the nyquist
           a_max = psp_get_a(X_a_max);

           PSP(n)=a/a_max;
        end

    end

end

%% Add in time to parameters
NAQ=[GCI(:)/fs NAQ(:)];
QOQ=[GCI(:)/fs QOQ(:)];
H1H2=[GCI(:)/fs H1H2(:)];
HRF=[GCI(:)/fs HRF(:)];
PSP=[GCI(:)/fs PSP(:)];

function [T1,T2] = findAmid_t(glot_adj,Amid,Tz)

% Function to find the start and stop positions of the quasi-open phase.

T1=0;
T2=0;
if Tz~=0
    n=Tz;

    while glot_adj(n) > Amid && n > 3
        n=n-1;
    end
    T1=n;
    n=Tz;

    while glot_adj(n) > Amid && n < length(glot_adj)-2
        n=n+1;
    end
    T2=n;
end

function a = psp_get_a(X)

% Function to calculate a coefficient, as per section 3.2.1 in Alku et al
% (1997)

%% Initial settings
NE_min = 0.01;
flag = 0;
N = 3; % An set to an initial value of 3, as per Alku et al (1997)
X=X(:);

% pre-allocate&compute
X2=X.^2;
k2 = ((0:100).^2)';

% init sum
k2_sum=sum(k2(1:N-1));
k4_sum=sum(k2(1:N-1).^2);
X_k1_sum = sum(X(1:N-1));
X_k1_k2_sum = sum(X((1:N-1)).*k2(1:N-1));
X2_k1_sum = sum(X2((1:N-1)));

%% Do processing
while flag == 0

    % re-allocate if necessary
    if N>numel(k2)
        k2 = ((0:numel(k2)*2).^2)';
    end

    % iteratively built sum
    k2_sum= k2_sum + k2(N);
    k4_sum= k4_sum + k2(N)*k2(N);
    X_k1_sum = X_k1_sum + X(N);
    X_k1_k2_sum = X_k1_k2_sum + X(N)*k2(N);
    X2_k1_sum = X2_k1_sum + X2(N);

    % Equation 5
    a = (N*X_k1_k2_sum-X_k1_sum.*k2_sum)/(N*k4_sum-k2_sum.^2);

    X_k1_a_k2 = X(1:N)-a*k2(1:N);

    % Equation 4
    b = 1/N*sum(X_k1_a_k2);
    NE = (sum((X_k1_a_k2-b).^2))/X2_k1_sum;

    % Increment or stop
    if NE < NE_min
        N = N+1;
    else
        flag = 1;
    end
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function cS = anorthite
% (Plagioclase feldspar)

cs_An = crystalSymmetry('-1',[8.1797 12.8748 14.1721], [93.13,115.89,91.24]*degree);
N = Miller({0,1,0},{0,0,1},{1,1,0},{1,-1,0},{1,1,-1},{1,-1,-1},{2,0,-1},cs_An);
dist = [0.65, 0.5, 1.25, 1.2, 1.2, 1.05, 0.95];
cS = crystalShape(N./dist);

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
% Exercise 4.47: Maximum determinant PSD matrix completion
% Boyd & Vandenberghe "Convex Optimization"
% Almir Mutapcic - Jan 2006
%
% Given a symmetric matrix A in R^(n-by-n) with some entries unspecified
% we find its completion such that A is positive semidefinite and
% it has a maximum determinant out of all possible completions.
% This problem can be formulated as a log det (and det_rootn) problem.
%
% This is a numerical instance of the specified book exercise.

% problem size
n = 4;

% create and solve the problem
cvx_begin
  % A is a PSD symmetric matrix (n-by-n)
  variable A(n,n) semidefinite;

  % constrained matrix entries.
  A(1,1) == 3; %#ok
  A(2,2) == 2; %#ok
  A(3,3) == 1; %#ok
  A(4,4) == 5; %#ok
  % Note that because A is symmetric, these off-diagonal
  % constraints affect the corresponding element on the
  % opposite side of the diagonal.
  A(1,2) == .5; %#ok
  A(1,4) == .25; %#ok
  A(2,3) == .75; %#ok

  % find the solution to the problem
  maximize( log_det( A ) )
  % maximize( det_rootn( A ) )
cvx_end

% display solution
fprintf('Matrix A with maximum determinant (%g) is:\n', det(A));
disp(A)
disp('Its eigenvalues are:')
disp(eig(A))

author: yu-jiang

repo: radpbook

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


text:
function output = callbaron(model)

% This sets up everything and more. Can be simplified significantly since
% baron handles its own computational tree etc
model = yalmip2nonlinearsolver(model);

% [Anonlinear*f(x) <= b;Anonlinear*f(x) == b]
cu = full([model.bnonlinineq;model.bnonlineq]);
cl = full([repmat(-inf,length(model.bnonlinineq),1);model.bnonlineq]);
Anonlinear = [ model.Anonlinineq; model.Anonlineq];

% Create string representing objective
obj = createmodelstring(model.c,model);
if nnz(model.Q)>0
    obj = [obj '+' createQstring(model.Q,model)];
end
if any(model.f) > 0
    obj = [obj '+' num2str(model.f)];
end
obj = baronify(obj);
if length(obj)>0
    obj = strrep(obj,'+-','-');
    obj = ['@(x) ' obj];
    obj = eval(obj);
else
    obj = @(x)0;
end

% Create string representing nonlinear consdbqtraints
if length(cu)>0
    con = '[';
    remove = [];
    for i = 1:length(cu)
        if isinf(cl(i)) & isinf(cu(i))
            remove = [remove;i];
        else
            con = [con createmodelstring(Anonlinear(i,:),model) ';'];
        end
    end
    cl(remove) = [];
    cu(remove) = [];
    con = [con ']'];
    con = baronify(con);
    con = ['@(x) ' con];
    con = eval(con);
else
    cu = [];
    cl = [];
    con = [];
end

% Linear constraints
ru = full([model.b;model.beq]);
rl = full([repmat(-inf,length(model.b),1);model.beq]);
A =  [model.A; model.Aeq];
if length(ru)>0
    remove = find(isinf(rl) & isinf(ru));
    A(remove,:)=[];
    rl(remove) = [];
    ru(remove) = [];
end

lb = model.lb;
ub = model.ub;

xtype = [];
xtype = repmat('C',length(lb),1);
xtype(model.binary_variables) = 'B';
xtype(model.integer_variables) = 'I';
x0 = model.x0;
opts = model.options.baron;
opts.prlevel = model.options.verbose;
if model.options.savedebug
    save barondebug obj con A ru rl cl cu lb ub x0 opts
end

solvertime = tic;
[x,fval,exitflag,info,allsol] = baron(obj,A,rl,ru,lb,ub,con,cl,cu,xtype,x0,opts);
solvertime = toc(solvertime);

% Check, currently not exhaustive...
switch exitflag
    case 1
        problem = 0;
    case 2
        problem = 1;
    case 3
        problem = 2;
    case 4
        problem = 3;
    case 5
        problem = 11;
    otherwise
        problem = 9;
end

% Save all data sent to solver?
if model.options.savesolverinput
    solverinput.obj = obj;
    solverinput.A = A;
    solverinput.rl = rl;
    solverinput.ru = ru;
    solverinput.lb = lb;
    solverinput.ub = ub;
    solverinput.con = con;
    solverinput.cl = cl;
    solverinput.cu = cu;
    solverinput.xtype = xtype;
    solverinput.opts = opts;
else
    solverinput = [];
end

% Save all data from the solver?
if model.options.savesolveroutput
    solveroutput.x = x;
    solveroutput.fval = fval;
    solveroutput.exitflag = exitflag;
    solveroutput.info=info;
    solveroutput.allsol=allsol;
else
    solveroutput = [];
end

if ~isempty(x)
    x = RecoverNonlinearSolverSolution(model,x);
end

% Standard interface
output = createOutputStructure(x,[],[],problem,'BARON',solverinput,solveroutput,solvertime);

function z = createOneterm(i,model)

[evalPos,pos] = ismember(i,model.evalVariables);
if pos
    % This is a nonlinear operator
    map = model.evalMap{pos};
    % We might hqave vectorized things to compute several expressions at the same time
    if strcmp(model.evalMap{pos}.fcn,'sum_square')
        z = ['('];
        for j = map.variableIndex
            jl = find(model.linearindicies == j);
            if isempty(jl)
                z =  [z '(' createmonomstring(model.monomtable(j,:),model) ')^2+'];
            else
                z =  [z 'x(' num2str(jl) ')^2+'];
            end
        end
        z = [z(1:end-1) ')'];
    elseif strcmp(model.evalMap{pos}.fcn,'slogfrac')
        j1 = find(model.linearindicies == map.variableIndex(1));
        j2 = find(model.linearindicies == map.variableIndex(2));
        z =  ['(log(1 + x(' num2str(j1) ')/x(' num2str(j2) ')))'];
    elseif strcmp(model.evalMap{pos}.fcn,'entropy')
        z = ['('];
        for j = map.variableIndex
            jl = find(model.linearindicies == j);
            if isempty(jl)
                z =  [z '(' createmonomstring(model.monomtable(j,:),model) ')^2+'];
            else
                z =  [z '-x(' num2str(jl) ')*log(x(' num2str(jl) '))+'];
            end
        end
        z = [z(1:end-1) ')'];
    elseif strcmp(model.evalMap{pos}.fcn,'logsumexp')
        z = ['('];
        for j = map.variableIndex
            jl = find(model.linearindicies == j);
            if isempty(jl)
                z =  [z '(' createmonomstring(model.monomtable(j,:),model) ')^2+'];
            else
                z =  [z 'exp(x(' num2str(jl) '))+'];
            end
        end
        z = ['log' z(1:end-1) ')'];
    elseif strcmp(model.evalMap{pos}.fcn,'power_internal1')
        j = find(map.computes == i);
        j = map.variableIndex(j);
        % we have f(x(j)), but have to map back to linear indicies
        jl = find(model.linearindicies == j);
        if isempty(jl)
            z =  [num2str(model.evalMap{pos}.arg{2}) '^(' createmonomstring(model.monomtable(j,:),model)  ')'];
        else
            z =  [num2str(model.evalMap{pos}.arg{2}) '^(x(' num2str(jl) '))'];
        end
    else
        j = find(map.computes == i);
        j = map.variableIndex(j);
        % we have f(x(j)), but have to map back to linear indicies
        jl = find(model.linearindicies == j);
        if isempty(jl)
            z =  [map.fcn '(' createmonomstring(model.monomtable(j,:),model)  ')'];
        else
            z =  [map.fcn '(x(' num2str(jl) '))'];
        end
    end
else
    i = find(model.linearindicies == i);
    z =  ['x(' num2str(i) ')'];
end

function monomstring = createmonomstring(v,model)

monomstring = '';
for i = find(v)
    z = createOneterm(i,model);
    if v(i)==1
        monomstring = [monomstring z '*'];
    else
        monomstring = [monomstring z '^' num2str(v(i),12) '*'];
    end
end
monomstring = monomstring(1:end-1);

function string = createmodelstring(row,model)
string = '';
index = find(row);
for i = index(:)'
    monomstring = createmonomstring(model.monomtable(i,:),model);
    string = [string  num2str(row(i),12) '*' monomstring '+'];
end
string = string(1:end-1);
string = strrep(string,'+-','-');
string = strrep(string,')-1*',')-');
string = strrep(string,')+1*',')+');

function string = createQstring(Q,model)
% Special code for the quadratic term
string = '';
[ii,jj,c]= find(triu(Q));
for i = 1:length(ii)
    if ii(i)==jj(i)
        % Quadratic term
        monomstring = createmonomstring(sparse(1,ii(i),1,1,size(Q,1)),model);
        string = [string  num2str(c(i),12) '*' monomstring '^2' '+'];
    else
        % Bilinear term
        monomstring1 = createmonomstring(sparse(1,ii(i),1,1,size(Q,1)),model);
        monomstring2 = createmonomstring(sparse(1,jj(i),1,1,size(Q,1)),model);
        string = [string  num2str(c(i),12) '*2*' monomstring1 '*' monomstring2 '+'];
    end
end
string = string(1:end-1);
string = strrep(string,'+-','-');
string = strrep(string,'-1*','-');
string = strrep(string,'+1*','+');

function s = baronify(s)
s = strrep(s,'+1*','+');
s = strrep(s,'-1*','-');
s = strrep(s,'sqrtm_internal','sqrt');
s = strrep(s,'cabs','abs');
% YALMIP catches log(1+z) and implements internal model slog(z)
% Replace with log
expression = 'slog((\w*)';
replace = 'log(1+$1';
s =  regexprep(s,expression,replace);
% YALMIP catches x^y and replaced with power1_internal
expression = 'slog((\w*)';
replace = 'log(1+$1';
s =  regexprep(s,expression,replace);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function value = r8_cotd ( degrees )

%*****************************************************************************80
%
%% R8_COTD returns the cotangent of an angle given in degrees.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real DEGREES, the angle in degrees.
%
%    Output, real VALUE, the cotangent of the angle.
%
  radians = pi * ( degrees / 180.0 );

  value = cos ( radians ) / sin ( radians );

  return
end


text:
%kVRML 'visualizes threshhold images using vrwave '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros VRML.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'First Input data object'
%
% Example: kVRML(i, {'i',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% VRML - visualizes threshhold images using vrwave
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1996-2003, Rainer Heintzmann,  All rights reserved.
%

function varargout = kVRML(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kVRML(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input'};
maxval={0};
minval={0};
istoggle=[0];
was_set=istoggle * 0;
paramtype={'InputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(0);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
callKhoros([w 'vrmlstart"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function edges_out = def_edges(X, nb, edges_in)

edges_out = edges_in;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function niftiWrite(ni,fName)
%  Writes a nifti structure to disk.
%
%   niftiWrite(ni,[fName])
%
% Wrapper to writeFileNifti mex file, which writes a VISTASOFT nifti
% structure to a NIFTI-1 file.
%
% INPUTS
%   ni - a nifti structure
%   fName - output file name  (defaults to ni.fname if empty).
%
% FILE FORMATS:
%   fName extension should be either '.nii' or 'nii.gz'
%
% EXAMPLE:
%   dFolder = mrtInstallSampleData('functional','mrBOLD_01');
%   pth =  fullfile(dFolder, 'Raw', 'T1_Inplane.nii.gz');
%   ni  = niftiRead(pth);
%   fName = fullfile(tempdir, 'myFile.nii.gz');
%   niftiWrite(ni,fName);
%
% See also:  niftiRead
%
% (C) Stanford VISTA, 2012-2015
%

if exist('fName', 'var') && ~isempty(fName)
    ni.fname = fName;
elseif ~isfield(ni,'fname') || isempty(ni.fname)
    error('A filename is required to save a NIFTI-1 structure to file.');
end

% EDIT GLU 2020-01-28
% We want to be able to write Nifti-2 files. If the the field version is not in
% the ni struct, or it is = 1, write a Nifti-1 file using the old mex code, to
% maximize backwards compatibility. If the version field is 2, write a Nifti-2
% using the nii_tool (https://github.com/xiangruili/dicm2nii)
if isfield(ni, 'version')
    if ni.version == 2
        % It will take file name from struct, and it will use the NIfti-2
        % from the struct as well.
        %
        % Modify it to format as nii_tool wants it
        ni2                = struct();
        ni2.img            = ni.data;
        ni                 = rmfield(ni, 'data');
        ni2.hdr            = ni;
        ni2.hdr.file_name  = ni2.hdr.fname;
        ni2.hdr            = rmfield(ni2.hdr, 'fname');
        % Use a function to setup the xyzt_units field
        ni2.hdr.xyzt_units = setSpaceTimeUnits(ni2.hdr.xyz_units, ni2.hdr.time_units);
        ni2.hdr            = rmfield(ni2.hdr, 'xyz_units');
        ni2.hdr            = rmfield(ni2.hdr, 'time_units');

        % Adapt back the rest of the fields using writeFileNifti code:
        ni2.hdr.pixdim              = zeros(1,1+numel(ni.pixdim));
        ni2.hdr.pixdim(2:1+numel(ni.pixdim)) = ni.pixdim;
        ni2.hdr.datatype            = ni.nifti_type;
        ni2.hdr                     = rmfield(ni2.hdr, 'nifti_type');
        ni2.hdr.pixdim(1)           = ni.qfac;
        ni2.hdr                     = rmfield(ni2.hdr, 'qfac');

        % When reading files we need to add 1 voxel to the header information to
        % complaint with matlab 1-based indexing.
        %
        % Here we undo such operation, we subtract one before writing the file to disk.
        if (~all(ni.qto_xyz(1:9)==0))
            ni.qto_ijk(1:3,4) = ni.qto_ijk(1:3,4) - 1;
            q_xyz = inv(ni.qto_ijk);
            q = matToQuat(q_xyz);
            ni2.hdr.qoffset_x = q.quatern_x;
            ni2.hdr.qoffset_y = q.quatern_y;
            ni2.hdr.qoffset_z = q.quatern_z;
        else
            ni2.hdr.qoffset_x = 0;
            ni2.hdr.qoffset_y = 0;
            ni2.hdr.qoffset_z = 0;
        end

        if (~all(ni.sto_xyz(1:9)==0))
            sto_ijk = ni.sto_ijk;
            sto_ijk(1:3,4) = ni.sto_ijk(1:3,4) - 1;
            sto_xyz = inv(sto_ijk);
            ni2.hdr.srow_x = sto_xyz(1,:);
            ni2.hdr.srow_y = sto_xyz(2,:);
            ni2.hdr.srow_z = sto_xyz(3,:);
        else
            ni2.hdr.srow_x = 0;
            ni2.hdr.srow_y = 0;
            ni2.hdr.srow_z = 0;
        end

        % --- nk fields
        % The fields freq_dim, phase_dim, slice_dim are all squished into the single
        % byte field dim_info (2 bits each, since the values for each field are
        % limited to the range 0..3):
        % ni.freq_dim =  bitand(uint8(3),          uint8(nii.hdr.hk.dim_info)    );
        % ni.phase_dim = bitand(uint8(3), bitshift(uint8(nii.hdr.hk.dim_info),-2));
        % ni.slice_dim = bitand(uint8(3), bitshift(uint8(nii.hdr.hk.dim_info),-4));
        % Here we undo the previous operations perfrmed during the file load.

        ni2.hdr.dim_info = ni2.hdr.freq_dim + 4*ni2.hdr.phase_dim + 16*ni2.hdr.slice_dim;

        % Save it.
        nii_tool('save', ni2);
    else
        writeFileNifti(ni);
    end
else
    writeFileNifti(ni);
end

return
end

% GLU 2020-01-28 added this function from Matlab's nifti implementation to
% convert back the spaceTime units to the nifti field required by
    function xyztCode = setSpaceTimeUnits(spaceUnitText, timeUnitText)
        %setSpaceTimeUnits: convert simplified space time units to standard
        %form.
        %    This is a helper function to convert simplified space and time
        %    units to the raw format as specified in the NIfTI header.

        spaceKey   = {0, 1, 2, 3};
        spaceValue = {'unknown', 'meter', 'mm', 'micron'};

        spaceMap = containers.Map(spaceValue, spaceKey);

        if isempty(find(strcmp(spaceValue, spaceUnitText), 1))
            error(message('images:nifti:spaceUnitNotSupported'));
        end

        spaceUnits = spaceMap(spaceUnitText);

        timeValue = {'unknown', 'sec', 'msec', 'usec', 'hz', 'ppm', 'rads'};
        timeKey = {0, 8, 16, 24, 32, 40, 48};

        timeMap = containers.Map(timeValue, timeKey);

        if isempty(find(strcmp(timeValue, timeUnitText), 1))
            error(message('images:nifti:timeUnitNotSupported'));
        end

        timeUnits = timeMap(timeUnitText);

        spaceUnitCode = bitand(uint8(spaceUnits),uint8(7));
        timeUnitCode  = bitand(uint8(timeUnits),uint8(56)); % 0x38

        xyztCode = bitor(spaceUnitCode, timeUnitCode);

    end

%% Old code
% The old code was not a wrapper to writeFileNifti, but instead converted
% the VISTASOFT nifti struct to a nifti-1 struct using the function
% niftiVista2ni, and then write the nifti-1 struct to file using the Shen
% function save_nii.
%
% We prefer the mex file writeFileNifti, which is fast which does the
% inverse conversion of the mex function readFileNifti (wrapped in the m
% file niftiRead).

% %% Set defaults
% compress = true;
% fname_supplied = false;
%
%
% %% Check inputs
% if nargin==0
%     help(mfilename)
%     return
% end
%
%
%
% if notDefined('fName')
%     fName = '';
% else
%     fname_supplied = true;
% end
%
%
% %% Check file extension
%
% % If no name was supplied check the nifti structure for the name
% if isempty(fName)
%     if isfield(ni,'fname') && ~isempty(ni.fname)
%         fName = ni.fname;
%     else
%         error('A filename is required to save a NIFTI-1 structure to file.');
%     end
% end
%
% % Handle the file extension (must be '.nii' for save_nii)
% [p, f, e] = fileparts(fName);
% switch e
%     case '.gz'
%         fName = fullfile(p,f);
%     case '.nii'
%         % If a filename was passed in and it was something'.nii' then we
%         % assume the user does not want the file compressed.
%         if fname_supplied; compress = false; end
%     otherwise
%         fName = fullfile(p,[f '.nii']);
% end
%
%
% %% Deal with old VISTASOFT nifti structures.
% if isfield(ni,'data')
%   % This is likely to be a old VISTASOFT nifti-1 structure, see
%   % niftiVista2ni.m
%   ni = niftiVista2ni(ni);
% end
%
%
% %% Save the file to disk using Shen's code.
% save_nii(ni,fName);
%
%
% %% GZip the file
% if compress
%     gzip(fName);
%
%     % Delete the unzipped file created by save_nii.m:
%     delete(fName);
% end
%
% return


text:
function geometry_test011 ( )

%*****************************************************************************80
%
%% TEST011 tests CIRCLE_DIA2IMP_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST011\n' );
  fprintf ( 1, '  CIRCLE_DIA2IMP_2D converts a diameter to an\n' );
  fprintf ( 1, '  implicit circle in 2D.\n' );

  theta = 2.0;

  p1(1,1) = 2.0 + 5.0 * cos ( theta );
  p1(2,1) = 3.0 + 5.0 * sin ( theta );

  p2(1,1) = 2.0 - 5.0 * cos ( theta );
  p2(2,1) = 3.0 - 5.0 * sin ( theta );

  r8vec_print ( dim_num, p1, '  P1:' )
  r8vec_print ( dim_num, p2, '  P2:' )

  [ r, center ] = circle_dia2imp_2d ( p1, p2 );

  circle_imp_print_2d ( r, center, '  The implicit circle:' );

  return
end


text:
function [merged_timestamps, merged_signals] = ...
    merge_signals(sig, offset)
    % Merge two signals sample-by-sample given the samples and offsets

    num_signals = length(sig);
    timestamps = cell(num_signals, 1);
    for i = 1:length(sig)
        timestamps{i} = offset(i):offset(i)+length(sig{i}) - 1;
    end

    if num_signals == 2
        [merged_timestamps, merged_signals] = ...
            merge_samples(timestamps{1}, sig{1}, ...
                          timestamps{2}, sig{2});
    elseif length(sig) == 4
        [merged_timestamps, merged_signals] = ...
            merge_samples4(timestamps{1}, sig{1}, ...
                           timestamps{2}, sig{2}, ...
                           timestamps{3}, sig{3}, ...
                           timestamps{4}, sig{4});
    end
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function [mu sigma2] = estimateGaussian(X)
%ESTIMATEGAUSSIAN This function estimates the parameters of a
%Gaussian distribution using the data in X
%   [mu sigma2] = estimateGaussian(X),
%   The input X is the dataset with each n-dimensional data point in one row
%   The output is an n-dimensional vector mu, the mean of the data set
%   and the variances sigma^2, an n x 1 vector
%

% Useful variables
[m, n] = size(X);

% You should return these values correctly
mu = zeros(n, 1);
sigma2 = zeros(n, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the mean of the data and the variances
%               In particular, mu(i) should contain the mean of
%               the data for the i-th feature and sigma2(i)
%               should contain variance of the i-th feature.
%

mu = sum(X) ./ m
sigma2 = sum((X - (mu' * ones(1, m))').^2) ./ m

% =============================================================

end

author: xjwhhh

repo: AndrewNgMachineLearning


text:
function c = struct2namevaluecell(s)
    % STRUCT2NAMEVALUECELL converts a struct into a cell of name-value pairs.eg. {name1, value1, name2, value2,...}
    % ex.
    %   origcell = {'name','fred','age',23}
    %   s = struct(origcell{:});
    %   c = struct2namevaluecell(s)
    %    % so that...
    %   isequal(origcell, c);   % TRUE
    %

    fn = fieldnames(s)';
    ns = numel(s);
    c = cell(numel(s), numel(fn)*2);
    c(:,1:2:end)=repmat(fn,ns,1);
    for r = 1:ns
        for f = 1:numel(fn)
        c(r,f*2) = {s(r).(fn{f})};
        end
    end
end

author: CelsoReyes


text:
function colorMap = diffMap(cMapSize)
% matRad difference colormap
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

colorMapData = ...
    [0         0    1.0000
    0.0323    0.0323    1.0000
    0.0645    0.0645    1.0000
    0.0968    0.0968    1.0000
    0.1290    0.1290    1.0000
    0.1613    0.1613    1.0000
    0.1935    0.1935    1.0000
    0.2258    0.2258    1.0000
    0.2581    0.2581    1.0000
    0.2903    0.2903    1.0000
    0.3226    0.3226    1.0000
    0.3548    0.3548    1.0000
    0.3871    0.3871    1.0000
    0.4194    0.4194    1.0000
    0.4516    0.4516    1.0000
    0.4839    0.4839    1.0000
    0.5161    0.5161    1.0000
    0.5484    0.5484    1.0000
    0.5806    0.5806    1.0000
    0.6129    0.6129    1.0000
    0.6452    0.6452    1.0000
    0.6774    0.6774    1.0000
    0.7097    0.7097    1.0000
    0.7419    0.7419    1.0000
    0.7742    0.7742    1.0000
    0.8065    0.8065    1.0000
    0.8387    0.8387    1.0000
    0.8710    0.8710    1.0000
    0.9032    0.9032    1.0000
    0.9355    0.9355    1.0000
    0.9677    0.9677    1.0000
    1.0000    1.0000    1.0000
    1.0000    1.0000    1.0000
    1.0000    0.9677    0.9677
    1.0000    0.9355    0.9355
    1.0000    0.9032    0.9032
    1.0000    0.8710    0.8710
    1.0000    0.8387    0.8387
    1.0000    0.8065    0.8065
    1.0000    0.7742    0.7742
    1.0000    0.7419    0.7419
    1.0000    0.7097    0.7097
    1.0000    0.6774    0.6774
    1.0000    0.6452    0.6452
    1.0000    0.6129    0.6129
    1.0000    0.5806    0.5806
    1.0000    0.5484    0.5484
    1.0000    0.5161    0.5161
    1.0000    0.4839    0.4839
    1.0000    0.4516    0.4516
    1.0000    0.4194    0.4194
    1.0000    0.3871    0.3871
    1.0000    0.3548    0.3548
    1.0000    0.3226    0.3226
    1.0000    0.2903    0.2903
    1.0000    0.2581    0.2581
    1.0000    0.2258    0.2258
    1.0000    0.1935    0.1935
    1.0000    0.1613    0.1613
    1.0000    0.1290    0.1290
    1.0000    0.0968    0.0968
    1.0000    0.0645    0.0645
    1.0000    0.0323    0.0323
    1.0000         0         0];

if nargin < 1
    colorMap = colorMapData;
elseif size(colorMapData,1) == cMapSize
    colorMap = colorMapData;
else
    %We have to interpolate the colormap
    newX = linspace(1,64,cMapSize);
    oldX = 1:64;
    colorMap = interp1(oldX,colorMapData,newX);
    %{
    %resample via HSV.. more color-true than above, but doesn't work with
    %every colormap
    hsv                        = rgb2hsv(cm);
    hsv(144:end,1)             = hsv(144:end,1)+1;
    ColorMap                   = interp1(linspace(0,1,size(cm,1)),hsv,linspace(0,1,cMapSize));
    ColorMap(cm(:,1)>1,1) = ColorMap(cm(:,1)>1,1)-1;
    ColorMap                   = hsv2rgb(cm);
    %}
end

end

repo: matRad


text:
%Simple function to test the cross method
d=10;
%elem_fun=@(x) sum(x); %Just sum of everything
p=0:d-1; p = 2.^p;
a=-5;b=5;
n=2^d;
h=(b-a)/(n-1);
%mv=@(x) x.^3;
mv=@(x) sqrt(x)+abs(x);
%elem_fun=@(x) 1.0./(dot((x-1),p)+1e-3); %Just sum of everything
%elem_fun=@(x) mv(1e-12+dot(x-1,p)*h);

%Compare functions of TT-tensors
x=tt_x(d,2); x=tt_tensor(x);
e=tt_ones(d,2);e=tt_tensor(e);
x=a*e+h*x; x=round(x,1e-13);
x1=kron(e,x); x2=kron(x,e);
rs=x1.^2+x2.^2; rs=round(rs,1e-13);
%fun=@(x) 1.0./sqrt(x);
%fun=@(x) 1.0./x;
fun=@(x) exp(-x.^4);
%fun = @(x) x.^2;
rs=x;
%elem_fun = @(ind) fun(rs(ind));
%elem_fun=@(ind) rs(ind);
%elem_fun=@(x) sqrt(x(1))+x(2);
elem_fun = @(ind) fun(rs(ind));
eps=1e-6;
D=ndims(rs);
%y=tt_rc2(2*d,2,elem_fun,1e-12);
y=tt_rc(D,2,elem_fun,1e-6,'nswp',40);
%v=tt_rand(size(y),ndims(y),2);
%y=y+v; y=round(y,1e-12);
%y1=tt_rc(D,2,elem_fun,1e-8,'nswp',40,'x0',y);

z=funcrs2(rs,fun,1e-6,rs,20);
z=round(z,1e-12);
z1=tt_rc(D,2,elem_fun,1e-8,'nswp',40,'x0',z);

author: oseledets

repo: TT-Toolbox


text:
classdef EpsilonSequence < IncrementalSequence

    methods (Access = protected)

        function generateAlphaSequence(obj)
            frac = 2;
            kmax = ceil(log10(obj.x0/obj.x1)/log10(frac));
            obj.alpha = obj.initialValue./frac.^(1:kmax);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function cmap = crameri(ColormapName,varargin)
% crameri returns perceptually-uniform scientific colormaps created
% by Fabio Crameri.
%
%% Syntax
%
%  crameri
%  cmap = crameri('ColormapName')
%  cmap = crameri('-ColormapName')
%  cmap = crameri(...,NLevels)
%  cmap = crameri(...,'pivot',PivotValue)
%  crameri(...)
%
%% Description
%
% crameri without any inputs displays the options for colormaps.
%
% cmap = crameri('ColormapName') returns a 256x3 colormap.  For a visual
% depiction of valid colormap names, type |crameri|.
%
% cmap = crameri('-ColormapName') a minus sign preceeding any ColormapName flips the
% order of the colormap.
%
% cmap = crameri(...,NLevels) specifies a number of levels in the colormap.  Default
% value is 256.
%
% cmap = crameri(...,'pivot',PivotValue) centers a diverging colormap such that white
% corresponds to a given value and maximum extents are set using current caxis limits.
% If no PivotValue is set, 0 is assumed.
%
% crameri(...) without any outputs sets the current colormap to the current axes.
%
%% Examples
% For examples, type:
%
%  showdemo crameri_documentation
%
%% Author Info
% This function was written by Chad A. Greene of the University of Texas
% Institute for Geophysics (UTIG), August 2018, using Fabio Crameri's
% scientific colormaps, version 4.0. http://www.fabiocrameri.ch/colourmaps.php
%
%% Citing this colormap:
% Please acknowledge the free use of these colormaps by citing
%
% Crameri, F. (2018). Scientific colour-maps. Zenodo. http://doi.org/10.5281/zenodo.1243862
%
% Crameri, F. (2018), Geodynamic diagnostics, scientific visualisation and
% StagLab 3.0, Geosci. Model Dev., 11, 2541-2562, doi:10.5194/gmd-11-2541-2018.
%
% For more on choosing effective and accurate colormaps for science, be sure
% to enjoy this fine beach reading:
%
% Thyng, K.M., C.A. Greene, R.D. Hetland, H.M. Zimmerle, and S.F. DiMarco. 2016. True
% colors of oceanography: Guidelines for effective and accurate colormap selection.
% Oceanography 29(3):9-13, http://dx.doi.org/10.5670/oceanog.2016.66.
%
% See also colormap and caxis.

%% Display colormap options:

if nargin==0
   figure('menubar','none','numbertitle','off','Name','crameri options:')

   if license('test','image_toolbox')
      imshow(imread('crameri7.0.png'));
   else
      axes('pos',[0 0 1 1])
      image(imread('crameri7.0.png'));
      axis image off
   end

   return
end

%% Error checks:

assert(isnumeric(ColormapName)==0,'Input error: ColormapName must be a string.')

%% Set defaults:

NLevels = 256;
autopivot = false;
PivotValue = 0;
InvertedColormap = false;

%% Parse inputs:

% Does user want to flip the colormap direction?
dash = strncmp(ColormapName,'-',1);
if any(dash)
   InvertedColormap = true;
   ColormapName(dash) = [];
end

% Standardize all colormap names to lowercase:
%ColormapName = lower(ColormapName);

% Oleron's too hard for me to remember, so I'm gonna use dem or topo.
if ismember(ColormapName,{'dem','topo'})
   ColormapName = 'oleron';
end

% Does the user want to center a diverging colormap on a specific value?
% This parsing support original 'zero' syntax and current 'pivot' syntax.
tmp = strncmpi(varargin,'pivot',3);
if any(tmp)
   autopivot = true;
   try
      if isscalar(varargin{find(tmp)+1})
         PivotValue = varargin{find(tmp)+1};
         tmp(find(tmp)+1) = 1;
      end
   end
   varargin = varargin(~tmp);
end

% Has user requested a specific number of levels?
tmp = isscalar(varargin);
if any(tmp)
   NLevels = varargin{tmp};
end

%% Load RGB values and interpolate to NLevels:

try
   S = load('CrameriColourMaps7.0.mat',ColormapName);
   cmap = S.(ColormapName);
catch
   error(['Unknown colormap name ''',ColormapName,'''. Try typing crameri with no inputs to check the options and try again.'])
end

% Interpolate if necessary:
if NLevels~=size(cmap,1)
   cmap = interp1(1:size(cmap,1), cmap, linspace(1,size(cmap,1),NLevels),'linear');
end

%% Invert the colormap if requested by user:

if InvertedColormap
   cmap = flipud(cmap);
end

%% Adjust values to current caxis limits?

if autopivot
   clim = caxis;
   maxval = max(abs(clim-PivotValue));
   cmap = interp1(linspace(-maxval,maxval,size(cmap,1))+PivotValue, cmap, linspace(clim(1),clim(2),size(cmap,1)),'linear');
end

%% Clean up

if nargout==0
   colormap(gca,cmap)
   clear cmap
end

%% Code to collect Fabio's data into a single .mat file:
% Unzip the latest folder, navigate to that filepath, and run this.
% Update the file list as needed.
%
% clear all
% f = {'acton','bam','bamO','bamako','batlow','batlowK','batlowW','berlin','bilbao','broc','brocO','buda','bukavu','cork',...
%    'corkO','davos','devon','fes','grayC','hawaii','imola','lajolla','lapaz','lisbon',...
%    'nuuk','oleron','oslo','roma','romaO','tofino','tokyo','turku','vanimo','vik','vikO'};
%
% for k = 1:length(f)
%    load([f{k},'/',f{k},'.mat'])
% end
%
% clear f k
% save('CrameriColourMaps7.0.mat')

author: opencobra

repo: cobratoolbox


text:
function  varargout = formulate_data(source_1, source_2, eI, mode)
% Copyright (c) 2014-present University of Illinois at Urbana-Champaign
% All rights reserved.
%
% Developed by:     Po-Sen Huang, Paris Smaragdis
%                   Department of Electrical and Computer Engineering
%                   Department of Computer Science
%
% eI.winSize: Size of window
% eI.seqLen: unique lengths (in ascending order)
%             files are chopped by these lengths. (ex: [1, 10, 100])
% eI.targetWhiten: Specify the path to the whitening data.
% data_ag: noisy data. cell array of different training lengths.
% target_ag: clean data. cell array of different training lengths.

% mode:
%     0: Training (noisy data and clean data)
%     1: Testing (just noisy data)
%     2: Error testing (noisy data and clean data, both loaded without chunking)

%% Testing code
input_fnames = {};
unique_lengths = [];

%% Set up. During testing, dont know the lengths so cant pre-allocate
eI.fs = 16000;

if mode,
    data_ag = {};
    target_ag = {};  % returns empty targets
    mixture_ag = {};
    item=0;
else % training -- chunk
    seqLenSizes = zeros(1,length(eI.seqLen));

    assert(numel(source_1) == numel(source_2));

    for ifile=1:numel(source_1)
%         if (strcmp(train_files(ifile).name,'.') || strcmp(train_files(ifile).name,'..')), continue; end

%         [train_wav,fs]=wavread([eI.DataPath,'train',filesep,train_files(ifile).name]);
%         eI.fs=fs;
%         train1 = train_wav(:,2); % singing
%         train2 = train_wav(:,1); % music
        train1 = source_1{ifile};
        train2 = source_2{ifile};

%         maxLength=max([length(train1), length(train2)]);
%         train1(end+1:maxLength)=eps;
%         train2(end+1:maxLength)=eps;

%         train1=train1./sqrt(sum(train1.^2));
%         train2=train2./sqrt(sum(train2.^2));

        shift_size=min(eI.circular_step, numel(train2)-1);
        for ioffset=1:shift_size:numel(train2)-shift_size % circle shift
            train2_shift = [train2(ioffset: end), train2(1: ioffset-1)];
            fprintf('%d %d, %d %d, %d\n', ioffset, numel(train2), 1, ioffset-1, numel(train2_shift));
            dmix = train1+ train2_shift;

            % input feature calculate
            [DATA, mixture_spectrum, eI]=compute_features_stft2(dmix, eI);

            [T, innfeat] = size(DATA');
            nfeat=size(mixture_spectrum,1);

            remainder = T;
            for i=length(eI.seqLen):-1:1
              num = floor(remainder/eI.seqLen(i));
              remainder = mod(remainder,eI.seqLen(i));
              seqLenSizes(i) = seqLenSizes(i)+num;
            end
        end

        data_ag = cell(1,length(eI.seqLen));
        target_ag = cell(1,length(eI.seqLen));
        mixture_ag = cell(1,length(eI.seqLen));
        for i=length(eI.seqLen):-1:1
            data_ag{i} = zeros(eI.inputDim*eI.seqLen(i),seqLenSizes(i));
            mixture_ag{i} = zeros(nfeat*eI.seqLen(i),seqLenSizes(i));
            if eI.cleanonly==1,
                target_ag{i} = zeros(nfeat*eI.seqLen(i),seqLenSizes(i));
            else
                target_ag{i} = zeros(2*nfeat*eI.seqLen(i),seqLenSizes(i));
            end

        end

    end
end

seqLenPositions = ones(1,length(eI.seqLen));
% winsize=eI.winsize;
nFFT= eI.nFFT; hop= eI.hop; scf=eI.scf;
% windows=sin(0:pi/winsize:pi-pi/winsize);
wn = sqrt( hann( eI.nFFT, 'periodic')); % hann window

 for ifile=1:numel(source_1)
%     if (strcmp(train_files(ifile).name,'.') || strcmp(train_files(ifile).name,'..')), continue; end

%     [train_wav,fs]=wavread([eI.DataPath,'train',filesep,train_files(ifile).name]);
%     eI.fs=fs;
%     train1 = train_wav(:,2); % singing
%     train2 = train_wav(:,1); % music
    train1 = source_1{ifile};
    train2 = source_2{ifile};

%     maxLength=max([length(train1), length(train2)]);
%     train1(end+1:maxLength)=eps;
%     train2(end+1:maxLength)=eps;

%     train1=train1./sqrt(sum(train1.^2));
%     train2=train2./sqrt(sum(train2.^2));

    shift_step=min(eI.circular_step,numel(train2)-1);
    for ioffset=1:shift_step:numel(train2)-shift_step % circle shift
        train2_shift = [train2(ioffset: end), train2(1: ioffset-1)];
        fprintf('%d %d, %d %d, %d\n', ioffset, numel(train2), 1, ioffset-1, numel(train2_shift));
        dmix = train1+ train2_shift;
        [DATA, mixture_spectrum, eI]=compute_features_stft2(dmix, eI);

%         spectrum.signal =scf* stft(train1, nFFT ,windows, hop, eI.fs);
%         spectrum.noise =scf* stft(train2_shift, nFFT ,windows, hop, eI.fs);

        if eI.MFCCorlogMelorSpectrum==2 || eI.MFCCorlogMelorSpectrum==3
            spectrum.signal = scf* stft2(train1, nFFT ,hop, 0, wn);
            spectrum.noise = scf* stft2(train2_shift, nFFT ,hop, 0, wn);
        else
            spectrum.signal = scf * stft3(train1, nFFT, hop, 0, wn, [], 0);
            spectrum.noise = scf* stft3(train2_shift, nFFT ,hop, 0, wn, [], 0);
        end

        multi_data = DATA;
        [nFeat,T] = size(multi_data);

        %% input normalize
        if eI.inputL1==1, % DATA (NUMCOFS x nSamp)
            % apply CMVN to targets (normalize such that the freq bean equal to zero mean var 1)
            cur_mean = mean(multi_data, 2);
            cur_std = std(multi_data, 0, 2);
            multi_data = bsxfun(@minus, multi_data, cur_mean);
            multi_data = bsxfun(@rdivide, multi_data, cur_std);
        elseif eI.inputL1==2,  % at each time frame, freq sum to 1
            l1norm = sum(multi_data,1)+eps;
            multi_data = bsxfun(@rdivide, multi_data, l1norm);
        end

        %% output normalize
         if eI.outputL1==1
            if eI.cleanonly==1,
                l1norm = sum(abs(spectrum.signal),1)+eps;
                spectrum.signal = bsxfun(@rdivide, spectrum.signal, l1norm);
            else
                l1norm = sum([abs(spectrum.noise); abs(spectrum.signal)],1)+eps;
                spectrum.noise = bsxfun(@rdivide, spectrum.noise, l1norm);
                spectrum.signal = bsxfun(@rdivide, spectrum.signal, l1norm);
            end
         end

         if eI.cleanonly==1,
            clean_data=abs(spectrum.signal);
         else
            clean_data=[abs(spectrum.noise); abs(spectrum.signal)];
         end
        %% zero pad
        if eI.num_contextwin > 1
            % winSize must be odd for padding to work
            if mod(eI.num_contextwin,2) ~= 1
                fprintf(1,'error! winSize must be odd!');
                return
            end;
            % pad with repeated frames on both sides so im2col data
            % aligns with output data
            nP = (eI.num_contextwin-1)/2;
            multi_data = [repmat(multi_data(:,1),1,nP), multi_data, ...
                repmat(multi_data(:,end),1,nP)];
        end;
        %% im2col puts winSize frames in each column
        multi_data_slid = im2col(multi_data,[nFeat, eI.num_contextwin],'sliding');
        % concatenate noise estimate to each input
         if mode == 1, % Testing
            c = find(unique_lengths == T);
            if isempty(c)
                % add new unique length if necessary
                data_ag = [data_ag, multi_data_slid(:)];
                unique_lengths = [unique_lengths, T];
            else
                data_ag{c} = [data_ag{c}, multi_data_slid(:)];
            end;
        elseif mode == 2, % Error analysis.
            c = find(unique_lengths == T);
            if isempty(c)
                % add new unique length if necessary
                data_ag = [data_ag, multi_data_slid(:)];
                target_ag = [target_ag, clean_data(:)];
                mixture_ag = [mixture_ag, mixture_spectrum(:)];
                unique_lengths = [unique_lengths, T];
            else
                data_ag{c} = [data_ag{c}, multi_data_slid(:)];
                target_ag{c} = [target_ag{c}, clean_data(:) ];
                mixture_ag{c} = [mixture_ag{c}, mixture_spectrum(:) ];
            end;
         elseif mode == 3 % formulate one data per cell
            item =item+1;
            % feadim x nframes
            data_ag{item} = multi_data_slid;
            target_ag{item} = clean_data;
            mixture_ag{item} = mixture_spectrum;
         else % training (mode == 0)
            %% put it in the correct cell area.
            while T > 0
                % assumes length in ascending order.
                % Finds longest length shorter than utterance
                c = find(eI.seqLen <= T, 1,'last');

                binLen = eI.seqLen(c);
                assert(~isempty(c),'could not find length bin for %d',T);
                % copy data for this chunk
                data_ag{c}(:,seqLenPositions(c))=reshape(multi_data_slid(:,1:binLen),[],1);
                target_ag{c}(:,seqLenPositions(c))=reshape(clean_data(:,1:binLen),[],1);
                mixture_ag{c}(:,seqLenPositions(c))=reshape(mixture_spectrum(:,1:binLen),[],1);

                seqLenPositions(c) = seqLenPositions(c)+1;
                % trim for next iteration
                T = T-binLen;
                if T > 0
                    multi_data_slid = multi_data_slid(:,(binLen+1):end);
                    clean_data = clean_data(:,(binLen+1):end);
                    mixture_spectrum = mixture_spectrum(:,(binLen+1):end);
                end;
            end;
        end;
    end;
 end

theoutputs = {data_ag, target_ag, mixture_ag};
varargout = theoutputs(1:nargout);

return;

%% Unit test
% % (TODO) add
% eI.DataPath=['.',filesep];
%
% eI.MFCCorlogMelorSpectrum=2; % 0- mfcc, 1- logmel, 2- spectrum
% eI.winsize = 1024;    eI.nFFT = 1024;    eI.hop =eI.winsize/2;    eI.scf=1;%scf = 2/3;
% eI.featDim =513;
% eI.num_contextwin=3;
% eI.inputDim = eI.featDim * eI.num_contextwin;
%
% train_files= dir( [eI.DataPath, 'train',filesep,'*wav']);
%
% eI.seqLen = [1 50 100];
% eI.inputL1=0;eI.outputL1=0;
% eI.circular_step=100000; %eI.hop = eI.winsize/2
% formulate_data(train_files, eI, 0)

%% Unit test
eI.DataPath=['.',filesep];

eI.MFCCorlogMelorSpectrum=2; % 0- mfcc, 1- logmel, 2- spectrum
eI.winsize = 1024;    eI.nFFT = 1024;    eI.hop =eI.winsize/2;    eI.scf=1;%scf = 2/3;
eI.featDim =513;
eI.num_contextwin=3;
eI.inputDim = eI.featDim * eI.num_contextwin;

train_files= dir( [eI.DataPath, 'train',filesep,'*wav']);

eI.seqLen = [1 50 100];
eI.inputL1=0;eI.outputL1=0;
eI.circular_step=100000; %eI.hop = eI.winsize/2
formulate_data(source_1, source_2, eI, 0)
end

author: posenhuang

repo: deeplearningsourceseparation


text:
% File:   itload.m
% Brief:  Load an IT++ itfile content to Matlab/Octave workspace
% Author: Tony Ottosson and Adam Piatyszek
%
% Usage: itload("fname.it")
%
% This functions loads all variables from an IT++ file format to the
% Matlab/Octave workspace.
%
% -------------------------------------------------------------------------
%
% Copyright (C) 1995-2010  (see AUTHORS file for a list of contributors)
%
% This file is part of IT++ - a C++ library of mathematical, signal
% processing, speech processing, and communications classes and functions.
%
% IT++ 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.
%
% IT++ 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 IT++.  If not, see <http://www.gnu.org/licenses/>.
%
% -------------------------------------------------------------------------

function itload(fname)

[fid, err_msg] = fopen(fname, 'rb', 'ieee-le');
if (fid == -1)
  fname = [fname '.it'];
  [fid, err_msg2] = fopen(fname, 'rb', 'ieee-le');
  if (fid == -1)
    error(err_msg);
  end
end

% Check file size
fseek(fid, 0, 'eof');
file_size = ftell(fid);
fseek(fid, 0, 'bof');

% Read "IT++" magic string
[d, n] = fread(fid, 5, 'char');
if (n ~= 5 | d(1:4) ~= [73 84 43 43]')
  error('Not an IT++ file!');
end

% Check the IT++ file version
if (d(5) ~= 3)
  error('Only IT++ file version 3 is supported by this function!');
end

while (1)
  p = ftell(fid); % Save current file position
  [d1, n] = fread(fid, 3, 'uint64'); % Read header, data, and total block sizes
  name = fgetstr(fid); % Read current variable name
  type = fgetstr(fid); % Read current variable type
  fseek(fid, p+d1(1), 'bof'); % Skip header bytes

  if (length(type) == 0) % A deleted entry -> skip it

  % --- bin ---
  elseif (strcmp(type, 'bin'))
    [d, n] = fread(fid, 1, 'char');
    assignin('caller', name, d);
  % --- int8 (char) ---
  elseif (strcmp(type, 'int8'))
    [d, n] = fread(fid, 1, 'int8');
    assignin('caller', name, d);
  % --- int16 (short) ---
  elseif (strcmp(type, 'int16'))
    [d, n] = fread(fid, 1, 'int16');
    assignin('caller', name, d);
  % --- int32 (int) ---
  elseif (strcmp(type, 'int32'))
    [d, n] = fread(fid, 1, 'int32');
    assignin('caller', name, d);
  % --- float32 (float) ---
  elseif (strcmp(type, 'float32'))
    [d, n] = fread(fid, 1, 'float32');
    assignin('caller', name, d);
  % --- float64 (double) ---
  elseif (strcmp(type, 'float64'))
    [d, n] = fread(fid, 1, 'float64');
    assignin('caller', name, d);
  % --- cfloat32 (complex<float>) ---
  elseif (strcmp(type, 'cfloat32'))
    [d, n] = fread(fid, 2, 'float32');
    d = complex(d(1), d(2));
    assignin('caller', name, d);
  % --- cfloat64 (complex<double>) ---
  elseif (strcmp(type, 'cfloat64'))
    [d, n] = fread(fid, 2, 'float64');
    d = complex(d(1), d(2));
    assignin('caller', name, d);

  % --- bvec ---
  elseif (strcmp(type, 'bvec'))
    [size, n] = fread(fid, 1, 'uint64');
    [d, n] = fread(fid, size, 'char');
    assignin('caller', name, d);
  % --- svec ---
  elseif (strcmp(type, 'svec'))
    [size, n] = fread(fid, 1, 'uint64');
    [d, n] = fread(fid, size, 'int16');
    assignin('caller', name, d);
  % --- ivec ---
  elseif (strcmp(type, 'ivec'))
    [size, n] = fread(fid, 1, 'uint64');
    [d, n] = fread(fid, size, 'int32');
    assignin('caller', name, d);
  % --- fvec ---
  elseif (strcmp(type, 'fvec'))
    [size, n] = fread(fid, 1, 'uint64');
    [d, n] = fread(fid, size, 'float32');
    assignin('caller', name, d);
  % --- dvec ---
  elseif (strcmp(type, 'dvec'))
    [size, n] = fread(fid, 1, 'uint64');
    [d, n] = fread(fid, size, 'float64');
    assignin('caller', name, d);
  % --- fcvec ---
  elseif (strcmp(type, 'fcvec'))
    [size, n] = fread(fid, 1, 'uint64');
    [d, n] = fread(fid, size*2, 'float32');
    d = complex(d(1:2:end), d(2:2:end));
    assignin('caller', name, d);
  % --- dcvec ---
  elseif (strcmp(type, 'dcvec'))
    [size, n] = fread(fid, 1, 'uint64');
    [d, n] = fread(fid, size*2, 'float64');
    d = complex(d(1:2:end), d(2:2:end));
    assignin('caller', name, d);
  % --- string ---
  elseif (strcmp(type, 'string'))
    [size, n] = fread(fid, 1, 'uint64');
    [d, n] = fread(fid, size, 'char');
    d = char(d);
    assignin('caller', name, d);

  % --- bmat ---
  elseif (strcmp(type, 'bmat'))
    [size, n] = fread(fid, 2, 'uint64');
    [d, n] = fread(fid, size', 'char');
    assignin('caller', name, d);
  % --- smat ---
  elseif (strcmp(type, 'smat'))
    [size, n] = fread(fid, 2, 'uint64');
    [d, n] = fread(fid, size', 'int16');
    assignin('caller', name, d);
  % --- imat ---
  elseif (strcmp(type, 'imat'))
    [size, n] = fread(fid, 2, 'uint64');
    [d, n] = fread(fid, size', 'int32');
    assignin('caller', name, d);
  % --- fmat ---
  elseif (strcmp(type, 'fmat'))
    [size, n] = fread(fid, 2, 'uint64');
    [d, n] = fread(fid, size', 'float32');
    assignin('caller', name, d);
  % --- dmat ---
  elseif (strcmp(type, 'dmat'))
    [size, n] = fread(fid, 2, 'uint64');
    [d, n] = fread(fid, size', 'float64');
    assignin('caller', name, d);
  % --- fcmat ---
  elseif (strcmp(type, 'fcmat'))
    [size, n] = fread(fid, 2, 'uint64');
    [d, n] = fread(fid, size(1)*size(2)*2, 'float32');
    d = reshape(complex(d(1:2:end), d(2:2:end)), size(1), size(2));
    assignin('caller', name, d);
  % --- dcmat ---
  elseif (strcmp(type, 'dcmat'))
    [size, n] = fread(fid, 2, 'uint64');
    [d, n] = fread(fid, size(1)*size(2)*2, 'float64');
    d = reshape(complex(d(1:2:end), d(2:2:end)), size(1), size(2));
    assignin('caller', name, d);

  % --- bArray ---
  elseif (strcmp(type, 'bArray'))
    [size, n] = fread(fid, 1, 'uint64');
    [d, n] = fread(fid, size, 'char');
    assignin('caller', name, d);
  % --- sArray ---
  elseif (strcmp(type, 'sArray'))
    [size, n] = fread(fid, 1, 'uint64');
    [d, n] = fread(fid, size, 'int16');
    assignin('caller', name, d);
  % --- iArray ---
  elseif (strcmp(type, 'iArray'))
    [size, n] = fread(fid, 1, 'uint64');
    [d, n] = fread(fid, size, 'int32');
    assignin('caller', name, d);
  % --- fArray ---
  elseif (strcmp(type, 'fArray'))
    [size, n] = fread(fid, 1, 'uint64');
    [d, n] = fread(fid, size, 'float32');
    assignin('caller', name, d);
  % --- dArray ---
  elseif (strcmp(type, 'dArray'))
    [size, n] = fread(fid, 1, 'uint64');
    [d, n] = fread(fid, size, 'float64');
    assignin('caller', name, d);
  % --- fcArray ---
  elseif (strcmp(type, 'fcArray'))
    [size, n] = fread(fid, 1, 'uint64');
    [d, n] = fread(fid, size*2, 'float32');
    d = complex(d(1:2:end), d(2:2:end));
    assignin('caller', name, d);
  % --- dcArray ---
  elseif (strcmp(type, 'dcArray'))
    [size, n] = fread(fid, 1, 'uint64');
    [d, n] = fread(fid, size*2, 'float64');
    d = complex(d(1:2:end), d(2:2:end));
    assignin('caller', name, d);

  % --- bvecArray ---
  elseif (strcmp(type, 'bvecArray'))
    [size, n] = fread(fid, 1, 'uint64');
    clear d2;
    for i=1:size;
      [size2, n] = fread(fid, 1, 'uint64');
      [d, n] = fread(fid, size2, 'char');
      d2{i} = d;
    end
    assignin('caller', name, d2);
  % --- svecArray ---
  elseif (strcmp(type, 'svecArray'))
    [size, n] = fread(fid, 1, 'uint64');
    clear d2;
    for i=1:size;
      [size2, n] = fread(fid, 1, 'uint64');
      [d, n] = fread(fid, size2, 'int16');
      d2{i} = d;
    end
    assignin('caller', name, d2);
  % --- ivecArray ---
  elseif (strcmp(type, 'ivecArray'))
    [size, n] = fread(fid, 1, 'uint64');
    clear d2;
    for i=1:size;
      [size2, n] = fread(fid, 1, 'uint64');
      [d, n] = fread(fid, size2, 'int32');
      d2{i} = d;
    end
    assignin('caller', name, d2);
  % --- vecArray ---
  elseif (strcmp(type, 'vecArray'))
    [size, n] = fread(fid, 1, 'uint64');
    clear d2;
    for i=1:size;
      [size2, n] = fread(fid, 1, 'uint64');
      [d, n] = fread(fid, size2, 'float64');
      d2{i} = d;
    end
    assignin('caller', name, d2);
  % --- cvecArray ---
  elseif (strcmp(type, 'cvecArray'))
    [size, n] = fread(fid, 1, 'uint64');
    clear d2;
    for i=1:size;
      [size2, n] = fread(fid, 1, 'uint64');
      [d, n] = fread(fid, size2*2, 'float64');
      d2{i} = complex(d(1:2:end), d(2:2:end));
    end
    assignin('caller', name, d2);
  % --- stringArray ---
  elseif (strcmp(type, 'stringArray'))
    [size, n] = fread(fid, 1, 'uint64');
    clear d2;
    for i=1:size;
      [size2, n] = fread(fid, 1, 'uint64');
      [d, n] = fread(fid, size2, 'char');
      d2{i} = char(d);
    end
    assignin('caller', name, d2);

  % --- bmatArray ---
  elseif (strcmp(type, 'bmatArray'))
    [size, n] = fread(fid, 1, 'uint64');
    clear d2;
    for i=1:size;
      [size2, n] = fread(fid, 2, 'uint64');
      [d, n] = fread(fid, size2', 'char');
      d2{i} = d;
    end
    assignin('caller', name, d2);
  % --- smatArray ---
  elseif (strcmp(type, 'smatArray'))
    [size, n] = fread(fid, 1, 'uint64');
    clear d2;
    for i=1:size;
      [size2, n] = fread(fid, 2, 'uint64');
      [d, n] = fread(fid, size2', 'int16');
      d2{i} = d;
    end
    assignin('caller', name, d2);
  % --- imatArray ---
  elseif (strcmp(type, 'imatArray'))
    [size, n] = fread(fid, 1, 'uint64');
    clear d2;
    for i=1:size;
      [size2, n] = fread(fid, 2, 'uint64');
      [d, n] = fread(fid, size2', 'int32');
      d2{i} = d;
    end
    assignin('caller', name, d2);
  % --- matArray ---
  elseif (strcmp(type, 'matArray'))
    [size, n] = fread(fid, 1, 'uint64');
    clear d2;
    for i=1:size;
      [size2, n] = fread(fid, 2, 'uint64');
      [d, n] = fread(fid, size2', 'float64');
      d2{i} = d;
    end
    assignin('caller', name, d2);
  % --- cmatArray ---
  elseif (strcmp(type, 'cmatArray'))
    [size, n] = fread(fid, 1, 'uint64');
    clear d2;
    for i=1:size;
      [size2, n] = fread(fid, 2, 'uint64');
      [d, n] = fread(fid, size2(1)*size2(2)*2, 'float64');
      d2{i} = reshape(complex(d(1:2:end), d(2:2:end)), size2(1), size2(2));
    end
    assignin('caller', name, d2);

  % --- else ---
  else
    warning(['Not supported type: ' type]);
  end

  if (p + d1(3) >= file_size)
    break;
  else
    fseek(fid, p+d1(3), 'bof');
  end

end

fclose(fid);

function str = fgetstr(fid)
str = '';
while (1)
  [d, n] = fread(fid, 1, 'char');
  if (d == 0)
    break;
  end
  str = [str char(d)];
end

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function [xsens , xsens_fields ] = xsens_read(file, row_lines, columns)
% xsens_read: reads .cvs files created by Xsens software.
%
% INPUT
%   fname: .cvs file name (string).
%   row_lines: number of rows in file.
%   columns: maximum column to be read.
%
% OUTPUT
%   xsens: data structure where each field is a column found in the file.
%   xsens_fields: cell vector with the names of each field found in the file.
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   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 Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% Reference:
%
%
% Version: 002
% Date:    2020/10/30
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

disp('xsens_read: processing data...')

%% OPEN FILE

xsens_f = fopen(file, 'r');
if xsens_f == -1
    error('xsens_read: ERROR: %s not found', file)
else
    fprintf('xsens_read: Loading %s...\n', file)
end

fprintf('xsens_read: %s has %d lines.\n', file, row_lines)

%% HEADER

header_rows = 13;

xsens_header = cell(header_rows,1);
for a = 1:header_rows
    xsens_header{a,1} = fgets(xsens_f);
end

%% FIELDS

xsens_line = fgets(xsens_f);
xsens_line = strrep(xsens_line,'[','_');    % Chance '[' for '_'
xsens_line = strrep(xsens_line,']','');     % Delete ']'

xsens_fields = textscan(xsens_line, '%s', 'Delimiter', ',' , 'EmptyValue', 0);

xsens_fields = xsens_fields{1};
xsens_fields(83) = {'GNSSLongitude'};
xsens_fields(84) = {'GNSSLatitude'};

[N,~] = size (xsens_fields);

fprintf('xsens_read: %s has %d columns.\n', file, N);

fprintf('xsens_read: only %d columns will be processed.\n', columns);

%% GET DATA

% Set pointer back to the beginning of file
% fseek(xsens_f,0,'bof');

% Pattern for textscan()
formatSpec = repmat('%f ', 1, N);

% Correction of numbers of rows
row_lines = row_lines - header_rows;

% Preallocation
xsens_imu = zeros(row_lines, columns);
xsens_gnss = zeros(floor(row_lines/2), columns);
xsens_ahrs = zeros(floor(row_lines/2), columns);
xsens_mag = zeros(floor(row_lines/2), columns);
xsens_att = zeros(floor(row_lines/2), columns);
xsens_baro = zeros(floor(row_lines/2), columns);

% Index for data allocation inside the while loop.
a = 1;  % index for reading IMU data
% b
c = 1;  % index for reading GNSS data
% d
e = 1;  % index for reading MAG data
% f
g = 1;  % index for reading AHRS data
% h
i = 1;  % index for reading ATTITUDE data
% j
k = 1;  % index for reading BARO data
% l

rows_block = 75000;		% Default vauel for rows_block

if rows_block > row_lines
    rows_block = floor(row_lines / 5);
end

fprintf('xsens_read: processing %s...\n', file)

loop_line = 1;
loop_ctr = 0;

while ~feof(xsens_f)

    fprintf('xsens_read: processing line %d...\n',  loop_line)
    loop_ctr = loop_ctr + 1;
    loop_line = loop_ctr * rows_block;

    data_c = textscan(xsens_f, formatSpec, rows_block, 'Delimiter',',');
    data = cell2mat(data_c(1:columns));

    % Delete rows where date is NaN for IMU data
    idl = ~isnan(data(: , 18));     % Column 18 is Acc_X
    data_imu = data(idl , :);

    % Delete rows where the difference between two adjacent times is negative
    idl = diff(data_imu(: , 2)) > 0;
    idl = [true; idl];
    data_imu = data_imu(idl , :);

    % Delete rows where date is NaN for GNSS data
    idl = ~isnan(data(: , 70));     % Column 70 is NUTimeOfWeek
    data_gnss = data(idl , :);

    % Delete rows where date is NaN for MAG data
    idl = ~isnan(data(: , 33));     % Column 33 is Mag_X
    data_mag = data(idl , :);

    % Delete rows where date is NaN for AHRS data
    idl = ~isnan(data(: , 24));     % Column 24 is AccHR_X
    data_ahrs = data(idl , :);

    % Delete rows where date is NaN for ATTITUDE data
    idl = ~isnan(data(: , 44));     % Column 44 is Quat_q0
    data_att = data(idl , :);

    % Delete rows where date is NaN for BARO data
    idl = ~isnan(data(: , 43));     % Column 43 is Pressure
    data_baro = data(idl , :);

    [M, ~] = size( data_imu );
    b = a + M - 1;
    xsens_imu(a:b , :) = data_imu;
    a = b + 1;

    [M, ~] = size( data_gnss );
    d = c + M - 1;
    xsens_gnss(c:d , :) = data_gnss;
    c = d + 1;

    [M, ~] = size( data_mag );
    f = e + M - 1;
    xsens_mag(e:f , :) = data_mag;
    e = f + 1;

    [M, ~] = size( data_ahrs );
    h = g + M - 1;
    xsens_ahrs(g:h , :) = data_ahrs;
    g = h + 1;

    [M, ~] = size( data_att);
    j = i + M - 1;
    xsens_att(i:j , :) = data_att;
    i = j + 1;

    [M, ~] = size( data_baro);
    l = k + M - 1;
    xsens_baro(k:l , :) = data_baro;
    k = l + 1;

    clear data_c data data_vld idl data_imu data_gnss data_mag data_ahrs

end

fclose(xsens_f);

% Delete epmty rows
idl = xsens_imu(: , 1) ~= 0.0 ;
xsens_imu = xsens_imu(idl , :);

idl = xsens_gnss(: , 1) ~= 0.0 ;
xsens_gnss = xsens_gnss(idl , :);

idl = xsens_mag(: , 1) ~= 0.0 ;
xsens_mag = xsens_mag(idl , :);

idl = xsens_ahrs(: , 1) ~= 0.0 ;
xsens_ahrs = xsens_ahrs(idl , :);

idl = xsens_att(: , 1) ~= 0.0 ;
xsens_att = xsens_att(idl , :);

idl = xsens_baro(: , 1) ~= 0.0 ;
xsens_baro = xsens_baro(idl , :);

%% MATRIX TO STRUCTURE

xsens = struct;

% IMU DATA
for a=1:23
    xsens.(xsens_fields{a}) = xsens_imu(:, a);
end

for a=27:29
    xsens.(xsens_fields{a}) = xsens_imu(:, a);
end

for a=36:42
    xsens.(xsens_fields{a}) = xsens_imu(:, a);
end

for a=60:66
    xsens.(xsens_fields{a}) = xsens_imu(:, a);
end

% AHRS DATA
for a=24:26
    xsens.(xsens_fields{a}) = xsens_ahrs(:, a);
end

for a=30:32
    xsens.(xsens_fields{a}) = xsens_ahrs(:, a);
end

% MAG DATA
for a=33:35
    xsens.(xsens_fields{a}) = xsens_mag(:, a);
end

% ATTITUDE DATA
for a=44:59
    xsens.(xsens_fields{a}) = xsens_att(:, a);
end

% BARO DATA
xsens.(xsens_fields{43}) = xsens_baro(:, 43);

% GNSS DATA
for a=70:columns
    xsens.(xsens_fields{a}) = xsens_gnss(:, a);
end

xsens.SampleTimeFine        = xsens_imu(:, 2);
xsens.AHRSSampleTimeFine    = xsens_ahrs(:, 2);
xsens.MagSampleTimeFine     = xsens_mag(:, 2);
xsens.AttitudeSampleTimeFine = xsens_att(:, 2);
xsens.BaroSampleTimeFine    = xsens_baro(:, 2);
xsens.GNSSSampleTimeFine    = xsens_gnss(:, 2);

% FREQUENCIES
dti = median(diff(xsens.SampleTimeFine)) * 10e-5;
xsens.freq_imu = (1/dti);

dta = median(diff(xsens.AHRSSampleTimeFine)) * 10e-5;
xsens.freq_ahrs = (1/dta);

dtm = median(diff(xsens.MagSampleTimeFine)) * 10e-5;
xsens.freq_mag = (1/dtm);

dtt = median(diff(xsens.AttitudeSampleTimeFine)) * 10e-5;
xsens.freq_att = (1/dtt);

dtb = median(diff(xsens.BaroSampleTimeFine)) * 10e-5;
xsens.freq_baro = (1/dtb);

dtg = median(diff(xsens.NUTimeOfWeek)) * 10e-4;
xsens.freq_gnss = (1/dtg);

% Header
xsens.header = xsens_header;

disp('xsens_read: data is ready.')

end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function scannerXform = getScannerXform(rawFile)
% Get the transform matrix from INPLANE coordinates (x,y,z, indices) to
% scanner coordinates (in mm).
%
%  scannerXform = getScannerXform([rawFile])
%
% The transform will yield:
%   scannerCoords = scannerXform * ipCoords;
%
% 8/2009: Siphoned off from RAS's code computeB0DirectionMap.
%
% That code is now shorter. And we can now use this transform for other things
% such as getting a 3-vector for the B0 direction in a scan
%
% Example: scannerXform = getScannerXform;

% Find the raw inplane dicoms
if notDefined('rawFile'),
    rawFile = fullpath('Raw/Anatomy/Inplane');
end

if exist(rawFile, 'dir')
	% allow rawFile to be a directory containing DICOM files
	pattern = fullfile(rawFile, '*.dcm');
	w = dir(pattern);
	if isempty(w)
		pattern = fullfile(rawFile, '*.DCM');
		w = dir(pattern);
	end

	if isempty(w)
		% still not found? We have a problem.
		error(['No DICOM files found in %s. Please provide a ' ...
			'file path if the raw file is not in DICOM format.'], ...
			rawFile);
	end

	% got here? then we have some files. Take the first one.
	rawFile = fullfile(rawFile, w(1).name);

elseif ~exist(rawFile, 'file')
	% no file or directory found with this name? Can't proceed.
	error('%s is not a file or directory.', rawFile);
end

% load the raw file and make sure we have an xform into scanner space
[p f ext] = fileparts(rawFile);
if strncmpi( ext, '.dcm', 4 )
	inplane = mrReadDicom(rawFile);
else
	inplane = mrLoad(rawFile);
end

if ~isfield(inplane, 'spaces') || isempty(inplane.spaces)
	error('No scanner-coordinate information found in the header fields');
end

% find the alignment to the scanner coordinate space
spaces = {inplane.spaces.name};
I = cellfind(spaces, 'Scanner');
if isempty(I)
	error('No scanner-coordinate information found in the header fields');
end

scannerXform = inplane.spaces(I).xform;

author: vistalab

repo: vistasoft


text:
function A = my_spdiags(v,diagonals,m,n)
%MY_SPDIAGS Construct a banded matrix from diagonal input
%
% A = my_spdiags(v,diagonals,m,n);
%
% Inputs:
%  v  the diagonals to be inserted in the matrix - each column of v is one
%     diagonal in diagonals
%  diagonals  the index of diagonals requested
%  m,n  the dimension of the requested matrix
% Outputs:
%  A  the requested banded sparse matrix
%

A = ...

end

author: odedstein

repo: sgi-introduction-course


text:
function myerode(mask_file,new_mask_file)
% This code is written by Liang Zhan (zhan.liang@ucla.edu)
% Input variables
% 1) make_file name (nifti format, no gz format)
% 2) new_mask_file name (this is output in nifti format)

nii=load_untouch_nii(mask_file);
mask=double(nii.img);
dim1=size(mask,1); dim2=size(mask,2); dim3=size(mask,3);
for zs=1:dim1
   mask(zs,:,:)=bwmorph(squeeze(mask(zs,:,:)),'erode');
end
for zs=1:dim2
   mask(:,zs,:)=bwmorph(squeeze(mask(:,zs,:)),'erode');
end
for zs=1:dim3
mask(:,:,zs)=bwmorph(mask(:,:,zs),'erode');
end
nii.img=mask;
save_untouch_nii(nii,new_mask_file);
end

author: HennyJie

repo: BrainGB

save_path: github-repos/MATLAB/HennyJie-BrainGB


text:
function [dt,t1] = dtiLoadDt6(fname, applyBrainMask)
%Load a dt6 (tensor) data file.
%
%  [dt,t1] = dtiLoadDt6(fname, [applyBrainMask=true]);
%
% Optionally, you can also load the T1 data file and apply a brain mask.
%
% The six entries are the diffusion tensor values, derived from the raw
% data.  The raw data can be in many different directions.  The diffusion
% tensor is a 3x3, symmetric, positive-definite matrix, D.  The six
% independent entries in the matrix are stored in a vector: (Dxx Dyy Dzz
% Dxy Dxz Dyz)
%
% The project-style dt6- contain no data, just pointers to data files.  The
% directory architecture of the data should be:
%
%
% HISTORY:
%  2007.10.29 RFD: wrote it, pulling parts from dtiLoadDT6. The dt6 loading
%  was a mess. Hopefully it is better now.
%
% (c) Stanford VISTA Team

%% Initialize the variables
dt = [];   % This is the main dt structure.
t1 = [];
if(~exist('fname','var') || isempty(fname))
    fname = [pwd filesep];
end
if(isdir(fname))
    [f,p] = uigetfile({'*.mat'}, 'Select a dt6 file...',fname);
    if(isequal(f,0)), disp('Load dt6 canceled.'); return; end
    fname = fullfile(p,f);
end
if(~exist('applyBrainMask','var') || isempty(applyBrainMask))
    applyBrainMask = true;
end

%% Load the dt6 information.
%  This information describes things like where the data files are.
[dt6Struct, dt.dataFile, defaultPath, dataDir] = dtiLoadDt6Info(fname);
[p,n] = fileparts(dataDir);
dt.subName = n;
[p,n] = fileparts(p);
if(~isempty(n))
    dt.subName = [n '_' dt.subName];
    [p,n] = fileparts(p);
    if(~isempty(n))
        dt.subName = [n '_' dt.subName];
    end
end

% Decide whether we will output the T1 data
if(nargout>1)
    loadT1 = true;
    t1.talairachScale = [];
else
    loadT1 = false;
end

%% Read the data
% We support both new NIFTI-based dt6 files and the old dt6 file format.
if(isfield(dt6Struct,'files'))   % Test whether new dt6 stye
    %
    % New project-style dt6- contains no data, just pointers to data files.
    %
    [dt.dt6, dt.xformToAcpc, dt.mmPerVoxel] = dtiLoadTensorsFromNifti(dt6Struct.files.tensors);
    dt.adcUnits = '';

    clear ni;
    ni = niftiRead(dt6Struct.files.b0);
    ni.data(isnan(ni.data)) = 0;
    dt.b0 = ni.data;

    clear ni;
    ni = niftiRead(dt6Struct.files.brainMask);
    dt.brainMask = ni.data;

    if(loadT1 && isfield(dt6Struct.files,'t1'))
        t1Fname = dt6Struct.files.t1;
        if(exist(t1Fname,'file'))
            ni = niftiRead(t1Fname);
            t1.img = ni.data;
            t1.xformToAcpc = ni.qto_xyz;
            t1.mmPerVoxel = ni.pixdim(1:3);
        else
            % try another level up
            t1Fname = dt6Struct.files.t1;
            if(exist(t1Fname,'file'))
                ni = niftiRead(t1Fname);
                t1.img = ni.data;
                t1.xformToAcpc = ni.qto_xyz;
                t1.mmPerVoxel = ni.pixdim(1:3);
            else
                disp(['Couldn''t load t1 file "' t1Fname '".']);
            end
        end
        % look for a t1 brain mask in there too
        t1MaskName = fullfile(fileparts(t1Fname),'t1_mask.nii.gz');
		if(~exist(t1MaskName,'file'))
		  t1MaskName = fullfile(fileparts(t1Fname),'t1_brain_mask.nii.gz');
		end
        if(exist(t1MaskName,'file'))
            ni = niftiRead(t1MaskName);
            t1.brainMask = ni.data;
            t1.brainMaskXform = ni.qto_xyz;
        end
    end
    dt.files = dt6Struct.files;
else
    %
    % It's an old dt6 file- the data are already loaded, but we need to
    % check things and pack the relevant structs.
    %
    if(~isfield(dt6Struct, 'xformToAcPc')), dt.xformToAcpc = eye(4);
    else dt.xformToAcpc = dt6Struct.xformToAcPc; end
    dt.dt6 = dt6Struct.dt6;
    dt.mmPerVoxel = dt6Struct.mmPerVox;

    dt6Struct.b0(isnan(dt6Struct.b0)) = 0;
    dt.b0 = dt6Struct.b0;

    if(isfield(dt6Struct,'brainMask'))
        dt.brainMask = dt6Struct.brainMask;
    elseif(isfield(dt6Struct,'dtBrainMask'))
        dt.brainMask = dt6Struct.dtBrainMask;
    else
        dt.brainMask = [];
    end

    if(isfield(dt6Struct, 'adcUnits'))
        dt.adcUnits = dt6Struct.adcUnits;
    else
        dt.adcUnits = '';
    end

    if(loadT1 && isfield(dt6Struct,'anat'))
        t1.img = double(dt6Struct.anat.img);
        t1.xformToAcpc = dt6Struct.anat.xformToAcPc;
        t1.mmPerVoxel = dt6Struct.anat.mmPerVox;
        if(isfield(dt6Struct.anat, 'brainMask'))
            t1.brainMask = dt6Struct.anat.brainMask;
            t1.brainMaskXform = dt6Struct.anat.xformToAcPc;
        end

        if(isfield(dt6Struct.anat, 'talScale'))
            t1.talairachScale = dt6Struct.anat.talScale;
        end
    end
    dt.files = [];
end

%% Set up the T1 data
if(loadT1 && isempty(t1.talairachScale) && isfield(t1,'brainMask') && ~isempty(t1.brainMask))
    % Approximate Talairach scales using brain mask
    disp('Computing Talairach scales from brain mask...');
    t1.talairachScale = mrAnatGetTalairachScalesFromMask(t1.brainMask, t1.brainMaskXform);
end

%% NaN's take forever on Pentium 4 CPUs!
% (eg. http://www.psychology.nottingham.ac.uk/staff/cr1/simd.html)
dt.dt6(isnan(dt.dt6)) = 0;

%% Apply the brain mask, if it exists
if(applyBrainMask && ~isempty(dt.brainMask))
    for ii=1:6
        tmp = dt.dt6(:,:,:,ii);
        tmp(~dt.brainMask) = 0;
        dt.dt6(:,:,:,ii) = tmp;
    end
end
if(isempty(dt.brainMask))
    dt.brainMask = dt.dt6(:,:,:,1)>0;
end

%% Convert adc units to our standard micrometers^2 / millisecond
if(isempty(dt.adcUnits))
    [dt.adcUnits,dt.adcScale,stdUnitStr] = dtiGuessDiffusivityUnits(dt.dt6);
    if(dt.adcScale~=1)
        dt.dt6 = dt.dt6*dt.adcScale;
    end
    dt.adcUnits = stdUnitStr;
end

%% Apply a spatial transformation
if(isfield(dt6Struct,'xformVAnatToAcpc'))
    dt.xformVAnatToAcpc = dt6Struct.xformVAnatToAcpc;
elseif(isfield(dt6Struct,'xformToMrVista') && exist('t1','var') && ~isempty(t1))
    % try to use an old vanatXform
    warning('mrDiffusion:deprecated','Using old xformToMrVista field- you may want to recompute this using Xform->Compute mrVista xform');
    dt.xformVAnatToAcpc = t1.xformToAcpc*dt6Struct.xformToMrVista;
end
if(loadT1 && isfield(dt6Struct,'t1NormParams'))
    t1.t1NormParams = dt6Struct.t1NormParams;
end

%% Set the default bounding box
dt.bb = sort(mrAnatXformCoords(dt.xformToAcpc, [1 1 1; size(dt.brainMask)]));
if(any(dt.mmPerVoxel<1))
    dt.renderMm = dt.mmPerVoxel;
else
    dt.renderMm = [1 1 1];
end

return;