text:
function [ idx, tri ] = nearest_neighbor ( varargin )

%NEAREST_NEIGHBOR    find nearest neighbors
%   IDX = NEAREST_NEIGHBOR(X) finds the nearest neighbor by Euclidean
%   distance to each point (column) in X from X. X is a matrix with points
%   as columns. IDX is a vector of indices into X, such that X(:, IDX) are
%   the nearest neighbors to X. e.g. the nearest neighbor to X(:, 2) is
%   X(:, IDX(2))
%
%   IDX = NEAREST_NEIGHBOR(P, X) finds the nearest neighbor by Euclidean
%   distance to each point in P from X. P and X are both matrices with the
%   same number of rows, and points are the columns of the matrices. Output
%   is a vector of indices into X such that X(:, IDX) are the nearest
%   neighbors to P
%
%   IDX = NEAREST_NEIGHBOR(I, X) where I is a logical vector or vector of
%   indices, and X has at least two rows, finds the nearest neighbor in X
%   to each of the points X(:, I).
%   I must be a row vector to distinguish it from a single point.
%   If X has only one row, the first input is treated as a set of 1D points
%   rather than a vector of indices
%
%   IDX = NEAREST_NEIGHBOR(..., Property, Value)
%   Calls NEAREST_NEIGHBOR with the indicated parameters set. Property
%   names can be supplied as just the first letters of the property name if
%   this is unambiguous, e.g. NEAREST_NEIGHBOR(..., 'num', 5) is equivalent
%   to NEAREST_NEIGHBOR(..., 'NumberOfNeighbors', 5). Properties are case
%   insensitive, and are as follows:
%      Property:                         Value:
%      ---------                         ------
%         NumberOfNeighbors             natural number, default 1
%            NEAREST_NEIGHBOR(..., 'NumberOfNeighbors', K) finds the closest
%            K points in ascending order to each point, rather than the
%            closest point. If Radius is specified and there are not
%            sufficient numbers, fewer than K neighbors may be returned
%
%         Radius                         positive, default +inf
%            NEAREST_NEIGHBOR(..., 'Radius', R) finds neighbors within
%            radius R. If NumberOfNeighbors is not set, it will find all
%            neighbors within R, otherwise it will find at most
%            NumberOfNeighbors. The IDX matrix is padded with zeros if not
%            all points have the same number of neighbors returned. Note
%            that specifying a radius means that the Delaunay method will
%            not be used.
%
%         DelaunayMode                   {'on', 'off', |'auto'|}
%            DelaunayMode being set to 'on' means NEAREST_NEIGHBOR uses the
%            a Delaunay triangulation with dsearchn to find the points, if
%            possible. Setting it to 'auto' means NEAREST_NEIGHBOR decides
%            whether to use the triangulation, based on efficiency. Note
%            that the Delaunay triangulation will not be used if a radius
%            is specified.
%
%         Triangulation                  Valid triangulation produced by
%                                        delaunay or delaunayn
%            If a triangulation is supplied, NEAREST_NEIGHBOR will attempt
%            to use it (in conjunction with dsearchn) to find the
%            neighbors.
%
%   [IDX, TRI] = NEAREST_NEIGHBOR( ... )
%   If the Delaunay Triangulation is used, TRI is the triangulation of X'.
%   Otherwise, TRI is an empty matrix
%
%   Example:
%
%     % Find the nearest neighbor in X to each column of X
%     x = rand(2, 10);
%     idx = nearest_neighbor(x);
%
%     % Find the nearest neighbors to each point in p
%     p = rand(2, 5);
%     x = rand(2, 20);
%     idx = nearest_neighbor(p, x)
%
%     % Find the five nearest neighbors to points x(:, [1 6 20]) in x
%     x = rand(4, 1000)
%     idx = nearest_neighbor([1 6 20], x, 'NumberOfNeighbors', 5)
%
%     % Find all neighbors within radius of 0.1 of the points in p
%     p = rand(2, 10);
%     x = rand(2, 100);
%     idx = nearest_neighbor(p, x, 'r', 0.1)
%
%     % Find at most 10 nearest neighbors to point p from x within a
%     % radius of 0.2
%     p = rand(1, 2);
%     x = rand(2, 30);
%     idx = nearest_neighbor(p, x, 'n', 10, 'r', 0.2)
%
%
%   See also DELAUNAYN, DSEARCHN, TSEARCH

%TODO    Allow other metrics than Euclidean distance
%TODO    Implement the Delaunay mode for multiple neighbors

% Copyright 2006 Richard Brown. This code may be freely used and
% distributed, so long as it maintains this copyright line
error(nargchk(1, Inf, nargin, 'struct'));

% Default parameters
userParams.NumberOfNeighbors = []    ; % Finds one
userParams.DelaunayMode       = 'auto'; % {'on', 'off', |'auto'|}
userParams.Triangulation      = []    ;
userParams.Radius             = inf   ;

% Parse inputs
[P, X, fIndexed, userParams] = parseinputs(userParams, varargin{:});

% Special case uses Delaunay triangulation for speed.

% Determine whether to use Delaunay - set fDelaunay true or false
nX  = size(X, 2);
nP  = size(P, 2);
dim = size(X, 1);

switch lower(userParams.DelaunayMode)
    case 'on'
        %TODO Delaunay can't currently be used for finding more than one
        %neighbor
        fDelaunay = userParams.NumberOfNeighbors == 1 && ...
            size(X, 2) > size(X, 1)                    && ...
            ~fIndexed                                  && ...
            userParams.Radius == inf;
    case 'off'
        fDelaunay = false;
    case 'auto'
        fDelaunay = userParams.NumberOfNeighbors == 1 && ...
            ~fIndexed                                  && ...
            size(X, 2) > size(X, 1)                    && ...
            userParams.Radius == inf                   && ...
            ( ~isempty(userParams.Triangulation) || delaunaytest(nX, nP, dim) );
end

% Try doing Delaunay, if fDelaunay.
fDone = false;
if fDelaunay
    tri = userParams.Triangulation;
    if isempty(tri)
        try
            tri   = delaunayn(X');
        catch
            msgId = 'Nearest_Neighbor:DelaunayFail';
            msg = ['Unable to compute delaunay triangulation, not using it. ',...
                'Set the DelaunayMode parameter to ''off'''];
            warning(msgId, msg);
        end
    end
    if ~isempty(tri)
        try
            idx = dsearchn(X', tri, P')';
            fDone = true;
        catch
            warning('Nearest_Neighbor:DSearchFail', ...
                'dsearchn failed on triangulation, not using Delaunay');
        end
    end
else % if fDelaunay
    tri = [];
end

% If it didn't use Delaunay triangulation, find the neighbors directly by
% finding minimum distances
if ~fDone
    idx = zeros(userParams.NumberOfNeighbors, size(P, 2));

    % Loop through the set of points P, finding the neighbors
    Y = zeros(size(X));
    for iPoint = 1:size(P, 2)
        x = P(:, iPoint);

        % This is the faster than using repmat based techniques such as
        % Y = X - repmat(x, 1, size(X, 2))
        for i = 1:size(Y, 1)
            Y(i, :) = X(i, :) - x(i);
        end

        % Find the closest points, and remove matches beneath a radius
        dSq = sum(abs(Y).^2, 1);
        iRad = find(dSq < userParams.Radius^2);
        if ~fIndexed
            iSorted = iRad(minn(dSq(iRad), userParams.NumberOfNeighbors));
        else
            iSorted = iRad(minn(dSq(iRad), userParams.NumberOfNeighbors + 1));
            iSorted = iSorted(2:end);
        end

        % Remove any bad ones
        idx(1:length(iSorted), iPoint) = iSorted';
    end
    %while ~isempty(idx) && isequal(idx(end, :), zeros(1, size(idx, 2)))
    %    idx(end, :) = [];
    %end
    idx( all(idx == 0, 2), :) = [];
end % if ~fDone
if isvector(idx)
    idx = idx(:)';
end
end % nearest_neighbor

%DELAUNAYTEST   Work out whether the combination of dimensions makes
%fastest to use a Delaunay triangulation in conjunction with dsearchn.
%These parameters have been determined empirically on a Pentium M 1.6G /
%WinXP / 512MB / Matlab R14SP3 platform. Their precision is not
%particularly important
function tf = delaunaytest(nx, np, dim)
switch dim
    case 2
        tf = np > min(1.5 * nx, 400);
    case 3
        tf = np > min(4 * nx  , 1200);
    case 4
        tf = np > min(40 * nx , 5000);

        % if the dimension is higher than 4, it is almost invariably better not
        % to try to use the Delaunay triangulation
    otherwise
        tf = false;
end % switch
end % delaunaytest

%MINN   find the n most negative elements in x, and return their indices
%  in ascending order
function I = minn(x, n)

% Make sure n is no larger than length(x)
n = min(n, length(x));

% Sort the first n
[xsn, I] = sort(x(1:n));

% Go through the rest of the entries, and insert them into the sorted block
% if they are negative enough
for i = (n+1):length(x)
    j = n;
    while j > 0 && x(i) < xsn(j)
        j = j - 1;
    end

    if j < n
        % x(i) should go into the (j+1) position
        xsn = [xsn(1:j), x(i), xsn((j+1):(n-1))];
        I   = [I(1:j), i, I((j+1):(n-1))];
    end
end

end %minn

%PARSEINPUTS    Support function for nearest_neighbor
function [P, X, fIndexed, userParams] = parseinputs(userParams, varargin)
if length(varargin) == 1 || ~isnumeric(varargin{2})
    P           = varargin{1};
    X           = varargin{1};
    fIndexed    = true;
    varargin(1) = [];
else
    P             = varargin{1};
    X             = varargin{2};
    varargin(1:2) = [];

    % Check the dimensions of X and P
    if size(X, 1) ~= 1
        % Check to see whether P is in fact a vector of indices
        if size(P, 1) == 1
            try
                P = X(:, P);
            catch
                error('Nearest_Neighbor:InvalidIndexVector', ...
                    'Unable to index matrix using index vector');
            end
            fIndexed = true;
        else
            fIndexed = false;
        end % if size(P, 1) == 1
    else % if size(X, 1) ~= 1
        fIndexed = false;
    end

    if ~fIndexed && size(P, 1) ~= size(X, 1)
        error('Nearest_Neighbor:DimensionMismatch', ...
            'No. of rows of input arrays doesn''t match');
    end
end
% Parse the Property/Value pairs
if rem(length(varargin), 2) ~= 0
    error('Nearest_Neighbor:propertyValueNotPair', ...
        'Additional arguments must take the form of Property/Value pairs');
end

propertyNames = {'numberofneighbors', 'delaunaymode', 'triangulation', ...
    'radius'};
while length(varargin) ~= 0
    property = varargin{1};
    value    = varargin{2};

    % If the property has been supplied in a shortened form, lengthen it
    iProperty = find(strncmpi(property, propertyNames, length(property)));
    if isempty(iProperty)
        error('Nearest_Neighbor:InvalidProperty', 'Invalid Property');
    elseif length(iProperty) > 1
        error('Nearest_Neighbor:AmbiguousProperty', ...
            'Supplied shortened property name is ambiguous');
    end
    property = propertyNames{iProperty};

    switch property
        case 'numberofneighbors'
            if rem(value, 1) ~= 0 || ...
                    value > length(X) - double(fIndexed) || ...
                    value < 1
                error('Nearest_Neighbor:InvalidNumberOfNeighbors', ...
                    'Number of Neighbors must be an integer, and smaller than the no. of points in X');
            end
            userParams.NumberOfNeighbors = value;

        case 'delaunaymode'
            fOn = strcmpi(value, 'on');
            if strcmpi(value, 'off')
                userParams.DelaunayMode = 'off';
            elseif fOn || strcmpi(value, 'auto')
                if userParams.NumberOfNeighbors ~= 1
                    if fOn
                        warning('Nearest_Neighbor:TooMuchForDelaunay', ...
                            'Delaunay Triangulation method works only for one neighbor');
                    end
                    userParams.DelaunayMode = 'off';
                elseif size(X, 2) < size(X, 1) + 1
                    if fOn
                        warning('Nearest_Neighbor:TooFewDelaunayPoints', ...
                            'Insufficient points to compute Delaunay triangulation');
                    end
                    userParams.DelaunayMode = 'off';

                elseif size(X, 1) == 1
                    if fOn
                        warning('Nearest_Neighbor:DelaunayDimensionOne', ...
                            'Cannot compute Delaunay triangulation for 1D input');
                    end
                    userParams.DelaunayMode = 'off';
                else
                    userParams.DelaunayMode = value;
                end
            else
                warning('Nearest_Neighbor:InvalidOption', ...
                    'Invalid Option');
            end % if strcmpi(value, 'off')

        case 'radius'
            if isscalar(value) && isnumeric(value) && isreal(value) && value > 0
                userParams.Radius = value;
                if isempty(userParams.NumberOfNeighbors)
                    userParams.NumberOfNeighbors = size(X, 2) - double(fIndexed);
                end
            else
                error('Nearest_Neighbor:InvalidRadius', ...
                    'Radius must be a positive real number');
            end

        case 'triangulation'
            if isnumeric(value) && size(value, 2) == size(X, 1) + 1 && ...
                    all(ismember(1:size(X, 2), value))
                userParams.Triangulation = value;
            else
                error('Nearest_Neighbor:InvalidTriangulation', ...
                    'Triangulation not a valid Delaunay Triangulation');
            end
    end % switch property

    varargin(1:2) = [];
end % while
if isempty(userParams.NumberOfNeighbors)
    userParams.NumberOfNeighbors = 1;
end
end %parseinputs

author: johannesgerer

repo: jburkardt-m


text:
function delay = PETHTransition(psth,timeBins,varargin)

%PETHTransition - Find a transition point in a peri-event time histogram.
%
% Find a transition point in a peri-stimulus time histogram using the maximum
% likelihood or mean square error algorithm of Friedman and Priebe (1998), or
% piecewise linear fit.
%
%  USAGE
%
%    delay = PETHTransition(peth,timeBins,<options>)
%
%    peth           PETH (computed using SyncHist)
%    timeBins       time bins
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'method'      either 'ml' (maximum likelihood, default), 'ls' (least
%                   squares) or 'pl' (piecewise linear)
%     'show'        either 'on' (default) or 'off'
%    =========================================================================
%
%  EXAMPLE
%
%    [raster,indices] = Sync(spikes,stimuli);     % compute spike raster data
%    figure;PlotSync(raster,indices);             % plot spike raster
%    [s,t] = SyncHist(raster,indices);            % compute PETH
%    delay = PETHTransition(s,t);                 % find transition point
%
%  SEE
%
%    See also Sync, SyncHist, SyncMap, PlotSync.

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

% Defaults
method = 'ml';
show = false;

% Check number of parameters
if nargin < 2 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help PETHTransition'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
  if ~ischar(varargin{i}),
    error(['Parameter ' num2str(i+1) ' is not a property (type ''help PETHTransition'' for details).']);
  end
  switch(lower(varargin{i})),

    case 'method',
		method = lower(varargin{i+1});
      if ~isstring_FMAT(method,'ml','ls','pl'),
        error('Incorrect value for property ''method'' (type ''help PETHTransition'' for details).');
      end

    case {'show','plot'},
    	show = lower(varargin{i+1});
    	if ~isstring_FMAT(show,'on','off'),
        error('Incorrect value for property ''show'' (type ''help PETHTransition'' for details).');
      end

    otherwise,
      error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help PETHTransition'' for details).']);

  end
