text:
function [y2,ntrimmed,spikes, yy] = splinetrim(y,varargin)
% :Usage:
% ::
%
%     function [y,ntrimmed,spikes, yfit] = splinetrim(y,[iqrmult],[knotrate],[X],['p'])
%
% Uses a robust measure of deviations in a timeseries gradient
% to find high-velocity 'spikes', presumed to be artifacts
%
% Uses spline interpolation to replace spikes with reasonable values.
%
% :Input:
%
%   **y:**
%        a timeseries
%
% :Optional Inputs:
%
%   **iqrmult:**
%        how many times the interquartile range above which velocities are
%           outliers, default is 1.5
%
%   **knotrate:**
%        sets knot points every k observations, default is 3
%
%   **X:**
%        matrix of session means or other linear regressors to remove
%
%   **p:**
%        plot the results
%
%   X and 'p' can be entered in any order, but after iqrmult and knotrate
%
% :Examples:
% ::
%
%    [y2,nt] = splinetrim(trialdat,3,5,'p'); nt
%
% ..
%    3/29/05, Tor Wager
% ..

% -------------------------------------------------
% * setup
% -------------------------------------------------

iqrmult = 3;  % how many times the interquartile range above which pts are outliers
knotrate = 5;   % in images
X = [];
plotme = 0;

if length(varargin) > 0, iqrmult = varargin{1};, end
if length(varargin) > 1, knotrate = varargin{2};, end

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}
        case 'p', plotme = 1;
        otherwise error('Unrecognized argument: acceptable is ''p'' for plot.')
        end
    else
        if ndims(varargin{i}) == 2 & sum(size(varargin{i}))>2
            X = varargin{i};
        end
    end
end
% filter y using X matrix; yf is residuals
if ~isempty(X),
    mfit = X * pinv(X) * y;
    y = y - mfit;
end

% ---------------------------------------------------------
% find spike regions
% ---------------------------------------------------------
% EXCLUDE FROM CONSIDERING AS KNOT POINTS BASED ON VELOCITY
veloc = gradient(y); veloc(end) = 0;    % clamp last value to 0
tmp = abs(veloc - median(veloc));

% threshold is mad + 1.5 times the interquartile range of abs. deviations
thr = median(tmp) + iqrmult * (prctile(tmp,75) - prctile(tmp,25));

spikes = tmp>thr;

% EXCLUDE FROM CONSIDERING AS KNOT POINTS BASED ON OUTLIER STATUS
d = abs(y - median(y));    % distance from mean
thr = median(d) + iqrmult * (prctile(d,75) - prctile(d,25));
sptmp = d>thr;
spikes = spikes + sptmp;

% smooth this some, so that low-deviation regions in the middle of spikes
% do not get counted
spikes2 = smooth_timeseries(spikes,4);  % excluded from being knot points

spikes = find(spikes>0);        % which points to interpolate in the end
spikes2 = find(spikes2>0);      % larger set of pts to not include as knot pts

% -------------------------------------------------
% * get spline fit
% -------------------------------------------------
bp = zeros(size(y));
bp(1:knotrate:length(y)) = 1;
bp(isnan(y)) = 0;                % ignore NaNs
%bp(end) = 1;                     % clamp last point to be an endpt
bp(spikes2) = 0;                 % do not put knot points on spikes

bp = find(bp);

nbp = length(bp);
bpy = zeros(nbp, 1);

% figure out medians for each segment -
% this is the interp knot point.
ytmp = y;
ytmp(spikes2) = NaN;  % get rid of spikes

for i = 1:nbp

    if i == 1
        st = 1;
    else
        st = bp(i) - round((bp(i) - bp(i-1)) ./ 2);
    end

    if i == nbp
        en = length(y);
    else
        en = bp(i) + round((bp(i+1) - bp(i)) ./ 2);
    end
    data = ytmp(st:en);
    bpy(i) = median(data);

end

% xx = 1:length(y);
yy = spline(bp, bpy, 1:length(y));
if ~iscol(yy), yy = yy'; end

% -------------------------------------------------
% * find pts that are really far from spline fit
% -------------------------------------------------
d = abs(y - yy);    % distance from spline fit
thr = median(d) + iqrmult * (prctile(d,75) - prctile(d,25));

spikes3 = d>thr;

spikes = unique([spikes; find(spikes3>0)]);        % final answer
% replace y data with spline fits where spikes occur
y2 = y;
y2(spikes) = yy(spikes);

ntrimmed = length(spikes);

if plotme
    figure; hold on;
    plot(y,'k','LineWidth',2);
    plot(y2,'g');
    plot(yy,'b');
    plot(bp, yy(bp), 'b.','MarkerSize', 10);
    legend({'Original' 'Adjusted' 'Spline fit'})
    plot(spikes,y2(spikes),'ro','MarkerFaceColor','r','MarkerSize',6);
end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function x = uminus( x )
    %UMINUS Unary minus.
    %   X = UMINUS(X) returns the negated TTeMPS block-mu tensor X. Only the supercore
    %   is touched.
    %
    %   See also UPLUS.

    %   TTeMPS Toolbox.
    %   Michael Steinlechner, 2013-2016
    %   Questions and contact: michael.steinlechner@epfl.ch
    %   BSD 2-clause license, see LICENSE.txt

    x.U{x.mu} = -x.U{x.mu};
end

author: NicolasBoumal

repo: manopt


