text:
% ST | GRASTA | Grassmannian Robust Adaptive Subspace Tracking Algorithm (He et al. 2012)
% process_video('ST', 'GRASTA', 'dataset/demo.avi', 'output/demo_ST-GRASTA.avi');

% clear, clc;
% load('dataset/trafficdb/traffic_patches.mat');
% V = im2double(imgdb{100});
% [M,m,n,p] = convert_video3d_to_2d(V);
% I = reshape(M(:,1),m,n); % imshow(I);

%%% predefine your video parameters
%VIDEO_ROWS  =  320; % Row of each frame
%VIDEO_COLS  =  240; % Column of each fram
%DIM         = VIDEO_ROWS * VIDEO_COLS;
DIM = size(M,1);

%%% GRASTA parameters
subsampling                 = 1;   % how much patial information will be used in your application
OPTIONS.RANK                = 1;     % the estimated rank
OPTIONS.rho                 = 1.8;
OPTIONS.MAX_MU              = 10000; % set max_mu large enough for initial subspace training
OPTIONS.MIN_MU              = 1;
OPTIONS.ITER_MAX            = 20;
OPTIONS.DIM_M               = DIM;   % your data's dimension
OPTIONS.USE_MEX             = 0;     % If you do not have the mex-version of Alg 2
% please set Use_mex = 0.
OPTS                        = struct(); % initiate a empty struct for OPTS
status.init                 = 0;        % status of grasta at each iteration
U_hat                       = zeros(1); % initiate a zero U_hat at beginning

%%% Initial subspace training [optional]
% You may use some frames to train the initial subspace.
OPTIONS.CONSTANT_STEP       = 0; % use adaptive step-size for initial subspace training
max_cycles                  = 30;
training_frames             = 10;
for outiter = 1:max_cycles
  frame_order = randperm(training_frames);
  for i = 1:training_frames
    % prepare the training frame
    %I = imread(fname);
    %I = double(rgb2gray(I));
    %I = I/max(max(I));
    I = M(:,frame_order(i));

    % random subsampling the frame I
    Z = round(subsampling * DIM);
    rp = randperm(DIM);
    idx = rp(1:Z)';
    I_Omega = I(idx);

    % training the background
    [U_hat,status,OPTS] = grasta_stream(I_Omega,idx,U_hat,status,OPTIONS,OPTS);
    %fprintf('Training %d/%d ...\n',outiter,max_cycles);
  end
end

%%% Real-time background/foreground separation
OPTIONS.CONSTANT_STEP = 1e-2; % use small constant step-size
nframes = size(M,2);
for i = 1:nframes
  % prepare the image I whether it is saved as file or caputured by
  % camera
  % I = imread(fname);
  %I = double(rgb2gray(I));
  %I = I/max(max(I));
  I = M(:,i);

  % random subsampling the frame I
  Z = round(subsampling * DIM);
  rp = randperm(DIM);
  idx = rp(1:Z)';
  I_Omega = I(idx);

  % tracking the background
  [U_hat,status,OPTS] = grasta_stream(I_Omega,idx,U_hat,status,OPTIONS,OPTS);

  % bg_img is the background
  L_hat = U_hat * status.w * status.SCALE;
  %bg_img = reshape(U_hat * status.w * status.SCALE, VIDEO_ROWS,VIDEO_COLS);

  % s_img is the separated foreground
  %s_hat = I(:) - U_hat * status.w * status.SCALE;
  %s_img = reshape(s_hat,VIDEO_ROWS,VIDEO_COLS);
  S_hat = I - L_hat;

  L(:,i) = L_hat;
  S(:,i) = S_hat;

  %fprintf('Processing %d/%d ...\n',i,nframes);
end

%%
% show_results(M,L,S,hard_threshold(S),p,m,n);

author: andrewssobral

repo: lrslibrary


text:
% OtStart Start using the Omnitrans interface.
%
% Usage:
%    OtStart
%
% When finished, release the Omnitrans resources with a call to OtStop().
%
% To be able to use this interface, ensure the Omnitrans directory
% (with OtInterface.dll and all those BPL files) are on the system path.
% Otherwise, OtStart will fail.

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% 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
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

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

% 2006-02-28 Created (Johan Meijdam, VORtech)

function OtStart()
   Ot('start');
end

author: Sable

repo: mcbench-benchmarks


text:
function F  = obj_find_valve_param_table_3way(x,Q_r)
% Objective function to find out the required area vs. spool displacement
% relationship for 3-way valve.
% Copyright 2010 MathWorks, Inc.

% x - vector of variable parameters. It is constructed of ten area values
% a_1 ... a_max - cross-sectional area of the valve at 10 successive
% spool position s_1 ... s_10. The positions are set by shifting the spool
% from zero to its maximum opening at constant speed.

% Q_r - required flow rate expressed in percentage with respect to its
% maximum flow and maximum displacement

% Q_max - maximum flow rate

assignin('base','a_1', x(1));
assignin('base','a_2', x(2));
assignin('base','a_3', x(3));
assignin('base','a_4', x(4));
assignin('base','a_5', x(5));
assignin('base','a_6', x(6));
assignin('base','a_7', x(7));
assignin('base','a_8', x(8));
assignin('base','a_9', x(9));
assignin('base','a_10', x(10));

model = 'valve_testrig_flow_char_3way';
load_system(model);
sim(model);

k = [1 1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2]; % Weight multipliers

% Computing objective function
F = 0;
for j = 1:11
    F = F + k(j) * (yout(j) - Q_r(j))^2;
end
end

% EOF


text:
function t = full(t)
%FULL Convert a ktensor to a (dense) tensor.
%
%   T = FULL(C) converts a ktensor to a (dense) tensor.
%
%   Examples
%   X = ktensor([3; 2], rand(4,2), rand(5,2), rand(3,2));
%   Y = full(A) %<-- equivalent dense tensor
%
%   See also KTENSOR, TENSOR, KTENSOR/DOUBLE.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

sz = size(t);
data = t.lambda' * khatrirao(t.u,'r')';
t = tensor(data,sz);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function rs = center(s)

%tstoolbox/@signal/center
%   Syntax:
%     * center(s)
%
%   Center signal by removing it's mean.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(1,1);
if nargin  < 1, help(mfilename); end

N = dlens(s,1);
points = data(s);
c = core(points - repmat(mean(points), N,1));
rs = signal(c,s);

rs = addhistory(rs, ['Centered signal around zero']);
rs = addcommandlines(rs, 's = center(s');

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% run the tests for Travis CI

import matlab.unittest.plugins.CodeCoveragePlugin
import matlab.unittest.plugins.codecoverage.CoberturaFormat
import matlab.unittest.TestRunner

suite = testsuite('IncludeSubfolders', false);
runner = TestRunner.withTextOutput;

% add a coverage report
reportFile = fullfile('..', 'coverage.xml');
reportFormat = CoberturaFormat(reportFile);
plugin = CodeCoveragePlugin.forFolder('..', 'Producing',reportFormat);
runner.addPlugin(plugin);

%% setup the path
addpath ..

%% Run all unit tests in my repository

% Run all unit tests in my repository.
results = runner.run(suite);

% Assert no tests failed.
assert(all(~[results.Failed]));

%% Build the toolbox distribution file
fprintf('---------------------------------- Build the MLTBX file ------------------------------------\n')
cd ..
matlab.addons.toolbox.packageToolbox('PackageToolbox.prj')

author: petercorke

repo: spatialmath-matlab


text:
function value = r8vec_norm_l0 ( n, a )

