text:
classdef (Abstract) BaseWidget < matlab.ui.componentcontainer.ComponentContainer & ...
        wt.mixin.ErrorHandling
    % Base class for a graphical widget

    % Copyright 2020-2021 The MathWorks Inc.

    %% Internal properties
    properties (AbortSet, Transient, NonCopyable, ...
            Access = {?wt.abstract.BaseWidget, ?wt.test.BaseWidgetTest} )

        % List of graphics controls that BackgroundColor should apply to
        BackgroundColorableComponents (:,1) matlab.graphics.Graphics

    end %properties

    properties (Transient, NonCopyable, Access = {?wt.test.BaseWidgetTest})

        % Used internally to indicate when setup has finished
        SetupFinished_I (1,1) logical = false

    end %properties

    properties (Transient, NonCopyable, UsedInUpdate, Access = {?wt.test.BaseWidgetTest})

        % Used internally to request update call
        RequestUpdate_I (1,1) logical = false

    end %properties

    %% Debugging Methods
    methods

        function forceUpdate(obj)
           % Forces update to run (For debugging only!)

            obj.update();

        end %function

    end %methods

    %% Setup
    properties ( Transient, NonCopyable, ...
            Access = {?wt.abstract.BaseWidget, ?wt.test.BaseWidgetTest} )

        % The internal grid to manage contents
        Grid matlab.ui.container.GridLayout

    end %properties

    methods (Access = protected)

        function setup(obj)
            % Configure the widget

            % Use CreateFcn to tell when setup finishes
            obj.CreateFcn = @(src,evt)obj.toggleSetupFinished_I();

            % Grid Layout to manage building blocks
            obj.Grid = uigridlayout(obj);
            obj.Grid.ColumnWidth = {'1x'};
            obj.Grid.RowHeight = {'1x'};
            obj.Grid.RowSpacing = 2;
            obj.Grid.ColumnSpacing = 2;
            obj.Grid.Padding = [0 0 0 0];

            % Listen to BackgroundColor changes
            addlistener(obj,'BackgroundColor','PostSet',...
                @(h,e)obj.updateBackgroundColorableComponents());

        end %function

    end %methods

    %% Protected Methods
    methods (Access = protected)

        function updateBackgroundColorableComponents(obj)
            % Update components that are affected by BackgroundColor

            obj.Grid.BackgroundColor = obj.BackgroundColor;
            hasProp = isprop(obj.BackgroundColorableComponents,'BackgroundColor');
            set(obj.BackgroundColorableComponents(hasProp),...
                "BackgroundColor",obj.BackgroundColor);

        end %function

        function requestUpdate(obj)
            % Request update method to run

            if ~obj.SetupFinished_I
                % Do nothing
            elseif verLessThan('matlab','9.10')
                % g228243 Force a call to update
                obj.update();
            else
                % Trigger set of a UsedInUpdate property to request update
                % during next drawnow. (for optimal efficiency)
                obj.RequestUpdate_I = true;
            end

        end %function

    end %methods

    %% Private Methods
    methods (Access = private)

        function toggleSetupFinished_I(obj)
            % Indicate setup is complete

            obj.SetupFinished_I = true;
            obj.CreateFcn = '';

        end %function

    end %methods

    %% Display customization
    methods (Hidden, Access = protected)

        function groups = getPropertyGroups(obj)
            % Customize how the properties are displayed

            % Ignore most superclass properties for default display
            persistent superProps
            if isempty(superProps)
                superProps = properties('matlab.ui.componentcontainer.ComponentContainer');
            end

            % Get the relevant properties
            allProps = properties(obj);
            propNames = setdiff(allProps, superProps, 'stable');

            % Remove props we don't need to see
            propNames( startsWith(propNames, "Font") ) = [];
            propNames( matches(propNames, "Enable") ) = [];

            % Split out the callbacks, fonts
            isCallback = endsWith(propNames, "Fcn");
            isColor = endsWith(propNames, "Color");
            normalProps = propNames(~isCallback & ~isColor);
            callbackProps = propNames(isCallback & ~isColor);

            % Define the groups
            groups = [
                matlab.mixin.util.PropertyGroup(callbackProps)
                matlab.mixin.util.PropertyGroup(normalProps)
                matlab.mixin.util.PropertyGroup(["Position", "Units"])
                ];

            % Ignore empty groups
            groups(~[groups.NumProperties]) = [];

        end %function

    end %methods

    %% Accessors
    methods

        function set.BackgroundColorableComponents(obj,value)
            obj.BackgroundColorableComponents = value;
            obj.updateBackgroundColorableComponents()
        end

    end %methods

end %classdef

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function UNew = fluidPeriodic3D(varargin)
% fluidPeriodic3D: solve fluid registraion in 3D with Periodic
%        boundary conditions
%
%
% author: Nathan D. Cahill
% email: nathan.cahill@rit.edu
% affiliation: Rochester Institute of Technology
% date: January 2014
% licence: GNU GPL v3
%
% Copyright Nathan D. Cahill
% Code available from https://github.com/tomdoel/npReg
%
%

% parse input arguments
[DU,F,mu,lambda,PixSize,M,N,P,RegularizerFactor,HX,HY,HZ] = parse_inputs(varargin{:});

% multiply F by adjoint of Navier-Lame equations
FNew = adjointNL(F/RegularizerFactor,mu,lambda,0,M,N,P);

% compute Fourier transform of new force field
FS = discreteFourierTransform(FNew,M,N,P);

% construct images of coordinates scaled by pi/(N or M or P)
[a,b,c] = ndgrid(pi*(0:(M-1))/(M-1),pi*(0:(N-1))/(N-1),pi*(0:(P-1))/(P-1));

% construct LHS factor
LHSfactor = mu.*(lambda+2*mu).*(2*cos(a) + 2*cos(b) + 2*cos(c) - 6).^2;

% if gamma is zero, set origin term to 1, as DC term does not matter
LHSfactor(1,1,1) = 1;

% solve for FFT of U
VS = cat(4,FS(:,:,:,1)./LHSfactor,FS(:,:,:,2)./LHSfactor,FS(:,:,:,3)./LHSfactor);

% perform inverse FFT
V = discreteFourierTransformInverse(VS,M,N,P);

% now perform Euler integration to construct new displacements
UNew = zeros(M,N,P,3);
UNew(:,:,:,1) = (1 - imfilter(V(:,:,:,1),HX,'replicate','same')).*V(:,:,:,1) - ...
    imfilter(V(:,:,:,2),HY,'replicate','same').*V(:,:,:,2) - ...
    imfilter(V(:,:,:,3),HZ,'replicate','same').*V(:,:,:,3);
UNew(:,:,:,2) = -imfilter(V(:,:,:,1),HY,'replicate','same').*V(:,:,:,1) + ...
    (1 - imfilter(V(:,:,:,2),HY,'replicate','same')).*V(:,:,:,2) - ...
    imfilter(V(:,:,:,3),HY,'replicate','same').*V(:,:,:,3);