end

if strcmp(method,'ml'),

	% Maximum likelihood

	f = psth;
	n = length(f);
	for i = 1:n,
		F(i) = sum(log(1:max([1 f(i)])));
	end
	for theta = 1:n-1,
		t1 = 1:theta;t1 = t1';
		t2 = theta+1:n;t2 = t2';
		lambda1 = sum(f(t1))/(theta+1);
		lambda2 = sum(f(t2))/(n-theta);
		likelihood(theta) =  ...
			- lambda1 * (theta+1) ...
			+ log(lambda1) * sum(f(t1)) ...
			- sum(F(t1)) ...
			- lambda2 * (n-theta) ...
			+ log(lambda2) * sum(f(t2)) ...
			- sum(F(t2));
	end
	likelihood(end+1) = likelihood(end); % dummy value
	d = find(likelihood == max(likelihood));
	d = d(1);
	delay = timeBins(d);

	if show,
		fig = figure;
		set(fig,'name','PETH Transition - Maximum Likelihood','number','off');
		subplot(2,1,1);hold on;
		bar(timeBins,psth);
		plot([delay delay],ylim,'k','linestyle','--');
		title(['delay = ' num2str(delay) ' s']);
		ylabel('Occurrences');
		subplot(2,1,2);hold on;
		plot(timeBins,likelihood,'Color',[1 0 0],'Marker','none','LineStyle','-');
		plot([delay delay],ylim,'k','linestyle','--');
		xlabel('Time (s)');
		ylabel('Likelihood');
	end

elseif strcmp(method,'ls'),

	% Least square error

	f0 = psth;
	n0 = length(f0);
	t0 = 1:n0;t0 = t0';
	for i = 1:n0,
		F0(i) = sum(f0(1:i-1));
	end
	F0 = F0';
	start = 1;
	stop = length(f0);
	f = f0(start:stop);
	n = length(f);
	for i = 1:length(f),
		F(i) = sum(f(1:i-1));
	end
	F = F';
	t = 1:n;t = t';
	squaredSum(1:n) = NaN;
	for theta = 2:n-1,
		t1 = 1:theta; t1 = t1';
		t2 = theta+1:n; t2 = t2';
		lambda1(theta) = ...
			(...
			- sum(t1 .* F(t1)) ...
			- theta * sum(F(t2)) ...
			+ theta * sum(theta-t2) * sum(F(t2) .* (theta-t2)) / sum((theta-t2).^2) ...
			) / (...
			+ theta^2 * sum(theta-t2) * sum(theta-t2) / sum((theta-t2).^2) ...
			- sum(t1.^2) ...
			- theta^2 * (n-theta) ...
			);
		lambda2(theta) = ...
			(...
			+ lambda1(theta) * theta * sum(theta-t2) ...
			- sum(F(t2).*(theta-t2)) ...
			) ...
			/ sum((theta-t2).^2);
		squaredSum(theta) =  ...
			+ sum((F(t1)-lambda1(theta)*t1).^2) ...
			+ sum((F(t2)-lambda1(theta)*theta+lambda2(theta)*(theta-t2)).^2);
	end
	d = find(squaredSum == min(squaredSum));
	d = d(1);
	delay = timeBins(d);

	if show,
		fig = figure;
		set(fig,'name','PETH Transition - Least Square Error','number','off');
		subplot(2,1,1);hold on;
		bar(timeBins,psth);
		plot([delay delay],ylim,'k','linestyle','--');
		title(['delay = ' num2str(delay) ' s']);
		ylabel('Occurrences');
		subplot(2,1,2);hold on;
		t1 = 1:d;t1 = t1';
		t2 = d+1:n;t2 = t2';
		plot(timeBins(t0),F0,'Color',[1 0 0],'Marker','none','LineStyle','-');
		plot(timeBins(t1+start-1),lambda1(d)*(t1)+F0(start),'Color',[0 0 0]);
		plot(timeBins(t2+start-1),lambda2(d)*(t2-d)+lambda1(d)*(d)+F0(start),'Color',[0 0 0]);
		plot([delay delay],ylim,'k','linestyle','--');
		xlabel('Time (s)');
		ylabel('Cumulative count');
	end

elseif strcmp(method,'pl')

	% Piecewise linear regression

	f = psth;

	n = length(f);
	t = 1:n;t = t';
	for theta = 2:n-2,
		t1 = 1:theta;t1 = t1';
		t2 = theta+1:n;t2 = t2';
		a1(theta) = LinearRegression(t1,f(t1));
		a2(theta) = LinearRegression(t2,f(t2),a1(theta)*theta);
		squaredSum(theta) = sum((f(t1)-a1(theta)*t1).^2) ...
			+ sum((f(t2)-a2(theta)*t2-a1(theta)*theta).^2);
	end
	squaredSum(1) = Inf;
	d = find(squaredSum == min(squaredSum));
	d = d(1);
	delay = timeBins(d);

	if show,
		fig = figure;
		set(fig,'name','PETH Transition - Piecewise Linear Fit','number','off');
		subplot(2,1,1);hold on;
		bar(timeBins,psth);
		plot([delay delay],ylim,'k','linestyle','--');
		title(['delay = ' num2str(delay) ' s']);
		ylabel('Occurrences');
		subplot(2,1,2);hold on;
		t1 = 1:d;t1 = t1';
		t2 = d+1:n;t2 = t2';
		plot(timeBins(t),f,'Color',[1 0 0],'Marker','o','LineStyle','none');
		plot(timeBins(t1),a1(d)*t1,'k');
		plot(timeBins(t2),a2(d)*t2+a1(d)*d,'k');
		plot([delay delay],ylim,'k','linestyle','--');
		xlabel('Time (s)');
		ylabel('Occurrences');
	end

end

function [a,b] = LinearRegression(x,y,b)

if nargin == 2, b = 0; end

y = y - b;
X = sum(x);
X2 = sum(x.^2);
Y = sum(y);
XY = sum(x.*y);
n = length(x);

b = (Y*X2-X*XY) / (n*X2-X.^2);
a = (n*XY-X*Y) / (n*X2-X.^2);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
classdef MaxCrossCorr < Algorithm

    methods (Access = public)

        function obj = MaxCrossCorr()
            obj.name = 'MaxCrossCorr';
            obj.inputPort = DataType.kSignal2;
            obj.outputPort = DataType.kFeature;
        end

        function result = compute(~,signal)
            result = max(xcorr(signal(:,1),signal(:,2)));
        end

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

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
laststream_range = 1+mod(laststream.smax:laststream.smax+size(laststream_chunk_clr,2)-1,laststream.buffer_len);
laststream.marker_pos(:,laststream_range) = 0;
laststream.buffer(:,laststream_range) = laststream_chunk_clr;
laststream.smax = laststream.smax + size(laststream_chunk_clr,2);

repo: REST


text:
function writeGrid(filename, grid_size, grid_spacing, pml_size, pml_alpha, Nt, dt, c_ref)
%WRITEGRID    Write grid and PML properties to a k-Wave HDF5 file.
%
% DESCRIPTION:
%       writeGrid creates and writes the wavenumber grids and PML variables
%       required by the k-Wave C++ code to the HDF5 file specified by the
%       user.
%
%       List of parameters that are written:
%           Nx
%           Ny
%           Nz
%           Nt
%           dt
%           dx
%           dy
%           dz
%           c_ref
%           ddx_k_shift_pos_r
%           ddx_k_shift_neg_r
%           ddy_k_shift_pos
%           ddy_k_shift_neg
%           ddz_k_shift_pos
%           ddz_k_shift_neg
%           x_shift_neg_r
%           y_shift_neg_r
%           z_shift_neg_r
%           pml_x_sgx
%           pml_y_sgy
%           pml_z_sgz
%           pml_x
%           pml_y
%           pml_z
%           pml_x_alpha
%           pml_y_alpha
%           pml_z_alpha
%           pml_x_size
%           pml_y_size
%           pml_z_size
%
% USAGE:
%       writeGrid(filename, grid_size, grid_spacing, pml_size, pml_alpha, Nt, dt, c_ref)
%
% INPUTS:
%       filename            - filename and location of the input HDF5 file
%       grid_size           - [Nx, Ny, Nz]
%       grid_spacing        - [dx, dy, dz]
%       pml_size            - [pml_x_size, pml_y_size, pml_z_size]
%       pml_alpha           - [pml_x_alpha, pml_y_alpha, pml_z_alpha]
%       Nt                  - number of time points
%       dt                  - time step
%       c_ref               - scalar sound speed used in the k-space
%                             operator and to define the pml variables
%
% ABOUT:
%       author              - Bradley Treeby
%       date                - 30th May 2013
%       last update         - 21st August 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also h5writeatt, writeAttributes, writeFlags, writeMatrix

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

%#ok<*INUSL>
%#ok<*NASGU>

% get literals
getH5Literals;

% unpack grid size inputs to make code easier to read
Nx          = grid_size(1);
Ny          = grid_size(2);
Nz          = grid_size(3);
dx          = grid_spacing(1);
dy          = grid_spacing(2);
dz          = grid_spacing(3);
pml_x_size  = pml_size(1);
pml_y_size  = pml_size(2);
pml_z_size  = pml_size(3);
pml_x_alpha = pml_alpha(1);
pml_y_alpha = pml_alpha(2);
pml_z_alpha = pml_alpha(3);

% =========================================================================
% CREATE WAVENUMBER AND PML VECTORS
% =========================================================================

% create the wavenumber grids (assuming Nx, Ny and Nz are even)
nx = ((-Nx/2:Nx/2-1)/Nx).';
nx(floor(Nx/2) + 1) = 0;
kx_vec = (2*pi/dx).*nx;

ny = ((-Ny/2:Ny/2-1)/Ny).';
ny(floor(Ny/2) + 1) = 0;
ky_vec = (2*pi/dy).*ny;

nz = ((-Nz/2:Nz/2-1)/Nz).';
nz(floor(Nz/2) + 1) = 0;
kz_vec = (2*pi/dz).*nz;

% force the vector operators be in the correct direction (Nx, 1, 1), (1, Ny, 1), (1, 1, Nz)
ky_vec = ky_vec.';
kz_vec = permute(kz_vec, [2 3 1]);