text:
function varargout = not_calculated(varargin)
% NOT_CALCULATED M-file for not_calculated.fig
%      NOT_CALCULATED by itself, creates a new NOT_CALCULATED or raises the
%      existing singleton*.
%
%      H = NOT_CALCULATED returns the handle to a new NOT_CALCULATED or the handle to
%      the existing singleton*.
%
%      NOT_CALCULATED('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in NOT_CALCULATED.M with the given input arguments.
%
%      NOT_CALCULATED('Property','Value',...) creates a new NOT_CALCULATED or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before not_calculated_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to not_calculated_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

% This file is part of the Matlab Toolbox for Dimensionality Reduction v0.7.2b.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, 2010
% University California, San Diego / Delft University of Technology

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

% Last Modified by GUIDE v2.5 30-Jul-2008 11:02:27

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @not_calculated_OpeningFcn, ...
                   'gui_OutputFcn',  @not_calculated_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 not_calculated is made visible.
function not_calculated_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 not_calculated (see VARARGIN)

% Choose default command line output for not_calculated
handles.output = 'Yes';

% Update handles structure
guidata(hObject, handles);

% Insert custom Title and Text if specified by the user
% Hint: when choosing keywords, be sure they are not easily confused
% with existing figure properties.  See the output of set(figure) for
% a list of figure properties.
if(nargin > 3)
    for index = 1:2:(nargin-3),
        if nargin-3==index, break, end
        switch lower(varargin{index})
         case 'title'
          set(hObject, 'Name', varargin{index+1});
         case 'string'
          set(handles.text1, 'String', varargin{index+1});
        end
    end
end

% Determine the position of the dialog - centered on the callback figure
% if available, else, centered on the screen
FigPos=get(0,'DefaultFigurePosition');
OldUnits = get(hObject, 'Units');
set(hObject, 'Units', 'pixels');
OldPos = get(hObject,'Position');
FigWidth = OldPos(3);
FigHeight = OldPos(4);
if isempty(gcbf)
    ScreenUnits=get(0,'Units');
    set(0,'Units','pixels');
    ScreenSize=get(0,'ScreenSize');
    set(0,'Units',ScreenUnits);

    FigPos(1)=1/2*(ScreenSize(3)-FigWidth);
    FigPos(2)=2/3*(ScreenSize(4)-FigHeight);
else
    GCBFOldUnits = get(gcbf,'Units');
    set(gcbf,'Units','pixels');
    GCBFPos = get(gcbf,'Position');
    set(gcbf,'Units',GCBFOldUnits);
    FigPos(1:2) = [(GCBFPos(1) + GCBFPos(3) / 2) - FigWidth / 2, ...
                   (GCBFPos(2) + GCBFPos(4) / 2) - FigHeight / 2];
end
FigPos(3:4)=[FigWidth FigHeight];
set(hObject, 'Position', FigPos);
set(hObject, 'Units', OldUnits);

% Show a question icon from dialogicons.mat - variables questIconData
% and questIconMap
load dialogicons.mat

IconData=questIconData;
questIconMap(256,:) = get(handles.figure1, 'Color');
IconCMap=questIconMap;

Img=image(IconData, 'Parent', handles.axes1);
set(handles.figure1, 'Colormap', IconCMap);

set(handles.axes1, ...
    'Visible', 'off', ...
    'YDir'   , 'reverse'       , ...
    'XLim'   , get(Img,'XData'), ...
    'YLim'   , get(Img,'YData')  ...
    );

% Make the GUI modal
set(handles.figure1,'WindowStyle','modal')

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

% --- Outputs from this function are returned to the command line.
function varargout = not_calculated_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;

% The figure can be deleted now
delete(handles.figure1);

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

handles.output = get(hObject,'String');

% Update handles structure
guidata(hObject, handles);

% Use UIRESUME instead of delete because the OutputFcn needs
% to get the updated handles structure.
uiresume(handles.figure1);

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

handles.output = get(hObject,'String');

% Update handles structure
guidata(hObject, handles);

% Use UIRESUME instead of delete because the OutputFcn needs
% to get the updated handles structure.
uiresume(handles.figure1);

% --- Executes when user attempts to close figure1.
function figure1_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isequal(get(handles.figure1, 'waitstatus'), 'waiting')
    % The GUI is still in UIWAIT, us UIRESUME
    uiresume(handles.figure1);
else
    % The GUI is no longer waiting, just close it
    delete(handles.figure1);
end

% --- Executes on key press over figure1 with no controls selected.
function figure1_KeyPressFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Check for "enter" or "escape"
if isequal(get(hObject,'CurrentKey'),'escape')
    % User said no by hitting escape
    handles.output = 'No';

    % Update handles structure
    guidata(hObject, handles);

    uiresume(handles.figure1);
end

if isequal(get(hObject,'CurrentKey'),'return')
    uiresume(handles.figure1);
end

repo: cluster


text:
clear all; close all;
I=imread('coins.png');
I=im2double(I);
V=zeros(size(I));
for i=1:size(V, 1)
    V(i,:)=0.02*i/size(V,1);
end
J=imnoise(I, 'localvar', V);
figure;
subplot(121);  imshow(I);
subplot(122);  imshow(J);

author: HuangCongQing


text:
%     MOD1(x,y) is like MOD(x,y), except that it returns a value between 1
%     and y, instead of between 0 and y-1.

function v = mod1(x,y)

v = mod(x-1,y)+1;

author: kristinbranson

repo: JAABA


text:
function pass = test_plotting()
% Check that the very basic plotting commands do not crash.

% Real-valued functions.
f1 = chebfun3(@(x,y,z) x.*y.*z, [-1 2 -1 2 -1 2]);

% A complex-valued function.
f2 = chebfun3(@(x,y,z) exp(cos(1i*x.*y.*z)));

% Obviously, we can't check if the plots are correct without human
% intervention, so all these tests are meant to do is make sure none of the
% plotting functions crash.

hfig = figure('Visible', 'off');

%% Real-valued function
% Non-GUI plots
pass(1) = doesNotCrash(@() plot(f1));
pass(2) = doesNotCrash(@() slice(f1, 'noslider'));
pass(3) = doesNotCrash(@() slice(f1, 0.5, -0.3, 0.9));
pass(4) = doesNotCrash(@() isosurface(f1, 'noslider'));
pass(5) = doesNotCrash(@() isosurface(f1, [0.5, -0.6]));
pass(6) = doesNotCrash(@() scan(f1));
pass(6) = doesNotCrash(@() scan(f1, 1));
pass(7) = doesNotCrash(@() scan(f1, 1, 'hold'));
pass(8) = doesNotCrash(@() scan(f1, 3));
% GUI-based plots:
pass(9) = doesNotCrash(@() isosurface(f1));
pass(10) = doesNotCrash(@() slice(f1));
pass(11) = doesNotCrash(@() surf(f1));

%% Complex-valued function
% Non-GUI plots
pass(12) = doesNotCrash(@() plot(f2));
pass(13) = doesNotCrash(@() slice(f2, 'noslider'));
pass(14) = doesNotCrash(@() slice(f2, 0.5, -0.3, 0.9));
pass(15) = doesNotCrash(@() isosurface(f2, 'noslider'));
pass(16) = doesNotCrash(@() isosurface(f2, [0.5, -0.6]));
pass(17) = doesNotCrash(@() scan(f2));
pass(18) = doesNotCrash(@() scan(f2, 'hold'));
pass(19) = doesNotCrash(@() scan(f2, 2, 'hold'));
% GUI-based plots:
pass(20) = doesNotCrash(@() isosurface(f2));
pass(21) = doesNotCrash(@() slice(f2));
pass(22) = doesNotCrash(@() surf(f2));

%% Test plotting on a non-standard domain
r = chebfun3(@(r,t,p) r, [0 1 0 2*pi pi/4 pi/2]);
t = chebfun3(@(r,t,p) t, [0 1 0 2*pi pi/4 pi/2]);
p = chebfun3(@(r,t,p) p, [0 1 0 2*pi pi/4 pi/2]);
x = r.*cos(t).*cos(p);
y = r.*sin(t).*cos(p);
z = r.*sin(p);
density = sin(10*t).*cos(10*r)+1;
pass(23) = doesNotCrash(@() plot(x,y,z,density));

close(hfig);

end

function pass = doesNotCrash(fn)
try
    fn();
    pass = true;
catch ME
    pass = false;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%	robot objects can be multiplied r1*r2 which is mechanically equivalent
%	to mounting robot r2 on the end of robot r1.
%

% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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.
%
% RTB 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 RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
function r2 = mtimes(r, l)
        %SerialLink.mtimes Join robots
        %
        % R = R1 * R2 is a robot object that is equivalent to mounting robot R2
        % on the end of robot R1.
            if isa(l, 'SerialLink')
                r2 = SerialLink(r);
				x = r2.links;
				for i = 1:length(x)
					new_links(i) = Link(x(i));
				end
				y = l.links;
				for i = 1:length(y)
					new_links(i + length(x)) = Link(y(i));
				end
				r2.links = new_links;
                r2.base = r.base;
                r2.n = length(r2.links);
            elseif isa(l, 'Link')
                r2 = SerialLink(r);
				x = r2.links;
				for i = 1:length(x)
					new_links(i) = Link(x(i));
				end
				for i = 1:length(l)
					new_links(i + length(x)) = Link(l(i));
				end

                r2.links = new_links;
                r2.n = length(r2.links);
            end
        end

author: petercorke

repo: robotics-toolbox-matlab


text:
% {make_libsvm} compiles Libsvm and the mex interface for Libsvm (C++).
%
%      make_libsvm
%
% Author: Stefano Melacci (2009)
%         mela@dii.unisi.it

delete *.obj
delete *.o
mex -c svm.cpp
if exist('svm.obj','file')>0
    mex mexGramSVMTrain.cpp svm.obj
elseif exist('svm.o','file')>0
    mex mexGramSVMTrain.cpp svm.o
end

author: viggin

repo: domain-adaptation-toolbox


text:
function [status, result] = lrsSpecifyRows(fileNameIn)
%read the number of rows of data in the file and replace the **** with the
%number of rows

if ~contains(fileNameIn,filesep)
    fileNameIn = [pwd filesep fileNameIn];
end

fid = fopen(fileNameIn,'r+');
if fid<0
    disp(fileNameOut)
    error('Could not open lrs output file.');
end

countRows = 0;
while 1
    tline = fgetl(fid);
    if countRows ~=0
        countRows = countRows + 1;
    end
    if strcmp(tline, 'begin')
        countRows = 1;
    elseif strcmp(tline, 'end')
        break
    elseif ~ischar(tline)
        error('Could not read lrs output file.');
    end
end
nRows = countRows -3;

%open the file and write the number of rows and columns
%The procedure to change the text in files under Linux/Unix using sed:
%Use Stream EDitor (sed) as follows:
%sed -i 's/old-text/new-text/g' input.txt
%The s is the substitute command of sed for find and replace
%It tells sed to find all occurrences of ‘old-text’ and replace with ‘new-text’ in a file named input.txt
[status, result] =  system(['sed -i ''s/\*\*\*\*\*/' int2str(nRows) '/g'' ' fileNameIn]);
if status == 1
    disp(result)
    error(['sed failed on file ' fileNameIn]);
end

end

author: opencobra

repo: cobratoolbox


text:
function [t,p1,p2] = intersectSegments(l1,l2)

% INTERSECTSEGMENTS  Intersect two 3d segments.
%   [T,P1,P2] = INTERSECTSEGMENTS(L1,L2) intersects the two segments
%   L1=[E1,F1] and L2=[E2,F2], and returns a vector of abscissas T=[T1;T2],
%   and two intersection points P1 and P2, so that
%
%       P1 = E1 + T1*V1
%       P2 = E2 + T2*V2
%
%   where V1=F1-E1 and V2=F2-E2 are director vectors.
%
%   Note that what we really intersect is the infinite line supporting the
%   segments. Some examples (intersection is at the cross '+'):
%
%       ----+----         +         |
%                                   |
%           |           / |     ----+---
%           |          /  |         |
%
%   If the lines do intersect in the 3D space, we have the intersection
%   point P = P1 = P2. Otherwise, P1 and P2 are the points of the lines
%   that are closest to the other line, that is, the line passing over P1
%   and P2 is orthogonal to both L1 and L2.
%
%   See also INTERSECTPVLINES, SEG2PVLIN.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

m1 = seg2pvLin(l1);
m2 = seg2pvLin(l2);

[t,p1,p2] = intersectPvLines(m1,m2);

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function isargstruct(varargin)
%ISARGSTRUCT tests if the given arg is a struct and returns an error otherwise
%
%   Usage: isargstruct(arg1,arg2,...)
%
%   Input options:
%       args        - list of args
%
%   ISARGSTRUCT(args) tests if all given args are a struct and returns
%   an error otherwise.
%
%   see also: isargchar

%*****************************************************************************
% Copyright (c) 2010-2012 Quality & Usability Lab                            *
%                         Deutsche Telekom Laboratories, TU Berlin           *
%                         Ernst-Reuter-Platz 7, 10587 Berlin, Germany        *
%                                                                            *
% This file is part of the Sound Field Synthesis-Toolbox (SFS).              *
%                                                                            *
% The SFS 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.                                                         *
%                                                                            *
% The SFS 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/>.            *
%                                                                            *
% The SFS is a toolbox for Matlab/Octave to  simulate and  investigate sound *
% field  synthesis  methods  like  wave  field  synthesis  or  higher  order *
% ambisonics.                                                                *
%                                                                            *
% http://dev.qu.tu-berlin.de/projects/sfs-toolbox       sfstoolbox@gmail.com *
%*****************************************************************************

% AUTHOR: Hagen Wierstorf
% $LastChangedDate: 2012-04-26 10:12:51 +0200 (Thu, 26 Apr 2012) $
% $LastChangedRevision: 710 $
% $LastChangedBy: wierstorf.hagen $

%% ===== Checking for struct =============================================
for ii = 1:nargin
    if ~isstruct(varargin{ii})
        error('%s need to be a struct.',inputname(ii));
    end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function g = mult_givens ( c, s, k, g )

%*****************************************************************************80
%
%% MULT_GIVENS applies a Givens rotation to two successive entries of a vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 March 2008
%
%  Author:
%
%    C original version by Lili Ju
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Richard Barrett, Michael Berry, Tony Chan, James Demmel,
%    June Donato, Jack Dongarra, Victor Eijkhout, Roidan Pozo,
%    Charles Romine, Henk van der Vorst,
%    Templates for the Solution of Linear Systems:
%    Building Blocks for Iterative Methods,
%    SIAM, 1994.
%    ISBN: 0898714710,
%    LC: QA297.8.T45.
%
%    Tim Kelley,
%    Iterative Methods for Linear and Nonlinear Equations,
%    SIAM, 2004,
%    ISBN: 0898713528,
%    LC: QA297.8.K45.
%
%    Yousef Saad,
%    Iterative Methods for Sparse Linear Systems,
%    Second Edition,
%    SIAM, 2003,
%    ISBN: 0898715342,
%    LC: QA188.S17.
%
%  Parameters:
%
%    Input, real C, S, the cosine and sine of a Givens
%    rotation.
%
%    Input, integer K, indicates the location of the first vector entry.
%
%    Input/output, real G(1:K+1), the vector to be modified.  On output,
%    the Givens rotation has been applied to entries G(K) and G(K+1).
%
  g1 = c * g(k) - s * g(k+1);
  g2 = s * g(k) + c * g(k+1);

  g(k)   = g1;
  g(k+1) = g2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function  [C] =  TWSC_ADMM( Y, D, S, W1, W2, Par )
% This routine solves the following trilateral weighted sparse coding problem
%
% min_{C,Z} |W1(Y-DSC)W2|_F,2 + |Z|_1 s.t.  C=Z
%
% inputs:
%        Y -- d*M data matrix, d is the data dimension, and M is the number
%             of image patches.
%        W1 -- d*d matrix of row weights
%        W2 -- M*M matrix of column weights
% outputs:
%        C -- d*M data matrix, sparse coding coefficient matrix
%        Z -- d*M data matrix, auxiliary variable, equal to C

tol = 1e-8;
Par.maxrho = 100;
Par.maxIter = 10;
Par.rho = 0.5;
Par.mu = 1.1;
Par.display = 0;
% Initializing optimization variables
C = zeros(size(S, 1), size(Y, 2));
Z = zeros(size(C));
U = zeros(size(C));
% Start main loop
iter = 0;
stopCZ = zeros(Par.maxIter, 1);
stopC = zeros(Par.maxIter, 1);
stopZ = zeros(Par.maxIter, 1);
while iter < Par.maxIter
    iter = iter + 1;
    Cpre = C;
    Zpre = Z;
    %% update C, fix Z and U
    % min_{C} ||W1 * (Y - DSC) * W2||_F^2 + 0.5 * rho * ||C - Z + 1/rho * U||_F^2
    % The solution is equal to solve A * X + X * B = E
    A = S' * D' * diag(W1.^2) * D * S;
    W2inv = diag(1./(W2.^2));
    B = 0.5 * Par.rho * W2inv;
    E = S' * D' * diag(W1.^2) * Y + 0.5 * (Par.rho * Z - U) * W2inv;
    C = sylvester(A, B, E);

    %     %% faster solution
    %     [Ua, Sa, ~] = svd(A);
    %     I1 = eye(size(A, 2));
    %     I2 = eye(size(B, 1));
    %     K = kron(I1, A) + kron(B', I2);
    %     invK = 1./diag(K);
    %     UTE = Ua'*E;
    %     vecUTE = UTE(:);
    %     vecUTC = invK .* vecUTE;
    %     MatvecUTC = reshape(vecUTC, [size(UTE, 1) size(UTE, 2)]);
    %     C = Ua*MatvecUTC;

    %% update Z, fix X and D
    % min_{Z} 0.5 * rho * ||Z - (C + 1/rho * U)||_F^2 + ||Z||_1
    Temp = C + U/Par.rho;
    Z = sign(Temp) .* max( abs(Temp) - 1/Par.rho, 0 );

    %% check the convergence conditions
    stopCZ(iter) = max(max(abs(C - Z)));
    stopC(iter) = max(max(abs(C - Cpre)));
    stopZ(iter) = max(max(abs(Z - Zpre)));
    if Par.display %&& (iter==1 || mod(iter,10)==0 || stopC<tol)
        disp(['iter ' num2str(iter) ', mu=' num2str(Par.mu,'%2.1e') ...
            ', max(||c-z||)=' num2str(stopCZ(iter),'%2.3e') ...
            ', max(||c-cpre||)=' num2str(stopC(iter),'%2.3e') ...
            ', max(||z-zpre||)=' num2str(stopZ(iter),'%2.3e')]);
    end
    if stopCZ(iter) < tol && stopC(iter) < tol && stopZ(iter) < tol
        break;
    else
        %% update the augmented multiplier D, fix Z and X
        U = U + Par.rho * (C - Z);
        Par.rho = min(Par.maxrho, Par.mu * Par.rho);
    end
end
return;

author: csjunxu


text:
function results = vl_test_twister(varargin)
% VL_TEST_TWISTER
vl_test_init ;

function test_illegal_args()
vl_assert_exception(@() vl_twister(-1), 'vl:invalidArgument') ;
vl_assert_exception(@() vl_twister(1, -1), 'vl:invalidArgument') ;
vl_assert_exception(@() vl_twister([1, -1]), 'vl:invalidArgument') ;

function test_seed_by_scalar()
rand('twister',1) ; a = rand ;
vl_twister('state',1) ; b = vl_twister ;
vl_assert_equal(a,b,'seed by scalar + VL_TWISTER()') ;

function test_get_set_state()
rand('twister',1) ; a = rand('twister') ;
vl_twister('state',1) ; b = vl_twister('state') ;
vl_assert_equal(a,b,'read state') ;

a(1) = a(1) + 1 ;
vl_twister('state',a) ; b = vl_twister('state') ;
vl_assert_equal(a,b,'set state') ;

function test_multi_dimensions()
b = rand('twister') ;
rand('twister',b) ;
vl_twister('state',b) ;
a=rand([1 2 3 4 5]) ;
b=vl_twister([1 2 3 4 5]) ;
vl_assert_equal(a,b,'VL_TWISTER([M N P ...])') ;

function test_multi_multi_args()
a=rand(1, 2, 3, 4, 5) ;
b=vl_twister(1, 2, 3, 4, 5) ;
vl_assert_equal(a,b,'VL_TWISTER(M, N, P, ...)') ;

function test_square()
a=rand(10) ;
b=vl_twister(10) ;
vl_assert_equal(a,b,'VL_TWISTER(N)') ;

author: Cloud-CV

repo: object-proposals


text:
function [aic, sbic] = aicsbic(errors,constant,p,q,X)
% Computes the Akaike and Schwartz/Bayes Information Criteria for an ARMA(P,Q) as parameterized in
% ARMAXFILTER
%
% USAGE:
%   [AIC] = aicsbic(ERRORS,CONSTANT,P,Q)
%   [AIC,SBIC] = aicsbic(ERRORS,CONSTANT,P,Q,X)
%
% INPUTS:
%   ERRORS   - A T by 1 length vector of errors from the regression
%   CONSTANT - Scalar variable: 1 to include a constant, 0 to exclude
%   P        - Non-negative integer vector representing the AR orders to include in the model.
%   Q        - Non-negative integer vector representing the MA orders to include in the model.
%   X        - [OPTIONAL]  a T by K  matrix of exogenous variables.
%
% OUTPUTS:
%   AIC       - The Akaike Information Criteria
%   SBIC      - The Schwartz/Bayes Information Criteria
%
% COMMENTS:
%   This is a helper for ARMAXFILTER and uses the same inputs, CONSTANT, P, Q and X.  ERRORS should
%   be the errors returned from a call to ARMAXFILTER with the same values of P, Q, etc.
%
% EXAMPLES:
%   Compute AIC and SBIC from an ARMA
%       [parameters, LL, errors] = armaxfilter(y, constant, p, q);
%       [aic,sbic] = aicsbic(errors,constant,p,q)
%
%  See also ARMAXFILTER, HETEROGENEOUSAR

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 10/19/2009

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<3 || nargin>5
    error('3 to 5 inputs required')
end
if nargin==3
    q = [];
    X=[];
elseif nargin==4
    X=[];
end
%%%%%%%%%%%%%%%
% y
%%%%%%%%%%%%%%%
if size(errors,2) > 1 || length(errors)==1
    error('ERRORS series must be a column vector.')
elseif isempty(errors)
    error('ERRORS is empty.')
end

%%%%%%%%%%%%%%%
% P
%%%%%%%%%%%%%%%
if size(p,2)>size(p,1)
    p=p';
end
if isempty(p)
    p=0;
end
if min(size(p))~=1
    error('P must be a column vector of included lags')
end
if  any(p<0) || any(floor(p)~=p)
    error('P must contain non-negative integers only')
end
if max(p)>=(length(errors)-max(p))
    error('Too many lags in the AR.  max(P)<T/2')
end
maxP=max(p);
if size(p,1)==1 && p==0
    p=[];
end
if length(unique(q))~=length(q)
    error('P must contain at most one of each lag')
end
%%%%%%%%%%%%%%%
% Q
%%%%%%%%%%%%%%%
if size(q,2)>size(q,1)
    q=q';
end
if isempty(q)
    q=0;
end
if min(size(q))~=1
    error('Q must be a column vector of included lags')
end
if  any(q<0) || any(floor(q)~=q)
    error('Q must contain non-negative integers only')
end
if max(q)>=length(errors)
    error('Too many lags in the AR.  max(Q)<T')
end
maxq=max(q);
if size(q,1)==1 && q==0
    q=[];
end
if length(unique(q))~=length(q)
    error('Q must contain at most one of each lag')
end
%%%%%%%%%%%%%%%
% Constant
%%%%%%%%%%%%%%%
if ~ismember(constant,[0 1])
    error('CONSTANT must be 0 or 1')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T = length(errors);
seregression=sqrt(errors'*errors/T);

lp = length(unique(p));
lq = length(unique(q));
K=constant+lp+lq+size(X,2);

aic = log(seregression^2) + 2*K/T;
sbic = log(seregression^2) + log(T)*K/T;

author: bashtage

repo: mfe-toolbox


text:
function [data] = spm2fieldtrip(D)

% SPM2FIELDTRIP converts an SPM8 meeg object into a FieldTrip raw data structure
%
% Use as
%   data = spm2fieldtrip(D)
% where D is the SPM meeg object which you can load in with SPM_EEG_LOAD
% and where data is a FieldTrip raw data structure as if it were returned
% by FT_PREPROCESSING.
%
% See also FT_PREPROCESSING, SPM_EEG_LOAD

if ~ft_hastoolbox('spm8up')
  % it should be version spm8 or higher, since spm99, spm2 and spm5 did not yet the "meeg" object
  ft_error('this requires SPM8 or later to be on your MATLAB path');
end

if ~isa(D, 'meeg')
  ft_error('this requires an SPM "meeg" object as input')
end

% this is how SPM8 represents it
data = D.ftraw(0);

clist      = D.condlist;
conditions = D.conditions;

data.trialinfo = zeros(D.ntrials,1);
for k = 1:numel(clist)
  fprintf('mapping condition label "%s" to condition code %d\n', clist{k}, k);
  sel=strcmp(clist{k}, conditions);
  data.trialinfo(sel) = k;
end

% FIXME the following is not correct
%
% data.sampleinfo = zeros(ntrials,2);
% for i=1:ntrials
%   data.sampleinfo(i,1) = D.indsample(i);
%   data.sampleinfo(i,2) = D.indsample(i) + D.nsamples;
% end

author: spm


text:
function [dist,score,points] = getNormGTJointDistRP_inference(keypointsAll,gt_annolist,jidx)

dist = cell(length(gt_annolist),1);
score = cell(length(gt_annolist),1);
points = cell(length(gt_annolist),1);

for imgidx = 1:length(gt_annolist)
    assert(strcmp(keypointsAll(imgidx).imgname,gt_annolist(imgidx).image.name) > 0);
    dist{imgidx} = inf(0,length(gt_annolist(imgidx).annorect));
end

for imgidx = 1:length(gt_annolist)

    figure(200); clf; imagesc(imread(gt_annolist(imgidx).image.name));
    hold on; axis equal;
    det = keypointsAll(imgidx).det{jidx+1,1};

    for detidx = 1:2
        [val,idx] = max(det(:,end));
        sc = det(idx,3);
        pp = det(idx,1:2);
        %             plot(pp(1,1),pp(1,2),'r*','Markersize',15);
        for ridx = 1:length(gt_annolist(imgidx).annorect)
            rect = gt_annolist(imgidx).annorect(ridx);
            assert(length(rect) == 1);
            refDist = util_get_head_size(rect);
            if (isfield(rect, 'annopoints') && isfield(rect.annopoints, 'point'))
                p = util_get_annopoint_by_id(rect.annopoints.point, jidx);
                if (~isempty(p))
                    if (p.x>=x1 && p.x<= x2 && p.y >= y1 && p.y <= y2)
                        d(1,ridx) = norm([p.x p.y] - pp)/refDist;
                    end
                    %                         plot(p.x,p.y,'y+','Markersize',15);
                end
            end
        end
        score{imgidx} = [score{imgidx};sc];
        points{imgidx} = [points{imgidx};pp];
        dist{imgidx} = [dist{imgidx};d];
    end
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function []= mdp_silent()

% mdp_silent Ask for running resolution functions of the MDP Toolbox
% in silent mode

% MDPtoolbox: Markov Decision Processes Toolbox
% Copyright (C) 2009  INRA
% 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 the <ORGANIZATION> 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 HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
% INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
% BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
% OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
% OF THE POSSIBILITY OF SUCH DAMAGE.

global mdp_VERBOSE;

mdp_VERBOSE = 0;

author: Sable

repo: mcbench-benchmarks


text:
function idx = getDealiasingIndexes(~, N, nVars)
%GETGDEALIASINGINDEXES  % Returns indexes for dealiasing procedure (2/3-rule).
%
% See also SPINOP3.

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

toOne = floor(N/2) + 1 - ceil(N/6):floor(N/2) + ceil(N/6);
idx = false(N, N, N);
idx(toOne, toOne, toOne) = 1;
idx = repmat(idx, nVars, 1);

end


text:
% sgwt_kernel_abspline3 : Monic polynomial / cubic spline / power law decay kernel
%
% function r = sgwt_kernel_abspline3(x,alpha,beta,t1,t2)
%
% defines function g(x) with g(x) = c1*x^alpha for 0<x<x1
% g(x) = c3/x^beta for x>t2
% cubic spline for t1<x<t2,
% Satisfying g(t1)=g(t2)=1
%
% Inputs :
% x : array of independent variable values
% alpha : exponent for region near origin
% beta : exponent decay
% t1, t2 : determine transition region
%
% Outputs :
% r - result (same size as x)

% This file is part of the SGWT toolbox (Spectral Graph Wavelet Transform toolbox)
% Copyright (C) 2010, David K. Hammond.
%
% The SGWT toolbox 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.
%
% The SGWT toolbox 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 the SGWT toolbox.  If not, see <http://www.gnu.org/licenses/>.

function r = sgwt_kernel_abspline3(x,alpha,beta,t1,t2)
  r=zeros(size(x));
  % compute spline coefficients
  % M a = v
  M=[[1 t1 t1^2 t1^3];...
     [1 t2 t2^2 t2^3];...
     [0 1 2*t1 3*t1^2];...
     [0 1 2*t2 3*t2^2]];
  %v=[t1^alpha ; t2^(-beta) ; alpha*t1^(alpha-1) ; -beta*t2^(-beta-1)];
  v=[1 ; 1 ; t1^(-alpha)*alpha*t1^(alpha-1) ; -beta*t2^(-beta-1)*t2^beta];
  a=M\v;

  r1=find(x>=0 & x<t1);
  r2=find(x>=t1 & x<t2);
  r3=find(x>=t2);
  r(r1)=x(r1).^alpha*t1^(-alpha);
  r(r3)=x(r3).^(-beta)*t2^(beta);

  x2=x(r2);
  r(r2)=a(1)+a(2)*x2+a(3)*x2.^2+a(4)*x2.^3;
%  tmp=polyval(flipud(a),x2);
%  keyboard

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [node,elem,T] = squarequadmesh(square,h)
%% SQUAREQUADMESH uniform rectangular mesh of a square
%
% [node,elem] = squarequadmesh([x0,x1,y0,y1],h) generates a uniform mesh of the
% rectangle [x0,x1]*[y0,y1] with mesh size h.
%
% [node,elem,T] = squarequadmesh([x0,x1,y0,y1],h) with additional output T
% contains auxiliary structure: edge, elem2edge, edge2elem, bdEdge
%
% Example
%
%   [node,elem] = squarequadmesh([0,1,0,2],0.5);
%   showmesh(node,elem);
%   findnode(node);
%   findquadelem(node,elem);
%
% See also: squaremesh, squaregradmeshquad, cubehexmesh
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Generate nodes
x0 = square(1); x1 = square(2);
y0 = square(3); y1 = square(4);
% [x,y] = meshgrid(x0:h:x1,y1:-h:y0);
[x,y] = ndgrid(x0:h:x1,y0:h:y1);
node = [x(:),y(:)];

%% Generate elements
ni = size(x,1); % number of rows
nj = size(x,2);
N = size(node,1);
nodeidx = reshape(1:N,ni,nj);
t2nidxMap = nodeidx(1:ni-1,1:nj-1);
k = t2nidxMap(:);
elem = [k k+ni k+ni+1 k+1];
NT = size(elem,1);
% 4 k+1 --- k+ni+1 3
%    |        |
% 1  k  ---  k+ni  2

%% Check if further structure is needed
if nargout<=2
    return;
end

%% Generate edges
Ne1 = (ni-1)*nj;
Ne2 = ni*(nj-1);
edge = zeros(Ne1+Ne2,2,'int32');
% vertical edges
e2nidxMap = nodeidx(1:ni-1,1:nj);
k = e2nidxMap(:);
edge(1:Ne1,:) = [k k+1];
% horizontal edges
e2nidxMap = nodeidx(1:ni,1:nj-1);
k = e2nidxMap(:);
edge(Ne1+(1:Ne2),:) = [k k+ni];

%% Generate index map between element and edges
elem2edge = zeros(NT,4,'int32');
edge2elem = zeros(Ne1+Ne2,2,'int32');
k = 1:NT;
veidx = k; % veidx = sub2ind([ni-1,nj],i,j);
heidx = Ne1 + (1:Ne2); % linear indexing
heidx = reshape(heidx,ni,nj-1); % matrix indexing
heidx = heidx(1:ni-1,:); % without top edges
elem2edge(k,1) = heidx(:);
elem2edge(k,3) = heidx(:)+1;
elem2edge(k,2) = veidx + (ni-1);
elem2edge(k,4) = veidx;
% vertical edges
edge2elem(veidx,2) = k;   % right element
edge2elem(veidx+ni-1,1) = k; % left element
% horizontal edges
edge2elem(heidx,1) = k;
edge2elem(heidx+1,2) = k;

%% Generate boundary edges
% for boundary edges, set the another element idx to zero
idx1 = edge2elem(:,1) == 0; % left and top
idx2 = edge2elem(:,2) == 0;
bdEdge = [edge(idx1,[2 1]); edge(idx2,:)];
% switch edge index on the left and top to have a consistent orientation

%% Auxstructure
T = struct('edge',edge,'elem2edge',elem2edge,'edge2elem',edge2elem,'bdEdge',bdEdge);
% neighbor and bdElem can be easily get from 2-D index of all elements

repo: ifem


text:
function cdf = birthday_cdf ( n )

%*****************************************************************************80
%
%% BIRTHDAY_CDF returns the Birthday Concurrence CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 August 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of people whose birthdays have been
%    disclosed.
%
%    Output, real CDF, the probability that at least
%    two of the N people have matching birthays.
%
  if ( n < 1 )
    cdf = 0.0;
    return
  elseif ( 365 < n )
    cdf = 1.0;
    return
  end
%
%  Compute the probability that N people have distinct birthdays.
%
  cdf = 1.0;
  for i = 1 : n
    cdf = cdf * ( 365 + 1 - i ) / 365.0;
  end
%
%  Compute the probability that it is NOT the case that N people
%  have distinct birthdays.  This is the cumulative probability
%  that person 2 matches person 1, or person 3 matches 1 or 2,
%  etc.
%
  cdf = 1.0 - cdf;

  return
end


text:
% Script demonstrating usage of the cbpdndlmd function.
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2017-04-29
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'Copyright' and 'License' files
% distributed with the library.

% Training images
S0 = zeros(512, 512, 2, 'single');
S0(:,:,1) = single(stdimage('lena.grey')) / 255;
S0(:,:,2) = single(stdimage('barbara.grey')) / 255;

%Reduce images size to speed up demo script
tmp = zeros(128, 128, 2, 'single');
for k = 1:size(S0,3),
  tmp(:,:,k) = imresize(S0(:,:,k), 0.25);
end
S0 = tmp;

% Filter input images and compute highpass images
npd = 16;
fltlmbd = 5;
[Sl, Sh] = lowpass(S0, fltlmbd, npd);

% Construct weight matrix and padded test image set
Shp = padarray(Sh, [7 7], 'post');
t = 0.5;
W = randn(size(Sh));
W(abs(W) > t) = 1;
W(abs(W) < t) = 0;
W = padarray(W, [7 7], 'post');
ShW = W .* Shp;

% Construct initial dictionary
D0 = zeros(8,8,32, 'single');
D0(3:6,3:6,:) = single(randn(4,4,32));

% Set up cbpdndl parameters
lambda = 0.05;
opt = [];
opt.Verbose = 1;
opt.MaxMainIter = 500;
opt.rho = 50*lambda + 0.5;
opt.sigma = size(Sh,3);
opt.AutoRho = 1;
opt.AutoRhoPeriod = 10;
opt.AutoSigma = 1;
opt.AutoSigmaPeriod = 10;
opt.XRelaxParam = 1.8;
opt.DRelaxParam = 1.8;

% Do standard dictionary learning and reconstruct
[D1, X1, optinf1] = cbpdndl(D0, ShW, lambda, opt);
DX1 = ifft2(bsxfun(@times, fft2(D1, size(X1,1), size(X1,2)), fft2(X1)), ...
           'symmetric');
Sr1 = squeeze(sum(DX1,3)) + padarray(Sl, [7 7], 'post');

% Do dictionary learning with mask decoupling and reconstruct
opt.W = W;
[D2, X2, optinf2] = cbpdndlmd(D0, ShW, lambda, opt);
DX2 = ifft2(bsxfun(@times, fft2(D2, size(X2,1), size(X2,2)), fft2(X2)), ...
           'symmetric');
Sr2 = squeeze(sum(DX2,3)) + padarray(Sl, [7 7], 'post');

% Display dictionaries
figure;
subplot(1,2,1);
imdisp(tiledict(D1));
title('Standard DL');
subplot(1,2,2);
imdisp(tiledict(D2));
title('DL with mask decoupling');

% Display reconstructions
figure;
subplot(2,2,1);
imdisp(Sr1(:,:,1));
title('Standard DL');
subplot(2,2,2);
imdisp(Sr2(:,:,1));
title('DL with mask decoupling');
subplot(2,2,3);
imdisp(Sr1(:,:,2));
title('Standard DL');
subplot(2,2,4);
imdisp(Sr2(:,:,2));
title('DL with mask decoupling');

% Plot functional value evolution
figure;
subplot(1,2,1);
semilogx(optinf1.itstat(:,2), 'LineWidth', 2);
ylim([15, 45]);
xlabel('Iterations');
ylabel('Functional value');
title('Standard DL');
subplot(1,2,2);
semilogx(optinf2.itstat(:,2), 'LineWidth', 2);
ylim([15, 45]);
xlabel('Iterations');
ylabel('Functional value');
title('DL with mask decoupling');

author: thfylsty


text:
function img_crop = imcrop_pad(img, bbox, padding, sz)
xi = linspace(-1, 1, sz(2));
yi = linspace(-1, 1, sz(1));
[xx,yy] = meshgrid(xi,yi);
yyxx = single([yy(:), xx(:)]') ; % 2xM
[im_h,im_w,im_c,~] = size(img);
if im_c == 1
    img = repmat(img,[1,1,3,1]);
end

target_crop_w = (1+padding)*(bbox(:,3)-bbox(:,1));
target_crop_h = (1+padding)*(bbox(:,4)-bbox(:,2));
target_crop_cx = (bbox(:,1)+bbox(:,3))/2;
target_crop_cy = (bbox(:,2)+bbox(:,4))/2;

cy_t = (target_crop_cy*2/(im_h-1))-1;
cx_t = (target_crop_cx*2/(im_w-1))-1;

h_s = target_crop_h/(im_h-1);
w_s = target_crop_w/(im_w-1);

s = reshape([h_s,w_s]', 2, 1, []); % x,y scaling
t = reshape([cy_t,cx_t]', 2, 1, []); % translation

g = bsxfun(@times, yyxx, s); % scale
g = bsxfun(@plus, g, t); % translate
g = reshape(g, 2, sz(1), sz(2), []);

img_crop = vl_nnbilinearsampler(img, g);

end

author: foolwood

repo: DCFNet

save_path: github-repos/MATLAB/foolwood-DCFNet


text:
function datapt = lvd_ThrottleTask(stateLogEntry, subTask, inFrame)
%lvd_ThrottleTask Summary of this function goes here
%   Detailed explanation goes here

    switch subTask
        case 'throttle'
            datapt = 100*stateLogEntry.throttle;

        case 't2w'
            throttle = stateLogEntry.throttle;
            tankStates = stateLogEntry.getAllActiveTankStates();

            tankMasses = zeros(size(tankStates));
            for(i=1:length(tankStates))
                tankMasses(i) = tankStates(i).tankMass;
            end

            powerStorageStates = stateLogEntry.getAllActivePwrStorageStates();
            storageSoCs = NaN(size(powerStorageStates));
            for(i=1:length(powerStorageStates))
                storageSoCs(i) = powerStorageStates(i).getStateOfCharge();
            end

            datapt = computeTWRatio(throttle, stateLogEntry.time, stateLogEntry.position, stateLogEntry.velocity, tankMasses, stateLogEntry.getTotalVehicleDryMass(), ...
                                    stateLogEntry.stageStates, stateLogEntry.lvState, tankStates, stateLogEntry.centralBody, storageSoCs, powerStorageStates);
        case 'totalthrust'
            [totalThrust] = getThrustParameters(stateLogEntry);
            datapt = totalThrust;

        case 'thrust_x'
            vehElemSet = stateLogEntry.getCartesianElementSetRepresentation();
            [~, ~, ~, rotMat] = inFrame.getOffsetsWrtInertialOrigin(stateLogEntry.time, vehElemSet);

            [~, thrustForceVector] = getThrustParameters(stateLogEntry);
            thrustForceVector = thrustForceVector(:);
            thrustForceVector = rotMat * thrustForceVector;

            datapt = thrustForceVector(1);
            if(isnan(datapt))
                datapt = 0;
            end

        case 'thrust_y'
            vehElemSet = stateLogEntry.getCartesianElementSetRepresentation();
            [~, ~, ~, rotMat] = inFrame.getOffsetsWrtInertialOrigin(stateLogEntry.time, vehElemSet);

            [~, thrustForceVector] = getThrustParameters(stateLogEntry);
            thrustForceVector = thrustForceVector(:);
            thrustForceVector = rotMat * thrustForceVector;

            datapt = thrustForceVector(2);
            if(isnan(datapt))
                datapt = 0;
            end

        case 'thrust_z'
            vehElemSet = stateLogEntry.getCartesianElementSetRepresentation();
            [~, ~, ~, rotMat] = inFrame.getOffsetsWrtInertialOrigin(stateLogEntry.time, vehElemSet);

            [~, thrustForceVector] = getThrustParameters(stateLogEntry);
            thrustForceVector = thrustForceVector(:);
            thrustForceVector = rotMat * thrustForceVector;

            datapt = thrustForceVector(3);
            if(isnan(datapt))
                datapt = 0;
            end

        case 'thrust_vector'
            [~, thrustForceVector] = getThrustParameters(stateLogEntry);

            datapt = thrustForceVector(:);
            if(any(isnan(datapt)))
                datapt = [0;0;0];
            end

        otherwise
            error('Unrecongized task: %s', subTask);

    end
end

function [totalThrust, thrustForceVector] = getThrustParameters(stateLogEntry)
    if(stateLogEntry.event.propagatorObj.canProduceThrust())
        ut = stateLogEntry.time;
        rVect = stateLogEntry.position;
        vVect = stateLogEntry.velocity;

        bodyInfo = stateLogEntry.centralBody;
        tankStates = stateLogEntry.getAllActiveTankStates();
        stageStates = stateLogEntry.stageStates;
        lvState = stateLogEntry.lvState;

        dryMass = stateLogEntry.getTotalVehicleDryMass();
        tankStatesMasses = [tankStates.tankMass]';

        throttleModel = stateLogEntry.throttleModel;
        steeringModel = stateLogEntry.steeringModel;

        altitude = norm(rVect) - bodyInfo.radius;
        pressure = getPressureAtAltitude(bodyInfo, altitude);

        powerStorageStates = stateLogEntry.getAllActivePwrStorageStates();
        storageSoCs = NaN(size(powerStorageStates));
        for(i=1:length(powerStorageStates)) %#ok<*NO4LP>
            storageSoCs(i) = powerStorageStates(i).getStateOfCharge();
        end

        attState = LaunchVehicleAttitudeState();
        attState.dcm = steeringModel.getBody2InertialDcmAtTime(ut, rVect, vVect, bodyInfo);

        throttle = throttleModel.getThrottleAtTime(ut, rVect, vVect, tankStatesMasses, dryMass, stageStates, lvState, tankStates, bodyInfo, storageSoCs, powerStorageStates);

        [~, totalThrust, thrustForceVector] = LaunchVehicleStateLogEntry.getTankMassFlowRatesDueToEngines(tankStates, tankStatesMasses, stageStates, throttle, lvState, pressure, ut, rVect, vVect, bodyInfo, steeringModel, storageSoCs, powerStorageStates, attState);
        thrustForceVector = 1000*thrustForceVector; %in order to recover kN
    else
        totalThrust = 0;
        thrustForceVector = [0;0;0];
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [vOut,data]=getNextMultipartition(data)
%%GETNEXTMULTIPARTITION Get the next multipartition of a set of elements
%           with possible repeats, or get the first multipartition in the
%           sequence. The multipartitions are given in decreasing
%           lexicographic order.
%
%INPUTS: data If the first multipartition is desired, then this is an mX1
%             or 1Xm vector containing how many times each element in the
%             set to be partitioned is repeated. For example, if one wishes
%             to partition the set (1,2,3,4), then this is [1;1;1;1]. For
%             the set [1;2;2;2;3;3;4], then this is [1;3;2;1]. Otherwise,
%             if one wishes to get the next multipartition in the sequence,
%             then this is the data variable returned from the previous
%             call to this function.
%
%OUTPUTS: vOut This is the current multipartition, or an empty vector if
%              one has passed the final multipartition in the sequence. The
%              multipartition is an mXnumParts matrix. The number of
%              columns is the number of parts of the partition. Each column
%              has the number of repeats of each item (corresponding to the
%              positions in the original data passed). See the examples
%              below for more explanation.
%         data This is a structure that can be passed to this function to
%              get the next multipartition in the sequence.
%
%This function implements Algorithm M of Chapter 7.2.1.5 of [1].
%
%As an example, there are nine multipartitions of the set (1,1,2,2). These
%are
%(1,1,2,2);  (1,1,2),(2);  (1,1),(2,2);  (1,1),(2),(2); (1,2,2),(1);
%(1,2),(1,2);  (1,2),(1),(2);  (1),(1),(2,2); (1),(1),(2),(2)
%However, represented in terms of the output vOut, we have 9 matrices:
% [2;  [2,0; [2,0; [2,0,0;  [1,1;  [1,1;  [1,1,0;  [1,1,0;  [1,1,0,0
%  2]   1,1]  0,2]  0,1,1]   2,0]   1,1]   1,0,1]   0,0,2]   0,0,1,1]
%The example below shows how to use this function and how to display the
%outputs in a more recognizable form as partitions.
%
%EXAMPLE:
%Here, we display all multipartitions of 1123, which means that n=[2;1;1].
% vTotal=[1;1;2;3];
% vUnique=unique(vTotal);
% n=[2;1;1];
% m=length(n);
% [vOut,data]=getNextMultipartition(n);
% while(~isempty(vOut))
%     %Display the current multipartition.
%     %We will build a string showing the partitions. of [1,1,2,3].
%     strDisp=[];
%
%     numParts=size(vOut,2);
%     for curPart=1:numParts
%         strDisp=[strDisp,sprintf('(')];
%         for k=1:m
%             for i=1:vOut(k,curPart)
%                 strDisp=[strDisp,sprintf('%i',vUnique(k))];
%             end
%         end
%         strDisp=[strDisp,sprintf(')')];
%     end
%     disp(strDisp)
%
%     %Get the next one.
%     [vOut,data]=getNextMultipartition(data);
% end
%
%REFERENCES:
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 4, Fascicle 3:
%    Generating all Combinations and Partitions, Upper Saddle River, NJ:
%    Addison-Wesley, 2009.
%
%February 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%If the first multipartition is desired.
if(~isstruct(data))
    n=data;
    data=[];

    m=length(n);
    nSum=sum(n);

    data.nSum=nSum;
    data.m=m;

    %Allocate space
    f=zeros(nSum+1,1);
    c=zeros(m*nSum+1,1);
    u=zeros(m*nSum+1,1);
    v=zeros(m*nSum+1,1);

    %Step M1
    c(1:m)=1:m;
    u(1:m)=n(1:m);
    v(1:m)=n(1:m);
    f(0+1)=0;
    a=0;
    l=0;

    f(1+1)=m;
    b=m;
else
    %If the algorithm is finished.
    if(data.algFinished)
       vOut=[];
       return;
    end

    f=data.f;
    c=data.c;
    u=data.u;
    v=data.v;
    a=data.a;
    b=data.b;
    l=data.l;
    m=data.m;
    nSum=data.nSum;
end

while(1)
    %Step M2
    j=a;
    k=b;
    x=0;

    while(j<b)
        u(k+1)=u(j+1)-v(j+1);

        if(u(k+1)==0)
            x=1;
            j=j+1;
        elseif(x==0)
            c(k+1)=c(j+1);
            v(k+1)=min(v(j+1),u(k+1));
            x=u(k+1)<v(j+1);
            k=k+1;
            j=j+1;
        else
            c(k+1)=c(j+1);
            v(k+1)=u(k+1);
            k=k+1;
            j=j+1;
        end
    end

    %Step M3
    if(k>b)
        a=b;
        b=k;
        l=l+1;
        f(l+1+1)=b;
        continue;%Return to M2
    end

    %Step M4, visit the partition
    numOut=0;
    vOut=zeros(m,nSum);%Allocate the maximum possible space needed.
    for k=0:l
        numOut=numOut+1;
        for j=f(k+1):(f(k+1+1)-1)
            vOut(c(j+1),numOut)=v(j+1);
        end
    end
    vOut=vOut(:,1:numOut);

    %Step M5
    while(1)

        j=b-1;
        while(v(j+1)==0)
            j=j-1;
        end

        if(~(j==a&&v(j+1)==1))
            v(j+1)=v(j+1)-1;

            for k=(j+1):(b-1)
               v(k+1)=u(k+1);
            end

            %Return to M2 to visit the next one, so return from this
            %function.
            data.f=f;
            data.c=c;
            data.u=u;
            data.v=v;
            data.a=a;
            data.b=b;
            data.l=l;
            data.algFinished=false;
            return
        end

        %Step M6
        if(l==0)%All multiset partitions have been found.
            %There is no need to update the other things in data, because
            %there are no more multipartitions to find.
            data.algFinished=true;
            return;
        end

        l=l-1;
        b=a;
        a=f(l+1);
    end
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [R,L] = match_query(D,SR)
% [R,L] = match_query(D,SR)
%     Match landmarks from an audio query against the database.
%     Rows of R are potential maxes, in format
%      songID  modalDTcount modalDT
%     i.e. there were <modalDTcount> occurrences of hashes
%     that occurred in the query and reference with a difference of
%     <modalDT> frames.
%     L returns the actual landmarks that this implies.
% 2008-12-29 Dan Ellis dpwe@ee.columbia.edu

%Rt = get_hash_hits(landmark2hash(find_landmarks(D,SR)));
Lq = find_landmarks(D,SR);
%Lq = fuzzify_landmarks(Lq);
% Augment with landmarks calculated half-a-window advanced too
landmarks_hopt = 0.032;
%Lq = [Lq;find_landmarks(D(round(landmarks_hopt/4*SR):end),SR)];
Lq = [Lq;find_landmarks(D(round(landmarks_hopt/2*SR):end),SR)];
%Lq = [Lq;find_landmarks(D(round(3*landmarks_hopt/4*SR):end),SR)];
% add in quarter-hop offsets too for even better recall

Hq = landmark2hash(Lq);
Rt = get_hash_hits(Hq);
nr = size(Rt,1);

if nr > 0

  % Find all the unique tracks referenced
  [utrks,xx] = unique(sort(Rt(:,1)));
  utrkcounts = diff(xx,nr);

  nutrks = length(utrks);

  R = zeros(nutrks,3);

  for i = 1:nutrks
    tkR = Rt(Rt(:,1)==utrks(i),:);
    % Find the most popular time offset
    [dts,xx] = unique(sort(tkR(:,2)),'first');
    dtcounts = 1+diff([xx',size(tkR,1)]);
    [vv,xx] = max(dtcounts);
    R(i,:) = [utrks(i),vv,dts(xx)];
  end

  % Sort by descending match count
  [vv,xx] = sort(R(:,2),'descend');
  R = R(xx,:);

  % Extract the actual landmarks
  H = Rt((Rt(:,1)==R(1,1)) & (Rt(:,2)==R(1,3)),:);
  % Restore the original times
  for i = 1:size(H,1)
    hix = find(Hq(:,3)==H(i,3));
    hix = hix(1);  % if more than one...
    H(i,2) = H(i,2)+Hq(hix,2);
    L(i,:) = hash2landmark(H(i,:));
  end

  % Return no more than 10 hits, and only down to half the #hits in
  % most popular
  if size(R,1) > 10
    R = R(1:10,:);
  end
  maxhits = R(1,2);
  nuffhits = R(:,2)>(maxhits/2);
  %R = R(nuffhits,:);

else
  R = [];
  disp('*** NO HITS FOUND ***');
end


text:
function t = set(t,varargin)
% T = SET(T.PROP,VAL); Set object properties
% Usage:
% t = set(t,property1,value1,property2,value2,...)
%	Sets values of properties 'property1', 'property2', etc. to values
%	'value1', value2', ... respectively

if nargin<3
    error('Sorry, Insufficient number of arguments');
elseif rem(nargin-1,2)
    error('Sorry, Property and value must be given in pairs');
else
    if ~isa(t,'acqparam')
        error('First input must belong to class acqparam');
    end;
    inargs = varargin;              % Input arguments
    taglist = fieldnames(struct(t)); % List of fields
    while length(inargs)>=2
        prop = inargs{1};           % One property at a time
        val  = inargs{2};
        inargs = inargs(3:end);
        tagindx = find(strcmpi(taglist(:),prop)); % Find which property?
        if ~isempty(tagindx)        % If regular field
            eval(['t.' taglist{tagindx} ' = val;']);
        else                        % If SampleInfo field
            sampinfotaglist = fieldnames(t.SampleInfo);
            subindx = find(strcmpi(sampinfotaglist(:),prop));
            if ~isempty(subindx)
                eval(['t.SampleInfo.' sampinfotaglist{subindx} '=val;']);
            else
                error(['Couldnt find requested field - ''',prop, ''' to update']);
            end;
        end;
    end;
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function K=commutationMatrix(m,n)
%%COMMUTATIONMATRIX Return a commutation matrix K for an mXn matrix. This
%          means that given an mXn matrix A that K*A(:)=B(:) where
%          B=transpose(A). As shown below, this matrix can also be used for
%          reversing the order of Kronecker products. As in [2], this
%          matrix is also known as a vec-permutation matrix.
%
%INPUTS: m,n The number of rows and the number of columns of the matrix
%            that is to be commutated by this commutation matrix.
%
%OUTPUTS: K An (m*n)X(m*n) commutation matrix.
%
%The matrix K is just a type of permutation matrix with 1's in the correct
%places.
%
%Properties of the commutation matrix are given in [1] and in [2]. For
%example, commutationMatrix(m,n)=commutationMatrix(n,m)' and also there is
%an inverse relation that
%commutationMatrix(m,n)*commutationMatrix(n,m)==eye(n*m) and
%Note that the Vec permutation matrix I_{m,n} as in [2] is actually
%commutationMatrix(n,m).
%
%EXAMPLE 1:
%Here, we demonstrate the property that the vec-permutation matrix
%(commutation matrix) allows one to get vec(A') from vec(A):
% m=12;
% n=17;
% A=randn(m,n);
% Imn=commutationMatrix(n,m);%The vec-permutation matrix as in [2].
% all(vec(A)==Imn*vec(A'))
%The result will be true, showing that
%vec(A)==commutationMatrix(n,m)*vec(A')
%which agrees with the identity given in the abstract to [2].
%
%EXAMPLE 2:
%Here, we give a demonstration of how to reverse the ordering of a
%Kronecker product of two matrices:
% m=12;
% n=17;
% A=randn(m,n);
% p=4;
% q=6;
% B=randn(p,q);
% Imp=commutationMatrix(p,m);
% Iqn=commutationMatrix(n,q);
% %Reversing the order of the matrix Kronecker product as in Equation 25 in
% %[2].
% all(vec(kron(B,A)==Imp*kron(A,B)*Iqn))
%The result will be true, showing that one can use two commutation matrices
%to reverse the ordering of Kronecker products of matrices.
%
%EXAMPLE 3:
%Now, we demonstrate that one can use a single commutation matrix (vec-
%permutation matrix) to reverse the order of the Kronecker product of two
%vectors.
% m=13;
% n=7;
% a=randn(m,1);
% b=randn(n,1);
% Imat=commutationMatrix(m,n);
% all(vec(kron(a,b)==Imat*kron(b,a)))
%The result will be true.
%
%REFERENCES:
%[1] J. R. Magnus and H. Neudecker, "The elimination matrix: Some lemmas
%    and applications," SIAM Journal on Algebraic Discrete Methods, vol. 1,
%    no. 4, pp. 422-449, 1980.
%[2] H. V. Henderson and S. R. Searle, "The vec-permutation matrix, the vec
%    operator and Kronecker products: A review," Linear and Multilinear
%    Algebra, vol. 9, no. 4, pp. 271-288, 1981.
%
%December 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

K=zeros(m*n,m*n);

idx=0;
col=0;

for i=1:n
    idx=idx+1;
    row=idx;

    for j=1:m
        col=col+1;
        K(row,col)=1;
        row=row+n;
    end
end

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
function x=modcent(x,r);
%MODCENT  Centered modulo
%   Usage:  y=modcent(x,r);
%
%   `modcent(x,r)` computes the modulo of *x* in the range $[-r/2,r/2[$.
%
%   As an example, to compute the modulo of *x* in the range $[-\pi,\pi[$ use
%   the call::
%
%     y = modcent(x,2*pi);

x=mod(x,r);
idx=x>r/2;
x(idx)=x(idx)-r;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
clear; clc;

%% DFT filter bank
N = 8;          % number of subbands
K = 16;         % overlapping factor
L = 2*K*N;      % length of analysis filters

hopt = fir1(L-1,1/N);       % prototype filter
[H,F] = make_bank_DFT(hopt,N);
H = sqrt(N)*H;                    % analysis section
F = sqrt(N)*F;                    % synthesis section

%% signals in the DFT filter bank
xn = sin(2*pi*0.01*(0:1000));
ITER = length(xn);
yn = zeros(ITER,1);
x = zeros(length(H),1);
b = zeros(length(F),1);

%% analysis and synthesis filter banks implementation
for n = 1:ITER
    x = [xn(n); x(1:end-1)];
    % block transformation is performed once for every N input samples
    if mod(n,N) == 0
        x_D = H.'*x;
        % =====================
        % insert subband processing here
        % =====================
        b = F*x_D + b;
    end
    yn(n) = real(b(1));
    b = [b(2:end); 0];
end

%% plot input and output signal
figure;
plot(0:length(xn)-1, xn); hold on;
plot(0:length(yn)-1, yn, 'r:');
xlabel('Time index, n'); ylabel('Amplitude');

%{
 %% plot frequency response
DFTpoint = 4096;
[Hz,w] = FreqResp(H, DFTpoint);
Hz = Hz.';
figure;
subplot(3,1,[1 2]); hold on;
for k = 1:N
    if mod(k,2) == 0
        plot(w/pi,20*log10(abs(Hz(k,:))));
    else
        plot(w/pi,20*log10(abs(Hz(k,:))),':');
    end
end
title('frequency response');
%}

author: CharlesThaCat

repo: acoustic-interference-cancellation

path: github-repos/MATLAB/CharlesThaCat-acoustic-interference-cancellation/acoustic-interference-cancellation-edb394499ea6f9c96445a3e9613bd64a854c289e/Subband processing/Book/After reading Chapter 2 -- Design Analysis & Synthesis Filter bank/DFTFilterBank.m


text:
function c = vandermonde_approx_2d_coef ( n, m, x, y, z )

%*****************************************************************************80
%
%% VANDERMONDE_APPROX_2D_COEF computes a 2D polynomial approximant.
%
%  Discussion:
%
%    We assume the approximating function has the form of a polynomial
%    in X and Y of total degree M.
%
%      p(x,y) = c00
%             + c10 * x                + c01 *  y
%             + c20 * x^2   + c11 * xy + c02 * y^2
%             + ...
%             + cm0 * x^(m) + ...      + c0m * y^m.
%
%    If we let T(K) = the K-th triangular number
%            = sum ( 1 <= I <= K ) I
%    then the number of coefficients in the above polynomial is T(M+1).
%
%    We have n data locations (x(i),y(i)) and values z(i) to approximate:
%
%      p(x(i),y(i)) = z(i)
%
%    This can be cast as an NxT(M+1) linear system for the polynomial
%    coefficients:
%
%      [ 1 x1 y1  x1^2 ... y1^m ] [ c00 ] = [  z1 ]
%      [ 1 x2 y2  x2^2 ... y2^m ] [ c10 ] = [  z2 ]
%      [ 1 x3 y3  x3^2 ... y3^m ] [ c01 ] = [  z3 ]
%      [ ...................... ] [ ... ] = [ ... ]
%      [ 1 xn yn  xn^2 ... yn^m ] [ c0m ] = [  zn ]
%
%    In the typical case, N is greater than T(M+1) (we have more data and
%    equations than degrees of freedom) and so a least squares solution is
%    appropriate, in which case the computed polynomial will be a least squares
%    approximant to the data.
%
%    The polynomial defined by the T(M+1) coefficients C could be evaluated
%    at the Nx2-vector x by the command
%
%      pval = r8poly_value_2d ( m, c, n, x )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of data points.
%
%    Input, integer M, the maximum degree of the polynomial.
%
%    Input, real X(N,1), Y(N,1), the data locations.
%
%    Input, real Z(N,1), the data values.
%
%    Output, real C(T(M+1),1), the coefficients of the approximating
%    polynomial.  C(1) is the constant term, and C(T(M+1)) multiplies Y^M.
%
  tm = triangle_num ( m + 1 );

  A = vandermonde_approx_2d_matrix ( n, m, tm, x, y );

  c = A \ z(:);

  return
end


text:
Network ncf_gemm {
Layer GEMM0 {
Type: CONV
Dimensions { K: 128, C: 2048, Y: 256, X: 1, R: 1, S: 1 }

}
Layer GEMM1 {
Type: CONV
Dimensions { K: 64, C: 2048, Y: 128, X: 1, R: 1, S: 1 }

}
Layer GEMM2 {
Type: CONV
Dimensions { K: 256, C: 2048, Y: 256, X: 1, R: 1, S: 1 }

}
Layer GEMM3 {
Type: CONV
Dimensions { K: 256, C: 256, Y: 2048, X: 1, R: 1, S: 1 }

}
Layer GEMM4 {
Type: CONV
Dimensions { K: 256, C: 256, Y: 2048, X: 1, R: 1, S: 1 }

}
Layer GEMM5 {
Type: CONV
Dimensions { K: 128, C: 256, Y: 2048, X: 1, R: 1, S: 1 }

}
Layer GEMM6 {
Type: CONV
Dimensions { K: 256, C: 128, Y: 2048, X: 1, R: 1, S: 1 }

}
Layer GEMM7 {
Type: CONV
Dimensions { K: 128, C: 64, Y: 2048, X: 1, R: 1, S: 1 }

}
Layer GEMM8 {
Type: CONV
Dimensions { K: 64, C: 128, Y: 2048, X: 1, R: 1, S: 1 }

}
Layer GEMM9 {
Type: CONV
Dimensions { K: 1, C: 2048, Y: 128, X: 1, R: 1, S: 1 }

}
Layer GEMM10 {
Type: CONV
Dimensions { K: 1, C: 128, Y: 2048, X: 1, R: 1, S: 1 }

}
Layer GEMM11 {
Type: CONV
Dimensions { K: 128, C: 1, Y: 2048, X: 1, R: 1, S: 1 }

}
}

author: maestro-project

repo: maestro


text:
function K = covCos(hyp, x, z, i)

% Stationary covariance function for a sinusoid with period p in 1d:
%
% k(x,z) = sf^2*cos(2*pi*(x-z)/p)
%
% where the hyperparameters are:
%
% hyp = [ log(p)
%         log(sf) ]
%
% Note that covPeriodicNoDC converges to covCos as ell goes to infinity.
%
% Copyright (c) by James Robert Lloyd, 2013-08-05.
%
% See also COVFUNCTIONS.M, COVPERIODICNODC.M.

if nargin<2, K = '2'; return; end                  % report number of parameters
if nargin<3, z = []; end                                   % make sure, z exists
xeqz = isempty(z); dg = strcmp(z,'diag');                       % determine mode

[n,D] = size(x);
if D>1, error('Covariance is defined for 1d data only.'), end
p   = exp(hyp(1));
sf2 = exp(2*hyp(2));

% precompute distances
if dg                                                               % vector kxx
  K = zeros(size(x,1),1);
else
  if xeqz                                                 % symmetric matrix Kxx
    K = repmat(x,1,n) - repmat(x',n,1);
  else                                                   % cross covariances Kxz
    K = repmat(x,1,size(z,1)) - repmat(z',n,1);
  end
end

K = 2*pi*K/p;
if nargin<4                                                        % covariances
    K = sf2*cos(K);
else                                                               % derivatives
  if i==1
    K = sf2*sin(K).*K;
  elseif i==2
    K = 2*sf2*cos(K);
  else
    error('Unknown hyperparameter')
  end
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [frames, descrs] = vl_phow(im, varargin)
% VL_PHOW  Extract PHOW features
%   [FRAMES, DESCRS] = VL_PHOW(IM) extracts PHOW features [1] from the
%   image IM. PHOW is simply dense SIFT applied at several resolutions. This function is a commodity interface to VL_DSIFT() and
%   VL_IMSMOOTH().
%
%   DESCRS has the same format of VL_SIFT() and VL_DSIFT(). FRAMES(1:2,:)
%   are the x,y coordinates of the center of each descriptor, FRAMES(3,:)
%   is the contrast of the descriptor, as returned by VL_DSIFT() (for
%   colour variant, contranst is computed on the intensity channel).
%   FRAMES(4,:) is the size of the bin of the descriptor.
%
%   By default,
%   VL_PHOW() computes the gray-scale variant of the descriptor.  The
%   COLOR option can be used to compute the color variant instead.
%
%   Verbose:: false
%     Set to true to turn on verbose output.
%
%   Sizes:: [4 6 8 10]
%     Scales at which the dense SIFT features are extracted. Each
%     value is used as bin size for the VL_DSIFT() function.
%
%   Fast:: true
%     Set to false to turn off the fast SIFT features computation by
%     VL_DSIFT().
%
%   Step:: 2
%     Step (in pixels) of the grid at which the dense SIFT features
%     are extracted.
%
%   Color:: 'gray'
%     Choose between 'gray' (PHOW-gray), 'rgb', 'hsv', and 'opponent'
%     (PHOW-color).
%
%   ContrastThreshold:: 0.005
%     Contrast threshold below which SIFT features are mapped to
%     zero. The input image is scaled to have intensity range in [0,1]
%     (rather than [0,255]) and this value is compared to the
%     descriptor norm as returned by VL_DSIFT().
%
%   WindowSize:: 1.5
%     Size of the Gaussian window in units of spatial bins.
%
%   Magnif:: 6
%     The image is smoothed by a Gaussian kernel of standard deviation
%     SIZE / MAGNIF. Note that, in the standard SIFT descriptor, the
%     magnification value is 3; here the default one is 6 as it seems
%     to perform better in applications.
%
%   FloatDescriptors:: false
%     If set to TRUE, the descriptors are returned in floating point
%     format.
%
%   See also: VL_DSIFT(), VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

% -------------------------------------------------------------------
%                                                 Parse the arguments
% -------------------------------------------------------------------

  opts.verbose = false ;
  opts.fast = true ;
  opts.sizes = [4 6 8 10] ;
  opts.step = 2 ;
  opts.color = 'gray' ;
  opts.floatdescriptors = false ;
  opts.magnif = 6 ;
  opts.windowsize = 1.5 ;
  opts.contrastthreshold = 0.005 ;
  opts = vl_argparse(opts,varargin) ;

  dsiftOpts = {'norm', 'windowsize', opts.windowsize} ;
  if opts.verbose, dsiftOpts{end+1} = 'verbose' ; end
  if opts.fast, dsiftOpts{end+1} = 'fast' ; end
  if opts.floatdescriptors, dsiftOpts{end+1} = 'floatdescriptors' ; end
  dsiftOpts(end+(1:2)) = {'step', opts.step} ;

% -------------------------------------------------------------------
%                                                Extract the features
% -------------------------------------------------------------------

  % standarize the image
  imageSize = [size(im,2) ; size(im,1)] ;
  if strcmp(lower(opts.color), 'gray')
    numChannels = 1 ;
    if size(im,3) > 1, im = rgb2gray(im) ; end
  else
    numChannels = 3 ;
    if size(im,3) == 1, im = cat(3, im, im, im) ; end
    switch lower(opts.color)
      case 'rgb'
      case 'opponent'
        % Note that the mean differs from the standard definition of opponent
        % space and is the regular intesity (for compatibility with
        % the contrast thresholding).
        %
        % Note also that the mean is added pack to the other two
        % components with a small multipliers for monochromatic
        % regions.
        mu = 0.3*im(:,:,1) + 0.59*im(:,:,2) + 0.11*im(:,:,3) ;
        alpha = 0.01 ;
        im = cat(3, mu, ...
                 (im(:,:,1) - im(:,:,2))/sqrt(2) + alpha*mu, ...
                 (im(:,:,1) + im(:,:,2) - 2*im(:,:,3))/sqrt(6) + alpha*mu) ;
      case 'hsv'
        im = rgb2hsv(im) ;
      otherwise
        opts.color = 'hsv' ;
        warning('Color space not recongized, defaulting to HSV color space.') ;
    end
  end

  if opts.verbose
    fprintf('%s: color space: %s\n', mfilename, opts.color) ;
    fprintf('%s: image size: %d x %d\n', mfilename, imageSize(1), imageSize(2)) ;
    fprintf('%s: sizes: [%s]\n', mfilename, sprintf(' %d', opts.sizes)) ;
  end

  for si = 1:length(opts.sizes)

    % Recall from VL_DSIFT() that the first descriptor for scale SIZE has
    % center located at XC = XMIN + 3/2 SIZE (the Y coordinate is
    % similar). It is convenient to align the descriptors at different
    % scales so that they have the same geometric centers. For the
    % maximum size we pick XMIN = 1 and we get centers starting from
    % XC = 1 + 3/2 MAX(OPTS.SIZES). For any other scale we pick XMIN so
    % that XMIN + 3/2 SIZE = 1 + 3/2 MAX(OPTS.SIZES).
    %
    % In pracrice, the offset must be integer ('bounds'), so the
    % alignment works properly only if all OPTS.SZES are even or odd.

    off = floor(1 + 3/2 * (max(opts.sizes) - opts.sizes(si))) ;

    % smooth the image to the appropriate scale based on the size
    % of the SIFT bins
    sigma = opts.sizes(si) / opts.magnif ;
    ims = vl_imsmooth(im, sigma) ;

    % extract dense SIFT features from all channels
    for k = 1:numChannels
      [f{k}, d{k}] = vl_dsift(...
        ims(:,:,k), ...
        dsiftOpts{:},  ...
        'size', opts.sizes(si), ...
        'bounds', [off off +inf +inf]) ;
    end

    % remove low contrast descriptors
    % note that for color descriptors the V component is
    % thresholded
    switch lower(opts.color)
      case {'gray', 'opponent'}
        contrast = f{1}(3,:) ;
      case 'rgb'
        contrast = mean([f{1}(3,:) ; f{2}(3,:) ; f{3}(3,:)],1) ;
      otherwise % hsv
        contrast = f{3}(3,:) ;
    end
    for k = 1:numChannels
      d{k}(:, contrast < opts.contrastthreshold) = 0 ;
    end

    % save only x,y, and the scale
    frames{si} = [f{1}(1:3, :) ; opts.sizes(si) * ones(1,size(f{1},2))] ;
    descrs{si} = cat(1, d{:}) ;
  end
  descrs = cell2mat(descrs) ;
  frames = cell2mat(frames) ;
end

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function [Dtf, Dtph] = spm_eeg_tf(S)
% Compute instantaneous power and phase in peri-stimulus time and frequency
% FORMAT [Dtf, Dtph] = spm_eeg_tf(S)
%
% S                     - input structure
%
% fields of S:
%   S.D                 - MEEG object or filename of M/EEG mat-file with
%
%   S.channels          - cell array of channel names. Can include generic
%                         wildcards: 'All', 'EEG', 'MEG' etc.
%
%   S.frequencies      - vector of frequencies of interest
%
%   S.timewin          - time window of interest in PST in ms.
%
%   S.method           - name for the spectral estimation to use. This
%                        corresponds to the name of a plug-in function that comes
%                        after 'spm_eeg_specest_' prefix.
%   S.settings         - plug-in specific settings
%
%   S.phase            - also save phase dataset (1) or not (0)
%                        phase dataset cannot be computed for some
%                        spectral estimation methods
%   S.prefix           - prefix added before the standard prefix (tf_ or tph_)
%
% Output:
% Dtf                   - M/EEG object with power (also written on disk)
% Dtph                  - M/EEG object with phase (also written on disk)
%__________________________________________________________________________
% This is a modular function for which plugins can be developed implementing
% specific spectral estimation methods. There are 3 basic plugins presently
% implemented and they can be used as templates for new plugins.
% The name of a plugin function should start with 'spm_eeg_specest_'
%
% morlet (spm_eeg_specest_morlet) - Morlet wavelet transform
%
% hilbert (spm_eeg_specest_hilbert) - filtering + Hilbert transform
%
% ft_mtmconvol (spm_eeg_specest_mtmconvol) - Fieldtrip implementation
%                                           of multi-taper spectral
%                                           analysis
%__________________________________________________________________________
% Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_tf.m 7449 2018-10-16 13:52:04Z vladimir $

SVNrev = '$Rev: 7449 $';

%-Startup
%--------------------------------------------------------------------------
spm('FnBanner', mfilename, SVNrev);
spm('FigName','M/EEG Time-Frequency'); spm('Pointer','Watch');

%-Configure the analysis
%--------------------------------------------------------------------------
if ~isfield(S, 'channels'),   S.channels = 'all';                 end
if ~isfield(S, 'timewin'),    S.timewin  = [-Inf Inf];            end
if ~isfield(S, 'phase'),      S.phase    = 0;                     end
if ~isfield(S, 'prefix'),     S.prefix   = '';                    end
if ~isfield(S, 'frequencies') || isempty(S.frequencies)
    S.frequencies = 1:48;
end
if ~isfield(S, 'method')
    S.method = 'morlet';
    S.settings = [];
end

D = spm_eeg_load(S.D);

chanind = D.selectchannels(S.channels);

if isempty(chanind)
    error('No channels selected.');
end

if isfield(S, 'settings')
    S1 = S.settings;
else
    S1 = [];
end

timeind = D.indsample(1e-3*min(S.timewin)):D.indsample(1e-3*max(S.timewin));

% remove uninetended non-uniformities in the frequency axis
frequencies = S.frequencies;

if length(frequencies)>1
    df = unique(diff(frequencies));
    if length(df) == 1 || (max(diff(df))/mean(df))<0.1
        df = mean(diff(frequencies));
        frequencies = (0:(length(frequencies)-1))*df + frequencies(1);
    end
end

S1.frequencies = frequencies;

if ~isequal(D.type, 'continuous')

    %-Run the analysis on all trials
    %--------------------------------------------------------------------------
    spm_progress_bar('Init', D.ntrials, 'trials done');
    if D.ntrials > 100, Ibar = floor(linspace(1, D.ntrials, 100));
    else Ibar = 1:D.ntrials; end

    for k = 1:D.ntrials
        spm_progress_bar('Set','ylabel','reading and computing...');
        trial = feval(['spm_eeg_specest_' S.method], S1, D(chanind, timeind, k), D.time(timeind));

        if k == 1
            spm_progress_bar('Set','ylabel','initialising output...');
            if isfield(trial, 'fourier')
                outdata = trial.fourier;
            else
                outdata = trial.pow;
            end

            Nchannels = size(outdata, 1);
            Nfrequencies = size(outdata, 2);
            Nsamples = size(outdata, 3);

            %-Generate output datasets
            %--------------------------------------------------------------------------
            Dtf = clone(D, [S.prefix 'tf_' D.fname], [Nchannels Nfrequencies Nsamples D.ntrials], 0, 1);
            Dtf = Dtf.frequencies(':', trial.freq);
            Dtf = timeonset(Dtf, trial.time(1));
            Dtf = fsample(Dtf, 1/diff(trial.time(1:2)));
            Dtf = transformtype(Dtf, 'TF');

            Dtf = chanlabels(Dtf, 1:Nchannels, D.chanlabels(chanind));
            Dtf = badchannels(Dtf, 1:Nchannels, D.badchannels(chanind));
            Dtf = chantype(Dtf, 1:Nchannels, D.chantype(chanind));
            Dtf = coor2D(Dtf, 1:Nchannels, coor2D(D,chanind));

            if S.phase && isfield(trial, 'fourier')
                Dtph = clone(Dtf, [S.prefix 'tph_' D.fname]);
                Dtph = transformtype(Dtph, 'TFphase');
            else
                if ~isfield(trial, 'fourier')
                    warning('Phase cannot be estimated with the requested method. Estimating power only.');
                end

                Dtph = [];
            end

        end

        spm_progress_bar('Set','ylabel','writing...');
        if isfield(trial, 'fourier')
            Dtf(:, :, :, k) = trial.fourier.*conj(trial.fourier);

            if S.phase
                Dtph(:, :, :, k) = angle(trial.fourier);
            end
        elseif isfield(trial, 'pow')
            Dtf(:, :, :, k) = trial.pow;
        else
            error('The plug-in returned unexpected output');
        end

        if ismember(k, Ibar), spm_progress_bar('Set', k); end
    end

    spm_progress_bar('Clear');

else % by channel for continuous data
    Nchannels = length(chanind);

    spm_progress_bar('Init', Nchannels , 'channels done');
    if Nchannels > 100, Ibar = floor(linspace(1, Nchannels, 100));
    else Ibar = 1:Nchannels; end

    for k = 1:Nchannels
        spm_progress_bar('Set','ylabel','reading and computing...');
        trial = feval(['spm_eeg_specest_' S.method], S1, D(chanind(k), timeind), D.time(timeind));

        if k == 1
            spm_progress_bar('Set','ylabel','initialising output...');
            if isfield(trial, 'fourier')
                outdata = trial.fourier;
            else
                outdata = trial.pow;
            end

            Nfrequencies = size(outdata, 2);
            Nsamples = size(outdata, 3);

            %-Generate output datasets
            %--------------------------------------------------------------------------
            Dtf = clone(D, [S.prefix 'tf_' D.fname], [Nchannels Nfrequencies Nsamples 1]);
            Dtf = Dtf.frequencies(':', trial.freq);
            Dtf = timeonset(Dtf, trial.time(1));
            Dtf = fsample(Dtf, 1/diff(trial.time(1:2)));
            Dtf = transformtype(Dtf, 'TF');

            Dtf = chanlabels(Dtf, 1:Nchannels, D.chanlabels(chanind));
            Dtf = badchannels(Dtf, 1:Nchannels, D.badchannels(chanind));
            Dtf = chantype(Dtf, 1:Nchannels, D.chantype(chanind));
            Dtf = coor2D(Dtf, 1:Nchannels, coor2D(D,chanind));

            ev  = Dtf.events;

            if ~isempty(ev)
                ev  = ev([ev.time]>=Dtf.time(1) & [ev.time]<=Dtf.time(end));
                Dtf = events(Dtf, 1, ev);
            end

            if S.phase && isfield(trial, 'fourier')
                Dtph = clone(Dtf, [S.prefix 'tph_' D.fname]);
                Dtph = transformtype(Dtph, 'TFphase');
            else
                if ~isfield(trial, 'fourier')
                    warning('Phase cannot be estimated with the requested method. Estimating power only.');
                end

                Dtph = [];
            end
        end

        spm_progress_bar('Set','ylabel','writing...');
        if isfield(trial, 'fourier')
            Dtf(k, :, :, 1) = trial.fourier.*conj(trial.fourier);

            if S.phase
                Dtph(k, :, :, 1) = angle(trial.fourier);
            end
        elseif isfield(trial, 'pow')
            Dtf(k, :, :, 1) = trial.pow;
        else
            error('The plug-in returned unexpected output');
        end

        if ismember(k, Ibar), spm_progress_bar('Set', k); end
    end

    spm_progress_bar('Clear');

end
%-Save new M/EEG dataset(s)
%--------------------------------------------------------------------------
Dtf = Dtf.history('spm_eeg_tf', S);
save(Dtf);
if ~isempty(Dtph)
    Dtph = Dtph.history('spm_eeg_tf', S);
    save(Dtph);
end

%-Cleanup
%--------------------------------------------------------------------------
fprintf('%-40s: %30s\n','Completed',spm('time'));                       %-#
spm('FigName','M/EEG Time Frequency: done'); spm('Pointer','Arrow');