UNew(:,:,:,3) = -imfilter(V(:,:,:,1),HZ,'replicate','same').*V(:,:,:,1) - ...
    imfilter(V(:,:,:,2),HZ,'replicate','same').*V(:,:,:,2) + ...
    (1 - imfilter(V(:,:,:,3),HZ,'replicate','same')).*V(:,:,:,3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function FS = discreteFourierTransformInverse(F,M,N,P);
% compute inverse DFT of 3-D vector field

% initialize resulting array
FS = F;

% first perform sine transform down columns
for p=1:P
    for n=1:N
        FS(:,n,p,:) = ifft(FS(:,n,p,:),M,1,'symmetric');
    end
end

% next perform sine transform across rows
for p=1:P
    for m=1:M
        FS(m,:,p,:) = ifft(FS(m,:,p,:),N,2,'symmetric');
    end
end

% finally perform sine transform across pages
for n=1:N
    for m=1:M
        FS(m,n,:,:) = ifft(FS(m,n,:,:),P,3,'symmetric');
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function FS = discreteFourierTransform(F,M,N,P);
% compute DFT of 3-D vector field

% initialize resulting array
FS = complex(F);

% first perform sine transform down columns
for p=1:P
    for n=1:N
        FS(:,n,p,:) = fft(FS(:,n,p,:),M,1);
    end
end

% next perform sine transform across rows
for p=1:P
    for m=1:M
        FS(m,:,p,:) = fft(FS(m,:,p,:),N,2);
    end
end

% finally perform sine transform across pages
for n=1:N
    for m=1:M
        FS(m,n,:,:) = fft(FS(m,n,:,:),P,3);
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function FNew = adjointNL(F,mu,lambda,gamma,M,N,P);
% multiply vector field F by adjoint Navier-Lame equations

% initialize FNew
FNew = zeros(M,N,P,3);

% construct filter that implements 3-D Laplacian
L = (lambda+2*mu)*cat(3,[0 0 0;0 1 0;0 0 0],[0 1 0;1 -6 1;0 1 0],[0 0 0;0 1 0;0 0 0]);

% we will need to use L to form two different filters
% L1 = -(lambda+2*mu)*L; L1(2,2,2) = gamma + L1(2,2,2);
% L2 = -mu*L; L2(2,2,2) = gamma + L2(2,2,2);

% construct grad div filters
GD11 = (lambda+mu)*cat(3,zeros(3,3),[0 1 0;0 -2 0;0 1 0],zeros(3,3));
GD22 = ipermute(GD11,[2 1 3]);
GD33 = ipermute(GD11,[3 2 1]);
GD23 = zeros(3,3,3);
GD23(2,1,1) = 1; GD23(2,3,3) = 1; GD23(2,1,3) = -1; GD23(2,3,1) = -1;
GD23 = GD23*(lambda+mu)/4;
GD12 = ipermute(GD23,[3 1 2]);
GD13 = ipermute(GD23,[2 3 1]);

% perform filtering
FNew(:,:,:,1) = imfilter(F(:,:,:,1),L-GD11,'replicate') + ...
    imfilter(F(:,:,:,2),-GD12,'replicate') + ...
    imfilter(F(:,:,:,3),-GD13,'replicate');
FNew(:,:,:,2) = imfilter(F(:,:,:,1),-GD12,'replicate') + ...
    imfilter(F(:,:,:,2),L-GD22,'replicate') + ...
    imfilter(F(:,:,:,3),-GD23,'replicate');
FNew(:,:,:,3) = imfilter(F(:,:,:,1),-GD13,'replicate') + ...
    imfilter(F(:,:,:,2),-GD23,'replicate') + ...
    imfilter(F(:,:,:,3),L-GD33,'replicate');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [DU,F,mu,lambda,PixSize,M,N,P,RegularizerFactor,HX,HY,HZ] = parse_inputs(varargin);

% get displacement field and check size
F = varargin{2};
PixSize = varargin{4}(1:3);
M = varargin{5};
N = varargin{6};
P = varargin{7};
mu = varargin{8};
lambda = varargin{9};
RegularizerFactor = varargin{10};
DU = varargin{11};
HX = varargin{12};
HY = varargin{13};
HZ = varargin{14};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: tomdoel

repo: pulmonarytoolkit


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Doyen Sahoo
% Contributors: Bin LI, Steven C.H. Hoi
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [] = displayTable( results )
 % Initial display
    cumReturns      = [results.benchmarks results.returns];
    dailyReturns    = [results.benchmarks_daily-1 results.returns_daily];

    % Compute the statsand display the important numbers in table using the
    % library functions
    addpath('../GUI/lib');

    % Get final Values
    [r c] = size(cumReturns);
    results.finalValues = cumReturns(r,:);

    % Get the mean returnfor every day - This is a simple average
    results.meanReturns = mean(dailyReturns);

    % Get annualised returns
    denominator = 252/results.dataFrequency;
    Y = r/denominator;
    results.annualisedReturns = results.finalValues.^(1/Y)-1;

    % Get standard deviation- a measureof risk
    results.standardDeviation = std(dailyReturns);

    % Get annualised standard deviation
    results.annualisedStandardDeviation = results.standardDeviation * sqrt(denominator);

    % Get sharpe ratios
    results.sharpeRatios = sharpe(dailyReturns, results.finalValues,results.dataFrequency);

    % Get Sortino ratios
    results.sortinoRatios = sortino(dailyReturns, 0);

    % Get Value risks at level 5%
    results.valueAtRisks = var5(dailyReturns);

    % Get Maximum draw down
    results.mdds = maxDD_general(cumReturns);

    % Get Calmar ratios
    results.calmars = calmar(cumReturns, results.mdds, results.dataFrequency);

    % Fill up the tables

    tableData   = [results.finalValues; results.meanReturns; results.annualisedReturns; results.standardDeviation; results.annualisedStandardDeviation; results.sharpeRatios; results.calmars; results.sortinoRatios; results.valueAtRisks; results.mdds];

    if ~exist ('OCTAVE_VERSION', 'builtin');
        Market      = tableData(:,1);
        Uniform     = tableData(:,2);
        BestStock   = tableData(:,3);
        BCRP        = tableData(:,4);
        Algorithm   = tableData(:,5);
        report      = {'Final Value','Mean Return for every period','Annualised Return','Standard Deviation','Annualised Standard Deviation','Sharpe Ratio','Calmar Ratio','Sortino Ratio','Value at Risk','Maximum Draw Down'};
        % tableData   = table(Market,Uniform,BestStock,BCRP, Algorithm,
        % 'RowNames', report); % - works only in Matlab 2014a onwards

        finalDisplay = report';
        finalDisplay(2:end+1) = finalDisplay;
        finalDisplay{1} = '';

        finalDisplay{1,2} = 'Market';
        finalDisplay{1,3} = 'Uniform';
        finalDisplay{1,4} = 'BestStock';
        finalDisplay{1,5} = 'BCRP';
        finalDisplay{1,6} = 'Algorithm';

        for i = 1:1:10
            for j = 1:1:5
                finalDisplay{i+1,j+1} = tableData(i,j);
            end
        end

        disp('Performance of the Algorithm compared to baselines based on several metrics');
        disp(finalDisplay);

    else
        tabFinal    = {'Sl. No.', 'Market', 'Uniform', 'BestStock', 'BCRP', 'Algorithm'};
        report      = {'Final Value','Mean Return for every period','Annualised Return','Standard Deviation','Annualised Standard Deviation','Sharpe Ratio','Calmar Ratio','Sortino Ratio','Value at Risk','Maximum Draw Down'};
%         for i = 2:1:11
%             tabFinal{i,1} = report{i-1};
%             for j = 2:1:6
%                 tabFinal{i,j} = tableData(i-1, j-1);
%             end
%         end
        disp('Performance of the Algorithm compared to baselines based on several metrics');

        % Do the display of column titles
        fprintf('\t');
        for j = 1:1:5
            fprintf(char(tabFinal{j}));
            fprintf('\t\t');
        end
        fprintf('\n');
        id = (1:10)';
        tableData = [id tableData];
        disp(tableData);
        disp('Metrics:');
        for i = 1:1:10
            fprintf(num2str(i));
            fprintf('. ');
            fprintf(report{i});
            fprintf('\n');
        end

    rmpath('../GUI/lib');
end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function rtk=udpos(rtk,tt)

global glc
VAR_POS=30^2; VAR_VEL=10^2; VAR_ACC=10^2; var=0;

% initialize position for first epoch
if norm(rtk.x(1:3))<=0
    for i=1:3,rtk=initx(rtk,rtk.sol.pos(i),VAR_POS,i);end
    if rtk.opt.dynamics==1
        for i=1:3,rtk=initx(rtk,rtk.sol.vel(i),VAR_VEL,i+3);end
        for i=1:3,rtk=initx(rtk,1e-6,VAR_ACC,i+6);end
    end
    return;
end

% static mode
if rtk.opt.mode==glc.PMODE_STATIC,return;end

% kinmatic mode without dynamics
if ~rtk.opt.dynamics
    for i=1:3,rtk=initx(rtk,rtk.sol.pos(i),VAR_POS,i);end
    return;
end

for i=1:3
    var=var+rtk.P(i,i);
end
var=var/3;

if var>VAR_POS
    for i=1:3,rtk=initx(rtk,rtk.sol.pos(i),VAR_POS,i);end
    for i=4:6,rtk=initx(rtk,rtk.sol.vel(i-3),VAR_VEL,i);end
    for i=7:9,rtk=initx(rtk,1e-6,VAR_ACC,i);end
    return;
end

nx=0;ix=zeros(rtk.nx,1);
for i=1:rtk.nx
    if rtk.x(i)~=0&&rtk.P(i,i)>0
        ix(nx+1)=i;
        nx=nx+1;
    end
end
ix(nx+1:end)=[];

if nx<9,return;end

% state transition of position/velocity/acceleration
F=eye(nx); x=zeros(nx,1); P=zeros(nx,nx);
for i=1:6,F(i,i+3)=tt;end
for i=1:3,F(i,i+6)=tt^2*0.5;end

for i=1:nx
    x(i)=rtk.x(ix(i));
    for j=1:nx
        P(i,j)=rtk.P(ix(i),ix(j));
    end
end

xp=F*x; P=F*P*F';

for i=1:nx
    rtk.x(ix(i))=xp(i);
    for j=1:nx
        rtk.P(ix(i),ix(j))=P(i,j);
    end
end

% process noise added to only acceleration
Q(1,1)=rtk.opt.prn(4)^2*abs(tt);
Q(2,2)=rtk.opt.prn(4)^2*abs(tt);
Q(3,3)=rtk.opt.prn(5)^2*abs(tt);
[~,Cne]=xyz2blh(x(1:3));
Qv=Cne'*Q*Cne;

for i=1:3
    for j=1:3
        rtk.P(i+6,j+6)=rtk.P(i+6,j+6)+Qv(i,j);
    end
end

return

repo: GINav


text:
%Set parameters: k and threshold of cof
k = 7;
thresholds = 0.9:0.01:1.10;

disp('Demo COF')

data=[(1:10)' rand(10,1)*0.1];
data=[data; 5 1.5;];
data=[data; [(1:1.5:10)' rand(7,1)*0.2 - 10]];

axis([0 10 -5 2]);
scatter(data(:, 1), data(:, 2));

for threshold=thresholds

    clf;
    subplot(121)
    [~, lof] = LOF(data, k);
    target_lof = data(lof>=threshold, :);
    normal_lof = data(lof<threshold, :);
    hold on
    scatter(normal_lof(:, 1), normal_lof(:, 2), 'b');
    scatter(target_lof(:, 1), target_lof(:, 2), 'rx')
    title(sprintf('LOF, threshold=%f', threshold));
    axis([0 10 -15 2]);
    hold off
    subplot(122)
    hold on
    [~, cof] = COF(data, k);
    target_cof = data(cof>=threshold, :);
    normal_cof = data(cof<threshold, :);

    scatter(normal_cof(:, 1), normal_cof(:, 2), 'b');
    scatter(target_cof(:, 1), target_cof(:, 2), 'rx')
    title(sprintf('COF, threshold=%f', threshold));
    axis([0 10 -15 2]);
    hold off
    pause
end

author: dsmi-lab-ntust

repo: AnomalyDetectionToolbox


text:
function h = plot_embedded_graph( vertex,edge,varargin)
i = edge(:,1);
j = edge(:,2);

% Create a long, NaN-separated list of line segments,
% rather than individual segments.

X = [ vertex(i,1) vertex(j,1) repmat(NaN,size(i))]';
Y = [ vertex(i,2) vertex(j,2) repmat(NaN,size(i))]';
Z = [ vertex(i,3) vertex(j,3) repmat(NaN,size(i))]';
X = X(:);
Y = Y(:);
Z = Z(:);

h = plot3(X, Y, Z, varargin{:} );

author: taiya

repo: cloudcontr

save_path: github-repos/MATLAB/taiya-cloudcontr


text:
clear;  clc;
%P2.13
nx = 0:3;       % Index for sequence x(n)
x = 1:4;        % Sequence x(n) = {1,2,3,4}
nh = 0:2;       % Index for impulse h(n)
h = 3:-1:1;     % Sequence h(n) = {3,2,1}
[y,ny] = conv_m(x,nx,h,nh); % Linear Convolution y(n) = h(n)*x(n)
htilde = [h';zeros(size(x)-[0,1])'];
%m = 0:length(htilde)-1;
Hcol = zeros(size(x)+size(h)-[1,1])';
H = htilde;
for k = 1:length(x)-1,
    Hcol=circshift(htilde,k);
    H=[H,Hcol];
end
ytilde = H*x';  %Performs convolution using Toeplitz matrix

author: Sable

repo: mcbench-benchmarks


text:
% Test file for chebtech/minus.m

function pass = test_minus(pref)

% Get preferences.
if ( nargin < 1 )
    pref = chebtech.techPref();
end

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

% A random number to use as an arbitrary additive constant.
alpha = randn() + 1i*randn();

for n = 1:2
    if ( n == 1 )
        testclass = chebtech1();
    else
        testclass = chebtech2();
    end

    %%
    % Check operation in the face of empty arguments.

    f = testclass.make();
    g = testclass.make(@(x) x, [], pref);
    pass(n, 1) = (isempty(f - f) && isempty(f - g) && isempty(g - f));

    %%
    % Check subtraction with scalars.

    f_op = @(x) sin(x);
    f = testclass.make(f_op, [], pref);
    pass(n, 2:3) = test_sub_function_and_scalar(f, f_op, alpha, x);

    %%
    % Check subtraction of two chebtech objects.

    f_op = @(x) zeros(size(x));
    f = testclass.make(f_op, [], pref);
    pass(n, 4:5) = test_sub_function_and_function(f, f_op, f, f_op, x);

    f_op = @(x) exp(x) - 1;
    f = testclass.make(f_op, [], pref);

    g_op = @(x) 1./(1 + x.^2);
    g = testclass.make(g_op, [], pref);
    pass(n, 6:7) = test_sub_function_and_function(f, f_op, g, g_op, x);

    g_op = @(x) cos(1e4*x);
    g = testclass.make(g_op, [], pref);
    pass(n, 8:9) = test_sub_function_and_function(f, f_op, g, g_op, x);

    g_op = @(t) sinh(t*exp(2*pi*1i/6));
    g = testclass.make(g_op, [], pref);
    pass(n, 10:11) = test_sub_function_and_function(f, f_op, g, g_op, x);

    %%
    % Check operation for array-valued chebtech objects.

    f_op = @(x) [zeros(size(x)) zeros(size(x)) zeros(size(x))];
    f = testclass.make(f_op, [], pref);
    pass(n, 12:13) = test_sub_function_and_function(f, f_op, f, f_op, x);

    f_op = @(x) [sin(x) cos(x) exp(x)];
    f = testclass.make(f_op, [], pref);
    pass(n, 14:15) = test_sub_function_and_scalar(f, f_op, alpha, x);

    g_op = @(x) [cosh(x) airy(1i*x) sinh(x)];
    g = testclass.make(g_op, [], pref);
    pass(n, 16:17) = test_sub_function_and_function(f, f_op, g, g_op, x);

    % This should fail with a dimension mismatch error.
    % (Actually, we can do this in R2016a anda above!)
    g_op = @(x) sin(x);
    g = testclass.make(g_op, [], pref);
    try
        h = f - g; %#ok<NASGU>
        pass(n, 18) = false;
        if ( verLessThan('matlab', '9.1') )
            pass(n, 18) = false;
        else
            pass(n, 18) = true;
        end
    catch ME
        if ( verLessThan('matlab', '9.1') )
            pass(n, 18) = strcmp(ME.identifier, 'MATLAB:dimagree');
        else
            pass(n, 18) = false;
        end
    end

    %%
    % Check that direct construction and MINUS give comparable results.

    tol = 10*eps;
    f = testclass.make(@(x) x, [], pref);
    g = testclass.make(@(x) cos(x) - 1, [], pref);
    h1 = f - g;
    h2 = testclass.make(@(x) x - (cos(x) - 1), [], pref);
    h3 = h1 - h2;
    pass(n, 19) = norm(h3.coeffs, inf) < tol;

    %%
    % Check that subtracting a CHEBTECH and an unhappy CHEBTECH gives an
    % unhappy result.

    f = testclass.make(@(x) cos(x+1));    % Happy
    g = testclass.make(@(x) sqrt(x+1));   % Unhappy
    h = f - g;  % Subtract unhappy from happy.
    pass(n, 20) = (~g.ishappy) && (~h.ishappy);
    h = g - f;  % Subtract happy from unhappy.
    pass(n, 21) = (~g.ishappy) && (~h.ishappy);
end

end

% Test the subtraction of a CHEBTECH F, specified by F_OP, to and from a scalar
% ALPHA using a grid of points X in [-1  1] for testing samples.
function result = test_sub_function_and_scalar(f, f_op, alpha, x)
    g1 = f - alpha;
    g2 = alpha - f;
    result(1) = isequal(g1, -g2);
    g_exact = @(x) f_op(x) - alpha;
    result(2) = norm(feval(g1, x) - g_exact(x), inf) <= ...
        10*max(vscale(g1)*eps);
end

% Test the subraction of two CHEBTECH objects F and G, specified by F_OP and
% G_OP, using a grid of points X in [-1  1] for testing samples.
function result = test_sub_function_and_function(f, f_op, g, g_op, x)
    h1 = f - g;
    h2 = g - f;
    result(1) = isequal(h1, -h2);
    h_exact = @(x) f_op(x) - g_op(x);
    norm(feval(h1, x) - h_exact(x), inf);
    result(2) = norm(feval(h1, x) - h_exact(x), inf) <= ...
        1e4*max(vscale(h1)*eps);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% SCRIPT TEST FOR THE 3 DOF spherical manipulator

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE 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.
%
% ARTE 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 Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.

close all

%general call to inverse dynamic function
%TAU = inversedynamic(param, Q, QD, QDD, GRAV, FEXT)

%three different poses
q1 = [0, pi/2, pi/2];
q2 = [0, pi/2, 0];

%load 3dof spherical robot
%the parameters are loaded from robots/example/3dofspherical/parameters.m
spherical=load_robot('example', '3dofspherical');

spherical.graphical.draw_transparent=1;

fprintf('\nTorques at each joint given position = [0 0 0], and velocity=[0 0 0] and standard gravity acting on Z0')
%Please note that the forces and moments are specified with respect to the
%X3 Y3 Z3 reference system.
%In this mechanism, forces acting on Fx3 or Fz3 do not account for tau, nor
%moments acting on My3 or Mx3
tau = inversedynamic(spherical, q1, [0 0 0], [0 0 0], [0 0 9.81]', [0 0 0 0 0 0]')

%draw the robot at this position
drawrobot3d(spherical,q1)
disp('press any key to continue')
pause

%Now, add a force on Fy, and observe the results on tau
tau = inversedynamic(spherical, q1, [0 0 0], [0 0 0], [0 0 9.81]', [0 1 0 0 0 0]')
%draw the robot at this position
drawrobot3d(spherical,q1)
disp('press any key to continue')
pause

%torques necessary to instantaneously bring the arm to the specified state
fprintf('\nTorques at each joint given position = [0 pi/2 0], and velocity=[1 1 1] and acceleration = [1 1 1] and gravity acting on Z0')
tau = inversedynamic(spherical, q2, [1 1 1], [1 1 1], [0 0 9.81]', [0 0 0 0 0 0]')

%draw the robot
drawrobot3d(spherical,q2)
disp('press any key to continue')
pause

%In this case, all torques are caused by gravity acting on the Center Of
%Mass of each link. Note that we have changed the direction of the gravity
%vector g
fprintf('\nTorques at each joint given position = [0 0 0], and velocity=[0 0 0] and acceleration = [0 0 0] and gravity acting on Y0')
tau = inversedynamic(spherical, q1, [0 0 0], [0 0 0], [0 9.81 0]', [0 0 0 0 0 0]')

%draw the robot
drawrobot3d(spherical,q1)

repo: ARTE


text:
classdef AffinityLoss < dagnn.Loss
    properties
        Alpha = 1;
    end

    methods
        function outputs = forward(obj, inputs, params)
            BatchSize = size(inputs{1}, 4);
            outputs{1} = mean(inputs{1}(:) .* (1-inputs{2}(:))) + obj.Alpha * mean(inputs{3}(:) .* inputs{2}(:));
            n = obj.numAveraged;
            m = n + BatchSize;
            obj.average = (n * obj.average + double(gather(outputs{1} * BatchSize))) / m ;
            obj.numAveraged = m;
        end

        function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
            [W,H,C,~] = size(inputs{1});
            NumSamples = W * H * C;
            derInputs{1} = derOutputs{1} * (1-inputs{2}) / NumSamples;
            derInputs{2} = derOutputs{1} * (-inputs{1} + obj.Alpha * inputs{3})/ NumSamples;
            derInputs{3} = derOutputs{1} * obj.Alpha * inputs{2} / NumSamples;
            derParams = {};
        end

        function reset(obj)
            obj.average = 0 ;
            obj.numAveraged = 0 ;
        end
        function obj = AffinityLoss(varargin)
            obj.load(varargin) ;
        end
    end
end

author: KuangJuiHsu


text:
function DEM_demo_EM
% Dual estimation of parameters and hyperparameters; under known causes:
% This demo focuses on conditional parameter estimation with DEM and
% provides a comparative evaluation using EM.  This proceeds by removing
% uncertainly about the input so that the D-step can be discounted.

% get basic convolution model
%==========================================================================
M       = spm_DEM_M('convolution model');

% free parameters
%--------------------------------------------------------------------------
P       = M(1).pE;                            % true parameters
ip      = [2 5];                              % free parameters
pE      = spm_vec(P);
pE(ip)  = 0;
np      = length(pE);
pE      = spm_unvec(pE,P);
pC      = sparse(ip,ip,exp(8),np,np);
M(1).pE = pE;
M(1).pC = pC;

% free hyperparameters
%--------------------------------------------------------------------------
M(1).Q  = {speye(M(1).l,M(1).l)};
M(1).R  = {speye(M(1).n,M(1).n)};

% level 2
%--------------------------------------------------------------------------
M(2).l  = 1;                                  % inputs
M(2).V  = exp(16);                            % very precise causes

% and generate data
%==========================================================================
N       = 32;                                 % length of data sequence
U       = exp(-([1:N] - 12).^2/(2.^2));       % this is the Gaussian cause
DEM     = spm_DEM_generate(M,U,{P},{8,32},{32});

% invert model
%==========================================================================
DEM.U   = U;
DEM     = spm_DEM(DEM);

% overlay true values
%--------------------------------------------------------------------------
spm_DEM_qU(DEM.qU,DEM.pU)

% EM: spm_nlsi_GN
%==========================================================================
G.f   =  inline('P.f*x + P.h*u','x','u','P','M');
G.g   =  inline('P.g*x','x','u','P','M');
G.m   =  DEM.M(1).m;
G.n   =  DEM.M(1).n;
G.l   =  DEM.M(1).l;
G.x   =  DEM.M(1).x;
G.pE  =  DEM.M(1).pE;
G.pC  =  DEM.M(1).pC;
G.hE  = -DEM.M(1).hE;

% exogenous inputs
%--------------------------------------------------------------------------
GU.u  = U';
GU.dt = 1;

% data and serial correlations
%--------------------------------------------------------------------------
t     = ((1:N) - 1);
K     = toeplitz(exp(-t.^2/(2*M(1).E.s^2)));
Q     = K*K';

GY.y  = DEM.Y';
GY.X0 = DEM.X';
GY.dt = 1;
GY.Q  = {kron(speye(G.l,G.l),Q)};

% EM with a Gauss-Newton-like optimization of free energy
%==========================================================================
[Ep,Cp,Eh,F] = spm_nlsi_GN(G,GU,GY);

% parameters
%--------------------------------------------------------------------------
ip    = [2 5];
qP    = spm_vec(DEM.qP.P);
qP    = qP(ip);
tP    = spm_vec(DEM.pP.P);
tP    = tP(ip);
pP    = spm_vec(DEM.M(1).pE);
pP    = pP(ip);
eP    = spm_vec(Ep);
eP    = eP(ip);

spm_figure('GetWin','DEM');
subplot(2,2,4)
bar([tP qP eP])
axis square
legend('true','DEM','EM')
title('parameters','FontSize',16)

cq    = 1.64*sqrt(diag(DEM.qP.C(ip,ip)));
ce    = 1.64*sqrt(diag(Cp(ip,ip)));
hold on
for i = 1:length(qP)
    plot([i i],       qP(i) + [-1 1]*cq(i),'LineWidth',8,'color','r')
    plot([i i] + 1/4, eP(i) + [-1 1]*ce(i),'LineWidth',8,'color','r')
end
hold off

return

% repeat for several realizations
%==========================================================================
clear QP EP QH EH
for i = 1:8

    % generate new data and DEM
    %----------------------------------------------------------------------
    DEM     = spm_DEM_generate(M,U,{P},{8,32},{32});
    DEM.U   = U;
    DEM     = spm_DEM(DEM);

    % EM
    %----------------------------------------------------------------------
    GY.y  = DEM.Y';
    [Ep,Cp,Eh,F] = spm_nlsi_GN(G,GU,GY);

    % retain parameter estimates
    %----------------------------------------------------------------------
    qP      = spm_vec(DEM.qP.P);
    qP      = qP(ip);
    eP      = spm_vec(Ep);
    eP      = eP(ip);

    QP(:,i) = qP;
    EP(:,i) = eP;

    QH(i) = DEM.qH.h{1}(1);
    EH(i) = Eh(1);
end

spm_figure('GetWin','Figure 1');

subplot(2,1,1)
bar(tP,'FaceColor',[1 1 1]*.9,'EdgeColor',[1 1 1]*.9)
hold on
plot([1 2] - 1/8,EP,'r.',[1 2] + 1/4,QP,'k.','Markersize',16)
hold off
axis square
set(gca,'XLim',[0 3])
legend('true','EM','DEM')
title('conditional estimates','FontSize',16)

author: spm


text:
function [Pro,FreeNodec] = nodeisoP2TansferUniform(elemc,elemf,FreeNode)
%  Only works for uniformrefinered.
%
%
% Created by Jie Zhou,Jan,17,2013.
% Talk with Long Chen
%% Data structure.
%% Data structure.
if ~exist('FreeNode','var'), FreeNode = []; end
[elem2dofc,edgec] = dofP2(elemc);
[elem2doff,edgef] = dofP2(elemf);

Nc = max(elemc(:));  Nf = max(elemf(:));
Ndofc = Nc + size(edgec,1);Ndoff = Nf + size(edgef,1);
NTc = size(elemc,1);  %NTc:coarse elements number.
%NTf = size(elemf,1);

elem2node2f = sparse(NTc,9); %very important, map from coarse space to fine space
%we just consider the middle points in fine edges.
elem2node2f(:,1)=elem2doff(1:NTc,4);
elem2node2f(:,2)=elem2doff(1:NTc,5);
elem2node2f(:,3)=elem2doff(1:NTc,6);

elem2node2f(:,4)=elem2doff(NTc+1:2*NTc,4);
elem2node2f(:,5)=elem2doff(NTc+1:2*NTc,5);
elem2node2f(:,6)=elem2doff(NTc+1:2*NTc,6);

elem2node2f(:,7)=elem2doff(2*NTc+1:3*NTc,4);
elem2node2f(:,8)=elem2doff(2*NTc+1:3*NTc,5);
elem2node2f(:,9)=elem2doff(2*NTc+1:3*NTc,6);

%% Assembel the transfer matrix
%        1       2         3        4       5       6         7           8          9
c0=1/2;
locPij=[ 0    1/2*c0    1/2*c0      0        0        0          0         0        0;...
         0    0         0         1/2*c0     0       1/2*c0      0         0        0;...
         0    0         0           0        0        0         1/2*c0     1/2*c0   0;...
         0    0         0         1/2*c0     1/2      0         1/2*c0     0        1/2; ...
         1/2  1/2*c0    0           0        0        0          0         1/2*c0   1/2; ...
         1/2  0         1/2*c0      0        1/2      1/2*c0     0         0        0];
%%Corase grid order          fine grid order
% 3                          *
% * *                        *  7
% *   *                      8    *
% *     *                    *      *
% 5       4                  ***9******
% *         *                *      *    4
% *           *              2  1   5      *
% *             *            *      *        *
% 1*****6********2           ***3********6*****

ii = zeros(40*NTc,1); jj = zeros(40*NTc,1); ss = zeros(40*NTc,1);
index = 0;
for i = 1:6
    for j = 1:9
        if locPij(i,j)~=0
            %Pij = locPij(i,j)*elem2edgeSignc(:,i).*elem2edgesignc2f(:,j);
            ii(index+1:index+NTc) =  elem2dofc(:,i);
            jj(index+1:index+NTc) =  elem2node2f(:,j);
            ss(index+1:index+NTc) =  locPij(i,j);
            index = index + NTc;
        end
    end
end
ii(ii==0)=[];
jj(jj==0)=[];
ss(ss==0)=[];

% Modification for boundary middle points.
s = accumarray([elem2doff(:,4);elem2doff(:,5);elem2doff(:,6)],1,[Ndoff 1]);
bdEdgeidxf = (s == 1);
bdEdgeidxfinjj = bdEdgeidxf(jj);
ss(bdEdgeidxfinjj) = 2*ss(bdEdgeidxfinjj);
ii = [(1:Ndofc)'; ii];
jj = [(1:Ndofc)'; jj];
ss = [ones(Ndofc,1); ss];
Pro     =  sparse(jj,ii,ss,Ndoff,Ndofc);
% Pro(1:Ndofc,1:Ndofc) = speye(Ndofc,Ndofc);
% The transfer for  is a identical matrix, so the local matrix is 6 by 9.

%% Truncated to free nodes
if ~isempty(FreeNode)
    FreeNodec = FreeNode(FreeNode<=Ndofc);
    Pro = Pro(FreeNode,FreeNodec);
end

repo: ifem


text:
classdef CutMeshProvisionalQuadrilater < CutMesh

    properties (Access = private)
       connec
       coord

       cutMesh
       subMesh
       subCutSubMesh

       subMesher

       fullSubCells
       cutSubCells
       levelSetSubMesh
    end

    properties (Access = private)
        lastNode
    end

    methods (Access = public)

        function obj = CutMeshProvisionalQuadrilater(cParams)
            obj.init(cParams);
            obj.lastNode = cParams.lastNode;
        end

        function compute(obj)
            obj.createSubMesher();
            obj.createSubMesh();
            obj.computeLevelSetInSubMesh();
            obj.classifyCells();
            obj.computeSubCutSubMesh();
            obj.computeXcoord();
            obj.computeCoord();
            obj.computeConnec();
            obj.computeCellContainingSubCell();
            obj.computeMesh();
            obj.computeBoundaryMesh();
            obj.computeBoundaryXCoordsIso();
            obj.computeBoundaryCellContainingSubCell();
            obj.computeInnerCutMesh();
            obj.computeBoundaryCutMesh();
        end

    end

    methods (Access = private)

        function createSubMesher(obj)
            s.mesh        = obj.backgroundMesh;
            s.lastNode    = obj.lastNode;
            obj.subMesher = SubMesher(s);
        end

        function createSubMesh(obj)
            obj.subMesh = obj.subMesher.subMesh;
        end

        function computeLevelSetInSubMesh(obj)
            ls = obj.levelSet;
            s.mesh   = obj.backgroundMesh;
            s.fValues = ls;
            f = P1Function(s);
            q = Quadrature.set(obj.backgroundMesh.type);
            q.computeQuadrature('CONSTANT');
            xV = q.posgp;
            lsSubMesh = squeeze(f.evaluate(xV));
            obj.levelSetSubMesh = [ls;lsSubMesh];
        end

       function classifyCells(obj)
            lsInElem = obj.computeLevelSetInElem();
            isFull  = all(lsInElem<0,2);
            isEmpty = all(lsInElem>0,2);
            isCut = ~isFull & ~isEmpty;
            obj.fullSubCells  = find(isFull);
            obj.cutSubCells   = find(isCut);
       end

        function lsElem = computeLevelSetInElem(obj)
            ls = obj.levelSetSubMesh;
            nodes = obj.subMesh.connec;
            nnode = size(nodes,2);
            nElem = size(nodes,1);
            lsElem = zeros(nElem,nnode);
            for inode = 1:nnode
                node = nodes(:,inode);
                lsElem(:,inode) = ls(node);
            end
        end

        function computeSubCutSubMesh(obj)
            s.backgroundMesh = obj.subMesh;
            s.cutCells       = obj.cutSubCells;
            s.levelSet       = obj.levelSetSubMesh;
            cMesh = CutMesh.create(s);
            cMesh.compute();
            obj.subCutSubMesh = cMesh;
        end

        function  computeCellContainingSubCell(obj)
            cellSubMesh = obj.subCutSubMesh.innerCutMesh.cellContainingSubcell;
            fCells      = obj.fullSubCells;
            cellSubMesh = [fCells;cellSubMesh];

            cell = obj.computeSubTriangleOfSubCell();
            obj.cellContainingSubcell = cell(cellSubMesh);
        end

        function cell = computeSubTriangleOfSubCell(obj)
            nnode  = size(obj.backgroundMesh.connec,2);
            cElems = transpose(obj.cutCells);
            cell = repmat(cElems,nnode,1);
            cell = cell(:);
        end

        function globalToLocal = computeGlobalToLocal(obj)
            bConnec = obj.backgroundMesh.connec;
            nnode   = size(bConnec,2);
            nElem   = size(bConnec,1);
            cell = repmat((1:nnode)',1,nElem);
            globalToLocal = cell(:);
        end

        function computeXcoord(obj)
            s.fullCells     = obj.fullSubCells;
            s.cutCells      = obj.subCutSubMesh.innerCutMesh.cellContainingSubcell;
            s.globalToLocal = obj.computeGlobalToLocal();
            s.localMesh     = obj.subMesher.localMesh;
            s.xIsoCutCoord  = obj.subCutSubMesh.innerCutMesh.xCoordsIso;
            xC = XcoordIsoComputer(s);
            obj.xCoordsIso = xC.compute();
        end

        function computeConnec(obj)
            connecCutInterior = obj.subCutSubMesh.innerCutMesh.mesh.connec;
            connecFull        = obj.subMesh.connec(obj.fullSubCells,:);
            obj.connec = [connecFull;connecCutInterior];
        end

        function computeCoord(obj)
            obj.coord  = obj.subCutSubMesh.innerCutMesh.mesh.coord;
        end

        function computeMesh(obj)
            sM.connec = obj.connec;
            sM.coord  = obj.coord;
            sM.kFace  = obj.backgroundMesh.kFace;
            obj.mesh = Mesh(sM);
        end

       function computeBoundaryMesh(obj)
            m = obj.subCutSubMesh.boundaryCutMesh.mesh;
            obj.boundaryMesh = m;
        end

        function computeBoundaryXCoordsIso(obj)
            xCutIso = obj.subCutSubMesh.xCoordsIsoBoundary;
            s.fullCells     = obj.fullSubCells;
            s.cutCells      = obj.cutSubCells;
            s.globalToLocal = obj.computeGlobalToLocal();
            s.localMesh     = obj.subMesher.localMesh;
            s.xIsoCutCoord  = xCutIso;
            xC = XcoordIsoComputer(s);
            xCutG = xC.computeXSubCut();
            obj.xCoordsIsoBoundary = xCutG;
        end

        function cellCont = computeBoundaryCellContainingSubCell(obj)
            cutC = obj.cutSubCells;
            cell = obj.computeSubTriangleOfSubCell();
            cellCont = cell(cutC);
            obj.cellContainingSubCellBoundary = cellCont;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [p, count] = analyze (A, mode, k)				    %#ok
%ANALYZE order and analyze a matrix using CHOLMOD's best-effort ordering.
%
%   Example:
%   [p count] = analyze (A)         orders A, using just tril(A)
%   [p count] = analyze (A,'sym')   orders A, using just tril(A)
%   [p count] = analyze (A,'row')   orders A*A'
%   [p count] = analyze (A,'col')   orders A'*A
%
%   an optional 3rd parameter modifies the ordering strategy:
%
%   [p count] = analyze (A,'sym',k) orders A, using just tril(A)
%   [p count] = analyze (A,'row',k) orders A*A'
%   [p count] = analyze (A,'col',k) orders A'*A
%
%   Returns a permutation and the count of the number of nonzeros in each
%   column of L for the permuted matrix A.  That is, count is returned as:
%
%      count = symbfact2 (A (p,p))       if ordering A
%      count = symbfact2 (A (p,:),'row') if ordering A*A'
%      count = symbfact2 (A (:,p),'col') if ordering A'*A
%
%   CHOLMOD uses the following ordering strategy:
%
%       k = 0:  Try AMD.  If that ordering gives a flop count >= 500 * nnz(L)
%          and a fill-in of nnz(L) >= 5*nnz(C), then try METIS_NodeND (where
%          C = A, A*A', or A'*A is the matrix being ordered.  Selects the best
%          ordering tried.  This is the default.
%
%       if k > 0, then multiple orderings are attempted.
%
%       k = 1 or 2: just try AMD
%       k = 3: also try METIS_NodeND
%       k = 4: also try NESDIS, CHOLMOD's nested dissection (NESDIS), with
%            default parameters.  Uses METIS's node bisector and CCOLAMD.
%       k = 5: also try the natural ordering (p = 1:n)
%       k = 6: also try NESDIS with large leaves of the separator tree
%       k = 7: also try NESDIS with tiny leaves and no CCOLAMD ordering
%       k = 8: also try NESDIS with no dense-node removal
%       k = 9: also try COLAMD if ordering A'*A or A*A', (AMD if ordering A).
%       k > 9 is treated as k = 9
%
%       k = -1: just use AMD
%       k = -2: just use METIS
%       k = -3: just use NESDIS
%
%       The method returning the smallest nnz(L) is used for p and count.
%       k = 4 takes much longer than (say) k = 0, but it can reduce nnz(L) by
%       a typical 5% to 10%.  k = 5 to 9 is getting extreme, but if you have
%       lots of time and want to find the best ordering possible, set k = 9.
%
%   If METIS is not installed for use in CHOLMOD, then the strategy is
%   different:
%
%       k = 1 to 4: just try AMD
%       k = 5 to 8: also try the natural ordering (p = 1:n)
%       k = 9: also try COLAMD if ordering A'*A or A*A', (AMD if ordering A).
%       k > 9 is treated as k = 9
%
%   See also METIS, NESDIS, BISECT, SYMBFACT, AMD

%   Copyright 2006-2007, Timothy A. Davis
%   http://www.cise.ufl.edu/research/sparse

error ('analyze mexFunction not found') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function savgenas() % autogenerated function wrapper
% Matlab script to write output from genas to a file.
% writes two files: one for results for magnitudes and below
% another for magnitudes and above.
%
 % turned into function by Celso G Reyes 2017

ZG=ZmapGlobal.Data; % used by get_zmap_globals

report_this_filefun(mfilename('fullpath'));

figure(mess);
clf ;
set(mess,'Name','Messages');
set(gca,'visible','off');
set(mess,'pos',[ 0.02  0.9 0.3 0.2])
format short

[tbin,zmag,zval] = find(ZABO);           % deal with sparse matrix results
xtz = t0b + (tbin*days(ZG.bin_days));
zmag = minmg+(zmag-1)*magstep;
[~,l] = sort(xtz);                     % sort in time
xtz = xtz(l);
zmag = zmag(l);
zval = zval(l);
tbin = tbin(l);
Z = [tbin'; xtz'; zmag'; zval'];

[newmatfile, newpath] = uiputfile(hodi  , 'Above -Save As'); %Syntax change Matlab Version 7, no window positioning on macs
fid = fopen(fullfile(newpath,newmatfile),'w');
fprintf(fid,'%3.0f %4.2f  %3.2f+  %6.4f\n',Z);

[tbin,zmag,zval] = find(ZBEL);
xtz = t0b + (tbin*days(ZG.bin_days));
zmag = minmg+(zmag-1)*magstep;
[xx,l] = sort(xtz);                     % sort in time
xtz = xtz(l);
zmag = zmag(l);
zval = zval(l);
tbin = tbin(l);
Z = [tbin'; xtz'; zmag'; zval'];

[newmatfile, newpath] = uiputfile(hodi, 'Below -Save As'); %Syntax change Matlab Version 7, no window positioning on macs
fid = fopen(fullfile(newpath,newmatfile),'w');
fprintf(fid,'%3.0f %4.2f  %3.2f-  %6.4f\n',Z);

disp('Output was saved in files \newline \newlinebelow.out and above.out,\newline \newlinePlease rename them if desired. ');

pause(5.0);
zmap_message_center();

end

author: CelsoReyes


text:
function [X,Y,Z]=cylinder2(x,y,z,r,N,o)
%
%   [X,Y,Z]=cylinder2(x,y,z,r,N,o)
%
%   Vectors x, y, z define the central line of cylindrical surface
%   and vector r defines the radius of cylindrical surface, and has
%   the same length as x, y, z.
%
%   N is the number of points around the circumference. The default value
%   is N = 100.
%
%   Matrix o defines ovality of the directrix (circle) of cylindrical
%   surface. Its elements take values between 0 and 1 and apply to the
%   primary and secondary axis of directrix. The defeault values are
%   o=ones(length(x),2).
%
%   Matrices X, Y, Z define the cylindrical surface, surf(X,Y,Z) displays
%   the cylindrical surface.
%
%   Example 1: Horn
%
%              f=linspace(0,4*pi,100);
%              x=cos(f); y=sin(f); z=f; r=-1/99*([1:100]-1)+1;
%              [X,Y,Z]=cylinder2(x,y,z,r);
%
%   Example 2: Pot
%
%              x=[0 0 1 1 .5 .5]; y=[0 0 3 3 1.5 1.5];
%              z=[0 0 2 2 1 1];   r=[0 1 1.5 1.2 1 0];
%              [X,Y,Z]=cylinder2(x,y,z,r,6); view([-150 25])
%
%   Example 3: Plumbing
%
%              x=[0 0 0 0 0 0  1 1 1 1 1 1]; y=[0 .2 .2 1 1 .5 .5 1 1 .2 .2 0];
%              z=[1 1 1 1 0 0 0 0 1 1 1 1]; r=[.2 .2 .1 .1 .1 .1 .1 .1 .1 .1 .2 .2];
%              o=[1 1 .5 .5 .5 1 1 .5 .5 .5 1 1; 1 1 1 1 1 .5 .5 1 1 1 1 1];
%              [X,Y,Z]=cylinder2(x,y,z,r,8,o);
%
%
%   Example 4: Triangle , closed central line
%
%              f=linspace(0,2*pi,4);
%              x=cos(f-pi/6); y=sin(f-pi/6); z=zeros(4,1); r=0.2*ones(4,1);
%              [X,Y,Z]=cylinder3(x,y,z,r); view([0 90])
%
%
%
%   Copyright (c) 2011, Version 2.2
%   Avni Pllana <avniu66@hotmail.com>

X=[]; Y=[]; Z=[];

if nargin<5
    N=100;
end

x=x(:); y=y(:); z=z(:); r=r(:);

nx=length(x);
ny=length(y);
nz=length(z);
nr=length(r);

an=[nx ny nz nr];
if ~ismember(diff(an),[0 0 0],'rows')
    disp(' ')
    disp('x, y, z, r must have the same length!')
else
    if nx<2
        disp(' ')
        disp('The length of x must be greater than 1!')
    else

      if nargin<6
        o=ones(nx,2);
      else
        [o1,o2]=size(o);
        if o2>o1
            o=o';
        end
      end

      C=[x y z];
      D=diff(C);

      for i=1:nx-1
         if norm(D(i,:))>0
            E(i,:)=D(i,:)/norm(D(i,:));
         else
            E(i,:)=[0 0 0];
         end

      end

      if norm(C(1,:)-C(nx,:))<1e-10 && norm(cross(E(1,:),E(nx-1,:)))>0
         Closed=1;
         nxx=nx+1;
         E=[E;E(1,:)];
      else
         Closed=0;
         nxx=nx;
      end

      f=linspace(0,2*pi,N+1);
      for j=1:N+1
          xcc(j)=cos(f(j));
          ycc(j)=sin(f(j));
          zcc(j)=0;
      end

      for i=1:nx

          if i==1

              if ~Closed

                  d=[];
                  k=i;

                  while k<nx-1
                      if norm(cross(E(k,:),E(k+1,:)))==0
                          k=k+1;
                          d=[d k];
                      else
                          break
                      end

                  end

                  if norm(E(1,:))==0;
                      vz=E(2,:);
                  else
                      vz=E(1,:);
                  end

                  if k==nx-1
                     if norm(cross([0 0 1],vz))~=0
                         va=cross([0 0 1],vz);
                     else
                         va=cross([1 1 1],vz);
                     end
                  else
                     va=cross(E(k+1,:),vz);
                  end

                  vy=va/norm(va);
                  vx=cross(vy,vz);
                  V=[vx;vy;vz];
                  Vp=V;
                  Vpp=V;
                  alf=0;

              else

                  d=[];
                  vb=(E(1,:)+E(nx-1,:))/2;
                  vz=vb/norm(vb);
                  va=cross(E(1,:),E(nx-1,:));
                  vy=va/norm(va);
                  vx=cross(vy,vz);
                  V=[vx;vy;vz];
                  Vpp=V;

                  alf=acos(E(1,:)*E(nx-1,:)');

                  if abs(alf>3*pi/4)
                     alf=0;
                  end

              end

          else

             if ismember(i,d)
                 V=Vp;
                 alf=0;
             else
                    d=[];
                    k=i;

                    while k<nx-1
                      if norm(cross(E(k,:),E(k+1,:)))==0
                          k=k+1;
                          d=[d k];
                      else
                          break
                      end

                    end

                    if ~isempty(d) || i==nx
                          if norm(E(i-1,:))==0;
                              vz=E(i-2,:);
                          else
                              if i==nx
                                vz=E(i-1,:);
                              else
                                vz=E(i,:);
                              end
                          end

                          va=Vpp(2,:);
                          vy=va/norm(va);
                          vx=cross(vy,vz);
                          V=[vx;vy;vz];

                          Vp=V;
                          alf=0;

                    end

                    if i<nxx

                          vb=(E(i,:)+E(i-1,:))/2;
                          vz=vb/norm(vb);
                          va=cross(E(i,:),E(i-1,:));
                          vy=va/norm(va);
                          vx=cross(vy,vz);
                          V=[vx;vy;vz];

                          %%%%%%%%%%%%%
                          ang=acos(vy*Vpp(2,:)');
                          Vang=cross(vy,Vpp(2,:));
                          if Vang*vz'>0
                              ang=-ang;
                          end

                          A=[cos(ang) sin(ang); -sin(ang) cos(ang)];
                          xy=A*[xcc;ycc];
                          xcc=xy(1,:);
                          ycc=xy(2,:);
                          %%%%%%%%%%%%
                          Vpp=V;

                          alf=acos(E(i,:)*E(i-1,:)');

                         if abs(alf>3*pi/4)
                            alf=0;
                         end
                     end
             end

          end

          xc=o(i,1)*r(i)/cos(alf/2)*xcc;
          yc=o(i,2)*r(i)*ycc;
          xyz=[xc;yc;zcc];

          M=pinv(V);

          xyz1=M*xyz+repmat([x(i);y(i);z(i)],1,N+1);

          X(:,i)=xyz1(1,:)';
          Y(:,i)=xyz1(2,:)';
          Z(:,i)=xyz1(3,:)';

      end

       figure
       surf(X,Y,Z)
       shading interp
       camlight
       grid on
       axis equal

    end
end


text:
function test_example_ft_realtime_hilbert

% MEM 4gb
% WALLTIME 00:10:00

%
%% Realtime neurofeedback application based on Hilbert phase estimation
%
% This page describes the neurofeedback application that was developed by Marco Rotonda for his MSc thesis **"Controllo volontario della sincronizzazione di fase intraemisferica nella banda gamma eeg mediante neurofeedback"** (Voluntary control of intrahemispheric phase sychronization in the EEG gamma band using neurofeedback).
%
% The thesis (in Italian) which describes all details and the results is available [here](/assets/pdf/example/rt_hilbert/rotonda_thesis_msc.pdf).
%
% The signal analysis details are available [here](/assets/pdf/example/gamma_analisys.pdf).
%
% Neurofeedback, as the word suggest, is feeding back to the subject it's neuronal activity analyzed by it's EEG activity in (almost) real time. In this training we give back to the subjest a visual stimulation (vertical bars) when it's gamma phase synchronization between F3-P3 and F4-P4 increased over the baseline (bars up) or decreased (bars down).
% Here is the feedback
%
%
% If you scroll the script you can see that this image is a particular point of view of a 3D visualization. I you change the POV you could see all the data from the last 3 minutes.
%
%% # Procedure
%
% Get the FieldTripBufferDemo from the workshop_bci2000 folder on the FieldTrip FTP server and . Start BCI2000 and modifiy the config. The script works originally with 40 channels (see the variable 'lab' in the realtime_hilbert script), the samplingrate used is 1000Hz. Use a blocksize of 1000 samples. Ensure BCI2000 runs for more than 2 minutes in the application tab. Set config and start BCI2000.
% Start MATLAB and in the shell type realtime_hilbert.
%
% Now it will start the function realtime_baseline that for 2 minutes will record the subject baseline and store the results for realtime_hilbert.
% Basically there are 12 channels: 4 for EEG (I used F3-P3-F4-P4), 4 for the EOG, 4 for the EMG (front and neck). In this script, those are selected from the forty originally recorded channels you can find in the variable 'lab'.
% The algorithm will care to take away both EOG and EMG artifacts.
%
%% # MATLAB code for ft_realtime_hilbert
%
function ft_realtime_hilbert()

% FT_REALTIME_HILBERT is a neurofeedback application based on Hilbert phase estimation.
%
% Use as
%   ft_realtime_hilbert()
% with the following configuration options that are coded inside the function
%   cfg.channel    = cell-array, see FT_CHANNELSELECTION (default = 'gui')
%   cfg.foilim     = [Flow Fhigh] (default = [0 120])
%   cfg.blocksize  = number, size of the blocks/chuncks that are processed (default = 1 second)
%   cfg.bufferdata = whether to start on the 'first or 'last' data that is available (default = 'first')
%
% The source of the data is configured as
%   cfg.dataset       = string
% or alternatively to obtain more low-level control as
%   cfg.datafile      = string ( es:  cfg.datafile='buffer://localhost:1972';)
%   cfg.headerfile    = string ( es:  cfg.headerfile='buffer://localhost:1972';)
%   cfg.filename      = string ( es:  cfg.filename = 'buffer://localhost:1972';)
%   cfg.eventfile     = string ( es:  cfg.eventfile = 'buffer://localhost:1972';)
%   cfg.dataformat    = string, default is determined automatic
%   cfg.headerformat  = string, default is determined automatic
%   cfg.eventformat   = string, default is determined automatic
%
% To stop the realtime function, you have to press Ctrl-C

% Take off the warning message to avoid problems with ATAN2 and hilbert.
warning off all;

% This is to save subject data
starttime= DATESTR(now, 30);
subjname = input('Insert the subject name.>>>>>', 's');
trialdata=strcat(starttime,subjname);

realtime_baseline(); %This it will launch the baseline script

load means;  % get means from file created by realtime_baseline

cfg = [];
cfg.datafile = 'buffer://localhost:1972';
cfg.headerfile = 'buffer://localhost:1972';
cfg.filename = 'buffer://localhost:1972';

% set the default configuration options
if ~isfield(cfg, 'dataformat'),     cfg.dataformat = [];      end % default is detected automatically
if ~isfield(cfg, 'headerformat'),   cfg.headerformat = [];    end % default is detected automatically
if ~isfield(cfg, 'eventformat'),    cfg.eventformat = [];     end % default is detected automatically
if ~isfield(cfg, 'blocksize'),      cfg.blocksize = 1;        end % in seconds
if ~isfield(cfg, 'overlap'),        cfg.overlap = 0.5;        end % in seconds
if ~isfield(cfg, 'channel'),        cfg.channel = {'all', '-Fp1', '-Fp2',...
                                    '-F7','-Fz','-F8','-FT7', '-FC3', '-FCz',...
                                    '-FC4', '-FT8', '-T3', '-C3', '-Cz', '-C4','-T4','-TP7',...
                                    '-CP3', '-CPz', '-CP4','-TP8','-A1', '-T5', '-Pz',...
                                    '-T6','-A2', '-O1', '-Oz', '-O2'};
end                                                               % This will select F3 F4 P3 P4 for EEG
                                                                  % X1 X2 X3 X4 for EMG
                                                                  % X5 X6 X7 X8 for EOG
if ~isfield(cfg, 'foilim'),         cfg.foilim = [0.1 100];   end
if ~isfield(cfg, 'bufferdata'),     cfg.bufferdata='last';   end

% translate dataset into datafile+headerfile
cfg = ft_checkconfig(cfg, 'dataset2files', 'yes');
cfg = ft_checkconfig(cfg, 'required', {'datafile' 'headerfile'});

% ensure that the persistent variables related to caching are cleared
clear ft_read_header

% start by reading the header from the realtime buffer
hdr = ft_read_header(cfg.headerfile, 'cache', true);

% define a subset of channels for reading
lab=['X1 '; 'X2 '; 'Fp1'; 'Fp2'; 'X3 '; 'X4 '; 'F7 '; 'F3 '; 'Fz '; 'F4 ';...
    'F8 ';'FT7'; 'FC3'; 'FCz'; 'FC4'; 'FT8'; 'T3 '; 'C3 '; 'Cz '; 'C4 ';...
    'T4 ';'TP7'; 'CP3'; 'CPz'; 'CP4'; 'TP8'; 'A1 '; 'T5 '; 'P3 '; 'Pz ';...
    'P4 '; 'T6 ';'A2 '; 'O1 '; 'Oz '; 'O2 '; 'X5 '; 'X6 '; 'X7 '; 'X8 '];
label=cellstr(lab);
cfg.channel = ft_channelselection(cfg.channel, label);
chanindx    = match_str(label, cfg.channel);
nchan       = length(chanindx);
if nchan==0
  error('no channels were selected');
end

% determine the size of blocks to process
blocksize = cfg.blocksize * hdr.Fs;
overlap   = cfg.overlap * hdr.Fs;
prevSample  = 0;
count       = 0;

% Create arrays that contains the rhos found
vectordim = 360;         % vector dimention for the lasts 3 minutes of data
vl=zeros(vectordim,1);  % the visual vector of data for left synchrony
vr=zeros(vectordim,1);  % the visual vector of data for right synchrony
vl1=vl(1);              % the first element of the vector for the visual feedback
vr1=vr(1);              % the first element of the vector for the visual feedback

% This is the ratio above which accept the coherence
meanl = baselinel;           % It has to be calculated from the baseline
meanr = baseliner;           % It has to be calculated from the baseline

% This is the mean above which no EMG feedback has given
meanf = baselinef+(stf*2);             % It has to be calculated from the baseline
meann = baselinen+(stn*2);             % It has to be calculated from the baseline

% This is used to plot the feedback step
fullscreen = get(0,'ScreenSize');
fig1 = figure('NumberTitle','off', ...
    'MenuBar','none', ...
    'Units','pixels', ...
    'Position',[0 0 fullscreen(3) fullscreen(4)]);

% plot the feedback on the second monitor
% set(gcf,'position',[1025,1,1024,768]);

% fig2=figure;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is the general BCI loop where realtime incoming data is handled
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while true

    % Create a matrix that contains the rhos found
    M = [vl vr]; %#ok`<NASGU>`

    % Create a matrix for the visual feedback
    K = [vl1 vr1];

    % determine number of samples available in buffer
    hdr = ft_read_header(cfg.headerfile, 'cache', true);

    % see whether new samples are available
    newsamples = (hdr.nSamples*hdr.nTrials-prevSample);

    if newsamples>=blocksize

        % determine the samples to process
        if strcmp(cfg.bufferdata, 'last') && count==0
            begsample  = hdr.nSamples*hdr.nTrials - blocksize + 1;
            endsample  = hdr.nSamples*hdr.nTrials;
        elseif strcmp(cfg.bufferdata, 'last')
          begsample  = prevSample+1;
          endsample  = prevSample+blocksize ;
        elseif strcmp(cfg.bufferdata, 'first')
          begsample  = prevSample+1;
          endsample  = prevSample+blocksize ;
        else
            error('unsupported value for cfg.bufferdata');
        end

        % this allows overlapping data segments
        if overlap && (begsample>overlap) %#ok`<BDLGI>`
          begsample = begsample - overlap;
          endsample = endsample - overlap;
        end

        % remember up to where the data was read
        prevSample  = endsample;
        count       = count + 1;
        % fprintf('processing segment %d from sample %d to %d\n', count, begsample, endsample);

        % read data segment from buffer
        dat = ft_read_data(cfg.datafile, 'header', hdr, 'begsample', begsample,...
            'endsample', endsample, 'chanindx', chanindx, 'checkboundary', false);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % from here onward it is specific to the hilbert phase sinchronisation from the data %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % put the data in a fieldtrip-like raw structure
        data.trial{1} = dat;
        data.time{1}  = offset2time(begsample, hdr.Fs, endsample-begsample+1);
        data.label    = hdr.label(chanindx);
        data.hdr      = hdr;
        data.fsample  = hdr.Fs;

        % correction of EOG based on algoritm fro
        % Author: German Gomez-Herrero
        %         german.gomezherrero@ieee.org
        %         http://www.cs.tut.fi/~gomezher/index.htm
        %         Institute of Signal Processing
        %         Tampere University of Technology, 2007
        % Reference
        % [1] P. He et al., Med. Biol. Comput. 42 (2004), 407-412
        % [2] S. Haykin. Adaptive Filter Theory, (1996), Prentice Hall

        data.trial{2}(9,:)=data.trial{1}(9,:)-data.trial{1}(10,:);
        data.trial{2}(11,:)=data.trial{1}(11,:)-data.trial{1}(12,:);
        opt.refdata=[data.trial{2}(9,:);data.trial{2}(11,:)];
        [data.trial{3}] = crls_regression(data.trial{1}(5:8,:), opt);

        % Build a FIR filter for EMG correction
        N    = 150;      % Order
        gammaband = [35 45];
        emgband = [60 80];
        emgfnband = [60 499];
        flag = 'scale';  % Sampling Flag
        Beta = 0.9;      % Window Parameter
        win = kaiser(N+1, Beta);

        % Correction between EMG and EEG based on Sheer D.E. "Biofeedback training
        % of 40-Hz eeg and behavior", pp. 325-362, on Behavior and
        % brain electrical activity (1975), Plenum Press. New York

        gammafilter  = fir1(N, gammaband/(hdr.Fs/2), 'bandpass', win, flag);
        datfiltgamma1 = filtfilt(gammafilter,1,data.trial{1}(5,:));
        datfiltgamma2 = filtfilt(gammafilter,1,data.trial{1}(6,:));
        datfiltgamma3 = filtfilt(gammafilter,1,data.trial{1}(7,:));
        datfiltgamma4 = filtfilt(gammafilter,1,data.trial{1}(8,:));
        datfiltgamma = [datfiltgamma1; datfiltgamma2; datfiltgamma3; datfiltgamma4];

        emgfilter = fir1(N, emgband/(hdr.Fs/2), 'bandpass', win, flag);
        datfiltemg1 = filtfilt(emgfilter,1,data.trial{1}(5,:));
        datfiltemg2 = filtfilt(emgfilter,1,data.trial{1}(6,:));
        datfiltemg3 = filtfilt(emgfilter,1,data.trial{1}(7,:));
        datfiltemg4 = filtfilt(emgfilter,1,data.trial{1}(8,:));
        datfiltemg = [datfiltemg1; datfiltemg2; datfiltemg3; datfiltemg4];

        datfiltemgsqr = datfiltemg.^2;
        datfiltgammasqr = datfiltgamma.^2;
        datfiltcrossqr = (datfiltemg.*datfiltgamma).^2;

        correction = datfiltgammasqr-(datfiltcrossqr./datfiltemgsqr);

        data.trial{3}(5:8,:) = datfiltgamma - correction;

        % Find the EMG on forehead and neck
        data.trial{2}(1,:)=data.trial{1}(1,:)-data.trial{1}(2,:);   % Frontal electrods
        data.trial{2}(3,:)=data.trial{1}(3,:)-data.trial{1}(4,:);   % Neck electrods
        emgfnfilter = fir1(N, emgfnband/(hdr.Fs/2), 'bandpass', win, flag);
        datfiltemgf = filter(emgfnfilter,1,data.trial{2}(1,:));
        datfiltemgn = filter(emgfnfilter,1,data.trial{2}(3,:));
        datfiltemgf = abs(datfiltemgf);
        datfiltemgn = abs(datfiltemgn);
        extrMaxValuef = datfiltemgf(find(diff(sign(diff(datfiltemgf)))==-2)+1);
        extrMaxValuen = datfiltemgn(find(diff(sign(diff(datfiltemgn)))==-2)+1);
        extrMaxIndexf =   find(diff(sign(diff(datfiltemgf)))==-2)+1;
        extrMaxIndexn =   find(diff(sign(diff(datfiltemgn)))==-2)+1;
        upf = extrMaxValuef;
        upn = extrMaxValuen;
        upf_t = data.time{1}(extrMaxIndexf);
        upn_t = data.time{1}(extrMaxIndexn);
        upf = interp1(upf_t,upf,data.time{1},'linear');
        upn = interp1(upn_t,upn,data.time{1},'linear');
        emgfmean = nanmean (upf'); %#ok`<UDIM>`
        emgnmean = nanmean (upn'); %#ok`<UDIM>`
        % plot(data.time{1},upf,'r')

%
        % Istantaneous (proto)phase difference found via Hilbert
        % Based on Pikovsky, A. R. (2001). Synchronization. A Universal
        % Concept In Nonlinear Sciences. Cambridge: Cambridge University
        % Press, pag. 368 A2.7

        % crate the data needed for phase coherence index
        chan1=data.trial{3}(5,:); % F3
        chan2=data.trial{3}(6,:); % F4
        chan3=data.trial{3}(7,:); % P3
        chan4=data.trial{3}(8,:); % P4
        chan1h = hilbert(chan1);
        chan2h = hilbert(chan2);
        chan3h = hilbert(chan3);
        chan4h = hilbert(chan4);
        chan1hi = imag(chan1h);
        chan2hi = imag(chan2h);
        chan3hi = imag(chan3h);
        chan4hi = imag(chan4h);

        % find the istantaneous left hemisphere (proto)phase difference
        phil = atan2(((chan1hi .* chan3)-(chan1 .* chan3hi)),...
                     ((chan1 .* chan3)+(chan1hi .* chan3hi)));

        % find the istantaneous right hemisphere (proto)phase difference
        phir = atan2(((chan2hi .* chan4)-(chan2 .* chan4hi)),...
                     ((chan2 .* chan4)+(chan2hi .* chan4hi)));

        % find the right hemisphere synchronization index
        sumsinr = sum(sin(phir))/blocksize;
        sumcosr = sum(cos(phir))/blocksize;
        rhor = sqrt(sumsinr.^2 + sumcosr.^2);

        % find the left hemisphere synchronization index
        sumsinl = sum(sin(phil))/blocksize;
        sumcosl = sum(cos(phil))/blocksize;
        rhol = sqrt(sumsinl.^2 + sumcosl.^2);

        % Give the visual feedback
        if (meanf>emgfmean) && (meann>emgnmean)
            clf;
            visual = bar3(K,0.3);
            view([-90 0]);
            grid off;
            shading interp;
            for i = 1:length(visual)
                zdata = get(visual(i),'Zdata');
                set(visual(i),'Cdata',zdata,'EdgeColor','none')
                colormap hot;
            end
            set(gca,'ZColor',[0.8 0.8 0.8],'Zlim',[-1 1],'YColor',[0.8 0.8 0.8],...
                'Ylim',[0 3],'XColor',[0.8 0.8 0.8],'Xlim',[0.85 1.15],...
                'Color',[0.8 0.8 0.8],'CLim', [-1 1]);

            % Create colorbar
            colorbar([0.5 0.148 0.02 0.73],'ZColor',[0.8 0.8 0.8],'YTick',[],...
                'YColor',[0.8 0.8 0.8],'XColor',[0.8 0.8 0.8]);

            % Create textboxes
            annotation('textbox',[0.496 0.89 0.03 0.04],'String',{'+'},...
                'HorizontalAlignment','center','FontSize',20,'FitBoxToText','off','EdgeColor','none');
            annotation('textbox',[0.50 0.11 0.02 0.04],'String',{'-'},...
                'HorizontalAlignment','center','FontSize',20,'FitBoxToText','off','EdgeColor','none');
            annotation('textbox',[0.25 0.49 0.05 0.07],'String',{'Emisfero','sinistro'},...
                'HorizontalAlignment','center','FontSize',14,'FitBoxToText','off','EdgeColor','none');
            annotation('textbox',[0.735 0.49 0.05 0.07],'String',{'Emisfero','destro'},...
                'HorizontalAlignment','center','FontSize',14,'FitBoxToText','off','EdgeColor','none');

            %create a copy of the plot for the experimenter
            % h1=gcf;
            % h2=figure;
            % objects=allchild(h1);
            % fig2=copyobj(get(h1,'children'),h2);

            % force MATLAB to update the figure
            drawnow ;

            % add the step to the array for the feedback and upgrade the M
            % for the final performance plot

            vl = vl([end 1:end-1]);
            vl(1)=rhol-meanl;
            vl1=vl(1);

            vr = vr([end 1:end-1]);
            vr(1)=rhor-meanr;
            vr1=vr(1);

        elseif  emgfmean>meanf
            clf;
            set(gca,'ZColor',[0.8 0.8 0.8],'Zlim',[0 1],'YColor',[0.8 0.8 0.8],...
                'Ylim',[0 3],'XColor',[0.8 0.8 0.8],'Xlim',[0.85 1.15],...
                'OuterPosition', [-0.0175 0.185 1 0.605],...
                'Color',[0.8 0.8 0.8],'CLim', [-50 50]);
            annotation(fig1,'textbox',[0.20 0.35 0.6292 0.1929],...
                'String',{'Rilassa i muscoli della fronte'},...
                'HorizontalAlignment','center','FontSize',20,'FitBoxToText','off','EdgeColor','none',...
                'Color',[1 0 0]);
            drawnow ;

        elseif  emgnmean>meann
            clf;
            set(gca,'ZColor',[0.8 0.8 0.8],'Zlim',[0 1],'YColor',[0.8 0.8 0.8],...
                'Ylim',[0 3],'XColor',[0.8 0.8 0.8],'Xlim',[0.85 1.15],...
                'OuterPosition', [-0.0175 0.185 1 0.605],...
                'Color',[0.8 0.8 0.8],'CLim', [-1 1]);
            annotation(fig1,'textbox',[0.20 0.35 0.6292 0.1929],...
                'String',{'Rilassa i muscoli del collo'},...
                'HorizontalAlignment','center','FontSize',20,'FitBoxToText','off','EdgeColor','none',...
                'Color',[1 0 0]);
            drawnow ;

        end % end feedbacks

    end % if enough new samples

    if count > vectordim
        cd 'C:\TEST\hilbert\risultati';
        save (trialdata);
        save hilbert.mat;
        close all hidden;
        cd 'C:\TEST';
        break;
    end

end % while true

%% # MATLAB code for ft_realtime_baseline
%
function ft_realtime_baseline()

% FT_REALTIME_BASELINE is a neurofeedback application based on Hilbert phase estimation.
%
% Use as
%   ft_realtime_baseline()
% with the following configuration options that are coded inside the function
%   cfg.channel    = cell-array, see FT_CHANNELSELECTION (default = 'gui')
%   cfg.foilim     = [Flow Fhigh] (default = [0 120])
%   cfg.blocksize  = number, size of the blocks/chuncks that are processed (default = 1 second)
%   cfg.bufferdata = whether to start on the 'first or 'last' data that is available (default = 'first')
%
% The source of the data is configured as
%   cfg.dataset       = string
% or alternatively to obtain more low-level control as
%   cfg.datafile      = string ( es:  cfg.datafile='buffer://localhost:1972';)
%   cfg.headerfile    = string ( es:  cfg.headerfile='buffer://localhost:1972';)
%   cfg.filename      = string ( es:  cfg.filename = 'buffer://localhost:1972';)
%   cfg.eventfile     = string ( es:  cfg.eventfile = 'buffer://localhost:1972';)
%   cfg.dataformat    = string, default is determined automatic
%   cfg.headerformat  = string, default is determined automatic
%   cfg.eventformat   = string, default is determined automatic
%
% To stop the realtime function, you have to press Ctrl-C

%Take off the warning message to avoid problems with ATAN2 and hilbert.
warning off all;

cfg = [];
cfg.datafile = 'buffer://localhost:1972';
cfg.headerfile = 'buffer://localhost:1972';
cfg.filename = 'buffer://localhost:1972';

% set the default configuration options
if ~isfield(cfg, 'dataformat'),     cfg.dataformat = [];      end % default is detected automatically
if ~isfield(cfg, 'headerformat'),   cfg.headerformat = [];    end % default is detected automatically
if ~isfield(cfg, 'eventformat'),    cfg.eventformat = [];     end % default is detected automatically
if ~isfield(cfg, 'blocksize'),      cfg.blocksize = 1;        end % in seconds
if ~isfield(cfg, 'overlap'),        cfg.overlap = 0.5;        end % in seconds
if ~isfield(cfg, 'channel'),        cfg.channel = {'all', '-Fp1', '-Fp2',...
                                    '-F7','-Fz','-F8','-FT7', '-FC3', '-FCz',...
                                    '-FC4', '-FT8', '-T3', '-C3', '-Cz', '-C4','-T4','-TP7',...
                                    '-CP3', '-CPz', '-CP4','-TP8','-A1', '-T5', '-Pz',...
                                    '-T6','-A2', '-O1', '-Oz', '-O2'};
end                                                               % This will select F3 F4 P3 P4 for EEG
                                                                  % X1 X2 X3 X4 for EMG
                                                                  % X5 X6 X7 X8 for EOG
if ~isfield(cfg, 'foilim'),         cfg.foilim = [0.1 100];   end
if ~isfield(cfg, 'bufferdata'),     cfg.bufferdata='last';   end

% translate dataset into datafile+headerfile
cfg = checkconfig(cfg, 'dataset2files', 'yes');
cfg = checkconfig(cfg, 'required', {'datafile' 'headerfile'});

% ensure that the persistent variables related to caching are cleared
clear ft_read_header

% start by reading the header from the realtime buffer
hdr = ft_read_header(cfg.headerfile, 'cache', true);

% define a subset of channels for reading
lab=['X1 '; 'X2 '; 'Fp1'; 'Fp2'; 'X3 '; 'X4 '; 'F7 '; 'F3 '; 'Fz '; 'F4 ';...
    'F8 ';'FT7'; 'FC3'; 'FCz'; 'FC4'; 'FT8'; 'T3 '; 'C3 '; 'Cz '; 'C4 ';...
    'T4 ';'TP7'; 'CP3'; 'CPz'; 'CP4'; 'TP8'; 'A1 '; 'T5 '; 'P3 '; 'Pz ';...
    'P4 '; 'T6 ';'A2 '; 'O1 '; 'Oz '; 'O2 '; 'X5 '; 'X6 '; 'X7 '; 'X8 '];
label=cellstr(lab);
cfg.channel = ft_channelselection(cfg.channel, label);
chanindx    = match_str(label, cfg.channel);
nchan       = length(chanindx);
if nchan==0
  error('no channels were selected');
end

% determine the size of blocks to process
blocksize = cfg.blocksize * hdr.Fs;
overlap   = cfg.overlap * hdr.Fs;
prevSample  = 0;
count       = 0;

% Create arrays that contains the rhos found
nummaxarray = 240; % the step is 500ms so 240 are 2 minutes
i=1;
vlmean=zeros(1,nummaxarray);
vrmean=zeros(1,nummaxarray);
frontmeans=zeros(1,nummaxarray);
neckmeans=zeros(1,nummaxarray);

% This is used to plot the screen for the subject
fullscreen = get(0,'ScreenSize');
fig1 = figure('NumberTitle','off', ...
    'MenuBar','none', ...
    'Units','pixels', ...
    'Position',[0 0 fullscreen(3) fullscreen(4)]);

% plot the feedback on the second monitor
% set(gcf,'position',[1025,1,1024,768]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is the general BCI loop where realtime incoming data is handled
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while true

    % determine number of samples available in buffer
    hdr = ft_read_header(cfg.headerfile, 'cache', true);

    % see whether new samples are available
    newsamples = (hdr.nSamples*hdr.nTrials-prevSample);

    if newsamples>=blocksize

        % determine the samples to process
        if strcmp(cfg.bufferdata, 'last') && count==0
            begsample  = hdr.nSamples*hdr.nTrials - blocksize + 1;
            endsample  = hdr.nSamples*hdr.nTrials;
        elseif strcmp(cfg.bufferdata, 'last')
          begsample  = prevSample+1;
          endsample  = prevSample+blocksize ;
        elseif strcmp(cfg.bufferdata, 'first')
          begsample  = prevSample+1;
          endsample  = prevSample+blocksize ;
        else
            error('unsupported value for cfg.bufferdata');
        end

        % this allows overlapping data segments
        if overlap && (begsample>overlap) %#ok`<BDLGI>`
          begsample = begsample - overlap;
          endsample = endsample - overlap;
        end

        % remember up to where the data was read
        prevSample  = endsample;
        count       = count + 1;
        % fprintf('processing segment %d from sample %d to %d\n', count, begsample, endsample);

        % read data segment from buffer
        dat = ft_read_data(cfg.datafile, 'header', hdr, 'begsample', begsample, 'endsample', endsample, 'chanindx', chanindx, 'checkboundary', false);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % from here onward it is specific to the hilbert phase sinchronisation from the data %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % put the data in a fieldtrip-like raw structure
        data.trial{1} = dat;
        data.time{1}  = offset2time(begsample, hdr.Fs, endsample-begsample+1);
        data.label    = hdr.label(chanindx);
        data.hdr      = hdr;
        data.fsample  = hdr.Fs;

        % correction of EOG based on algoritm fro
        % Author: German Gomez-Herrero
        %         german.gomezherrero@ieee.org
        %         http://www.cs.tut.fi/~gomezher/index.htm
        %         Institute of Signal Processing
        %         Tampere University of Technology, 2007
        % Reference
        % [1] P. He et al., Med. Biol. Comput. 42 (2004), 407-412
        % [2] S. Haykin. Adaptive Filter Theory, (1996), Prentice Hall

        data.trial{2}(9,:)=data.trial{1}(9,:)-data.trial{1}(10,:);
        data.trial{2}(11,:)=data.trial{1}(11,:)-data.trial{1}(12,:);
        opt.refdata=[data.trial{2}(9,:);data.trial{2}(11,:)];
        [data.trial{3}] = crls_regression(data.trial{1}(5:8,:), opt);

        % Build a FIR filter for EMG correction
        N    = 150;      % Order
        gammaband = [35 45];
        emgband = [60 80];
        emgfnband = [60 499];
        flag = 'scale';  % Sampling Flag
        Beta = 0.9;      % Window Parameter
        win = kaiser(N+1, Beta);

        % Correction between EMG and EEG based on Sheer D.E. "Biofeedback training
        % of 40-Hz eeg and behavior", pp. 325-362, on Behavior and
        % brain electrical activity (1975), Plenum Press. New York

        gammafilter  = fir1(N, gammaband/(hdr.Fs/2), 'bandpass', win, flag);
        datfiltgamma1 = filtfilt(gammafilter,1,data.trial{1}(5,:));
        datfiltgamma2 = filtfilt(gammafilter,1,data.trial{1}(6,:));
        datfiltgamma3 = filtfilt(gammafilter,1,data.trial{1}(7,:));
        datfiltgamma4 = filtfilt(gammafilter,1,data.trial{1}(8,:));
        datfiltgamma = [datfiltgamma1; datfiltgamma2; datfiltgamma3; datfiltgamma4];

        emgfilter = fir1(N, emgband/(hdr.Fs/2), 'bandpass', win, flag);
        datfiltemg1 = filtfilt(emgfilter,1,data.trial{1}(5,:));
        datfiltemg2 = filtfilt(emgfilter,1,data.trial{1}(6,:));
        datfiltemg3 = filtfilt(emgfilter,1,data.trial{1}(7,:));
        datfiltemg4 = filtfilt(emgfilter,1,data.trial{1}(8,:));
        datfiltemg = [datfiltemg1; datfiltemg2; datfiltemg3; datfiltemg4];

        datfiltemgsqr = datfiltemg.^2;
        datfiltgammasqr = datfiltgamma.^2;
        datfiltcrossqr = (datfiltemg.*datfiltgamma).^2;

        correction = datfiltgammasqr-(datfiltcrossqr./datfiltemgsqr);

        data.trial{3}(5:8,:) = datfiltgamma - correction;

        % Find the EMG on forehead and neck
        data.trial{2}(1,:)=data.trial{1}(1,:)-data.trial{1}(2,:);   % Frontal electrods
        data.trial{2}(3,:)=data.trial{1}(3,:)-data.trial{1}(4,:);   % Neck electrods
        emgfnfilter = fir1(N, emgfnband/(hdr.Fs/2), 'bandpass', win, flag);
        datfiltemgf = filter(emgfnfilter,1,data.trial{2}(1,:));
        datfiltemgn = filter(emgfnfilter,1,data.trial{2}(3,:));
        datfiltemgf = abs(datfiltemgf);
        datfiltemgn = abs(datfiltemgn);
        extrMaxValuef = datfiltemgf(find(diff(sign(diff(datfiltemgf)))==-2)+1);
        extrMaxValuen = datfiltemgn(find(diff(sign(diff(datfiltemgn)))==-2)+1);
        extrMaxIndexf =   find(diff(sign(diff(datfiltemgf)))==-2)+1;
        extrMaxIndexn =   find(diff(sign(diff(datfiltemgn)))==-2)+1;
        upf = extrMaxValuef;
        upn = extrMaxValuen;
        upf_t = data.time{1}(extrMaxIndexf);
        upn_t = data.time{1}(extrMaxIndexn);
        upf = interp1(upf_t,upf,data.time{1},'linear');
        upn = interp1(upn_t,upn,data.time{1},'linear');
        emgfmean = nanmean (upf'); %#ok`<UDIM>`
        emgnmean = nanmean (upn'); %#ok`<UDIM>`
        % plot(data.time{1},upf,'r')

%
        % Istantaneous (proto)phase difference found via Hilbert
        % Based on Pikovsky, A. R. (2001). Synchronization. A Universal
        % Concept In Nonlinear Sciences. Cambridge: Cambridge University
        % Press, pag. 368 A2.7

        % crate the data needed for phase coherence index
        chan1=data.trial{3}(5,:); % F3
        chan2=data.trial{3}(6,:); % F4
        chan3=data.trial{3}(7,:); % P3
        chan4=data.trial{3}(8,:); % P4
        chan1h = hilbert(chan1);
        chan2h = hilbert(chan2);
        chan3h = hilbert(chan3);
        chan4h = hilbert(chan4);
        chan1hi = imag(chan1h);
        chan2hi = imag(chan2h);
        chan3hi = imag(chan3h);
        chan4hi = imag(chan4h);

        % find the istantaneous left hemisphere (proto)phase difference
        phil = atan2(((chan1hi .* chan3)-(chan1 .* chan3hi)),...
                     ((chan1 .* chan3)+(chan1hi .* chan3hi)));

        % find the istantaneous right hemisphere (proto)phase difference
        phir = atan2(((chan2hi .* chan4)-(chan2 .* chan4hi)),...
                     ((chan2 .* chan4)+(chan2hi .* chan4hi)));

        % find the right hemisphere synchronization index
        sumsinr = sum(sin(phir))/blocksize;
        sumcosr = sum(cos(phir))/blocksize;
        rhor = sqrt(sumsinr.^2 + sumcosr.^2);

        % find the left hemisphere synchronization index
        sumsinl = sum(sin(phil))/blocksize;
        sumcosl = sum(cos(phil))/blocksize;
        rhol = sqrt(sumsinl.^2 + sumcosl.^2);

        % update the array for the mean
        vlmean(i)=rhol;
        vrmean(i)=rhor;
        frontmeans(i)=emgfmean;
        neckmeans(i)=emgnmean;

        i=i+1;

    end % if enough new samples

    % screen for the subject
    clf;
    set(gca,'ZColor',[0.8 0.8 0.8],'Zlim',[-50 50],'YColor',[0.8 0.8 0.8],...
        'Ylim',[0 3],'XColor',[0.8 0.8 0.8],'Xlim',[0.85 1.15],...
        'DataAspectRatio',[0.2 0.1 3],'OuterPosition', [-0.0175 0.185 1 0.605],...
        'Color',[0.8 0.8 0.8],'CLim', [-50 50]);
    annotation(fig1,'textbox',[0.18 0.35 0.6292 0.1929],...
        'String',{'Rilassati mantenendo gli occhi aperti'},...
        'HorizontalAlignment','center','FontSize',20,'FitBoxToText','off','EdgeColor','none',...
        'Color',[1 0 0]);
    drawnow ;

    if count > (nummaxarray-1)
        close all hidden;
        break;
    end

end % while true

% print the mean baseline
baselinel = mean (vlmean);
fprintf ('mean baseline left= %d \n', baselinel);
baseliner = mean (vrmean);
fprintf ('mean baseline right= %d \n', baseliner);
baselinef = mean (frontmeans);
fprintf ('mean baseline forehead= %d \n', baselinef);
baselinen = mean (neckmeans);
fprintf ('mean baseline neck= %d \n', baselinen);
stf = std(frontmeans);
fprintf ('Std baseline front= %d \n', stf);
stn = std(neckmeans);
fprintf ('Std baseline neck= %d \n', stn);

save means.mat;   % save to the working folder - (edit jonaweber 2010)


text:
function zeta_values_test ( )

%*****************************************************************************80
%
%% ZETA_VALUES_TEST demonstrates the use of ZETA_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ZETA_VALUES_TEST:\n' );
  fprintf ( 1, '  ZETA_VALUES returns values of \n' );
  fprintf ( 1, '  the Riemann Zeta function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N        Zeta(N)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, zeta ] = zeta_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '%6d      %24.16f\n', n, zeta );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ vX, mX ] = SolveLsL2Prox( mA, vB, paramLambda, numIterations )
% ----------------------------------------------------------------------------------------------- %
%[ vX, mX ] = SolveLsL1Prox( mA, vB, lambdaFctr, numIterations )
% Solve L2 Regularized Least Squares Using Proximal Gradient (PGM) Method.
% Input:
%   - mA                -   Input Matirx.
%                           The model matrix.
%                           Structure: Matrix (m X n).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - vB                -   input Vector.
%                           The model known data.
%                           Structure: Vector (m X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - paramLambda       -   Parameter Lambda.
%                           The L1 Regularization parameter.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: (0, inf).
%   - numIterations     -   Number of Iterations.
%                           Number of iterations of the algorithm.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range {1, 2, ...}.
% Output:
%   - vX                -   Output Vector.
%                           Structure: Vector (n X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References
%   1.  Wikipedia PGM - https://en.wikipedia.org/wiki/Proximal_gradient_method.
% Remarks:
%   1.  Using vanilla PGM.
% Known Issues:
%   1.  A
% TODO:
%   1.  B
% Release Notes:
%   -   1.0.000     23/08/2017
%       *   First realease version.
% ----------------------------------------------------------------------------------------------- %

mAA = mA.' * mA;
vAb = mA.' * vB;
vX  = pinv(mA) * vB; %<! Dealing with "Fat Matrix"

stepSize = 1 / (2 * (norm(mA, 2) ^ 2));
% stepSize = 1 / sum(mA(:) .^ 2); %<! Faster to calculate, conservative (Hence slower)

mX = zeros([size(vX, 1), numIterations]);
mX(:, 1) = vX;

for ii = 2:numIterations

    vG = (mAA * vX) - vAb;
    vX = ProxL2(vX - (stepSize * vG), stepSize * paramLambda);

    mX(:, ii) = vX;

end

end

function [ vX ] = ProxL2( vX, lambdaFactor )

% Soft Thresholding L2
% See https://math.stackexchange.com/a/2094088/33
vX = (1 - (lambdaFactor / max(norm(vX, 2), lambdaFactor))) * vX;

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function varargout = spectrumLabel(hAxesTarget)
%spectrumLabel   Add a spectrum to the X-Axis of a plot.
%
%    H = SPECTRUMLABEL(HTARGET) adds a spectrum colorbar and wavelength
%    labels to the X axis of the axes object HTARGET.  For plots whose X
%    values are in the range 300 to 830, the colorbar will show the ROYGBIV
%    spectrum based on the 1964 CIE 10-degree observer.  Values outside the
%    range will be shown as black.  H contains a handle to the axes object
%    containing the spectrum.
%
%    Examples:
%    ---------
%    % Show the D65 illuminant.
%    [lambda, D65] = illuminant('d65');
%    figure;
%    plot(lambda, D65)
%    title('D_{65} illuminant')
%    spectrumLabel(gca)
%
%    % Show the CIE 1964 10-degree observer color matching functions.
%    [lambda, x, y, z] = colorMatchFcn('1964_full');
%    figure;
%    plot(lambda, x, 'r', lambda, y, 'g', lambda, z, 'b')
%    title('CIE 1964 10-degree observer')
%    spectrumLabel(gca)
%
%    See also colorMatchFcn, illuminant.

%    Copyright 1993-2005 The MathWorks, Inc.

hFig = ancestor(hAxesTarget, 'figure');

% Add a new axes to the same figure as the target axes.
figure(hFig);
hAxesSpectrum = axes;
set(hAxesSpectrum, 'visible', 'off')

% Areas outside the visible spectrum are black.
set(hAxesSpectrum, 'color', [0 0 0])

% Position the axes as appropriate.
targetPosition = get(hAxesTarget, 'position');
spectrumPosition = [targetPosition(1), ...
                    targetPosition(2), targetPosition(3), 1];
set(hAxesSpectrum, 'position', spectrumPosition)
set(hAxesSpectrum, 'units', 'pixels')

spectrumPosition = get(hAxesSpectrum, 'position');
set(hAxesSpectrum, 'position', [spectrumPosition(1), ...
                                spectrumPosition(2) - 20, ...
                                spectrumPosition(3), ...
                                20])

% Line the X limits of the two axes up and display the spectrum.
set(hAxesSpectrum, 'xtick', get(hAxesTarget, 'xtick'))

[lambda, RGB] = createSpectrum('1964_full');
axes(hAxesSpectrum);
image(lambda, 1:size(RGB,1), RGB);

% Turn off the unneeded axes labels.
set(hAxesTarget, 'xtick', []);
set(hAxesSpectrum, 'ytick', []);

% Make the figure visible.
set(hAxesSpectrum, 'units', 'normalized')
set(hAxesSpectrum, 'visible', 'on')

% Return the handle if requested.
if (nargout == 1)

    varargout{1} = hAxesSpectrum;

end

% Add a callback to update the label if the X limits of the target changes.
% TO DO.

repo: DeSCI


text:
function [data,units] = compute_dmean_wing_area(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  data{i} = diff(trx(fly).mean_wing_area) ./ trx(fly).dt;

end
units = parseunits('mm^2/s');

author: kristinbranson

repo: JAABA