% create vector derivative and shift variables
ddx_k_shift_pos = ifftshift( 1i*kx_vec .* exp( 1i*kx_vec*dx/2), 1);
ddx_k_shift_neg = ifftshift( 1i*kx_vec .* exp(-1i*kx_vec*dx/2), 1);
ddy_k_shift_pos = ifftshift( 1i*ky_vec .* exp( 1i*ky_vec*dy/2), 2);
ddy_k_shift_neg = ifftshift( 1i*ky_vec .* exp(-1i*ky_vec*dy/2), 2);
ddz_k_shift_pos = ifftshift( 1i*kz_vec .* exp( 1i*kz_vec*dz/2), 3);
ddz_k_shift_neg = ifftshift( 1i*kz_vec .* exp(-1i*kz_vec*dz/2), 3);

% create vector shift operators
x_shift_neg = ifftshift( exp(-1i*kx_vec*dx/2), 1);
y_shift_neg = ifftshift( exp(-1i*ky_vec*dy/2), 2);
z_shift_neg = ifftshift( exp(-1i*kz_vec*dz/2), 3);

% create reduced variables for use with real-to-complex FFT
Nx_r                    = floor(Nx/2) + 1;
Ny_r                    = floor(Ny/2) + 1;
Nz_r                    = floor(Nz/2) + 1;
ddx_k_shift_pos_r       = ddx_k_shift_pos(1:Nx_r);
ddx_k_shift_neg_r       = ddx_k_shift_neg(1:Nx_r);
x_shift_neg_r           = x_shift_neg(1:Nx_r);
y_shift_neg_r           = y_shift_neg(1:Ny_r);
z_shift_neg_r           = z_shift_neg(1:Nz_r);

% create vector PML variables
pml_x       = getPML(Nx, dx, dt, c_ref, pml_x_size, pml_x_alpha, false, 1);
pml_x_sgx   = getPML(Nx, dx, dt, c_ref, pml_x_size, pml_x_alpha, true,  1);
pml_y       = getPML(Ny, dy, dt, c_ref, pml_y_size, pml_y_alpha, false, 2);
pml_y_sgy   = getPML(Ny, dy, dt, c_ref, pml_y_size, pml_y_alpha, true,  2);
pml_z       = getPML(Nz, dz, dt, c_ref, pml_z_size, pml_z_alpha, false, 3);
pml_z_sgz   = getPML(Nz, dz, dt, c_ref, pml_z_size, pml_z_alpha, true,  3);

% cleanup unused variables
clear ddx_k_shift_pos ddx_k_shift_neg x_shift_neg y_shift_neg z_shift_neg;

% =========================================================================
% STORE FLOATS
% =========================================================================

% list of variables stored as floats
variable_names = {...
    'dt', 'dx', 'dy', 'dz', ...
    'ddx_k_shift_pos_r', 'ddx_k_shift_neg_r', 'x_shift_neg_r', ...
    'ddy_k_shift_pos',   'ddy_k_shift_neg',   'y_shift_neg_r', ...
    'ddz_k_shift_pos',   'ddz_k_shift_neg',   'z_shift_neg_r', ...
    'pml_x_sgx',   'pml_y_sgy',   'pml_z_sgz', ...
    'pml_x',       'pml_y',       'pml_z', ...
    'pml_x_alpha', 'pml_y_alpha', 'pml_z_alpha', ...
    'c_ref'};

% change float variables to be in single precision (float in C++), then
% add to HDF5 file
for index = 1:length(variable_names)

    % cast matrix to single precision
    eval([variable_names{index} ' = ' MATRIX_DATA_TYPE_MATLAB '(' variable_names{index} ');']);

    % write to HDF5 file
    writeMatrix(filename, eval(variable_names{index}), variable_names{index});

end

% =========================================================================
% STORE INTEGERS
% =========================================================================

% integer variables
variable_names = {'Nx', 'Ny', 'Nz', 'Nt',...
    'pml_x_size' , 'pml_y_size' , 'pml_z_size'};

% change all the index variables to be in 64-bit unsigned integers (long in C++)
for index = 1:length(variable_names)

    % cast matrix to 64-bit unsigned integer
    eval([variable_names{index} ' = ' INTEGER_DATA_TYPE_MATLAB '(' variable_names{index} ');']);

    % write to HDF5 file
    writeMatrix(filename, eval(variable_names{index}), variable_names{index});

end

repo: Aero-matlab


text:
function value = circle_area_2d ( r )

%*****************************************************************************80
%
%% CIRCLE_AREA_2D returns the area of a circle in 2D.
%
%  Integration region:
%
%    Points (X,Y) such that
%
%      ( X - XC )**2 + ( Y - YC )**2 <= R**2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%
%    Output, real CIRCLE_AREA_2D, the area of the circle.
%
  value = pi * r * r;

  return
end


text:
% Script unit_test_main
% Performs the complete unit test for uniQC. All unit test specified in
% MrUnitTest are included.
%
%  unit_test_main
%
%
%   See also MrUnitTest

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2019-04-04
% Copyright (C) 2019 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

% create test suite
import matlab.unittest.TestSuite;
% run complete test suit
UTall = TestSuite.fromClass(?MrUnitTest);
resultsAll = run(UTall);
disp(table(resultsAll));

author: translationalneuromodeling

repo: tapas


text:
% This files demos creating a interactive html file from a triangular mesh
% and a second set of deformed vertices

% Written by Benjamin Irving 2013/03/25

%From example 1 a mesh can be plotted in matlab as follows:
%Creating the mesh
[x,y,z,v] = flow;
fv=isosurface(x,y,z,v,-3);

% Creating a second set of vertices from the first set... and setting the
% second and third dimensions of the vertices to 0 as one example of a
% mesh transformation
fv2=fv;
fv2.vertices(:,2)=0;
fv2.vertices(:,3)=0;

% running the html creating function
x3mesh_deform(fv.faces, fv.vertices, fv2.vertices, 'name', 'Example1')
disp('Mesh exported to html. File is in htmlfigs subfolder')

author: Sable

repo: mcbench-benchmarks


text:
%{
Copyright (c) 2015, Tom Mertens
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%}

% Contruction of Laplacian pyramid
%
% Arguments:
%   image 'I'
%   'nlev', number of levels in the pyramid (optional)
%
% tom.mertens@gmail.com, August 2007
%
%
% More information:
%   'The Laplacian Pyramid as a Compact Image Code'
%   Burt, P., and Adelson, E. H.,
%   IEEE Transactions on Communication, COM-31:532-540 (1983).
%

function pyr = laplacian_pyramid(I,nlev)

r = size(I,1);
c = size(I,2);

if ~exist('nlev')
    % compute the highest possible pyramid
    nlev = floor(log(min(r,c)) / log(2));
end

% recursively build pyramid
pyr = cell(nlev,1);
filter = pyramid_filter;
J = I;
for l = 1:nlev - 1
    % apply low pass filter, and downsample
    I = downsample(J,filter);
    odd = 2*size(I) - size(J);  % for each dimension, check if the upsampled version has to be odd
    % in each level, store difference between image and upsampled low pass version
    pyr{l} = J - upsample(I,odd,filter);
    J = I; % continue with low pass image
end
pyr{nlev} = J; % the coarest level contains the residual low pass image

author: Mericam

repo: exposure-fusion


text:
classdef ImageTile < handle
    % ImageTile This class holds a list of images.
    %
    % Upon request they can be placed into a single image in a given pattern.
    % Any non-color image will be converted to a color image.
    %
    % Copyright by Holger Caesar, 2014

    properties
        images
    end

    methods (Access = public)
        function[obj] = ImageTile()
            % [obj] = ImageTile()
            %
            % Initialize and empty image tile.

            obj.images = cell(0, 1);
        end

        function addImage(obj, image)
            % addImage(obj, image)
            %
            % Stores a copy of the image in the tile.

            obj.images{end+1, 1} = image;
        end

        function addTiling(obj, tiling)
            % addTiling(obj, tiling)
            %
            % Copies all images from that tiling to the current one.

            obj.images(end+1:end+numel(tiling.images), 1) = tiling.images;
        end

        function[totalX] = getTotalX(obj)
            totalX = numel(obj.images);
        end

        function[image] = getImage(obj, idx)
            image = obj.images{idx};
        end

        function[] = setImage(obj, idx, image)
            obj.images{idx} = image;
        end

        function[totalImage] = getTiling(obj, varargin)
            % [totalImage] = getTiling(obj, varargin)
            %
            % Compute an image that contains all images held by this tiling.

            % Get essential values
            numImages = numel(obj.images);
            sqNumImages = ceil(sqrt(numImages));

            % Parse input
            p = inputParser;
            addParameter(p, 'totalX', sqNumImages);
            addParameter(p, 'totalY', []);
            addParameter(p, 'delimiterPixels', 0);
            addParameter(p, 'minTotalSize', [0, 0]);
            addParameter(p, 'keepAspectRatio', true);
            addParameter(p, 'backgroundBlack', true);
            parse(p, varargin{:});

            totalX = p.Results.totalX;
            totalY = p.Results.totalY;
            delimiterPixels = p.Results.delimiterPixels;
            minTotalSize = p.Results.minTotalSize;
            keepAspectRatio = p.Results.keepAspectRatio;
            backgroundBlack = p.Results.backgroundBlack;

            % Early abort if there are no images
            if isempty(obj.images) || any(cellfun(@isempty, obj.images)),
                totalImage = [];
                return;
            end;

            % Compute tiling if only one size is known
            if isempty(totalY) && ~isempty(totalX),
                totalY = ceil(numImages / totalX);
            end;
            if isempty(totalX) && ~isempty(totalY),
                totalX = ceil(numImages / totalY);
            end;
            assert(totalY * totalX >= numImages);

            % Determine mean image size
            imSizes = cell2mat(cellfun(@(x) [size(x, 1), size(x, 2)], obj.images, 'UniformOutput', false));
            meanImSize = ceil(mean(imSizes, 1));

            % Determine total image size
            totalSize = meanImSize .* [totalY, totalX];

            % Correct mean image size if specified (and then meanImSize as
            % well)
            if totalSize(1) < minTotalSize(1),
                totalSize(1) = minTotalSize(1);
                meanImSize(1) = ceil(totalSize(1) / totalY);
            end;
            if totalSize(2) < minTotalSize(2),
                totalSize(2) = minTotalSize(2);
                meanImSize(2) = ceil(totalSize(2) / totalX);
            end;

            % Add delimiter pixels
            totalSize(1) = totalSize(1) + delimiterPixels * (totalY - 1);
            totalSize(2) = totalSize(2) + delimiterPixels * (totalX - 1);

            % Initialize image
            if backgroundBlack,
                totalImage = zeros(totalSize(1), totalSize(2), 3, 'uint8');
            else
                totalImage = ones(totalSize(1), totalSize(2), 3, 'uint8') * 255;
            end;

            for posIdx = 1 : numImages,
                % Get current image
                posImage = obj.images{posIdx};

                % Convert to color if necessary
                if size(posImage, 3) ~= 3,
                    posImage = repmat(posImage, [1, 1, 3]);
                end;

                % Convert to 8 bit if necessary
                posImage = im2uint8(posImage);

                if keepAspectRatio,
                    % Scale longer side to mean size and fill missing
                    % pixels with white color

                    % Resize longer side to mean size
                    imageSize = [size(posImage, 1), size(posImage, 2)];
                    imageRatio = imageSize(1) / imageSize(2);
                    meanImRatio = meanImSize(1) / meanImSize(2);
                    if imageRatio > meanImRatio,
                        % Image is too high
                        targetSize = [meanImSize(1), nan];
                    elseif imageRatio < meanImRatio,
                        % Image is too wide
                        targetSize = [nan, meanImSize(2)];
                    else
                        % Ratio is correct
                        targetSize = meanImSize;
                    end;
                    posImage = imresize(posImage, targetSize);

                    % Fill missing pixels with black color
                    fullPosImage = repmat(uint8(255), [meanImSize, 3]);
                    moveYStart = 1 + floor((meanImSize(1) - size(posImage, 1)) / 2);
                    moveYEnd = moveYStart + size(posImage, 1) - 1;
                    moveXStart = 1 + floor((meanImSize(2) - size(posImage, 2)) / 2);
                    moveXEnd = moveXStart + size(posImage, 2) - 1;
                    fullPosImage(moveYStart:moveYEnd, moveXStart:moveXEnd, :) = posImage;
                    posImage = fullPosImage;
                else
                    % Scale to mean size and destroy ratio
                    posImage = imresize(posImage, meanImSize);
                end;

                % Find out current position (go through x first!)
                posY = 1 + floor((posIdx - 1) / totalX);
                posX = posIdx - (posY-1) * totalX;

                % Insert current image into totalImage
                startY = 1 + meanImSize(1) * (posY-1) + delimiterPixels * (posY - 1);
                endY = meanImSize(1) * posY           + delimiterPixels * (posY - 1);

                startX = 1 + meanImSize(2) * (posX-1) + delimiterPixels * (posX - 1);
                endX = meanImSize(2) * posX           + delimiterPixels * (posX - 1);

                totalImage(startY:endY, startX:endX, :) = posImage;
            end;
        end
    end