%*****************************************************************************80
%
%% R8VEC_NORM_L0 returns the l0 "norm" of an R8VEC.
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%    The l0 "norm" simply counts the number of nonzero entries in the vector.
%    It is not a true norm, but has some similarities to one.  It is useful
%    in the study of compressive sensing.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, real A(N), the vector.
%
%    Output, real R8VEC_NORM_L0, the value of the norm.
%
  value = sum ( a(1:n) ~= 0.0 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = digamma(x)
%DIGAMMA   Digamma function.
% DIGAMMA(X) returns digamma(x) = d log(gamma(x)) / dx
% If X is a matrix, returns the digamma function evaluated at each element.

% Reference:
%
%    J Bernardo,
%    Psi ( Digamma ) Function,
%    Algorithm AS 103,
%    Applied Statistics,
%    Volume 25, Number 3, pages 315-317, 1976.
%
% From http://www.psc.edu/~burkardt/src/dirichlet/dirichlet.f

large = 9.5;
d1 = -0.5772156649015328606065121;  % digamma(1)
d2 = pi^2/6;
small = 1e-6;
s3 = 1/12;
s4 = 1/120;
s5 = 1/252;
s6 = 1/240;
s7 = 1/132;
s8 = 691/32760;
s9 = 1/12;
s10 = 3617/8160;

% Initialize
y = zeros(size(x));

% illegal arguments
i = find(x == -Inf | isnan(x));
if ~isempty(i)
  x(i) = NaN;
  y(i) = NaN;
end

% Negative values
i = find(x < 0);
if ~isempty(i)
  % Use the reflection formula (Jeffrey 11.1.6):
  % digamma(-x) = digamma(x+1) + pi*cot(pi*x)
  y(i) = digamma(-x(i)+1) + pi*cot(-pi*x(i));
  % This is related to the identity
  % digamma(-x) = digamma(x+1) - digamma(z) + digamma(1-z)
  % where z is the fractional part of x
  % For example:
  % digamma(-3.1) = 1/3.1 + 1/2.1 + 1/1.1 + 1/0.1 + digamma(1-0.1)
  %               = digamma(4.1) - digamma(0.1) + digamma(1-0.1)
  % Then we use
  % digamma(1-z) - digamma(z) = pi*cot(pi*z)
end

i = find(x == 0);
if ~isempty(i)
  y(i) = -Inf;
end

%  Use approximation if argument <= small.
i = find(x > 0 & x <= small);
if ~isempty(i)
  y(i) = y(i) + d1 - 1 ./ x(i) + d2*x(i);
end

%  Reduce to digamma(X + N) where (X + N) >= large.
while(1)
  i = find(x > small & x < large);
  if isempty(i)
    break
  end
  y(i) = y(i) - 1 ./ x(i);
  x(i) = x(i) + 1;
end

%  Use de Moivre's expansion if argument >= large.
% In maple: asympt(Psi(x), x);
i = find(x >= large);
if ~isempty(i)
  r = 1 ./ x(i);
  y(i) = y(i) + log(x(i)) - 0.5 * r;
  r = r .* r;
  y(i) = y(i) - r .* ( s3 - r .* ( s4 - r .* (s5 - r .* (s6 - r .* s7))));
end


text:
function idxs = get_unary_idxs_global_sort(unProb,max_detections)

unProbThresh = unProb;
scores = sort(reshape(unProb, size(unProb,1)*size(unProb,2),1));
j = 1;
while (length(find(sum(unProbThresh,2) > 0)) > max_detections && j <= length(scores))
    unProbThresh = unProb;
    for m = 1:size(unProb,2)
        unProbThresh(unProb(:,m) < scores(j),m) = 0;
    end
    j = j + 1;
end

min_det_score_found = scores(j);
fprintf('found min_det_score: %1.2f\n',min_det_score_found);

idxs = find(sum(unProbThresh,2) > 0);
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


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

%% Syntax
% |[C_smooth]=patchSmoothFaceMeasure(F,V,C,smoothPar);|

%% Description
%

%% Examples

%%
clear; close all; clc;

%%
% Plot settings
markerSize=50;

%% Demonstrating

smoothPar.lambda=0.5; %Lambda for smoothing
smoothPar.n=25; %Number of smooth iterations

for testCase=1:6
    [F,V]=graphicsModels(testCase);
    switch testCase
        case 1 %Cut patch data with unused points
            [F,V]=geoSphere(4,1);
            Z=V(:,3);
            F=F(mean(Z(F),2)>0,:);
        case 2 %Cut patch data
            [F,V]=geoSphere(4,1);
            Z=V(:,3);
            F=F(mean(Z(F),2)>0,:);
            [F,V]=patchCleanUnused(F,V);
        case 3 %Quadrilateral patch data
            [F,V]=quadSphere(4,1);
        case 4 %A larger mesh
            [F,V]=graphicsModels(1);
        case 5 %A tangled mesh
            [F,V]=graphicsModels(8);
        case 6
            r=1; %Sphere radius
            rc=2.5; %Central radius
            nr=16;
            nc=25;
            patchType='honey';
            [F,V]=patchTorus(r,nr,rc,nc,patchType);
    end

    %%
    Z=V(:,3);
    C_clean=mean(Z(F),2);
    C=C_clean+(max(C_clean)-min(C_clean))/5*randn(size(F,1),1);

    %%
    [C_smooth]=patchSmoothFaceMeasure(F,V,C,smoothPar);

    %%
    cFigure;
    subplot(1,2,1); hold on;
    title('Raw face data');
    gpatch(F,V,C,'none');
    colorbar(gca); caxis([min(C_clean) max(C_clean)]);
    axisGeom;
    camlight headlight;

    subplot(1,2,2); hold on;
    title('Smoothed face data');
    gpatch(F,V,C_smooth,'none');
    colorbar(gca); caxis([min(C_clean) max(C_clean)]);
    axisGeom;
    camlight headlight;

    drawnow;

end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [ n_data, x, fx ] = elliptic_ka_values ( n_data )

%*****************************************************************************80
%
%% ELLIPTIC_KA_VALUES returns values of the complete elliptic integral K(ALPHA).
%
%  Discussion:
%
%    This is one form of what is sometimes called the complete elliptic integral
%    of the first kind.
%
%    The function is defined by the formula:
%
%      K(ALPHA) = integral ( 0 <= T <= PI/2 )
%        dT / sqrt ( 1 - sin ( ALPHA )^2 * sin ( T )^2 )
%
%    In Mathematica, the function can be evaluated by:
%
%      EllipticK[(Sin[alpha*Pi/180])^2]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function, measured
%    in degrees.
%
%    Output, real FX, the value of the function.
%
  n_max = 18;

  fx_vec = [ ...
     0.1570796326794897E+01, ...
     0.1573792130924768E+01, ...
     0.1582842804338351E+01, ...
     0.1598142002112540E+01, ...
     0.1620025899124204E+01, ...
     0.1648995218478530E+01, ...
     0.1685750354812596E+01, ...
     0.1731245175657058E+01, ...
     0.1786769134885021E+01, ...
     0.1854074677301372E+01, ...
     0.1935581096004722E+01, ...
     0.2034715312185791E+01, ...
     0.2156515647499643E+01, ...
     0.2308786798167196E+01, ...
     0.2504550079001634E+01, ...
     0.2768063145368768E+01, ...
     0.3153385251887839E+01, ...
     0.3831741999784146E+01 ];

  x_vec = [ ...
      0.0E+00, ...
      5.0E+00, ...
     10.0E+00, ...
     15.0E+00, ...
     20.0E+00, ...
     25.0E+00, ...
     30.0E+00, ...
     35.0E+00, ...
     40.0E+00, ...
     45.0E+00, ...
     50.0E+00, ...
     55.0E+00, ...
     60.0E+00, ...
     65.0E+00, ...
     70.0E+00, ...
     75.0E+00, ...
     80.0E+00, ...
     85.0E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function J=imshift_3D(I,u,kernel)
% USAGE    : J=imshift_3D(I,u);
% FUNCTION : Shifts the 3D image I by u{1} along the lines and by u{2}
% along the columns and u{3} in the third dimension. 'u' is an arbitrary
% vector field structured as:
%                   u{1} = u1(x,y,z)
%                   u{2} = u2(x,y,z)
%                   u{3} = u2(x,y,z).
% If 'u' is a single vector (i.e. u1 = a, u2 = b and u3 = c) then all the
% pixels are shifted by the same amount. Note the structure of u is alwasys
% assumed to be a cell.
%
% This function uses cubic spline interpolation used in 'interp_3D.m' with
% half point symmetric image extensions (can be changed).
%
% DATE     : 23 November 2014 (updated 28 Jan 2015 by C Gilliam)
% AUTHOR   : Thierry Blu, mailto:thierry.blu@m4x.org

[M,N,P]=size(I);

if nargin == 2,
%     kernel= 'cubicspline';
    kernel='cubicOMOMS';
% kernel = 'bilinear';
end

if length(u{1})>1               % check to see if u is a single vector or a field
    [M0,N0,P0]=size(u{1});
    if M0==M&N0==N&P0==P

        % Generate x, y and z:
        x = repmat((1:M)',[1,N,P]);
        y = repmat((1:N), [M,1,P]);
        z = repmat(shiftdim((1:P),-1), [M,N,1]);

        % Interpolation:
        J=interp_3D(double(x-u{1}),double(y-u{2}),double(z-u{3}),double(I),kernel);  %single/double
        J=cast(J,class(I));
    else
        error('Input image and flow dimensions do not match!')
    end
else
    u=[u{:}];
    integer_shift=(norm(round(u)-u)<=1e-6);

    if integer_shift
        u=round(u);
        xshift=abs(u(1));
        yshift=abs(u(2));
        zshift=abs(u(3));
        switch sign(u(1))
            case -1
                I1=[I;I(end-(1:xshift)+1,:,:)];
            case 0
                I1=I;
            case 1
                I1=[I((xshift:-1:1),:,:);I];
        end
        switch sign(u(2))
            case -1
                I1=[I1 I1(:,end-(1:yshift)+1,:)];
            case 0
            case 1
                I1=[I1(:,(yshift:-1:1),:) I1];
        end

        switch sign(u(3))
            case -1
                I1=cat(3,I1,I1(:,:,end-(1:zshift)+1));
            case 0
            case 1
                I1=cat(3,I1(:,:,(zshift:-1:1)), I1);
        end

        J=I1((1:M)+max(0,-u(1)),(1:N)+max(0,-u(2)), (1:P)+max(0,-u(3)));
    else
        % Generate x, y and z:
        x = repmat((1:M)',[1,N,P]);
        y = repmat((1:N), [M,1,P]);
        z = repmat(shiftdim((1:P),-1), [M,N,1]);

        % Interpolation:
        J=interp_3D(x-u(1),y-u(2),z-u(3),double(I),kernel);
        J=cast(J,class(I));
    end
end

author: thomaskuestner


text:
% Fig. W1 Web Appendix W8    Feedback Control of Dynamic Systems, 6e
%                                 Franklin, Powell, Emami
%

clear all;
close all;

F=[0 1;0 0];
G=[0;1];
H=[1 0];
J=0;
T=1;
[Phi,Gam]=c2d(F,G,T);
j=sqrt(-1);
Pc=[.78+.18*j;.78-.18*j];
K=acker(Phi,Gam,Pc);
Pe=[.2+.2*j;.2-.2*j];
L=acker(Phi',H',Pe)';
[A,B,C,D]=dreg(Phi,Gam,H,J,K,L);
A=Phi-Gam*K-L*H;
B=L;
C=K;
D=0;
[Ac,Bc,Cc,Dc]=series(A,B,C,D,Phi,Gam,H,J);
[Acl,Bcl,Ccl,Dcl]=feedback(Ac,Bc,Cc,Dc,0,0,0,1);
tf=30;
N=tf/T+1;
td=0:1:tf;
yd=dstep(Acl,Bcl,Ccl,Dcl,1,N);
axis([0 30 0 1.5])
plot(td,yd,'-',td,yd,'*'),
xlabel('time (sec)')
ylabel('plant output  y(t)')
title('Fig. 8.20  Step response of the continuous and digital systems')
nicegrid;
hold on

% use command structure from section 7.8

Nx=[1;0];
A2=[Phi  -Gam*K;
   L*H  Phi-L*H-Gam*K];
B2=[Gam*K*Nx;Gam*K*Nx];
C2=[1 0 0 0];
D2=0;
y2=dstep(A2,B2,C2,D2,1,N);
plot(td,y2,'-',td,y2,'o')
text(6.5,.25,'o-----o-----o-----o  Command structure from Fig 7.48(b)')
text(6.5,.45,'*-----*-----*-----*  Command structure from Fig 7.15')
hold off


text:
function animatorz(action) % autogenerated function wrapper
    % turned into function by Celso G Reyes 2017

    animator(action, @slicemapz);
    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun(mfilename('fullpath'));

    global ps1 ps2 plin pli

    switch(action)
        case 'start'
            [ps1, ps2, plin, pli] = animator_start(@animatorz);% ButtonMotion, ButtonUp

        case 'move'
            animator_move(ps2, pli, plin)
        case 'stop'
            animator_stop(gcf);
            slicemapz('newslice');
    end
end

author: CelsoReyes


text:
function imgOut = ReinhardGaussianFilter(img, s, alpha_i)
%
%
%      imgOut = ReinhardGaussianFilter(img, s, alpha_i)
%
%
%      Input:
%          -img: an image to be filtered
%          -s: size of the filter in pixel
%          -alpha_i:
%
%      Output:
%          -imgOut: the filtered image
%
%     Copyright (C) 2010 Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

alpha_s_sq = (alpha_i * s)^2;

% %Kernel of the filter
% s2 = round(s * 5);
% [X,Y] = meshgrid(-s2:s2, -s2:s2);
% H = exp(-(X.^2 + Y.^2) / alpha_s_sq) / (pi * alpha_s_sq);
%
% %Filtering
% imgOut = imfilter(img, H, 'replicate');

[r, c, col] = size(img);

[x, y] = meshgrid(1:c, 1:r);
x = x - (c / 2);
y = y - (r / 2);
kernel = exp(-(x.^2 + y.^2) / alpha_s_sq) / (pi * alpha_s_sq);
kernel_f = fft2(kernel);

imgOut = zeros(size(img));

for i=1:col
    imgOut(:,:,i) = fftshift(ifft2(fft2(img(:,:,i)) .* kernel_f));
end

end

author: banterle


text:
function varargout = GUIEditSegyTraceHeader(varargin)
% GUIPlotXY Application M-file for GUIEditSegyTraceHeader.fig
%    FIG = GUIPLOTXY launch GUIEditSegyTraceHeader GUI.
%    GUIPLOTXY('callback_name', ...) invoke the named callback.
%
% OBSOLETE ??
%

% Last Modified by GUIDE v2.0 16-Jun-2002 11:41:23

if (nargin == 0)|(isstruct(varargin{1}))  % LAUNCH GUI

	fig = openfig(mfilename,'reuse');

	% Use system color scheme for figure:
	set(fig,'Color',get(0,'defaultUicontrolBackgroundColor'));

	% Generate a structure of handles to pass to callbacks, and store it.
	handles = guihandles(fig);
	guidata(fig, handles);

    if nargin>0,
       data = guihandles(fig); % initialize it to contain handles
       data.H=varargin{1};
       data.Hname=fieldnames(data.H(1));

       data.NShowHeaders=7; % NUMBER OF HEADER VALUES TO LIST
       data.ShowNTraces=10; % NUMBER OF TRACES TO LIST
       data.StartTrace=1;

       set(handles.popX,'String',data.Hname);
       set(handles.popY,'String',data.Hname);

       set(handles.popX,'Value',2);
       set(handles.popY,'Value',4);

       guidata(fig,data);

       set(fig,'HandleVisibility','On')
       GUIPlotXY('actionPlotData',fig,handles)
       set(fig,'HandleVisibility','CallBack')

    end

	% Wait for callbacks to run and window to be dismissed:
	uiwait(fig);

	if nargout > 0
	  data=guidata(fig);
	  if isstruct(data.H)
	    varargout{1}=data.H
	  else
	    varargout{1} = fig;
	  end
    end

    close(handles.figure1);
    try;set(handles.figure1,'Visible','off');end

elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK

	try
		if (nargout)
			[varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
		else
			feval(varargin{:}); % FEVAL switchyard
		end
	catch
		disp(lasterr);
	end

end

%| ABOUT CALLBACKS:
%| GUIDE automatically appends subfunction prototypes to this file, and
%| sets objects' callback properties to call them through the FEVAL
%| switchyard above. This comment describes that mechanism.
%|
%| Each callback subfunction declaration has the following form:
%| <SUBFUNCTION_NAME>(H, EVENTDATA, HANDLES, VARARGIN)
%|
%| The subfunction name is composed using the object's Tag and the
%| callback type separated by '_', e.g. 'slider2_Callback',
%| 'figure1_CloseRequestFcn', 'axis1_ButtondownFcn'.
%|
%| H is the callback object's handle (obtained using GCBO).
%|
%| EVENTDATA is empty, but reserved for future use.
%|
%| HANDLES is a structure containing handles of components in GUI using
%| tags as fieldnames, e.g. handles.figure1, handles.slider2. This
%| structure is created at GUI startup using GUIHANDLES and stored in
%| the figure's application data using GUIDATA. A copy of the structure
%| is passed to each callback.  You can store additional information in
%| this structure at GUI startup, and you can change the structure
%| during callbacks.  Call guidata(h, handles) after changing your
%| copy to replace the stored original so that subsequent callbacks see
%| the updates. Type "help guihandles" and "help guidata" for more
%| information.
%|
%| VARARGIN contains any extra arguments you have passed to the
%| callback. Specify the extra arguments by editing the callback
%| property in the inspector. By default, GUIDE sets the property to:
%| <MFILENAME>('<SUBFUNCTION_NAME>', gcbo, [], guidata(gcbo))
%| Add any extra arguments after the last argument, before the final
%| closing parenthesis.

% --------------------------------------------------------------------

% --------------------------------------------------------------------
function varargout = popX_Callback(h, eventdata, handles, varargin)
actionPlotData(h,handles)

% --------------------------------------------------------------------
function varargout = popY_Callback(h, eventdata, handles, varargin)
actionPlotData(h,handles)

% --------------------------------------------------------------------
function varargout = figure1_ResizeFcn(h, eventdata, handles, varargin)

% --------------------------------------------------------------------
function varargout = pbDone_Callback(h, eventdata, handles, varargin)
uiresume;

function actionPlotData(h,handles)
data=guidata(h);
Xv=get(handles.popX,'value');
Yv=get(handles.popY,'value');

Xheader=data.Hname(Xv);
Yheader=data.Hname(Yv);

X=[data.H.(char(Xheader))];
Y=[data.H.(char(Yheader))];

axes(handles.axMain);
plot(X,Y,'k.');
xlabel(char(Xheader))
ylabel(char(Yheader))
zoom on

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
function ma_setAnimatorSliderBounds(handles, stateLog)
    if(size(stateLog,1) > 1)
        warpRate = handles.warpRateCombo.UserData;

        if(isempty(warpRate))
            lowerStep = 1/(size(stateLog,1)-1);
        else
            lowerStep = abs(warpRate / (stateLog(end,1) - stateLog(1,1)));
        end

        if(lowerStep <= 0.1)
            upperStep = 0.1;
        else
            upperStep = lowerStep;
        end

        set(handles.movieFrameSlider,'Min',stateLog(1,1),'Max',stateLog(end,1),'Value',stateLog(1,1));
        set(handles.movieFrameSlider,'SliderStep',[lowerStep,upperStep]);
        set(handles.movieFrameSlider,'TooltipString', sprintf('Current step size: %0.3f s', warpRate));
    else
        set(handles.movieFrameSlider,'Min',stateLog(1,1),'Max',stateLog(1,1)+1,'Value',stateLog(1,1));
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function grad = B_log(future_layers, input,curr_layer)
future_grad = GetFutureGrad(future_layers, curr_layer);

grad = 1./(input+curr_layer.const).*future_grad;

end

author: singaxiong

repo: SignalGraph


text:
function dl=lpcaa2dl(aa)
%LPCAA2DL LPC: Convert area coefficients to dct of log area DL=(AA)

% note: we do not correct for sinc distortion; perhaps we should multiply by
% k=1:p-1;s=[sqrt(0.5)/p 2*sin(pi*k/(2*p))./(pi*k)];

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: lpcaa2dl.m,v 1.4 2007/05/04 07:01:38 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,p2]=size(aa);
dl=rdct(log(aa(:,2:p2-1)./aa(:,p2*ones(1,p2-2))).').';

author: decouples


text:
function aedat = OffsetTime(aedat, offset)

%{
Offset all data timestamps by the(positive going) offset parameter.
No trimming at zero - this can result in negative timestamps.
Recalculate first and last timestamps
%}

dbstop if error

if ~isfield(aedat, 'data')
    disp('No data found')
    return
end

% Are we working with 32 or 64 bit timestamps? It depends on the fileFormat
% from which we imported; default is 32 bit
if ~isfield(aedat.info, 'fileFormat') || aedat.info.fileFormat == 2
    % Assume offset param is given in secs (AedatTools standard)
    offset = int32(offset * 1e6);
    % Special
    if isfield(aedat.data, 'special')
        aedat.data.special.timeStamp = uint32(int32(aedat.data.special.timeStamp) + offset);
    end
    % Polarity
    if isfield(aedat.data, 'polarity')
        aedat.data.polarity.timeStamp = uint32(int32(aedat.data.polarity.timeStamp) + offset);
    end

    % Frames
    % This assumes that timestamps have been simplified to aedat2 standard, if
    % they came from aedat3 file
    if isfield(aedat.data, 'frame')
        aedat.data.frame.timeStampStart = uint32(int32(aedat.data.frame.timeStampStart)   + offset);
        aedat.data.frame.timeStampEnd   = uint32(int32(aedat.data.frame.timeStampEnd)     + offset);
    end
    % Imu6
    if isfield(aedat.data, 'imu6')
        aedat.data.imu6.timeStamp = uint32(int32(aedat.data.imu6.timeStamp) + offset);
    end

    % Sample
    if isfield(aedat.data, 'sample')
        aedat.data.sample.timeStamp = uint32(int32(aedat.data.sample.timeStamp) + offset);
    end

    % Ear
    if isfield(aedat.data, 'ear')
        aedat.data.ear.timeStamp = uint32(int32(aedat.data.ear.timeStamp) + offset);
    end

    % Point1D
    if isfield(aedat.data, 'point1D')
        aedat.data.point1D.timeStamp = uint32(int32(aedat.data.point1D.timeStamp) + offset);
    end

    %Point2D
    if isfield(aedat.data, 'point2D')
        aedat.data.point2D.timeStamp = uint32(int32(aedat.data.point2D.timeStamp) + offset);
    end

    %Point3D
    if isfield(aedat.data, 'point3D')
        aedat.data.point3D.timeStamp = uint32(int32(aedat.data.point3D.timeStamp) + offset);
    end

    %Point4D
    if isfield(aedat.data, 'point4D')
        aedat.data.point4D.timeStamp = uint32(int32(aedat.data.point4D.timeStamp) + offset);
    end
else

    % Assume offset param is given in secs (AedatTools standard)
    offset = int64(offset * 1e6);
    % Special
    if isfield(aedat.data, 'special')
        aedat.data.special.timeStamp = uint64(int64(aedat.data.special.timeStamp) + offset);
    end
    % Polarity
    if isfield(aedat.data, 'polarity')
        aedat.data.polarity.timeStamp = uint64(int64(aedat.data.polarity.timeStamp) + offset);
    end

    % Frames
    % This assumes that timestamps have been simplified to aedat2 standard, if
    % they came from aedat3 file
    if isfield(aedat.data, 'frame')
        aedat.data.frame.timeStampStart = uint64(int64(aedat.data.frame.timeStampStart)   + offset);
        aedat.data.frame.timeStampEnd   = uint64(int64(aedat.data.frame.timeStampEnd)     + offset);
    end
    % Imu6
    if isfield(aedat.data, 'imu6')
        aedat.data.imu6.timeStamp = uint64(int64(aedat.data.imu6.timeStamp) + offset);
    end

    % Sample
    if isfield(aedat.data, 'sample')
        aedat.data.sample.timeStamp = uint64(int64(aedat.data.sample.timeStamp) + offset);
    end

    % Ear
    if isfield(aedat.data, 'ear')
        aedat.data.ear.timeStamp = uint64(int64(aedat.data.ear.timeStamp) + offset);
    end

    % Point1D
    if isfield(aedat.data, 'point1D')
        aedat.data.point1D.timeStamp = uint64(int64(aedat.data.point1D.timeStamp) + offset);
    end

    %Point2D
    if isfield(aedat.data, 'point2D')
        aedat.data.point2D.timeStamp = uint64(int64(aedat.data.point2D.timeStamp) + offset);
    end

    %Point3D
    if isfield(aedat.data, 'point3D')
        aedat.data.point3D.timeStamp = uint64(int64(aedat.data.point3D.timeStamp) + offset);
    end

    %Point4D
    if isfield(aedat.data, 'point4D')
        aedat.data.point4D.timeStamp = uint64(int64(aedat.data.point4D.timeStamp) + offset);
    end
end

aedat = FindFirstAndLastTimeStamps(aedat);

author: panpanfei


text:
%% Simulation Data
simu = simulationClass();                   % Initialize Simulation Class
simu.simMechanicsFile = 'RM3.slx';          % Specify Simulink Model File
simu.rampTime = 100;                        % Wave Ramp Time [s]
simu.endTime=200;                           % Simulation End Time [s]
simu.dt = 0.1;                              % Simulation Time-Step [s]
simu.explorer = 'off';
simu.cicEndTime = 20;

%% Wave Information
% Irregular Waves using BS Spectrum with State Space Calculation
waves = waveClass('irregular');             % Initialize Wave Class and Specify Type
waves.height = 2.5;                              % Significant Wave Height [m]
waves.period = 8;                                % Peak Period [s]
waves.spectrumType = 'PM';                  % Specify Wave Spectrum Type
waves.phaseSeed = 1;

%% Body Data
% Float
body(1) = bodyClass('../../../../examples/RM3/hydroData/rm3.h5');
body(1).geometryFile = '../../../../examples/RM3/geometry/float.stl';
body(1).mass = 'equilibrium';
body(1).inertia = [20907301 21306090.66 37085481.11];

% Spar/Plate
body(2) = bodyClass('../../../../examples/RM3/hydroData/rm3.h5');
body(2).geometryFile = '../../../../examples/RM3/geometry/plate.stl';
body(2).mass = 'equilibrium';
body(2).inertia = [94419614.57 94407091.24 28542224.82];

%% PTO and Constraint Parameters
% Floating (3DOF) Joint
constraint(1) = constraintClass('Constraint1'); % Initialize Constraint Class for Constraint1
constraint(1).location = [0 0 0];                    % Constraint Location [m]

% Translational PTO
pto(1) = ptoClass('PTO1');                      % Initialize PTO Class for PTO1
pto(1).stiffness = 0;                                   % PTO Stiffness [N/m]
pto(1).damping = 1200000;                             % PTO Damping [N/(m/s)]
pto(1).location = [0 0 0];                           % PTO Location [m]

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function [secInMin, secInHr, secInDay, secInYear] = getSecondsInVariousTimeUnits()
%getSecondsInVariousTimeUnits Summary of this function goes here
%   Detailed explanation goes here
    global options_UseEarthTimeSystem;

    secInMin = 60;
    secInHr = 60*secInMin;
	if(options_UseEarthTimeSystem == true)
        secInDay = 24*secInHr;
        secInYear = 365*secInDay;
    else
        secInDay = 6*secInHr;
        secInYear = 426*secInDay;
	end
end


text:
function param = robThreeDynamicsExtractParam(model)

% ROBTHREEDYNAMICSEXTRACTPARAM Extract parameters from the robot three dynamics model.

% FGPLVM

param = [];

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function varargout = ge_html(varargin)
% GE_HTML M-file for ge_html.fig
%      GE_HTML, by itself, creates a new GE_HTML or raises the existing
%      singleton*.
%
%      H = GE_HTML returns the handle to a new GE_HTML or the handle to
%      the existing singleton*.
%
%      GE_HTML('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in GE_HTML.M with the given input arguments.
%
%      GE_HTML('Property','Value',...) creates a new GE_HTML or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before ge_html_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to ge_html_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help ge_html

% Last Modified by GUIDE v2.5 03-Aug-2011 16:39:30

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @ge_html_OpeningFcn, ...
                   'gui_OutputFcn',  @ge_html_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before ge_html is made visible.
function ge_html_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to ge_html (see VARARGIN)

% Open ActiveX-Control from Microsoft-Internet-Explorer
handles.ie = actxcontrol('Shell.Explorer.2',[2 2 996 700 ]);

% Choose default command line output for ge_html
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes ge_html wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = ge_html_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Navigate to the html;
% Navigate(handles.ie,'C:\Users\sk\Documents\MATLAB\ge_ma.html');

% locate to the ge_ma.html
[FileName,PathName,FilterIndex] = uigetfile('*.html','open GE-html ') ;
if isequal(FileName,0)
   disp('User selected Cancel')
else
   Navigate(handles.ie,[PathName, FileName]);
end;

% Wait for initialising
pause(4);   % if pause is too short so that Matlab produce error, take a little more time

% Get the Document
handles.myDoc = handles.ie.Document;

% Update handles structure
guidata(hObject, handles);

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% show/ hide borders
handles.myDoc.parentWindow.execScript('toggleBorders()', 'Jscript');

% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.myDoc.parentWindow.execScript('toggleRotation()', 'Jscript');

% --- Executes on button press in pushbutton4.
function pushbutton4_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.myDoc.parentWindow.execScript('toggleNav()', 'Jscript');

% --- Executes on button press in pushbutton5.
function pushbutton5_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.myDoc.parentWindow.execScript('toggleRoads()', 'Jscript');

% --- Executes on button press in pushbutton6.
function pushbutton6_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.myDoc.parentWindow.execScript('showBuildings()', 'Jscript');

% --- Executes on button press in pushbutton7.
function pushbutton7_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton7 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.myDoc.parentWindow.execScript('hideBuildings()', 'Jscript');

function edit1_Callback(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit1 as text
%        str2double(get(hObject,'String')) returns contents of edit1 as a double

% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in pushbutton8.
function pushbutton8_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton8 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

searchStr = get(handles.edit1,'String');
handles.myDoc.parentWindow.execScript(['findmyaddress(''',searchStr ,''')' ], 'Jscript');

author: Sable


text:
function rgb = createRGBStack(img1, img2, img3)
%CREATERGBSTACK Concatenate 2 or 3 grayscale stacks to form a color stack
%
%   RGB = createRGBStack(RED, GREEN)
%   RGB = createRGBStack(RED, GREEN, BLUE)
%   Create a new 3D image containing RGB values from 2 or 3 grayscale
%   images. All input images must have the same size M-by-N-ny-P.
%   The resulting RGB image has size M-by-N-by-3-by-P.
%
%   It is also possible to specify empty red or green channel:
%   RGB = createRGBStack([], GREEN, BLUE)
%   RGB = createRGBStack(RED, [], BLUE)
%
%   Example
%     % Colorize in red the bright parts of a 3D image
%     metadata = analyze75info('brainMRI.hdr');
%     I = analyze75read(metadata) * 3;
%     I2 = I;
%     I2(I > 200) = 0;
%     rgb = createRGBStack(I, I2, I2);
%     slicer(rgb)
%
%   See also
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-10-26,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

%% Process input arguments

% ensure initialisation of variable
if nargin < 3
    img3 = [];
end

% choose one of the 3 inputs as reference
ref = img1;
if isempty(img1)
    ref = img2;
    if isempty(img2)
        ref = img3;
        if isempty(ref)
            error('At least one channel must be specified');
        end
    end
end

%% Initialize result

% size of reference and new image
dim = size(ref);
newDim = [dim(1:2) 3 dim(3)];

% create empty result image
if islogical(ref)
    rgb = false(newDim);
else
    rgb = zeros(newDim, class(ref));
end

%% Fill up result with data

% Red channel
if ~isempty(img1)
    rgb(:,:,1,:) = img1;
end

% green channel
if ~isempty(img2)
    if sum(size(img2) ~= dim) > 0
        error('Red and Green images must have the same size');
    end
    rgb(:,:,2,:) = img2;
end

% blue channel
if ~isempty(img3)
    if sum(size(img3) ~= dim) > 0
        error('Blue and reference images must have the same size');
    end
    rgb(:,:,3,:) = img3;
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function files_multiple ( )

%*****************************************************************************80
%
%% FILES_MULTIPLE demonstrates how to work with multiple files.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 December 2012
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FILES_MULTIPLE:\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Demonstrate how to work with multiple files.\n' );

  files_multiple_test01 ( );
  files_multiple_test02 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FILES_MULTIPLE:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function files_multiple_test01 ( )

%*****************************************************************************80
%
%% TEST01 writes data to four files which are open simultaneously.
%
%  Discussion:
%
%    The files contain the integers which are divisible by 2, 3, 5, and 7.
%
%    divisor1.txt  divisor2.txt  divisor3.txt  divisor4.txt
%        0             0             0             0
%        2             3             5             7
%        4             6            10            14
%        6             9            15            21
%        8            12            20            28
%       ...
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2012
%
%  Author:
%
%    John Burkardt
%

%
%  FILENUM can be increased, and the program will still work the same way.
%
  filenum = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FILES_MULTIPLE_TEST01\n' );
  fprintf ( 1, '  Create 4 files, and count to 100.\n' );
  fprintf ( 1, '  "divisor1.txt" will contain multiples of 2;\n' );
  fprintf ( 1, '  "divisor2.txt" will contain multiples of 3;\n' );
  fprintf ( 1, '  "divisor3.txt" will contain multiples of 5;\n' );
  fprintf ( 1, '  "divisor4.txt" will contain multiples of 7;\n' );
%
%  Set the divisors.
%  PRIME is a function which can return the I-th prime number.
%
  for i = 1 : filenum
    divisor(i) = prime ( i );
  end
%
%  Define the file names for an arbitrary number of files
%  by incrementing the numeric information in a "template".
%
  template = 'divisor0.txt';

  for i = 1 : filenum
    template = filename_inc ( template );
    filename(i,:) = template(:);
  end
%
%  Associate a unit number with each file.
%
  for i = 1 : filenum
    fileunit(i) = fopen ( filename(i,:), 'wt' );
    fprintf ( fileunit(i), 'Multiples of %d\n', divisor(i) );
  end

  for j = 0 : 100
    for i = 1 : filenum
      if ( mod ( j, divisor(i) ) == 0 )
        fprintf ( fileunit(i), '%d\n', j );
      end
    end
  end
%
%  Close the files.
%
  for i = 1 : filenum
    fclose ( fileunit(i) );
  end

  return
end
function files_multiple_test02 ( )

%*****************************************************************************80
%
%% FILES_MULTIPLE_TEST02 writes selected data to four files which are open simultaneously.
%
%  Discussion:
%
%    The vector X contains 100 random numbers, and we update X 20 times.
%    We want to keep track of X(10), X(25), X(64) and X(81) on each step.
%
%              x1.txt        x2.txt        x3.txt        x4.txt
%           X(10) values  X(25) values  X(64) values  X(81) values
%     1         2.1           1.7           5.3           9.8
%     2         2.0           3.6           6.2           7.3
%     3         2.5           4.1           5.9           6.8
%    ..
%    10         5.0           5.2           5.3           5.6
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 December 2012
%
%  Author:
%
%    John Burkardt
%

%
%  FILENUM can be increased, but more entries in K would have to be defined
%  by the user.
%
  filenum = 4;
  n = 100;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FILES_MULTIPLE_TEST02\n' );
  fprintf ( 1, '  Create 4 files.\n' );
  fprintf ( 1, '  The vector X contains 100 values,\n' );
  fprintf ( 1, '  and is set 20 times.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Write X(10), X(25), X(64) and X(81) to the files\n' );
  fprintf ( 1, '  "x1.txt", "x2.txt", "x3.txt" and "x4.txt"\n' );
%
%  Define the selected X indices.
%
  k(1:4) = [ 10, 25, 64, 81 ];
%
%  Define the file names for an arbitrary number of files
%  by incrementing the numeric information in a "template".
%
  template = 'x0.txt';

  for i = 1 : filenum
    template = filename_inc ( template );
    filename(i,:) = template(:);
  end
%
%  Associate a unit number with each file
%
  for i = 1 : filenum
    fileunit(i) = fopen ( filename(i,:), 'wt' );
    fprintf ( fileunit(i), 'Values of X(%d)\n', k(i) );
  end
%
%  Initial X vector is random.
%
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_01 ( n, seed );
%
%  Write selected data to file, then update data.
%
  for j = 0 : 20

    for i = 1 : filenum
      fprintf ( fileunit(i), '  %2d  %14.6g\n', j, x(k(i)) );
    end

    if ( j < 20 )
      y(1) = ( x(1) + x(2) ) / 2.0;
      y(2:99) = ( x(1:98) + x(2:99) + x(3:100) ) / 3.0;
      y(100) = ( x(99) + x(100) ) / 2.0;

      x(1:100) = y(1:100);
    end

  end
%
%  Close the files.
%
  for i = 1 : filenum
    fclose ( fileunit(i) );
  end

  return
end
function filename = filename_inc ( filename )

%*****************************************************************************80
%
%% FILENAME_INC increments a partially numeric filename.
%
%  Discussion:
%
%    It is assumed that the digits in the name, whether scattered or
%    connected, represent a number that is to be increased by 1 on
%    each call.  If this number is all 9's on input, the output number
%    is all 0's.  Non-numeric letters of the name are unaffected..
%
%    If the name is empty, then the routine stops.
%
%    If the name contains no digits, the empty string is returned.
%
%  Example:
%
%      Input            Output
%      -----            ------
%      'a7to11.txt'     'a7to12.txt'  (typical case.  Last digit incremented)
%      'a7to99.txt'     'a8to00.txt'  (last digit incremented, with carry.)
%      'a9to99.txt'     'a0to00.txt'  (wrap around)
%      'cat.txt'        ' '           (no digits in input name.)
%      ' '              STOP%         (error.)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILENAME, the string to be incremented.
%
%    Output, string FILENAME, the incremented string.
%
  lens = length ( filename );

  if ( lens <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILENAME_INC - Fatal error%\n' );
    fprintf ( 1, '  The input filename is empty.\n' );
    error ( 'FILENAME_INC - Fatal error%' );
  end

  change = 0;

  for i = lens : -1 : 1

    c = filename(i);

    if ( '0' <= c && c <= '8' )

      c = c + 1;

      filename(i) = c;

      return

    elseif ( c == '9' )

      change = change + 1;

      c = '0';

      filename(i) = c;

    end

  end
%
%  No digits were found.  Return blank.
%
  if ( change == 0 )
    filename = ' ';
  end

  return
end
function p = prime ( n )

%*****************************************************************************80
%
%% PRIME returns returns any of the first PRIME_MAX prime numbers.
%
%  Discussion:
%
%    PRIME_MAX is 1600, and the largest prime stored is 13499.
%
%    Thanks to Bart Vandewoestyne for pointing out a typo, 18 February 2005.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964, pages 870-873.
%
%    Daniel Zwillinger,
%    CRC Standard Mathematical Tables and Formulae,
%    30th Edition,
%    CRC Press, 1996, pages 95-98.
%
%  Parameters:
%
%    Input, integer N, the index of the desired prime number.
%    In general, is should be true that 0 <= N <= PRIME_MAX.
%    N = -1 returns PRIME_MAX, the index of the largest prime available.
%    N = 0 is legal, returning PRIME = 1.
%
%    Output, integer P, the N-th prime.  If N is out of range, P
%    is returned as -1.
%
  prime_max = 1600;

  prime_vector(1:1600) = [
        2,    3,    5,    7,   11,   13,   17,   19,   23,   29, ...
       31,   37,   41,   43,   47,   53,   59,   61,   67,   71, ...
       73,   79,   83,   89,   97,  101,  103,  107,  109,  113, ...
      127,  131,  137,  139,  149,  151,  157,  163,  167,  173, ...
      179,  181,  191,  193,  197,  199,  211,  223,  227,  229, ...
      233,  239,  241,  251,  257,  263,  269,  271,  277,  281, ...
      283,  293,  307,  311,  313,  317,  331,  337,  347,  349, ...
      353,  359,  367,  373,  379,  383,  389,  397,  401,  409, ...
      419,  421,  431,  433,  439,  443,  449,  457,  461,  463, ...
      467,  479,  487,  491,  499,  503,  509,  521,  523,  541, ...
      547,  557,  563,  569,  571,  577,  587,  593,  599,  601, ...
      607,  613,  617,  619,  631,  641,  643,  647,  653,  659, ...
      661,  673,  677,  683,  691,  701,  709,  719,  727,  733, ...
      739,  743,  751,  757,  761,  769,  773,  787,  797,  809, ...
      811,  821,  823,  827,  829,  839,  853,  857,  859,  863, ...
      877,  881,  883,  887,  907,  911,  919,  929,  937,  941, ...
      947,  953,  967,  971,  977,  983,  991,  997, 1009, 1013, ...
     1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, ...
     1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, ...
     1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, ...
     1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, ...
     1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, ...
     1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, ...
     1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, ...
     1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, ...
     1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, ...
     1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, ...
     1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, ...
     1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, ...
     1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, ...
     1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, ...
     2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, ...
     2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, ...
     2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, ...
     2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, ...
     2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, ...
     2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, ...
     2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, ...
     2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, ...
     2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, ...
     2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, ...
     2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, ...
     2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, ...
     3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, ...
     3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, ...
     3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, ...
     3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, ...
     3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, ...
     3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, ...
     3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, ...
     3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, ...
     3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, ...
     3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, ...
     3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, ...
     3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, ...
     4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, ...
     4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, ...
     4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, ...
     4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, ...
     4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, ...
     4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, ...
     4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, ...
     4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657, ...
     4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, ...
     4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, ...
     4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, ...
     4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, ...
     5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, ...
     5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, ...
     5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279, ...
     5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, ...
     5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, ...
     5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, ...
     5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639, ...
     5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, ...
     5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, ...
     5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, ...
     5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, ...
     5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, ...
     6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133, ...
     6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, ...
     6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, ...
     6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, ...
     6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, ...
     6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, ...
     6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673, ...
     6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761, ...
     6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, ...
     6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, ...
     6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997, ...
     7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, ...
     7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, ...
     7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, ...
     7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, ...
     7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499, ...
     7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, ...
     7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, ...
     7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723, ...
     7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829, ...
     7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919, ...
     7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017, ...
     8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, ...
     8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191, 8209, 8219, ...
     8221, 8231, 8233, 8237, 8243, 8263, 8269, 8273, 8287, 8291, ...
     8293, 8297, 8311, 8317, 8329, 8353, 8363, 8369, 8377, 8387, ...
     8389, 8419, 8423, 8429, 8431, 8443, 8447, 8461, 8467, 8501, ...
     8513, 8521, 8527, 8537, 8539, 8543, 8563, 8573, 8581, 8597, ...
     8599, 8609, 8623, 8627, 8629, 8641, 8647, 8663, 8669, 8677, ...
     8681, 8689, 8693, 8699, 8707, 8713, 8719, 8731, 8737, 8741, ...
     8747, 8753, 8761, 8779, 8783, 8803, 8807, 8819, 8821, 8831, ...
     8837, 8839, 8849, 8861, 8863, 8867, 8887, 8893, 8923, 8929, ...
     8933, 8941, 8951, 8963, 8969, 8971, 8999, 9001, 9007, 9011, ...
     9013, 9029, 9041, 9043, 9049, 9059, 9067, 9091, 9103, 9109, ...
     9127, 9133, 9137, 9151, 9157, 9161, 9173, 9181, 9187, 9199, ...
     9203, 9209, 9221, 9227, 9239, 9241, 9257, 9277, 9281, 9283, ...
     9293, 9311, 9319, 9323, 9337, 9341, 9343, 9349, 9371, 9377, ...
     9391, 9397, 9403, 9413, 9419, 9421, 9431, 9433, 9437, 9439, ...
     9461, 9463, 9467, 9473, 9479, 9491, 9497, 9511, 9521, 9533, ...
     9539, 9547, 9551, 9587, 9601, 9613, 9619, 9623, 9629, 9631, ...
     9643, 9649, 9661, 9677, 9679, 9689, 9697, 9719, 9721, 9733, ...
     9739, 9743, 9749, 9767, 9769, 9781, 9787, 9791, 9803, 9811, ...
     9817, 9829, 9833, 9839, 9851, 9857, 9859, 9871, 9883, 9887, ...
     9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973,10007, ...
    10009,10037,10039,10061,10067,10069,10079,10091,10093,10099, ...
    10103,10111,10133,10139,10141,10151,10159,10163,10169,10177, ...
    10181,10193,10211,10223,10243,10247,10253,10259,10267,10271, ...
    10273,10289,10301,10303,10313,10321,10331,10333,10337,10343, ...
    10357,10369,10391,10399,10427,10429,10433,10453,10457,10459, ...
    10463,10477,10487,10499,10501,10513,10529,10531,10559,10567, ...
    10589,10597,10601,10607,10613,10627,10631,10639,10651,10657, ...
    10663,10667,10687,10691,10709,10711,10723,10729,10733,10739, ...
    10753,10771,10781,10789,10799,10831,10837,10847,10853,10859, ...
    10861,10867,10883,10889,10891,10903,10909,10937,10939,10949, ...
    10957,10973,10979,10987,10993,11003,11027,11047,11057,11059, ...
    11069,11071,11083,11087,11093,11113,11117,11119,11131,11149, ...
    11159,11161,11171,11173,11177,11197,11213,11239,11243,11251, ...
    11257,11261,11273,11279,11287,11299,11311,11317,11321,11329, ...
    11351,11353,11369,11383,11393,11399,11411,11423,11437,11443, ...
    11447,11467,11471,11483,11489,11491,11497,11503,11519,11527, ...
    11549,11551,11579,11587,11593,11597,11617,11621,11633,11657, ...
    11677,11681,11689,11699,11701,11717,11719,11731,11743,11777, ...
    11779,11783,11789,11801,11807,11813,11821,11827,11831,11833, ...
    11839,11863,11867,11887,11897,11903,11909,11923,11927,11933, ...
    11939,11941,11953,11959,11969,11971,11981,11987,12007,12011, ...
    12037,12041,12043,12049,12071,12073,12097,12101,12107,12109, ...
    12113,12119,12143,12149,12157,12161,12163,12197,12203,12211, ...
    12227,12239,12241,12251,12253,12263,12269,12277,12281,12289, ...
    12301,12323,12329,12343,12347,12373,12377,12379,12391,12401, ...
    12409,12413,12421,12433,12437,12451,12457,12473,12479,12487, ...
    12491,12497,12503,12511,12517,12527,12539,12541,12547,12553, ...
    12569,12577,12583,12589,12601,12611,12613,12619,12637,12641, ...
    12647,12653,12659,12671,12689,12697,12703,12713,12721,12739, ...
    12743,12757,12763,12781,12791,12799,12809,12821,12823,12829, ...
    12841,12853,12889,12893,12899,12907,12911,12917,12919,12923, ...
    12941,12953,12959,12967,12973,12979,12983,13001,13003,13007, ...
    13009,13033,13037,13043,13049,13063,13093,13099,13103,13109, ...
    13121,13127,13147,13151,13159,13163,13171,13177,13183,13187, ...
    13217,13219,13229,13241,13249,13259,13267,13291,13297,13309, ...
    13313,13327,13331,13337,13339,13367,13381,13397,13399,13411, ...
    13417,13421,13441,13451,13457,13463,13469,13477,13487,13499 ];

  if ( n == -1 )
    p = prime_max;
  elseif ( n == 0 )
    p = 1;
  elseif ( n <= prime_max )
    p = prime_vector(n);
  else
    p = -1;
  end

  return
end
function [ r, seed ] = r8vec_uniform_01 ( n, seed )

%*****************************************************************************80
%
%% R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real R(N), the vector of pseudorandom values.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  r = zeros ( n, 1 );

  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_UNIFORM_01 - Fatal error%\n' );
    fprintf ( 1, '  Input SEED = 0%\n' );
    error ( 'R8VEC_UNIFORM_01 - Fatal error%' );
  end

  for i = 1 : n

    k = floor ( seed / 127773 );

    seed = 16807 * ( seed - k * 127773 ) - k * 2836;

    if ( seed < 0 )
      seed = seed + 2147483647;
    end

    r(i) = seed * 4.656612875E-10;

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
function h = conv(f, g, varargin)
%CONV   Convolution of CHEBFUN objects.
%   H = CONV(F, G) produces the convolution of CHEBFUN objects F and G:
%                     -
%                    /
%           H(x) =   |    F(t) G(x-t) dt,  x in [a + c, b + d]
%                    /
%                   -
%   where domain(F) is [a, b] and domain(G) is [c, d]. The integral is taken
%   over all t for which the integrand is defined: max(a, x - d) <= t <= min(b,
%   x - c).  The breakpoints of H are all pairwise sums of the breakpoints of F
%   and G.
%
%   H = CONV(F, G, 'same') will truncate the domain of H so that it is the same
%   as F. This is useful when F and G represent rapidly decaying functions on
%   large but finite intervals which are used to approximate infinity.
%
%   If F and G are simple, in the sense that their FUNS are CHEBTECH objects, a
%   fast algorithm due to Hale and Townsend is used [1]. Otherwise, the integral
%   is computed by brute force. CONV(F, G, 'old') forces the brute force
%   approach, even when the fast algorithm may be used.
%
%   Note that CONV only supports piecewise-smooth functions on bounded domains.
%
% Examples:
%     cheb.x; f = 0.8 - abs(x-0.2);
%     phi = @(t) chebfun(@(x) exp(-x^2/(4*t))/sqrt(4*pi*t));
%     fsmooth = conv(f,phi(1e-2),'same');
%     plot(f,'b',fsmooth','r')
%
%     f = chebfun(1/2); g = f;
%     subplot(2, 2, 1), plot(f)
%     for j = 2:4, g = conv(f, g); subplot(2, 2, j), plot(g), end
%     figure, for j = 1:4, subplot(2,2,j), plot(g), g = diff(g); end
%
% Reference:
%   [1] N. Hale and A. Townsend, "An algorithm for the convolution of Legendre
%   series", SIAM Journal on Scientific Computing, Vol. 36, No. 3,
%   pp. A1207-A1220, 2014.

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

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

% Check transpose state:
if ( xor(f(1).isTransposed, g(1).isTransposed) )
    error('CHEBFUN:CHEBFUN:conv:transposed', ...
        'CHEBFUN dimensions do not agree.');
end
transState = f(1).isTransposed;

% Support for quasimatrices:
nf = numColumns(f);
ng = numColumns(g);
if ( nf > 1 || ng > 1 )
    if ( nf == ng )
        f = mat2cell(f);
        g = mat2cell(g);
        h = f;
        for k = 1:nf
            h{k} = conv(f{k}, g{k}, varargin{:});
        end
    elseif ( nf == 1 )
        g = mat2cell(g);
        h = g;
        for k = 1:ng
            h{k} = conv(f, g{k}, varargin{:});
        end
    elseif ( ng == 1 )
        f = mat2cell(f);
        h = f;
        for k = 1:nf
            h{k} = conv(f{k}, g, varargin{:});
        end
    else
        error('CHEBFUN:CHEBFUN:conv:dimagree', 'CHEBFUN dimensions must agree.');
    end
    if ( ~transState )
        h = horzcat(h{:});
    else
        h = vertcat(h{:});
    end
    return
end

% Return empty for an empty input:
if ( isempty(f) || isempty(g) )
    h = chebfun();
    return
end

% Parse inputs:
oldMethod = false;
same = false;
for k = 1:numel(varargin)
    vk = varargin{k};
    if ( strcmpi(vk, 'old') )
        oldMethod = true;
    elseif ( strcmpi(vk, 'same') )
        same = true;
    elseif ( strcmpi(vk, 'full') )
        % Do nothing.
    elseif ( strcmpi(vk, 'valid') )
        % TODO: Supoprt 'valid'. Presumably where domains of f and g overlap?
        error('CHEBFUN:CHEBFUN:conv:validFlag', '''valid'' is not yet supprted.');
    else
        error('CHEBFUN:CHEBFUN:conv:badInput', 'Unknown input option %s.', vk);
    end
end

% Return a warning if F and G have too many pieces (the computation is probably
% going to be very slow):
if ( ( numel(f.funs) + numel(g.funs) ) > 50 )
    % Give a warning and proceed.
   warning('CHEBFUN:CHEBFUN:conv:piecewise',...
       ['Convolving CHEBFUNs with many pieces can be very slow.\n', ...
        'Try calling MERGE() on the inputs before calling CONV().']);
end

% Extract the domain:
[a, b] = domain(f);
[c, d] = domain(g);

% No support for unbounded domains:
if ( any(isinf([a b c d])) )
    error('CHEBFUN:CHEBFUN:conv:bounded', ...
        'CONV only supports CHEBFUN objects on bounded domains.');
end

if ( oldMethod || issing(f) || issing(g) )
    % Call the old (and slow) version of CONV if we are not based on CHEBTECHS.
    h = oldConv(f, g);

else

    % Ensure g is the signal (i.e., on the larger domain) and f is the filter:
    if ( (b - a) > (d - c) )
        h = conv(g, f);
        return
    end

    % Initialize the output:
    h = chebfun(0, [a + c, b + d]);
    % Deal with piecewise CHEBFUN objects by looping over each interaction:
    for j = 1:numel(f.funs)
        for k = 1:numel(g.funs)
            % Compute the contribution of jth fun of f with kth fun of g:
            hjk = conv(f.funs{j}, g.funs{k});
            % Add this contribution:
            h = myplus(h, chebfun(hjk));
        end
    end

    % Make sure that point values are not added twice:
    dom = domain(h);
    intDom = dom(2:end-1);
    h.pointValues(2:end-1) = 1/2*(feval(h, intDom.', 'left') + ...
        feval(h, intDom.', 'right'));

end

% Truncate:
if ( same )
    h = restrict(h, [a, b]);
end

% Transpose:
if ( transState )
    h = h.';
end

% Simplify the result.
h = simplify(h);

end

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

function h = myplus(f, g)
% Modified PLUS() which pads with zeros to fulfil domain requirements.
%  Note f is always on the largest possible domain. g is on a subdomain of f

% Tidy the domains:
[f, g] = tweakDomain(f, g);
[c, d] = domain(g);

fTmp = restrict(f, [c, d]); % f{c, d}
hTmp = fTmp + g;            % h{c, d} = f{c, d} + g;
h = defineInterval(f, [c, d], hTmp);

end

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

function h = oldConv(f, g)
% The old convolution algorithm based on quadrature. See section 2 of [1] for a
% description of the algorithm.

% Find all breakpoints in the convolution:
[A, B] = meshgrid(f.domain, g.domain);
dom = unique(A(:) + B(:)).';

% Coalesce breaks that are close due to roundoff:
dom(diff(dom) < 10*eps*max(abs(dom([1,end])))) = [];
dom(isnan(dom)) = [];

% Combine vertical and horizontal scales:
hs = max(hscale(f), hscale(g));
vs = 2*max([vscale(f), vscale(g)]);

% Set preferences:
%
% TODO:  CHEBFUN is not supposed to set the refinementFunction preference
% because it doesn't belong to the list of "abstract" preferences required of
% all techs.  Do we really need to alter it here?
p = chebfunpref();
p.splitting = false;
p.blowup = false;
p.techPrefs.extrapolate = true;
p.techPrefs.refinementFunction = 'nested';
p.techPrefs.sampleTest = false;

% Construct FUNS:
funs = cell(1, length(dom)-1);
for k = 1:length(dom)-1
    data.domain = dom(k:k+1);
    data.vscale = vs;
    data.hscale = hs;
    newFun = bndfun(@(x) convIntegral(x, f, g), data, p);
    vs = max(get(newFun, 'vscale'), vs);
    funs{k} = newFun;
end

% Construct CHEBFUN:
h = chebfun(funs);
h.isTransposed = f.isTransposed;

end

function out = convIntegral(x, f, g)
%CONVINTEGRAL   Evaluate convolution integral.
%   Y = CONVINTEGRAL(X, F, G) evaluates the convolution of the CHEBFUNs F and G
%   at the points X.

a = f.domain(1);
b = f.domain(end);
c = g.domain(1);
d = g.domain(end);

out = 0*x;
for k = 1:length(x)
    A = max(a, x(k) - d);
    B = min(b, x(k) - c);
    if ( A < B )
        ends = union(x(k) - g.domain, f.domain);
        dom = [A, ends((A < ends) & (ends < B)), B];
        for j = 1:length(dom)-1
            % INTEGRAL is not available in versions of MATLAB prior to R2012a,
            % so if we're running on an older version, fall back to QUADGK.
            integrand = @(t) feval(f, t).*feval(g, x(k) - t);
            if ( verLessThan('matlab', '7.14') )
                out(k) = out(k) + quadgk(integrand, dom(j), dom(j+1), ...
                    'AbsTol', 1e-15, 'RelTol', 100*eps);
            else
                out(k) = out(k) + integral(integrand, dom(j), dom(j+1), ...
                    'AbsTol', 1e-15, 'RelTol', 1e-15);
            end
        end
    end
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function str_out = structToString(data, show, level, str)
% Display the content of a structure

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       Andrea Gatti
%  Contributors:     Andrea Gatti
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

    if (nargin  < 2)
        show = true;
    end
    if (nargin  < 3)
        level = 0;
    end
    if (nargin  < 4)
        str = '';
    end

    key = fieldnames(data);
    for k = 1 : numel(key)
        if isstruct(data.(key{k}))
            str = sprintf('%s%*c -> %s\n', str, 4 * level, char(32), key{k});
            str = structToString(data.(key{k}), level+1, false, str);
        else
            tmp = toIniString(key{k}, data.(key{k}));
            str = sprintf('%s%*c -> %s\n', str, 4 * level, char(32), tmp{1});

        end
    end

    if show
        fprintf(str);
    end

    if (nargout == 1)
        str_out = str;
    end

end

% toIniString -----------------------------------------------------
function cell_str = toIniString(variable_name, value, format, cell_str)
% Convert any variable to ini string format
% SYNTAX:
%   cell_str = toIniString(variable_name, value)
%   cell_str = toIniString(variable_name, value, format)
%   cell_str = toIniString(variable_name, value, cell_str)
%   cell_str = toIniString(variable_name, value, format, cell_str)
switch nargin
    case 1
        error('Error in Ini_Manager.toIniString, too few parameters');
    case 2

        format = '';
        cell_str = {};
    case 3
        if iscellstr(format)
            cell_str = format;
            format = '%g';
        else
            cell_str = {};
        end
    case 4
        % ok
end
if ~isdeployed
    variable_name = ['<strong>' variable_name '</strong>'];
end
if ischar(value) % is string
    cell_str{numel(cell_str) + 1} = [variable_name ' = "' value(:)' '"'];
elseif isnumeric(value)
    if isempty(format)
        format = '%g';
    else
        format = strtrim(format);
    end
    if numel(value) > 1 % it's an array of values
        format = [format ' '];
        tmp = sprintf(format,value);
        cell_str{numel(cell_str) + 1} = [variable_name ' = [' tmp(1:end-1) ']'];
        clear tmp;
    else
        cell_str{numel(cell_str) + 1} = [variable_name ' = ' sprintf(format,value)];
    end
else % generic converter (may not work properly)
    toString = @(var) strtrim(regexprep(evalc('disp(var)'), '\n', ''));
    if iscell(value)
        if ~isempty(value) && ischar(value{1})
            cell_str{numel(cell_str) + 1} = [variable_name ' = [' Ini_Manager.strCell2Str(value) ']'];
        else
            cell_str{numel(cell_str) + 1} = [variable_name ' = [' toString(value) ']'];
        end
    else
        cell_str{numel(cell_str) + 1} = [variable_name ' = ' toString(value)];
    end
end

% I want a column array
if size(cell_str,1) < size(cell_str,2)
    cell_str = cell_str';
end
end

author: goGPS-Project


text:
% Fig. 6.47   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%

clear all;
close all;

num=[20 1];
den=[1 0 0];
w=logspace(-3,1,100);
[m,p]=bode(num,den,w);
loglog(w,m);
grid;
xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Fig. 6.47 Compensated open-loop transfer function.');