end

author: nightrome

repo: matconvnet-calvin


text:
function s=xticksi(ah)
%XTIXKSI labels the x-axis of a plot using SI multipliers S=(AH)
%
%  Inputs:  AH       axis handle [default: current axes]
%
% Outputs:  S        optional global SI multiplier (see usage below)
%
% Usage:   (1) plot(...);
%              xticksi;
%
%          (2) plot(...);
%              xlabel(['Frequency (' xticksi 'Hz)']);
%
% The first form will label the tick marks on the x-axis of the current plot
% using SI multipliers where appropriate. This is particularly useful for log
% plots which MATLAB does not label very well.
% The second form will, if possible, use a single SI multiplier for all the tick
% marks; this global multiplier can be incorporated into the axis label as shown.

%	   Copyright (C) Mike Brookes 2009
%      Version: $Id: xticksi.m,v 1.2 2010/01/02 20:02:58 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~nargin
    ah=gca;
end
if nargout
s=xyzticksi(1,ah);
else
    xyzticksi(1,ah);
end

author: decouples


text:
function geometry_test179 ( )

%*****************************************************************************80
%
%% TEST179 tests SOCCER_SIZE_3D and SOCCER_SHAPE_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2005
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST179\n' );
  fprintf ( 1, '  For the truncated icosahedron, or soccer ball,\n' );
  fprintf ( 1, '  SOCCER_SIZE_3D returns dimension information;\n' );
  fprintf ( 1, '  SOCCER_SHAPE_3D returns face and order information.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We will use this information to compute the\n' );
  fprintf ( 1, '  areas and centers of each face.\n' );

  [ point_num, edge_num, face_num, face_order_max ] = soccer_size_3d ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of points =   %d\n', point_num );
  fprintf ( 1, '  Number of edges =    %d\n', edge_num );
  fprintf ( 1, '  Number of faces =    %d\n', face_num );
  fprintf ( 1, '  Maximum face order = %d\n', face_order_max );

  [ point_coord, face_order, face_point ] = soccer_shape_3d ( ...
    point_num, face_num, face_order_max );
%
%  Compute the area of each face.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Face  Order  Area\n' );
  fprintf ( 1, '\n' );

  for i = 1 : face_num

    for j = 1 : face_order(i)
      k = face_point(j,i);
      v(1:dim_num,j) = point_coord(1:dim_num,k);
    end

    [ area, normal ] = polygon_area_3d ( face_order(i), v );

    fprintf ( 1, '  %6d  %5d  %8f\n', i, face_order(i), area );

  end
%
%  Find the center of each face.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Face  Center\n' );
  fprintf ( 1, '\n' );

  for i = 1 : face_num

    vave(1:dim_num) = 0.0;

    for j = 1 : face_order(i)
      k = face_point(j,i);
      vave(1:dim_num) = vave(1:dim_num) + point_coord(1:dim_num,k)';
    end

    vave(1:dim_num) = vave(1:dim_num) / face_order(i);

    fprintf ( 1, '  %6d  %10f  %10f  %10f\n', i, vave(1:dim_num) );

  end

  return
end


text:
function y=ar_med_filter(signal,Fs)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function AR_MED_FILTER takes input SIGNAL with Sampling Frequency, %
% Fs, and applies the Yule Walker method based AR filter. The order iof the
% filter is found by Maximum kurtosis. After the application od AR filter,
% the signal is passed through Minimum Entropy Deconvolution. This combined
% AR+MED method brings out the Bearing faults hidden in Noise.
%
% The function plots two figures for AR alone and another for AR+MED
%
% Example:
%   load('s4.mat');
%   signal=s4;
%   Fs=12000;
%   ar_med_filter(signal,Fs);

% The File 's4.mat' is the vibration signal recorded from a OR faulty
% bearing with a sampling frequency of 12000Hz. The Fault frequency is 161
% Hz and is brought out.
%
% This program isa based on the paper:
% Sawalhi N, Randall RB and Endo H (2007) The enhancement of fault detection
% and diagnosis in rolling element bearings using minimum entropy
% deconvolution combined with spectral kurtosis. Mechanical Systems and
% Signal Processing. 21:2616-2633
%
% This function is basically written for Bearing fault diagnosis from
% Vibration signal.
%
%Dont forget to rate or comment on the matlab central site
%http://www.mathworks.in/matlabcentral/fileexchange/authors/258518
%
%Author:Santhana Raj.A
%https://sites.google.com/site/santhanarajarunachalam/
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sig=signal;

clear y y_e A

for order= 1:100

[A,E]=aryule(sig,order);
y(:,order)=filter(1,A,sig);
kurt(:, order)=kurtosis(y(:,order));

end

[~,index]=sort(kurt,1,'descend');
output=y(:,index(1));

figure();
N=4*2048;T=N/Fs;sig_f=abs(fft(output(1:N),N));
sig_n=sig_f/(norm(sig_f));freq_s=(0:N-1)/T;
plot(freq_s(3:350),sig_n(3:350)); title('AR alone');

%MED
[ar_med f_armed kurt_armed] = med2d(output,30,[],0.01,0);

figure();
N=4*2048;T=N/Fs;sig_f=abs(fft(ar_med(1:N),N));
sig_n=sig_f/(norm(sig_f));freq_s=(0:N-1)/T;
plot(freq_s(3:350),sig_n(3:350)); title('AR +MED');

y=ar_med;

end


text:
function [Q, R] = qr(A, econ)
%QR   QR factorization of an array-valued CHEBFUN.
%   [Q, R] = QR(A) or QR(A, 0), where A is a column CHEBFUN with n columns,
%   produces a column CHEBFUN Q with n orthonormal columns and an n x n upper
%   triangular matrix R such that A = Q*R.
%
% See also SVD, MRDIVIDE, RANK.

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

% Check inputs:
if ( (nargin == 2) && (econ ~= 0) )
    error('CHEBFUN:CHEBFUN:qr:twoargs',...
      ['Use qr(A) or qr(A, 0) for QR decomposition of an array-valued ' ...
      'CHEBFUN or quasimatrix']);
elseif ( A(1).isTransposed )
    error('CHEBFUN:CHEBFUN:qr:transpose',...
        'CHEBFUN QR works only for column CHEBFUN objects.')
elseif ( ~all(isfinite(A(1).domain)) )
    error('CHEBFUN:CHEBFUN:qr:infdomain', ...
        'CHEBFUN QR does not support unbounded domains.');
end

numCols = numColumns(A);
if ( numCols == 1 )
    % Trivial case: If A has only one column we simply scale it.
    R = sqrt(innerProduct(A, A));
    if ( R ~= 0 )
       Q = A./R;
    else
       Q = 1./sqrt(diff(A.domain)) + 0*A;
    end
    return
end

% Attempt to convert to an array-valued CHEBFUN:
[A, isArrayValued] = quasi2cheb(A);

if ( isArrayValued && (numel(A.funs) == 1) )
    % Array-valued CHEBFUN with a single FUN.

    % Call QR at the FUN level:
    [Q, R] = qr(A.funs{1});
    Q = chebfun({Q});

elseif ( isArrayValued )
    % Array-valued CHEBFUN with multiple FUNS.

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Developer note:
    %   Here we essentially use a panel-factored QR which allows us to do a QR
    %   factorization on each fun individually and then combine the result.
    %   Here's an example of this in a 2-FUN case:
    %    [A1] = [Q1*R^1] = [Q1 0][R^1] = [Q1 0][Q^1 ~][R] = [Q1*Q^1]R
    %    [A2] 1 [Q2*R^2]   [0 Q2][R^2] 2 [0 Q2][Q^2 ~][0] 3 [Q2*Q^2]
    %                                  ^
    %                               here [Q^:=Qhat, R] = qr(Rhat:=[R^1;R^2])
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    numFuns = numel(A.funs);

    % Step 1: Perform QR on each piece.
    Q = cell(numFuns, 1); Rhat = Q;
    for k = 1:numFuns
        [Q{k}, Rhat{k}] = qr(A.funs{k});
    end

    % Step 2: Compute [Qhat, R] = qr(Rhat),
    [Qhat, R] = qr(cell2mat(Rhat));
    R = R(1:numCols,:);       % Extract first block row.
    Qhat = Qhat(:,1:numCols); % Extract first block column.

    % Step 2b: Ensure the diagonal is non-negative. (A = QR = (Q*S)*(S*R))
    s = sign(diag(R)); s(~s) = 1;
    S = spdiags(s, 0, numCols, numCols);
    Qhat = Qhat*S;
    R = S*R;

    % Step 2c: Separate the segments of Qhat back into a cell.
    m = cellfun(@(v) size(v, 1), Rhat); % m(k) = length of A.FUN{k}.
    Qhat = mat2cell(Qhat, m, numCols);

    % Step 3: Fold Qhat back in to Q.
    Q = cellfun(@mtimes, Q, Qhat, 'UniformOutput', false);

    % Construct a new CHEBFUN from the computed FUNS:
    Q = chebfun(Q);

else
    % Quasimatrix case (tricky/slow):

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Developer note:
    %   Currently (5th Feb 2016) the only way this is reachable is in the
    %   case of a quasimatrix consisting of SINGFUN or DELTAFUN objects,
    %   neither of which return anything sensible when we attempt to compute
    %   a QR factorization. Try for example
    %    x = chebfun('x', [0 1]);
    %    qr([1 x sqrt(x)])
    %
    %   This case is not tested (which is OK)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Legendre-Vandermonde matrix:
    L = legpoly(0:numCols-1, domain(A), 'norm', 1);
    % Convert so that L is also quasimatrix:
    L = cheb2quasi(L);
    % Call abstract QR:
    [Q, R] = abstractQR(A, L, @innerProduct, @normest);

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function ori = transformReferenceFrame(ori,cs1,cs2)
% change reference frame of an orientation
%
% Orientations are always described with respect to a cartesian reference
% frame x, y, z aligned in a fixed way with the crystal coordiante system
% a, b, c. Typical alignments are x||a and z||c* or x||a* and z||c. This
% function allows to change the aligment of the reference frame while NOT
% changing the orientation.
%
% Syntax
%   ori = ori.transformReferenceFrame(cs)
%   mori = mori.transformReferenceFrame(cs1,cs2)
%
% Input
%  ori - @orientation
%  mori - misorientation
%  cs, cs1, cs2 - @crystalSymmetry
%

% only applicable for crystal symmetry
if ~isa(cs1,'crystalSymmetry')
  warning('Symmetry missmatch!')
  ori.CS = cs1;
  return
end

% basis transformation into reference frame
M = transformationMatrix(ori.CS,cs1);

% check symmetries are compatible
if ori.CS.id ~= cs1.id || norm(eye(3)-M*M.')>0.01 || ...
    all(norm(ori.CS.axes - cs1.axes)./norm(cs1.axes)<10^-2) || ...
    (~isempty(ori.CS.mineral) && ~isempty(cs1.mineral) && ~strcmpi(ori.CS.mineral,cs1.mineral))
  warning('Symmetry missmatch! The following crystal frames seem to be different\n\n  %s\n  %s \n',char(ori.CS,'verbose'),char(cs1,'verbose'));
end

if det(M)>10*eps
  ori = times(ori, rotation('matrix',M^-1), 0);
  ori.CS = cs1;
end

% do the same for the second symmetry
if nargin == 3, ori = inv(transformReferenceFrame(inv(ori),cs2)); end

% this is some testing code
% cs1 = crystalSymmetry('triclinic',[1 2 3],[70 80 120]*degree,'Z||a*')
% cs2 = crystalSymmetry('triclinic',[1 2 3],[70 80 120]*degree,'Z||b','X||a*')
% o = orientation.byEuler(30*degree,50*degree,120*degree,cs1)
% o * Miller(1,0,0,cs1)
% o2 = transformReferenceFrame(o,cs2)
% o2 * Miller(1,0,0,cs2)

author: mtex-toolbox

repo: mtex

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


text:
function ebsdNew = interp(ebsd, varargin)
% interpolate at arbitrary points (x,y)
%
% Syntax
%
%   ebsdNew = interp(ebsd,xNew,yNew)
%
%   ebsdNew = interp(ebsd,xNew,yNew,'method','invDist')
%
% Input
%  ebsd - @EBSD
%  xNew, yNew - new x,y coordinates
%
% Output
%  ebsdNew - @EBSD with coordinates (xNew,yNew)
%
% Options
%  method - 'invDist', 'nearest'
%
% See also
%

ebsdNew = interp(ebsd.gridify,varargin{:});


text:
function [names,isdirs] = glob(pattern,prefix)
%GLOB   Filename expansion via wildcards.
% GLOB(PATTERN) returns a cell array of file/directory names which match the
% PATTERN.
% [NAMES,ISDIRS] = GLOB(PATTERN) also returns a logical vector indicating
% which are directories.
%
% Two types of wildcards are supported:
% * matches zero or more characters, besides /
% ** matches zero or more characters, including /, ending with a /
% *** is interpreted as ** followed by *, which means it matches zero or
% more characters, including /
%
% For example, 'a*b' matches 'ab','acb','acdb', but not 'a/b'.
% 'a**b' matches 'ab','a/b','ac/b','ac/d/b', but not 'acb' or 'a/cb'.
% 'a***b' matches all of the above plus 'a/cb','ac/d/cb',etc.
%
% 'a//b' is not considered a valid filename, so 'a/*/b' will not return
% 'a//b' or 'a/b'.
%
% Examples:
%   % if 'work' is a subdirectory, this returns only 'work', not its contents:
%   glob('work')
%   % returns 'work/fun.m' (not 'fun.m'):
%   glob('work/fun.m')
%   % all m-files in 'work', prefixed with 'work/':
%   glob('work/*.m')
%   % all files named 'fun.m' in 'work' or any subdirectory of 'work':
%   glob('work/**fun.m')
%   % all m-files in 'work' or any subdirectory of 'work':
%   glob('work/***.m')
%   % all files named 'fun.m' any subdirectory of 'work' (but not 'work'):
%   glob('work/**/fun.m')
%   % all files named 'fun.m' in any subdirectory of '.':
%   glob('**fun.m')
%   % all files in all subdirectories:
%   glob('***')
%
% See also globstrings.

% Written by Tom Minka, 28-Apr-2004
% (c) Microsoft Corporation. All rights reserved.

if nargin < 2
  prefix = '';
end

names = {};
isdirs = [];
if isempty(pattern)
  return
end

% break the pattern into path components
[first,rest] = strtok(pattern,'/');

i = strfind(first,'**');
if ~isempty(i)
  % double-star pattern
  i = i(1); % process first occurrence
  rest = fullfile(first((i+2):end),rest);
  first = first(1:(i-1));
  new_pattern = [first rest];
  % if the pattern was 'a**b/c', new_pattern is 'ab/c'
  [names,isdirs] = glob(new_pattern,prefix);
  first = [first '*'];
  rest = ['**' rest];
  % if the pattern was 'a**b/c', it is now 'a*/**b/c'
end

% expand the first component
fullfirst = fullfile(prefix,first);
if ~iswild(fullfirst) & isdir(fullfirst)
  first_files = struct('name',first,'isdir',1);
else
  first_files = stripdots(dir(fullfirst));
end
% for each match, add it to the results or recurse on the rest of the pattern
for i = 1:length(first_files)
  new_prefix = fullfile(prefix,first_files(i).name);
  if isempty(rest)
    names{end+1} = new_prefix;
    isdirs(end+1) = first_files(i).isdir;
  elseif first_files(i).isdir
    [new_names, new_isdirs] = glob(rest,new_prefix);
    names = cellcat(names, new_names);
    isdirs = [isdirs; new_isdirs];
  end
end
names = names(:);
isdirs = isdirs(:);

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

function files = stripdots(files)
% omit . and .. from the results of DIR

names = {files.name};
ok = (~strcmp(names,'.') & ~strcmp(names,'..'));
files = files(ok);

function c = cellcat(c,c2)

c = {c{:} c2{:}};

function tf = iswild(pattern)

tf = ~isempty(strfind(pattern,'*'));

function s = regexp_quote(s)

regexprep(s,'[!#]^','#^');
regexprep(s,'[!#]$','#$');

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
%% Set Parameters
% Preferences
global ggamma rrho
ggamma = 2;        % coefficient of relative risk aversion
rrho = 0.01;       % rate of time preference

% Production function
global ddelta aalpha
ddelta = .025;     % capital depreciation
aalpha = 1 / 3;    % capital share

% Aggregate shock
global ssigmaTFP rrhoTFP
ssigmaTFP = .007;  % standard deviation of TFP shock
rrhoTFP = .95;     % quarterly autocorrelation of TFP shock

% Idiosyncratic shocks
global z
zz1 = 0;           % unemployed
zz2 = 1;           % employed
z = [zz1,zz2];

% Transition probabilities
global lla
llambda1 = 1 / 2;  % expected duration of unemployment is 2 quarters
llambda2 = (llambda1 / (zz2 * .93 - zz1))*(zz2 - zz2 * .93); % unemployment rate 7%
lla = [llambda1,llambda2];

% Tax system
global mmu ttau
mmu = .15;        % UI replacement rate 15%
ttau = (mmu / zz2) * (lla(2) / lla(1));	     % labor income tax

% Labor supply
global zAvg
zAvg = (lla(1) * z(2) + lla(2) * z(1)) / (lla(1) + lla(2));

%% Approximation Parameters
% Wealth grid
global I amin amax a da aa aaa
I = 100;           % can't be too fine for some reason
amin = 0;
amax = 100;
a = linspace(amin,amax,I)';
da = (amax - amin) / (I - 1);
aa = [a,a];
aaa = reshape(aa,2*I,1);

% Labor productivity
global zz zzz
zz = ones(I,1) * z;
zzz = reshape(zz,2*I,1);

% Idiosyncratic shocks for income
global Aswitch
Aswitch = [-speye(I) * lla(1),speye(I) * lla(1);speye(I) * lla(2),-speye(I) * lla(2)];

% Steady state computations
global rmin rmax r0 maxit crit Delta Ir crit_S
rmin = .0001;       % lower bound for steady state interest rate
rmax = rrho;        % upper bound for steady state interest rate
r0 = .005;          % initial guess for steady state interest rate
maxit = 100;        % maximum iterations on steady state HJB
crit = 1e-6;        % error criterion for steady state value function convergence
Delta = 1e4;        % update size for implicit scheme on steady state HJB
Ir = 100;           % maximum iterations on steady state interest rate
crit_S = 1e-5;      % error criterion for steady state interest rate

% Number of variables in the system
global nVars nEErrors n_v n_g n_p
n_v = 2 * I;
n_g = 2 * I-1 + 1;
n_p = 6;
nVars = n_v + n_g + n_p;
nEErrors = 2 * I;

% Spline parameters
global n_splined
n_knots = 12;
c_power = 7;
x = a';
n_post = 2;	        % This is from two income states
n_prior = 1;
n_splined = n_prior*n_knots*n_post;

author: gregkaplan

repo: phact

save_path: github-repos/MATLAB/gregkaplan-phact


text:
function [V,F,Q] = torus(n,m,r,varargin)
  % TORUS Construct a triangle mesh of a unit torus.
  %
  % [V,F] = torus(n,m,r)
  % [V,F] = torus(n,m,r,'ParameterName',ParameterValue, ...)
  %
  % Inputs:
  %   n  number of vertices around inner ring
  %   m  number of vertices around outer ring
  %   r  radius of the inner ring
  %   Optional:
  %     'R'  followed by outer ring radius {1}
  % Outputs:
  %   V  #V by 3 list of mesh vertex positions
  %   F  #F by 3 list of triangle mesh indices
  %
  % Example:
  %   % Roughly even shaped triangles
  %   n = 40;
  %   r = 0.4;
  %   [V,F] = torus(n,round(r*n),r);
  R = 1;
  params_to_variables = containers.Map( ...
    {'R'},{'R'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  [V,F] = create_regular_grid(n,m,true,true);

  V = V*2*pi;
  th = V(:,2);
  phi = V(:,1);
  V = [cos(phi).*(R+r*cos(th)) sin(phi).*(R+r*cos(th)) r*sin(th)];
  Q = [F(1:2:end-1,[1 2]) F(2:2:end,[2 3])];

end

author: alecjacobson

repo: gptoolbox


text:
function yval = spline_beta_val ( beta1, beta2, ndata, tdata, ydata, tval )

%*****************************************************************************80
%
%% SPLINE_BETA_VAL evaluates a cubic beta spline approximant.
%
%  Discussion:
%
%    The cubic beta spline will approximate the data, but is not
%    designed to interpolate it.
%
%    If BETA1 = 1 and BETA2 = 0, the cubic beta spline will be the
%    same as the cubic B spline approximant.
%
%    With BETA1 = 1 and BETA2 large, the beta spline becomes more like
%    a linear spline.
%
%    In effect, two "phantom" data values are appended to the data,
%    so that the spline will interpolate the first and last data values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real BETA1, the skew or bias parameter.
%    BETA1 = 1 for no skew or bias.
%
%    Input, real BETA2, the tension parameter.
%    BETA2 = 0 for no tension.
%
%    Input, integer NDATA, the number of data values.
%
%    Input, real TDATA(NDATA), the abscissas of the data.
%
%    Input, real YDATA(NDATA), the data values.
%
%    Input, real TVAL, a point at which the spline is to be evaluated.
%
%    Output, real YVAL, the value of the function at TVAL.
%

%
%  Find the nearest interval [ TDATA(LEFT), TDATA(RIGHT) ] to TVAL.
%
  [ left, right ] = r8vec_bracket ( ndata, tdata, tval );
%
%  Evaluate the 5 nonzero beta spline basis functions in the interval,
%  weighted by their corresponding data values.
%
  u = ( tval - tdata(left) ) / ( tdata(right) - tdata(left) );

  delta = ( ( 2.0E+00   ...
    * beta1 + 4.0E+00 ) ...
    * beta1 + 4.0E+00 ) ...
    * beta1 + 2.0E+00 + beta2;

  yval = 0.0E+00;
%
%  Beta function associated with node LEFT - 1, (or "phantom node"),
%  evaluated in its 4th interval.
%
  bval = 2.0E+00 * ( beta1 * ( 1.0E+00 - u ) )^3  / delta;

  if ( 0 < left-1 )
    yval = yval + ydata(left-1) * bval;
  else
    yval = yval + ( 2.0E+00 * ydata(1) - ydata(2) ) * bval;
  end
%
%  Beta function associated with node LEFT,
%  evaluated in its third interval.
%
  a = beta2 + ( 4.0E+00 + 4.0E+00 * beta1 ) * beta1;

  b = - 6.0E+00 * beta1 * ( 1.0E+00 - beta1 ) * ( 1.0E+00 + beta1 );

  c = ( (     - 6.0E+00   ...
      * beta1 - 6.0E+00 ) ...
      * beta1 + 0.0E+00 ) ...
      * beta1 - 3.0E+00 * beta2;

  d = ( (     + 2.0E+00   ...
      * beta1 + 2.0E+00 ) ...
      * beta1 + 2.0E+00 ) ...
      * beta1 + 2.0E+00 * beta2;

  bval = ( a + u * ( b + u * ( c + u * d ) ) ) / delta;

  yval = yval + ydata(left) * bval;
%
%  Beta function associated with node RIGHT,
%  evaluated in its second interval.
%
  a = 2.0E+00;

  b = 6.0E+00 * beta1;

  c = 3.0E+00 * beta2 + 6.0E+00 * beta1 * beta1;

  d = - 2.0E+00 * ( 1.0E+00 + beta2 + beta1 + beta1 * beta1 );

  bval = ( a + u * ( b + u * ( c + u * d ) ) ) / delta;

  yval = yval + ydata(right) * bval;
%
%  Beta function associated with node RIGHT+1, (or "phantom node"),
%  evaluated in its first interval.
%
  bval = 2.0E+00 * u^3 / delta;

  if ( right+1 <= ndata )
    yval = yval + ydata(right+1) * bval;
  else
    yval = yval + ( 2.0E+00 * ydata(ndata) - ydata(ndata-1) ) * bval;
  end

  return
end


text:
function cvx_end

%CVX_END  Completes a cvx specification.
%   CVX_BEGIN marks the end of a new cvx model, and instructs cvx to
%   complete its processing. For standard, complete models, cvx will send
%   a transformed version of the problem to a solver to obtain numeric
%   results, and replace references to cvx variables with numeric values.

try
    evalin( 'caller', 'cvx_verify' );
    evalin( 'caller', 'cvx_finish' );
    evalin( 'caller', 'cvx_pop' );
catch exc
    if ~isequal( exc.identifier, 'CVX:IncompatibleSolver' ),
        evalin( 'caller', 'cvx_pop' );
    end
    rethrow( exc );
end

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

author: yu-jiang

repo: radpbook

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


text:
% The Hurst exponent
%--------------------------------------------------------------------------
% The first 20 lines of code are a small test driver.
% You can delete or comment out this part when you are done validating the
% function to your satisfaction.
%
% Bill Davidson, quellen@yahoo.com
% 13 Nov 2005

function []=hurst_exponent()
disp('testing Hurst calculation');

n=100;
data=rand(1,n);
plot(data);

hurst=estimate_hurst_exponent(data);

[s,err]=sprintf('Hurst exponent = %.2f',hurst);disp(s);

%--------------------------------------------------------------------------
% This function does dispersional analysis on a data series, then does a
% Matlab polyfit to a log-log plot to estimate the Hurst exponent of the
% series.
%
% This algorithm is far faster than a full-blown implementation of Hurst's
% algorithm.  I got the idea from a 2000 PhD dissertation by Hendrik J
% Blok, and I make no guarantees whatsoever about the rigor of this approach
% or the accuracy of results.  Use it at your own risk.
%
% Bill Davidson
% 21 Oct 2003

function [hurst] = estimate_hurst_exponent(data0)   % data set

data=data0;         % make a local copy

[M,npoints]=size(data0);

yvals=zeros(1,npoints);
xvals=zeros(1,npoints);
data2=zeros(1,npoints);

index=0;
binsize=1;

while npoints>4

    y=std(data);
    index=index+1;
    xvals(index)=binsize;
    yvals(index)=binsize*y;

    npoints=fix(npoints/2);
    binsize=binsize*2;
    for ipoints=1:npoints % average adjacent points in pairs
        data2(ipoints)=(data(2*ipoints)+data((2*ipoints)-1))*0.5;
    end
    data=data2(1:npoints);

end % while

xvals=xvals(1:index);
yvals=yvals(1:index);

logx=log(xvals);
logy=log(yvals);

p2=polyfit(logx,logy,1);
hurst=p2(1); % Hurst exponent is the slope of the linear fit of log-log plot

return;


text:
classdef PropertyGetter < Algorithm

    properties (Access = public)
        property; % a string representing the property to access
    end

    methods (Access = public)
        function obj = PropertyGetter(property)
            if nargin > 0
                obj.property = property;
            end

            obj.name = 'getter';
            obj.inputPort = DataType.kAny;
            obj.outputPort = DataType.kAny;
        end

        function computedSignal = compute(obj,data)
            computedSignal = data.(obj.property);
        end

        function str = toString(obj)
            str = sprintf('%s_%s',obj.name,obj.property);
        end

        function editableProperties = getEditableProperties(obj)
            editableProperties = Property('property',obj.property);
        end
    end
end


text:
function pco_errdisp(txt,errorcode)
%if errorcode display error text with errornumber in HEX
%
err=pco_uint32err(errorcode);
if(err)
 disp([txt,' failed with error 0x',num2str(err,'%08X')]);
end
end

%if(errorcode<0)
%disp([txt,' failed with error 0x',num2str(4294967296+errorcode,'%08X')]);
%elseif(errorcode~=0)
% disp([txt,' failed with error 0x',num2str(errorcode,'%08X')]);
%end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
% LFHelperBuild4DFreq - Helper function used to construct 4D frequency-domain filters
%
% Much of the complexity in constructing MD frequency-domain filters, especially around including
% aliased components and controlling the filter rolloff, is common between filter shapes.  This
% function wraps much of this complexity.
%
% This gets called by the LFBuild4DFreq* functions.

% Copyright (c) 2013-2020 Donald G. Dansereau

function [H, FiltOptions] = LFHelperBuild4DFreq( LFSize, BW, FiltOptions, DistFunc )

FiltOptions = LFDefaultField('FiltOptions', 'Precision', 'single');
FiltOptions = LFDefaultField('FiltOptions', 'Rolloff', 'Gaussian'); %Gaussian or Butter
FiltOptions = LFDefaultField('FiltOptions', 'Aspect4D', 1);
FiltOptions = LFDefaultField('FiltOptions', 'Window', false);
FiltOptions = LFDefaultField('FiltOptions', 'Extent4D', 1.0);
FiltOptions = LFDefaultField('FiltOptions', 'IncludeAliased', false);

% avoid rounding error during comparisons
FiltOptions.Extent4D = cast(FiltOptions.Extent4D, FiltOptions.Precision);
FiltOptions.Aspect4D = cast(FiltOptions.Aspect4D, FiltOptions.Precision);

if( length(LFSize) == 1 )
	LFSize = LFSize .* [1,1,1,1];
end
if( length(LFSize) == 5 )
	LFSize = LFSize(1:4);
end
if( length(FiltOptions.Extent4D) == 1 )
	FiltOptions.Extent4D = FiltOptions.Extent4D .* [1,1,1,1];
end
if( length(FiltOptions.Aspect4D) == 1 )
	FiltOptions.Aspect4D = FiltOptions.Aspect4D .* [1,1,1,1];
end
ExtentWithAspect = FiltOptions.Extent4D.*FiltOptions.Aspect4D / 2;

t = LFNormalizedFreqAxis( LFSize(1), FiltOptions.Precision ) .* FiltOptions.Aspect4D(1);
s = LFNormalizedFreqAxis( LFSize(2), FiltOptions.Precision ) .* FiltOptions.Aspect4D(2);
v = LFNormalizedFreqAxis( LFSize(3), FiltOptions.Precision ) .* FiltOptions.Aspect4D(3);
u = LFNormalizedFreqAxis( LFSize(4), FiltOptions.Precision ) .* FiltOptions.Aspect4D(4);
[tt,ss,vv,uu] = ndgrid(cast(t,FiltOptions.Precision),cast(s,FiltOptions.Precision), ...
	cast(v,FiltOptions.Precision),cast(u,FiltOptions.Precision));
P = [tt(:),ss(:),vv(:),uu(:)]';
clear s t u v ss tt uu vv

if( ~FiltOptions.IncludeAliased )
	Tiles = [0,0,0,0];
else
	Tiles = ceil(ExtentWithAspect) % todo[optimization]: optimization possible for large extents
end

if( FiltOptions.IncludeAliased )
	AADist = inf(LFSize, FiltOptions.Precision);
end
WinDist = 0;

% Tiling proceeds only in positive directions, symmetry is enforced afterward, saving some computation overall
for( TTile = 0:Tiles(1) )
	for( STile = 0:Tiles(2) )
		for( VTile = 0:Tiles(3) )
			for( UTile = 0:Tiles(4) )
				Dist = inf(LFSize, FiltOptions.Precision);

				if( FiltOptions.Window )
					ValidIdx = ':';
					WinDist = bsxfun(@minus, abs(P)', ExtentWithAspect);
					WinDist = sum(max(0,WinDist).^2, 2);
					WinDist = reshape(WinDist, LFSize);
				else
					ValidIdx = find(all(bsxfun(@le, abs(P), ExtentWithAspect')));
				end

				Dist(ValidIdx) = DistFunc(P(:,ValidIdx), FiltOptions);
				Dist = Dist + WinDist;

				if( FiltOptions.IncludeAliased )
					AADist(:) = min(AADist(:), Dist(:));
				end

				P(4,:) = P(4,:) + FiltOptions.Aspect4D(4);
			end
			P(4,:) = P(4,:) - (Tiles(4)+1) .* FiltOptions.Aspect4D(4);
			P(3,:) = P(3,:) + FiltOptions.Aspect4D(3);
		end
		P(3,:) = P(3,:) - (Tiles(3)+1) .* FiltOptions.Aspect4D(3);
		P(2,:) = P(2,:) + FiltOptions.Aspect4D(2);
	end
	P(2,:) = P(2,:) - (Tiles(2)+1) .* FiltOptions.Aspect4D(2);
	P(1,:) = P(1,:) + FiltOptions.Aspect4D(1);
end

if( FiltOptions.IncludeAliased )
	Dist = AADist;
	clear AADist
end

H = zeros(LFSize, FiltOptions.Precision);

switch lower(FiltOptions.Rolloff)

	case 'gaussian'
		BW = BW.^2 / log(sqrt(2));
		H(:) = exp( -Dist / BW );  % Gaussian rolloff

	case 'butter'
		FiltOptions = LFDefaultField('FiltOptions', 'Order', 3);
		Dist = sqrt(Dist) ./ BW;
		H(:) = sqrt( 1.0 ./ (1.0 + Dist.^(2*FiltOptions.Order)) ); % Butterworth-like rolloff

	case 'sinc'
		Dist = sqrt(Dist) ./ BW;
		H(:) = sinc(Dist); % sinc-shaped rolloff

	otherwise
		error('unrecognized rolloff method');

end

H = ifftshift(H);

% force symmetric -- needed to visualize in 4D to get this one right
H = max(H, H(mod(LFSize(1):-1:1,LFSize(1))+1, mod(LFSize(2):-1:1,LFSize(2))+1, ...
	mod(LFSize(3):-1:1,LFSize(3))+1, mod(LFSize(4):-1:1,LFSize(4))+1));
H = max(H, H(:, mod(LFSize(2):-1:1,LFSize(2))+1, :, mod(LFSize(4):-1:1,LFSize(4))+1));
H = max(H, H(mod(LFSize(1):-1:1,LFSize(1))+1, :, mod(LFSize(3):-1:1,LFSize(3))+1, :));
% todo: confirm this is sufficient

repo: LFToolbox


text:
function [colors, msh, roiMask] = meshDrawROIs2(msh, rois, nodes, colors, update, edges);
% 'Draw' ROIs on a mesh as colors over different nodes,
% returning the updated mesh colors. Version for mrVista2.
%
% colors = meshDrawROIs2(msh, rois, <oldColors=mesh curvature>, ...
%                          <perim=1>, <update=0>);
%
% INPUTS:
% msh: mesh structure on which to display ROIs.
%
% rois: struct array of ROI objects. NOTE: the .coords field of each ROI
% should be relative to the gray nodes (e.g., I|P|R space or
% VOLUME{1}.coords).
%
% If you provide an ROI with the name 'mask', rather than rendering the
% ROI, it will be used as a mask to show only part of the data on the mesh.
%
% nodes: gray nodes matrix (as in VOLUME{1}.nodes). These can be obtained
% from mrGray .gray graphs (and readGrayGraph) or using mrgGrowGray.
% Apparently the latter fixes an old bug in mapping nodes.
%
% colors: existing mesh colors over which to superimpose the ROIs. If
% omitted, gets from mesh.
%
%
% update: flag to update the mesh colors with the new colors. Default is
% 0, don't update, just return colors.
%
% OUTPUTS:
%   colors: 4 x nVertices set of colors for the mesh, w/ ROI colors added.
%
%   msh: updated mesh structure, including any vertexGrayMaps and
%   connection matrices needed for the computation.
%
%   roiMask: logical mask of size 1 x nVertices indicating where ROIs are.
%
% ras, 07/06.
if isempty(rois), return; end
if ~exist('colors','var') | isempty(colors), colors = meshGet(msh, 'colors'); end
if ~exist('update','var') | isempty(update), update = 0;                      end

% get the ROI mapping mode: may want to make this
% part of mrmPreferences, but maybe this is enough:
prefs = mrmPreferences;
if isequal(prefs.layerMapMode, 'layer1')
    roiMapMode = 'layer1';
else
    roiMapMode = 'any';
end

%%%%%params
nNodes = size(nodes, 2);
nVertices = size(msh.initVertices, 2);
v2g = msh.vertexGrayMap;
if isempty(v2g)
    myErrorDlg('Need Vertex/Gray Map.');
elseif isequal(unique(v2g(:)), 0)
    myErrorDlg('Mesh doesn''t map to any data.');
elseif isequal(roiMapMode, 'any') & size(v2g, 1)==1
	% we only have the mapping to layer 1 -- we need the other layers
	try, vs = msh.mmPerVox; end   % voxel size
	if length(vs) < 3, vs = [1 1 1]; end    % back-compatibility
	msh.vertexGrayMap = mrmMapVerticesToGray(msh.initVertices, nodes, vs, edges);

end

% We need the connection matrix to compute clustering and smoothing, below.
conMat = meshGet(msh,'connectionmatrix');
if isempty(conMat)
    msh = meshSet(msh,'connectionmatrix',1);
    conMat = meshGet(msh,'connectionmatrix');
end

% initialize ROI mask if it's requested as an output:
if nargout>=3
    roiMask = logical(zeros(1, nVertices));
end

%%%%%%Main part:
% Substitute the ROI color for appropriate nearest neighbor vertices
for ii = 1:length(rois)

	%% find indices of those mesh vertices belonging to this ROI:
	roiVertInds = findROIVertices(rois(ii), nodes, v2g, roiMapMode);

	%% adjust the set of ROI vertices to show perimeter / full ROI
	roiVertInds = adjustPerimeter(roiVertInds, msh, prefs, rois(ii).fillMode);

    % update the colors to reflect the ROIs
    if(strcmp(rois(ii).name,'mask'))
        % Don't show the ROI- just use it to mask the data.
        oldColors = mrmGet(msh, 'colors');
        tmp = logical(ones(size(dataMask)));
        tmp(roiVertInds) = 0;
        colors(1:3,tmp) = oldColors(1:3,tmp);
        clear tmp;
    else
        if (ischar(rois(ii).color))
            [stdColorValues, stdColorLabels]  = meshStandardColors;
            jj = findstr(stdColorLabels, rois(ii).color);
            colVal = stdColorValues(:,jj);
        else
            % Check to see if we have a 3x1 vector.
            colVal = rois(ii).color(:);
            if (length(colVal)~=3)
                error('Bad ROI color');
            end
        end
        colors(roiVertInds,1:3) = repmat(colVal'*255, [length(roiVertInds) 1]);

        % build ROI mask if requested
        if nargout>=3
            roiMask(roiVertInds) = true;
        end
    end
end

% round and clip
colors(colors>255) = 255;
colors(colors<0) = 0;
colors = round(colors);

return
% /--------------------------------------------------------------------/ %

% /--------------------------------------------------------------------/ %
function roiVertInds = findROIVertices(roi, nodes, v2g, roiMapMode);
% Returns a list of indices to each vertex belonging to the current ROI.

% find the nodes in the segmentation mapping to this ROI
[x nodeInds] = ismember(round(roi.coords'), nodes([2 1 3],:)', 'rows');
nodeInds = nodeInds(nodeInds>0);

switch lower(roiMapMode)
	case 'layer1'
		% The following will give us *a* vertex index for each gray node.
		% However, the same gray node may map to multiple vertices, and we
		% want them all. (Otherwise, the color overlay will have holes.)
		% So, we loop until we've got them all.
		[nodesToMap roiVertInds] = ismember(nodeInds, v2g(1,:));
		roiVertInds = roiVertInds(roiVertInds>0);
		while(any(nodesToMap))
			v2g(1,roiVertInds) = 0;
			[nodesToMap I] = ismember(nodeInds, v2g(1,:));
			roiVertInds = [roiVertInds; I(I>0)];
		end

	case 'any'
		% if any of the nodes mapping to a given vertex are in the
		% ROI, include that vertex for drawing the ROI
		[I vertInds] = ismember(v2g, nodeInds);
		roiVertInds = find(sum(I)>0); % find columns w/ at least 1 member

	case 'data'
		% take ROI value from the same nodes as the data mapping,
		% rounding up (e.g., for 'mean' data mapping, will behave
		% like 'any'

	otherwise
		error('Invalid ROI Draw Mode preference.')
end

return
% /--------------------------------------------------------------------/ %

% /--------------------------------------------------------------------/ %
function roiVertInds = adjustPerimeter(roiVertInds, msh, prefs, fillMode);
% Sometimes we want to dilate the ROIs a bit. The following does this
% by finding all the neighboring vertices for all the ROI vertices (now
% stored in roiVertInds). The sparse connection matrix comes in handy
% once again, making this a very fast operation.

% We need the connection matrix
conMat = meshGet(msh,'connectionmatrix');
if isempty(conMat)
    msh = meshSet(msh,'connectionmatrix',1);
    conMat = meshGet(msh,'connectionmatrix');
end

%% draw perimeter or fill in the ROI?
if isequal(fillMode, 'perimeter')
	% Draw a perimeter only
	origROI = roiVertInds;

	perimThickness = max(0, prefs.roiDilateIterations);

	if perimThickness > 0
		for t = 1:perimThickness
			neighbors = conMat(:,roiVertInds);
			[roiVertInds cols] = find([neighbors]);
			roiVertInds = unique(roiVertInds);
		end

		% Subtract the original from the dilated version.
		roiVertInds = setdiff(roiVertInds, origROI);

	else    % any way to make the perimeter thinner?
		neighbors = conMat(:,roiVertInds);
		[roiVertInds cols] = find([neighbors]);
		roiVertInds = setdiff(unique(roiVertInds), origROI);

	end

else
	% We can dilate the area if we are not also rendering the perimeter...
	for t = 1:prefs.roiDilateIterations
		neighbors = conMat(:,roiVertInds);
		[roiVertInds cols] = find([neighbors]);
		roiVertInds = unique(roiVertInds);
	end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%%***********************************************************
%% lmiexamp3: generate SDP data for the following LMI problem
%%
%%  max  eta
%%  s.t. [A*P+P*A'  P*B'] +  eta*[0    0]  <= [-G  0]
%%       [B*P       0   ]        [0    I]     [ 0  0]
%%        P >= 0
%%  P and eta are variables, P symmetric.
%%
%%  Ref: Body et al, Linear matrix inequalities in system and
%%       control theory,  p. 10.
%%***********************************************************
%% Here is an example on how to use this function to
%% find an optimal P.
%%
%% A = [-1  0  0; 0 -2  0; 1  1 -1];
%% B = [1  3  5; 2 4 6];
%% G = ones(3,3);
%%
%% [blk,Avec,C,b] = lmiexamp3(A,B,G);
%% [obj,X,y,Z] = sqlp(blk,Avec,C,b);
%% n = size(A,2); N = n*(n+1)/2;
%% blktmp{1,1} = 's'; blktmp{1,2} = n;
%% P = smat(blktmp,y(1:N));
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 20 Apr 02
%%***********************************************************

   function [blk,Avec,C,b] = lmiexamp3(A,B,G);
%%
   [m,n] = size(A);
   [m2,n2] = size(B);
   if (n ~= n2); error('lmiexamp3: A, B not compatible'); end;
%%
   blk{1,1} = 's'; blk{1,2} = m + m2;
   I = speye(n);
   Avec(1,1) = lmifun2(A,I,I,B);
   tmp =  [sparse(m,n+m2);  sparse(m2,n) speye(m2,m2)];
   Avec{1} = [Avec{1} svec(blk,tmp,1)];
%%
   C{1,1} = [-G  sparse(m,m2); sparse(m2,n+m2)];
%%
   N = n*(n+1)/2;
   b = [zeros(N,1); 1];
%%**********************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% make a x-section plus topography...

report_this_filefun(mfilename('fullpath'));

% make a track
lis1 = linspace(lat1,lat2,1000);
lis2 = linspace(lon1,lon2,1000);

tr = [lis2 ; lis1]; tr = tr';
z = [];
% get the topo at each point

for i = 1:length(tr)
    x = find(abs(vlon - tr(i,1)) == min(abs(vlon - tr(i,1))) );
    y = find(abs(vlat - tr(i,2)) == min(abs(vlat - tr(i,2))) );
    z = [z tmap(y,x)  ];
end

figure
axes('pos',[0.1 0.4 0.75 0.3])
plot(xsecx,-xsecy,'or');
hold on

axes('pos',[0.1 0.7 0.75 0.15])
pl = plot(z,'k'); hold on
l = z >= 0; l2 = find(z >= 0);
pl = plot(l2, z(l),'ks');
set(pl,'Markersize',6,'markerfacecolor','k')
l = z < 0; l2 = find(z < 0);
pl = plot(l2, z(l),'bs');
set(pl,'Markersize',6,'markerfacecolor','b')

set(pl,'Linewidth',2)
set(gca,'XTick',[])
grid

axes('pos',[0.1 0.1 0.75 0.3])
pcolor(gx,gy,re3);

author: CelsoReyes


text:
function dr=tide_pl(eu,rp,GMp,pos)

global glc
GME=3.986004415E+14; %earth gravitational constant
H3=0.292; L3=0.015;
dr=zeros(3,1);

r=norm(rp);
if r<=0,return;end

ep=rp/r;

K2=GMp/GME*glc.RE_WGS84^4/r^3;
K3=K2*glc.RE_WGS84/r;
latp=asin(ep(3)); lonp=atan2(ep(2),ep(1));
cosp=cos(latp); sinl=sin(pos(1)); cosl=cos(pos(1));

%step1 in phase (degree2)
p=(3*sinl*sinl-1)/2;
H2=0.6078-0.0006*p;
L2=0.0847+0.0002*p;
a=dot(ep,eu);
dp=K2*3*L2*a;
du=K2*(H2*(1.5*a^2-0.5)-3*L2*a^2);

%step1 in phase (degree3)
dp=dp+K3*L3*(7.5*a^2-1.5);
du=du+K3*(H3*(2.5*a^3-1.5*a)-L3*(7.5*a^2-1.5)*a);

%step1 out-of-phase (only radial)
du=du+3.0/4.0*0.0025*K2*sin(2.0*latp)*sin(2.0*pos(1))*sin(pos(2)-lonp);
du=du+3.0/4.0*0.0022*K2*cosp*cosp*cosl*cosl*sin(2.0*(pos(2)-lonp));

dr(1)=dp*ep(1)+du*eu(1);
dr(2)=dp*ep(2)+du*eu(2);
dr(3)=dp*ep(3)+du*eu(3);

return

repo: GINav


text:
clear;
disp('Program started');
% control mode (0: impedance control
%               1: admittance control(velocity command)
%               2: admittance control(position command))
control_mode = 0;

% reference trajectory mode (0: horizontal line
%                            1: circle)
trajectory_mode = 0;

% dyanmic params(the same as those used in vrep)
I1 = 1;
I2 = 1;
l1 = 0.3;
l2 = 0.3;
m1 = 0.5;
m2 = 0.5;
a1 = 0.15;
a2 = 0.15;
g = 9.81;
% impedance model params
kpx = 150;
kpy = 150;
kdx = 70;
kdy = 70;
mx = 1;
my = 1;

% for moving avarage force filter
windowSize = 5;
b = (1/windowSize)*ones(1,windowSize);
a = 1;
raw_externalforces = zeros(2,windowSize);
filtered_externalforce = zeros(2,1);

vrep=remApi('remoteApi'); % using the prototype file (remoteApiProto.m)
vrep.simxFinish(-1);      % just in case, close all opened connections
id=vrep.simxStart('127.0.0.1',19997,true,true,5000,5); % connect to vrep server
if (id>-1)
    disp('Connected to remote API server');
    % timestep
    dt = 0.01;
    vrep.simxSetFloatingParameter(id,vrep.sim_floatparam_simulation_time_step,dt,vrep.simx_opmode_oneshot_wait);
    % set sychronous mode
    vrep.simxSynchronous(id,true);
    % start the simulation:
    vrep.simxStartSimulation(id,vrep.simx_opmode_oneshot_wait);
    % trigger one step simulation and pause a short time to clear dirty
    % data in vrep buffer
    vrep.simxSynchronousTrigger(id);
%     pause(0.3);

    handles = two_link_workcell_init(vrep,id);
    armJoints = handles.armJoints;
    forcesensor = handles.forcesensor;
    x0 = [0;0;0;0];
    tau = [0;0];

    % set joint control mode in vrep
    if control_mode == 2
        for i=1:2
          vrep.simxSetObjectIntParameter(id,armJoints(i),2001,1,vrep.simx_opmode_oneshot_wait);
        end
    else
        for i=1:2
          vrep.simxSetObjectIntParameter(id,armJoints(i),2001,0,vrep.simx_opmode_oneshot_wait);
        end
    end
    % get init joint states
    for i=1:2
        [res,x0(i)] = vrep.simxGetJointPosition(id,armJoints(i),vrep.simx_opmode_buffer);vrchk(vrep, res, true);
        [res,x0(i+2)] = vrep.simxGetObjectFloatParameter(id,armJoints(i),2012,vrep.simx_opmode_buffer);vrchk(vrep, res, true);
        [res,tau(i)] = vrep.simxGetJointForce(id,armJoints(i),vrep.simx_opmode_buffer);vrchk(vrep, res, true);
    end
    % get external force in force sensor frame
    [res,state,fexternal,tauexternal] = vrep.simxReadForceSensor(id,forcesensor,vrep.simx_opmode_buffer);vrchk(vrep, res, true);
    % transform external force into the robot base frame(here the robot base frame coincides with the world frame in vrep)
    [res,eulerAngles]=vrep.simxGetObjectOrientation(id,forcesensor,-1,vrep.simx_opmode_buffer);vrchk(vrep, res, true);
    forcesensorrotation = rotx(double(eulerAngles(1)*180/pi))*roty(double(eulerAngles(2)*180/pi))*rotz(double(eulerAngles(3)*180/pi));
    fexternal_inertial = forcesensorrotation*[fexternal(1);fexternal(2);fexternal(3)];
    raw_externalforces(1,:) = fexternal_inertial(1)*ones(1,windowSize);
    raw_externalforces(2,:) = fexternal_inertial(3)*ones(1,windowSize);
    % record data
    recordData.x = [];
    recordData.tau = [];
    recordData.fmea = [];
    recordData.u = [];
    recordData.t = [];
    % max simulation duration
    maxsimtime = 20;
    % current time
    current_time = 0;
    % reference joint trajectory
    qref = x0(1:2);
    qdotref = [0;0];
    qdotdotref = [0;0];
    % reference ee trajectory
    eeposref = Direct_Kinematics(x0(1:2),[l1;l2]);
    eevref = [0;0];
    eearef = [0;0];
    % for impedance model
    xe_admittance = [0;0];
    xedot_admittance = [0;0];
    xedotdot_admittance = [0;0];
    q_admittance = [0;0];
    qdot_admittance = [0;0];
    qdotdot_admittance = [0;0];

    while vrep.simxGetConnectionId(id)~=-1
        % get robot joint states
        for i=1:2
            [res,x0(i)] = vrep.simxGetJointPosition(id,armJoints(i),vrep.simx_opmode_buffer);vrchk(vrep, res, true);
            [res,x0(i+2)] = vrep.simxGetObjectFloatParameter(id,armJoints(i),2012,vrep.simx_opmode_buffer);vrchk(vrep, res, true);
            [res,tau(i)] = vrep.simxGetJointForce(id,armJoints(i),vrep.simx_opmode_buffer);vrchk(vrep, res, true);
        end
        % get external force
        [res,state,fexternal,tauexternal] = vrep.simxReadForceSensor(id,forcesensor,vrep.simx_opmode_buffer);vrchk(vrep, res, true);

        % transfor external force form force sensor frame to robot base frame
        [res,eulerAngles]=vrep.simxGetObjectOrientation(id,forcesensor,-1,vrep.simx_opmode_buffer);vrchk(vrep, res, true);
        forcesensorrotation = rotx(double(eulerAngles(1)*180/pi))*roty(double(eulerAngles(2)*180/pi))*rotz(double(eulerAngles(3)*180/pi));
        fexternal_inertial = forcesensorrotation*[fexternal(1);fexternal(2);fexternal(3)];

        % shift raw external forces array
        for i=1:windowSize-1
            raw_externalforces(:,i) = raw_externalforces(:,i+1);
        end
        raw_externalforces(1,windowSize) = fexternal_inertial(1);
        raw_externalforces(2,windowSize) = fexternal_inertial(3);
        % force filtering
        filtered_externalforces_x = filter(b,a,raw_externalforces(1,:));
        filtered_externalforces_y = filter(b,a,raw_externalforces(2,:));
        filtered_externalforce = [filtered_externalforces_x(windowSize);filtered_externalforces_y(windowSize)];

        % feedbacked joint force in v-rep is the force sufferred by joint,
        % so change the sign to get the joint drive torque
        tau = -tau;

        if (current_time > maxsimtime)
            break;
        end

        % get reference trajectory point
        [eeposref,eevref,eearef] = ReferenceTrajectory(current_time,trajectory_mode);

       %% impedance control
        if control_mode == 0
            J =Geometric_Jacobian(x0(1:2),[l1;l2]);
            dJ_dt = Geometric_Jacobian_Derivative(x0,[l1;l2]);
            qdotdotref = J\(eearef - dJ_dt*x0(3:4));
            Mass = Mass_Matrix(x0(1:2),[I1 I2 l1 l2 m1 m2 a1 a2 g]');
            u = Impedance_Controller(x0,[l1;l2],[kpx;kpy;kdx;kdy],eeposref,eevref)...
                + Mass*qdotdotref...
                + Coriolis_Centrifugal_Torque(x0,[l1;l2;m1;m2;a1;a2;g])...
                + Gravity_Torque(x0(1:2),[l1;l2;m1;m2;a1;a2;g]);
            % limit joint tourque
            for i=1:2
                if u(i) > 25
                    u(i) = 25;
                elseif u(i) < -25
                    u(i) = -25;
                end
            end
            % send torque command to vrep
            for i=1:2
                if u(i) < 0
                    res = vrep.simxSetJointTargetVelocity(id,armJoints(i),-99999,vrep.simx_opmode_oneshot);vrchk(vrep, res, true);
                    res = vrep.simxSetJointForce(id,armJoints(i),-u(i),vrep.simx_opmode_oneshot);vrchk(vrep, res, true);
                else
                    res = vrep.simxSetJointTargetVelocity(id,armJoints(i),99999,vrep.simx_opmode_oneshot);vrchk(vrep, res, true);
                    res = vrep.simxSetJointForce(id,armJoints(i),u(i),vrep.simx_opmode_oneshot);vrchk(vrep, res, true);
                end
            end
       %% admittance control
        elseif control_mode == 1
            % --- velocity command --- (qdot_d + qdot_admittance)
            % compute qdotdot corresponds to impedance model
            qdotdot_admittance =  Admittance_Controller(x0,filtered_externalforce,[I1 I2 l1 l2 m1 m2 a1 a2 g]',[mx;my;kpx;kpy;kdx;kdy],eeposref,eevref);
            J =Geometric_Jacobian(x0(1:2),[l1;l2]);
            dJ_dt = Geometric_Jacobian_Derivative(x0,[l1;l2]);
            qdotdot_tracking =  J\(eearef - dJ_dt*x0(3:4));
            qdot_admittance = qdot_admittance + (qdotdot_admittance)*dt;
            qdotref  = J\eevref + qdot_admittance;
            u = qdotref;
            % send velocity command to vrep
            for i=1:2
                vrep.simxSetJointTargetVelocity(id,armJoints(i),qdotref(i),vrep.simx_opmode_oneshot);vrchk(vrep, res, true);
            end
        elseif control_mode == 2
            % --- position command ---  qd + q_admittance
            %compute qdotdot corresponds to impedance model
            qdotdot_admittance =  Admittance_Controller(x0,filtered_externalforce,[I1 I2 l1 l2 m1 m2 a1 a2 g]',[mx;my;kpx;kpy;kdx;kdy],eeposref-eevref*dt,eevref-eearef*dt);%-eevref*dt -eearef*dt
            J =Geometric_Jacobian(x0(1:2),[l1;l2]);
            dJ_dt = Geometric_Jacobian_Derivative(x0,[l1;l2]);
            qdotdot_tracking =  J\(eearef - dJ_dt*x0(3:4));
            q_admittance = q_admittance + qdot_admittance*dt +  1/2*qdotdot_admittance*dt^2;
            qdot_admittance = qdot_admittance + qdotdot_admittance*dt;
            qref = Inverse_Kinematics(x0(1:2),eeposref,[l1 l2]) + q_admittance;
            u = qref;
            % send position command to vrep
            for i=1:2
                vrep.simxSetJointTargetPosition(id,armJoints(i),u(i),vrep.simx_opmode_oneshot);vrchk(vrep, res, true);
            end
        end

        %% record data
        recordData.x = [recordData.x,x0];
        recordData.tau = [recordData.tau,tau];
        recordData.fmea = [recordData.fmea,raw_externalforces(:,end)];
        recordData.u = [recordData.u,u];
        recordData.t = [recordData.t,current_time];

        % update current_time
        current_time = current_time + dt;
        vrep.simxSynchronousTrigger(id);
    end
    % display result
    titles = {'q1','q2','dq1','dq2'};
    figure('name','state');
    for i=1:4
        subplot(2,2,i);
        plot(recordData.t,recordData.x(i,:));
        title(titles{i});
        xlabel('t/s');
    end

    titles = {'u1','u2'};
    figure('name','input');
    for i=1:2
        subplot(1,2,i);
        plot(recordData.t,recordData.u(i,:));
        title(titles{i});
        xlabel('t/s');
    end

    figure('name','external force');
    subplot(1,2,1);
    plot(recordData.t,recordData.fmea(1,:));
    title('fextx');
    xlabel('t/s');
    hold off

    subplot(1,2,2);
    plot(recordData.t,recordData.fmea(2,:));
    title('fexty');
    xlabel('t/s');
    hold off
else
    disp('Failed connecting to remote API server1');
end
vrep.simxStopSimulation(id,vrep.simx_opmode_oneshot_wait);
vrep.simxFinish(id);
vrep.delete();

author: xuhuairuogu

repo: V-REP-Simulation-Projects


text:
function ft_sourcewrite(cfg, source)

% FT_SOURCEWRITE exports source-reconstructed results to gifti or nifti format file.
% The appropriate output file depends on whether the source locations are described by
% on a cortically constrained sheet (gifti) or by a regular 3D lattice (nifti).
%
% Use as
%  ft_sourcewrite(cfg, source)
% where source is a source structure obtained from FT_SOURCEANALYSIS and
% cfg is a structure that should contain
%
%  cfg.filename  = string, filename without the extension
%  cfg.filetype  = string, can be 'nifti', 'gifti' or 'cifti' (default is automatic)
%  cfg.parameter = string, functional parameter to be written to file
%  cfg.precision = string, can be 'single', 'double', etc.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
% If you specify this the input data will be read from a *.mat
% file on disk. This mat file should contain only a single variable,
% corresponding with the input data structure.
%
% See also FT_SOURCEANALYSIS, FT_SOURCEDESCRIPTIVES, FT_VOLUMEWRITE

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar source
ft_preamble provenance source
ft_preamble trackconfig

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% ensure that the required options are present
cfg = ft_checkconfig(cfg, 'required', {'parameter', 'filename'});

% get the options
cfg.parameter      = ft_getopt(cfg, 'parameter');
cfg.filename       = ft_getopt(cfg, 'filename');
cfg.filetype       = ft_getopt(cfg, 'filetype');        % the default is determined further down
cfg.brainstructure = ft_getopt(cfg, 'brainstructure');  % is used for cifti
cfg.parcellation   = ft_getopt(cfg, 'parcellation');    % is used for cifti
cfg.precision      = ft_getopt(cfg, 'precision');       % is used for cifti

% check if the input data is valid for this function
if strcmp(cfg.filetype, 'cifti')

  % keep the transformation matrix
  if isfield(source, 'transform')
    transform = source.transform;
  elseif isfield(source, 'brainordinate') && isfield(source.brainordinate, 'transform')
    transform = source.brainordinate.transform;
  else
    transform = [];
  end

  if isfield(source, 'brainordinate')
    % it is a parcellated source representation, i.e. the main structure one channel for each parcel
    brainordinate = source.brainordinate;
    source = rmfield(source, 'brainordinate');

    % split them and check individually
    source        = ft_checkdata(source, 'datatype', {'timelock', 'freq', 'chan'}, 'feedback', 'yes');
    brainordinate = ft_checkdata(brainordinate, 'datatype', 'parcellation', 'parcellationstyle', 'indexed', 'hasunit', 'yes');

    % merge them again
    source = copyfields(brainordinate, source, setdiff(fieldnames(brainordinate), {'cfg'}));
  else
    source = ft_checkdata(source, 'datatype', 'source', 'hasunit', true, 'feedback', 'yes');
  end

  % keep the transformation matrix
  if ~isempty(transform)
    source.transform = transform;
  end

end % if cifti

if isempty(cfg.filetype)
  if isfield(source, 'dim')
    % source positions are on a regular 3D lattice, save as nifti
    cfg.filetype = 'nifti';
  elseif isfield(source, 'tri')
    % there is a specification of a 2D cortical sheet, save as gifti
    cfg.filetype = 'gifti';
  else
    ft_error('the input data does not look like a 2D sheet, nor as a 3D regular volume');
  end
end

switch (cfg.filetype)
  case 'nifti'
    if numel(cfg.filename)<=4 || ~strcmp(cfg.filename(end-3:end), '.nii');
      cfg.filename = cat(2, cfg.filename, '.nii');
    end

    % convert functional data into 4D
    dat = getsubfield(source, cfg.parameter);
    dat = reshape(dat, source.dim(1), source.dim(2), source.dim(3), []);

    if ~isfield(source, 'transform')
      source.transform = pos2transform(source.pos);
    end

    % this is a bit of a kludge, but ft_write_mri can write 3D and 4D nifti
    ft_write_mri(cfg.filename, dat, 'dataformat', 'nifti', 'transform', source.transform);

  case 'gifti'
    if numel(cfg.filename)<=4 || ~strcmp(cfg.filename(end-3:end), '.gii');
      cfg.filename = cat(2, cfg.filename, '.gii');
    end

    % this is a bit of a kludge, but ft_write_headshape can write gifti including data
    ft_write_headshape(cfg.filename, source, 'data', getsubfield(source, cfg.parameter), 'format', 'gifti');

  case 'cifti'
    % brainstructure should represent the global anatomical structure, such as CortexLeft, Thalamus, etc.
    % parcellation should represent the detailled parcellation, such as BA1, BA2, BA3, etc.
    ft_write_cifti(cfg.filename, source, 'parameter', cfg.parameter, 'brainstructure', cfg.brainstructure, 'parcellation', cfg.parcellation, 'precision', cfg.precision);

  otherwise
    ft_error('unsupported output format "%s"', cfg.filetype);
end % switch filetype

ft_postamble debug
ft_postamble trackconfig
ft_postamble previous source
ft_postamble provenance

author: spm