text:
classdef TestDft
    %TestDft

    methods (Static)
        function test_1d_real
            % 1x100 real signal
            x = sin(linspace(0, 6*pi, 100));

            % 1D DFT
            y = cv.dft(x);
            validateattributes(y, {class(x)}, {'vector', 'numel',numel(x)});

            % inverse transform
            xx = cv.dft(y, 'Inverse',true, 'Scale',true);
            validateattributes(xx, {class(x)}, {'vector', 'numel',numel(x)});

            err = norm(x - xx);
        end

        function test_1d_complex
            % 1x100 complex signal
            x = sin(linspace(0, 6*pi, 100)); % real
            x(:,:,2) = 0;                    % imag

            % 1D DFT
            y = cv.dft(x, 'ComplexInput',true);
            validateattributes(y, {class(x)}, {'size',size(x)});

            % inverse transform
            xx = cv.dft(y, 'Inverse',true, 'Scale',true);
            validateattributes(xx, {class(x)}, {'size',size(x)});

            err = norm(x(:) - xx(:));
        end

        function test_2d
            % image (256x256)
            img = cv.imread(fullfile(mexopencv.root(), 'test', 'blox.jpg'), ...
                'Grayscale',true);
            I = double(img);

            % 2D DFT
            J = cv.dft(I, 'ComplexOutput',true);
            validateattributes(J, {class(I)}, {'size',[size(I) 2]});

            % 2D DFT is simply the 1D DFT, performed along the rows and then
            % along the columns: fft2(I) = fft(fft(I).').'
            J2 = cv.dft(I, 'Rows',true, 'ComplexOutput',true);
            J2 = permute(J2, [2 1 3]);
            J2 = cv.dft(J2, 'Rows',true, 'ComplexOutput',true);
            J2 = permute(J2, [2 1 3]);
            %assert(norm(J(:) - J2(:)) < 1e-9);

            % inverse transform
            II = cv.dft(J, 'Inverse',true, 'Scale',true, 'RealOutput',true);
            validateattributes(II, {class(I)}, {'2d', 'size',size(I)});
            err = norm(I - II);
        end

        function test_ccs_1d
            % real 1D vector
            if true
                N = 100; % even case
            else
                N = 101; % odd case
            end
            x = sin(linspace(0, 6*pi, N));
            x = x + randn(size(x));

            % 1D DFT (packed CCS format)
            y = cv.dft(x);

            % un-pack CCS into full complex array (both even/odd length cases)
            if rem(N, 2) == 0
                dc = y(1);           % DC component
                pos = y(2:2:end-1);  % positive frequencies (real part)
                neg = y(3:2:end-1);  % negative frequencies (imaginary part)
                nyq = y(end);        % nyquist frequency
                yy = cat(3, ...
                    [dc, pos, nyq, fliplr(pos)], ...
                    [0, neg, 0, fliplr(-neg)]);
            else
                dc = y(1);
                pos = y(2:2:end);
                neg = y(3:2:end);
                nyq = [];  % no nyquist
                yy = cat(3, ...
                    [dc, pos, fliplr(pos)], ...
                    [0, neg, fliplr(-neg)]);
            end

            % compare
            y2 = cv.dft(x, 'ComplexOutput',true);
            assert(isequal(y2, yy));

            % compare against builtin FFT function
            yy = complex(yy(:,:,1), yy(:,:,2));
            y3 = fft(x);
            %assert(norm(y3 - yy) < 1e-9);
        end

        function test_ccs_2d
            %TODO
            error('mexopencv:testskip', 'todo');
        end

        function test_error_argnum
            try
                cv.dft();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function M=find_motif34(m,n)
%FIND_MOTIF34       Motif legend
%
%   Motif_matrices = find_motif34(Motif_id,Motif_class);
%   Motif_id = find_motif34(Motif_matrix);
%
%   This function returns all motif isomorphs for a given motif id and
%   class (3 or 4). The function also returns the motif id for a given
%   motif matrix
%
%   1. Input:       Motif_id,           e.g. 1 to 13, if class is 3
%                   Motif_class,        number of nodes, 3 or 4.
%
%      Output:      Motif_matrices,     all isomorphs for the given motif
%
%   2. Input:       Motif_matrix        e.g. [0 1 0; 0 0 1; 1 0 0]
%
%      Output       Motif_id            e.g. 1 to 13, if class is 3
%
%
%Mika Rubinov, UNSW, 2007-2008

persistent M3 ID3 M4 ID4

if isscalar(m)
    if n==3
        if isempty(ID3);
            load motif34lib M3 ID3;
        end
        ind=find(ID3==m).';
        M=zeros(3,3,length(ind));
        for i=1:length(ind)
            M(:,:,i)=reshape([0 M3(ind(i),1:3) 0 ...
                M3(ind(i),4:6) 0],3,3);
        end
    elseif n==4
        if isempty(ID4);
            load motif34lib M4 ID4;
        end
        ind=find(ID4==m).';
        M=zeros(4,4,length(ind));
        for i=1:length(ind)
            M(:,:,i)=reshape([0 M4(ind(i),1:4) 0 ...
                M4(ind(i),5:8) 0 M4(ind(i),9:12) 0],4,4);
        end
    end
else
    n=size(m,1);
    M=eval(['find(motif' int2str(n) 'struct_bin(m))']);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ xdk, ddk ] = dif_derivk_table ( nd, xd, dd, k )

%*****************************************************************************80
%
%% DIF_DERIVK_TABLE computes the divided difference table for K-th derivative.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carl deBoor,
%    A Practical Guide to Splines,
%    Springer, 2001,
%    ISBN: 0387953663,
%    LC: QA1.A647.v27.
%
%  Parameters:
%
%    Input, integer ND, the size of the input table.
%
%    Input, real XD(ND), the abscissas for the divided
%    difference table.
%
%    Input, real DD(ND), the divided difference table.
%
%    Input, integer K, the index of the derivative.
%    0 <= K
%
%    Input, real XDK(ND-K), the abscissas for the divided
%    difference table for the derivative.
%
%    Output, real DDK(ND-K), the divided difference
%    table for the derivative.
%
  if ( k < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DIF_DERIVK_TABLE - Fatal error!\n' );
    fprintf ( 1, '  K < 0.\n' );
    error ( 'DIF_DERIVK_TABLE - Fatal error!' )
  end

  if ( nd <= k )
    ndk = 0;
    xdk = [];
    ydk = [];
    return
  end
%
%  Shift the abscissas to zero.
%
  ndk = nd;

  xd_temp(1:ndk) = xd(1:ndk);
  dd_temp(1:ndk) = dd(1:ndk);

  [ xd_temp, dd_temp ] = dif_shift_zero ( nd, xd_temp, dd_temp );
%
%  Repeatedly differentiate.
%
  for j = 1 : k

    ndk = ndk - 1;

    for i = 1 : ndk
      dd_temp(i) = i * dd_temp(i+1);
    end

  end

  ddk(1:ndk) = dd_temp(1:ndk);
  xdk(1:ndk) = 0.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function r=interleave(x,y)
% interleave two vectors, x and y

xc=num2cell(x);
yc=num2cell(y);
r=cell(1,numel(x)+numel(y));
r(1:2:2*numel(x))=xc;
r(2:2:2*numel(y))=yc;
r=[r{:}];

repo: REST


text:
function add_legend_entry(plot_obj,lh,legend_string,append_name)
% A private method that stores legend information for when
% 'display_legend' is called.
% Inputs:
%   lh: The handle for the curve.
%   legend_string: A string describing the curve.
%   append_name: A boolean indicating whether to prepend the system
%     name (set in 'set_system') to the 'legend_string'.

assert(ischar(legend_string))

if ~strcmp(legend_string,'')
    plot_obj.handles_vec = [plot_obj.handles_vec lh];
    if exist('append_name','var') && append_name
	plot_obj.legend_strings = {plot_obj.legend_strings{:},[plot_obj.sys_name ' ' legend_string]};
    else
	plot_obj.legend_strings = {plot_obj.legend_strings{:},legend_string};
    end
end
end

author: nesl


text:
function point_num = levels_index_size_own ( dim_num, level_max )

%*****************************************************************************80
%
%% LEVELS_INDEX_SIZE_OWN sizes a sparse grid made from OWN 1D rules.
%
%  Discussion:
%
%    The sparse grid is presumed to have been created from products
%    of OPEN WEAKLY NESTED 1D quadrature rules.
%
%    OWN rules include Gauss Hermite and Gauss Legendre.
%
%    The sparse grid is the logical sum of product grids with total LEVEL
%    between LEVEL_MIN and LEVEL_MAX.
%
%    Oddly enough, in order to count the number of points, we will
%    behave as though LEVEL_MIN was zero.  This is because our computation
%    concentrates on throwing away all points generated at lower levels,
%    but, in fact, if we start at a nonzero level, we need to include
%    on that level all the points that would have been generated on lower
%    levels.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 July 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Output, integer POINT_NUM, the number of points in the grid.
%

%  Special case.
%
  if ( level_max == 0 )
    point_num = 1;
    return
  end
%
%  The outer loop generates LEVELs from LEVEL_MIN to LEVEL_MAX.
%
%  The normal definition of LEVEL_MIN:
%
%   level_min = max ( 0, level_max + 1 - dim_num )
%
%  Our somewhat artificial temporary local definition of LEVEL_MIN:
%
  if ( dim_num == 1 )
    level_min = level_max;
    point_num = 1;
  else
    level_min = 0;
    point_num = 0;
  end

  for level = level_min : level_max
%
%  The middle loop generates the next partition that adds up to LEVEL.
%
    level_1d = [];
    more = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );
%
%  Transform each 1D level to a corresponding 1D order.
%
      order_1d = level_to_order_open ( dim_num, level_1d );

      for dim = 1 : dim_num
%
%  Account for the repetition of the center point.
%
        if ( 1 < order_1d(dim) )
          order_1d(dim) = order_1d(dim) - 1;
        end

      end

      point_num = point_num + prod ( order_1d(1:dim_num) );

      if ( ~more )
        break
      end

    end
  end

  return
end


text:
%%*******************************************************************
%%  Read in a problem in SDPA sparse format.
%%
%%  [blk,At,C,b] = read_sdpa(fname)
%%
%%  Input: fname = name of the file containing SDP data in
%%                 SDPA foramt.
%%  Important: the data is assumed to contain only
%%             semidefinite and linear blocks.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%******************************************************************

function [blk,At,C,b] = read_sdpa(fname)

%%
%%  Open the file for input
%%
compressed = 0;
if exist(fname,'file')
    fid = fopen(fname,'r');
elseif exist([fname,'.Z'],'file');
    compressed = 1;
    unix(['uncompress ',fname,'.Z']);
    fid = fopen(fname,'r');
elseif exist([fname,'.gz'],'file');
    compressed = 2;
    unix(['gunzip ',fname,'.gz']);
    fid = fopen(fname,'r');
else
    fprintf('*** Problem not found, please specify the correct path or problem name. \n');
    blk = []; At = []; C = []; b = [];
    return;
end
%%
%%  Clean up special characters and comments from the file
%%
[datavec,count] = fscanf(fid,'%c'); %#ok
linefeeds = strfind(datavec,char(10));
comment_chars = '*"=';
cumidx = [];
for i=1:length(comment_chars)
    idx = strfind(datavec,comment_chars(i));
    cumidx = [cumidx,idx]; %#ok
end
for j=length(cumidx):-1:1
    if (cumidx(j)==1) || (strcmp(datavec(cumidx(j)-1),char(10)))
        datavec(cumidx(j):linefeeds(find(cumidx(j)<linefeeds,1,'first')))='';
    else
        datavec(cumidx(j):linefeeds(find(cumidx(j)<linefeeds,1,'first'))-1)='';
    end
end
special_chars = ',{}()';
cumidx=[];
for i=1:length(special_chars)
    idx = strfind(datavec,special_chars(i));
    cumidx = [cumidx,idx]; %#ok
end
datavec(cumidx) = blanks(length(cumidx));
clear linefeeds;
%%
%% Close the file
%%
fclose('all');
if compressed==1; unix(['compress ',fname]); end;
if compressed==2; unix(['gzip ',fname]); end;
%%
%%  Next, read in basic problem size parameters.
%%
datavec = sscanf(datavec,'%f');
if size(datavec,1) < size(datavec,2); datavec = datavec'; end;
m = datavec(1);
numblk  = datavec(2);
blksize = datavec(3:numblk+2);
if size(blksize,1) > size(blksize,2); blksize = blksize'; end
%%
%% Get input  b.
%%
idxstrt = 2+numblk;
b = datavec(idxstrt+1:idxstrt+m);
idxstrt = idxstrt+m;
b = -b;
%%
%% Construct blk
%%
deblksize = 100;
spblkidxtmp = find( (blksize>1) & (blksize < deblksize) );
spblkidxtmp = sort(spblkidxtmp);
deblkidx = find( (blksize<=1) | (blksize >= deblksize) );
denumblk = length(deblkidx);
linblkidx = zeros(1,denumblk);
for p = 1:denumblk
    n = blksize(deblkidx(p));
    if (n > 1);
        blk{p,1} = 's'; blk{p,2} = n; %#ok
        n2 = n*(n+1)/2;
        At{p,1} = sparse(n2,m); %#ok
        C{p,1} = sparse(n,n); %#ok
    else
        linblkidx(p) = p;
        blk{p,1} = 'l'; blk{p,2} = abs(n); %#ok
        At{p,1} = sparse(abs(n),m); %#ok
        C{p,1} = sparse(abs(n),1); %#ok
    end
end
if ~isempty(spblkidxtmp)
    maxnumblk = 200;
    spnumblk = ceil(length(spblkidxtmp)/maxnumblk);
    for q = 1:spnumblk
        if (q < spnumblk)
            spblkidxall{q} = spblkidxtmp([(q-1)*maxnumblk+1: q*maxnumblk]); %#ok
        else
            spblkidxall{q} = spblkidxtmp([(q-1)*maxnumblk+1: length(spblkidxtmp)]); %#ok
        end
        tmp = blksize(spblkidxall{q});
        blk{denumblk+q,1} = 's'; %#ok
        blk{denumblk+q,2} = tmp; %#ok
        n2 = sum(tmp.*(tmp+1))/2;
        At{denumblk+q,1} = sparse(n2,m); %#ok
        C{denumblk+q,1} = sparse(sum(tmp),sum(tmp)); %#ok
    end
else
    spnumblk = 0;
end
linblkidx(denumblk+1:denumblk+spnumblk) = 0;
%%
%% Construct single blocks of A,C
%%
len = length(datavec);
Y = reshape(datavec(idxstrt+1:len),5,(len-idxstrt)/5)';
clear datavec;
Y = sortrows(Y,[1 2]);
matidx = [0; find(diff(Y(:,1)) ~= 0); size(Y,1)];
%%
for k = 1:length(matidx)-1
    idx = matidx(k)+1 : matidx(k+1);
    Ytmp  = Y(idx,1:5);
    matno = Ytmp(1,1);
    Ytmp2 = Ytmp(:,2);
    for p = 1:denumblk
        n  = blksize(deblkidx(p));
        idx = find(Ytmp2 == deblkidx(p));
        ii = Ytmp(idx,3); jj = Ytmp(idx,4); vv =Ytmp(idx,5);
        len = length(idx);
        if (n > 1)
            idxtmp = find(ii > jj);
            if ~isempty(idxtmp);
                tmp = jj(idxtmp);
                jj(idxtmp) = ii(idxtmp); ii(idxtmp) = tmp;
            end
            tmp = -sparse(ii,jj,vv,n,n);
            tmp = tmp + triu(tmp,1)';
        else
            tmp = -sparse(ii,ones(len,1),vv,abs(n),1);
        end
        if (matno == 0)
            C{p,1} = tmp; %#ok
        else
            if (n > 1)
                At{p,1}(:,matno) = svec(blk(p,:),tmp,1); %#ok
            else
                At{p,1}(:,matno) = tmp; %#ok
            end
        end
    end
end
%%
%% Construct big sparse block of A,C
%%
if (spnumblk > 0)
    Y1 = Y(:,1);
    diffY1 = find(diff([-1; Y1; inf]));
    for kk = 1:length(diffY1)-1
        idx = diffY1(kk) : diffY1(kk+1)-1;
        matno = Y1(diffY1(kk));
        Ytmp  = Y(idx,1:5);
        Ytmp2 = Ytmp(:,2);
        maxYtmp2  = Ytmp2(length(Ytmp2));
        minYtmp2  = Ytmp2(1);
        diffYtmp2 = Ytmp2(diff([-1; Ytmp2])~=0);
        for q = 1:spnumblk
            spblkidx    = spblkidxall{q};
            maxspblkidx = spblkidx(length(spblkidx));
            minspblkidx = spblkidx(1);
            count = 0;
            if (minYtmp2 <= maxspblkidx) && (maxYtmp2 >= minspblkidx)
                tmpblksize = blksize(spblkidx);
                n = sum(tmpblksize);
                cumspblksize = [0 cumsum(tmpblksize)];
                n2 = sum(tmpblksize.*(tmpblksize+1))/2;
                idx = zeros(n2,1); offset = zeros(n2,1);
                for t = 1:length(diffYtmp2)
                    p = find(spblkidx == diffYtmp2(t));
                    if ~isempty(p)
                        idxtmp = find(Ytmp2 == spblkidx(p));
                        len = length(idxtmp);
                        idx(count+1:count+len) = idxtmp;
                        offset(count+1:count+len) = cumspblksize(p);
                        count = count + len;
                    end
                end
                idx = idx(1:count); offset = offset(1:count);
                ii = Ytmp(idx,3)+offset; jj = Ytmp(idx,4)+offset; vv = Ytmp(idx,5);
                idxtmp = find(ii > jj);
                if ~isempty(idxtmp);
                    tmp = jj(idxtmp);
                    jj(idxtmp) = ii(idxtmp); ii(idxtmp) = tmp;
                end
                idxeq = find(ii==jj);
                tmp = spconvert([ii jj -vv; jj ii -vv; n n 0]) ...
                    + spconvert([ii(idxeq) jj(idxeq) vv(idxeq); n n 0]);
                if (matno == 0)
                    C{denumblk+q,1} = tmp; %#ok
                else
                    At{denumblk+q,1}(:,matno) = svec(blk(denumblk+q,:),tmp,1); %#ok
                end
            end
        end
    end
end
%%
%% put all linear blocks together as a single linear block
%%
idx = find(linblkidx);
if (length(idx) > 1)
    sdpidx = find(linblkidx==0);
    blktmp = 0; Atmp = []; Ctmp = [];
    for k = 1:length(idx)
        tmp = linblkidx(idx(k));
        blktmp = blktmp+blk{tmp,2};
        Atmp = [Atmp; At{tmp}]; %#ok
        Ctmp = [Ctmp; C{tmp}]; %#ok
    end
    At = At(sdpidx); C = C(sdpidx); blk = blk(sdpidx,:);
    len = length(sdpidx);
    blk(2:len+1,:) = blk;
    blk{1,1} = 'l'; blk{1,2} = blktmp;
    At(2:len+1,1) = At; C(2:len+1,1) = C;
    At{1,1} = Atmp; C{1,1} = Ctmp;
end
%%******************************************************************

author: yu-jiang

repo: radpbook

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


text:
% +XYZFUN
%
% Contains Functions that operate on a grid in 3D (XYZ) space
%
%
% Files
%   bgrid3dB - This subroutine assigns creates a 3D grid with
%   zgrid3d  -

author: CelsoReyes


text:
function [B]=gen_CoriolisMatrix(q,Pcii,Icii,mcii,dq)
 %% This function is used to calculate Coriolis matrix of the KUKA iiwa 7 R 800 manipulator

% Arreguments:
%--------------------
% q: is 1x7 vector, joint nagles vector of the manipulator
% dq: is 1x7 vector, joint angular velocity vector
% Pcii: is 3X7 matrix while each column represents the local coordinates
% of the center of mass of each link.
% Icii: is (3x3x7) matrix, each 3x3 matrix of which represnets the
% associated link inertial tensor represented in its local inertial frame
% mcii: is (1x7) vector, each element of which specifies a mass of one of
% the links

% Return value:
%--------------------
% M: 7x7 matrix, the mass matrix of the KUKA iiwa 7 R 800.

% Copyright: Mohammad SAFEEA, 9th-April-2018

T=kukaGetKenimaticModelAccelerated(q);
B=GetCoriolisMatrix1(T,Pcii,Icii,mcii,dq);

repo: KST-Kuka-Sunrise-Toolbox


text:
function x = inv_posdef(A)
% INV_POSDEF        Invert positive definite matrix.
% INV_POSDEF(A) is like INV(A) but faster and more numerically stable.
% See test_inv_posdef for a timing test.

% Written by Tom Minka

U = cholproj(A);
iU = inv_triu(U);
x = iU*iU';

author: andrewssobral

repo: lrslibrary


text:
% Add the ImageGraph to path (in a folder named 'ImageGraphs'). Find it here:
% http://www.mathworks.com/matlabcentral/fileexchange/53614-image-graphs
addpath(fullfile(fileparts(mfilename('fullpath')), 'ImageGraphs'));

%% Read the image and features from the sample file
image = im2double(imread('docia.png'));
features = image(:, size(image, 2) / 2 + 1 : end, :);
image = image(:, 1 : size(image, 2) / 2, :);

% The eigendecomposition uses a lot of memory and may render the computer
% unresponsive, so better to test it first with a small image.
image = imresize(image, 0.5);
features = imresize(features, 0.5);

%% Semantic soft segmentation
% This function outputs many intermediate variables, if needed.
% The results may vary a bit from run to run, as there are 2 stages that use
% k-means for intialization & grouping.
sss = SemanticSoftSegmentation(image, features);

% To use the features generated using our network implementation,
% just feed them as the 'features' variable to the function. It will do
% the prepocessing described in the paper and give the processed
% features as an output.
% If you are dealing with many images, storing the features after
% preprocessing is recommended as raw hyperdimensional features
% take a lot of space. Check the 'preprocessFeatures.m' file.

% Visualize
figure; imshow([image features visualizeSoftSegments(sss)]);
title('Semantic soft segments');

% There's also an implementation of Spectral Matting included
sm = SpectralMatting(image);
% You can group the soft segments from Spectral Matting using
% semantic features, the way we presented our comparisons in the paper.
sm_gr = groupSegments(sm, features);
figure; imshow([image visualizeSoftSegments(sm) visualizeSoftSegments(sm_gr)]);
title('Matting components');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% set system limits
system = mr.opts('MaxGrad', 10, 'GradUnit', 'mT/m', ...
    'MaxSlew', 50, 'SlewUnit', 'T/m/s', 'rfRingdownTime', 20e-6, ...
    'rfDeadTime', 100e-6, 'adcDeadTime', 10e-6);

% high-level sequence parameters
alpha=90;
sliceThickness=50e-3;
Nx = 2048;
Nrep = 1;
TE=30e-3;
TR = 6000e-3;

rf_90_duration=2.6e-3; % excitation pulse duration
rf_180_duration= 4.6e-3; % refocusing pulse duration

% create and populate the sequence
seq = mr.Sequence(system);              % Create a new sequence object

% Create excitation pulse and gradient
rf_90= mr.makeSLRpulse(pi/2,'duration',rf_90_duration,'timeBwProduct',7.88,'dwell',rf_90_duration/500,'passbandRipple',1,'stopbandRipple',1e-2,'filterType','ms','system',system);

%Create refocusing pulse and gradient
rf_180_1 = mr.makeAdiabaticPulse('wurst','duration',rf_180_duration,'bandwidth',6000,'dwell',rf_180_duration/500,'n_fac',20,'use','refocusing','system',system);

timebwproduct_90=8;   % "experimental" value for the used for SLR excitation
timebwproduct_180=22; % "experimental" value for the used refocusing pulse

grad_amplitude_90 = (timebwproduct_90/2.6e-3)/sliceThickness; %BW/thickness
grad_amplitude_180 = (timebwproduct_180/4.6e-3)/sliceThickness; %BW/thickness

% we cannot use mr.calcDuration(rf_90) and mr.calcDuration(rf_180_1) here because this already includes delays and RF dead times
gx = mr.makeTrapezoid('x','flatTime',rf_90_duration,'amplitude',grad_amplitude_90,'system',system);
gy = mr.makeTrapezoid('y','flatTime',rf_180_duration,'amplitude',grad_amplitude_180,'system',system);
gz = mr.makeTrapezoid('z','flatTime',rf_180_duration,'amplitude',grad_amplitude_180,'system',system);
gz_2 = mr.makeTrapezoid('z','flatTime',rf_180_duration,'amplitude',grad_amplitude_180,'system',system);

% refocusing area for the 90-degree pulse (needs to be remembered in the "refocusing budget")
area_90_toRefocus=gx.amplitude*(0.5*gx.fallTime+gx.flatTime-mr.calcRfCenter(rf_90));

% gradient spoiling
spoilMoment=10/sliceThickness;
gzSpoil=mr.makeTrapezoid('z','area',spoilMoment,'system',system);
gxSpoil=mr.makeTrapezoid('x','duration',2*mr.calcDuration(gzSpoil),'area',spoilMoment,'system',system);
gzSpoil_2=mr.makeTrapezoid('z','duration',2*mr.calcDuration(gzSpoil),'area',spoilMoment,'system',system);

% split x gradient at the end of the plato and calculate the rf_90 delay
gx_parts=mr.splitGradientAt(gx,mr.calcDuration(gx)-gx.fallTime+system.rfRingdownTime);
[gx_p1,rf_90,~]=mr.align('right',gx_parts(1),rf_90,mr.makeDelay(system.rfDeadTime+rf_90_duration+system.rfRingdownTime));
gx_parts(2).delay=0; % set the delay to 0 because it will be used in a separate block

% split y gradient at the end of the plato and calculate the rf_180_1 delay
gy_parts=mr.splitGradientAt(gy,mr.calcDuration(gy)-gy.fallTime+system.rfRingdownTime);
rf_180_1.delay=max(mr.calcDuration(gy_parts(1))-rf_180_duration-system.rfRingdownTime, mr.calcDuration(gzSpoil));
assert(rf_180_1.delay>=system.rfDeadTime);
gy_p1=gy_parts(1);
gy_p1.delay=rf_180_1.delay+rf_180_duration+system.rfRingdownTime-mr.calcDuration(gy_p1);
gy_parts(2).delay=0; % set the delay to 0 because it will be used in a separate block

% combine y gradient with delays to a new gradient the rf_180_2 delay
rf_180_2=rf_180_1;
rf_180_2.delay=max(mr.calcDuration(gxSpoil),mr.calcDuration(gzSpoil_2));
gy_tmp=mr.splitGradientAt(gy,mr.calcDuration(gy)-gy.fallTime);
gy_tmp(1).delay=rf_180_2.delay+rf_180_duration-mr.calcDuration(gy_tmp(1));
gySpoil=mr.makeExtendedTrapezoidArea('y',gy.amplitude,0,spoilMoment+0.5*gy.amplitude*gy.fallTime,system);
gySpoil.delay=mr.calcDuration(gy_tmp(1));
gy_comb=mr.addGradients({gy_parts(2),gy_tmp(1),gySpoil},'system',system);
gy_comb_parts=mr.splitGradientAt(gy_comb,rf_180_2.delay+rf_180_duration+system.rfRingdownTime);%
gy_comb_parts(2).delay=0;

%
gxSpoil_2=mr.makeTrapezoid('x','area',spoilMoment,'system',system);
rf_180_3=rf_180_1;
rf_180_3.delay=mr.calcDuration(gxSpoil_2);
gz.delay=rf_180_3.delay-gz.riseTime;

%Additional Spoiler gradients
gzSpoil_semiFinal=mr.makeExtendedTrapezoidArea('z',0,gz.amplitude,spoilMoment,system);
gxSpoil_semiFinal=mr.makeTrapezoid('x','area',spoilMoment+area_90_toRefocus,'system',system);
gySpoil_semiFinal=mr.makeTrapezoid('y','area',2*spoilMoment,'system',system);
[gxSpoil_semiFinal,gySpoil_semiFinal,gzSpoil_semiFinal]=mr.align('right',gxSpoil_semiFinal,gySpoil_semiFinal,gzSpoil_semiFinal);

rf_180_4=rf_180_1;
rf_180_4.delay=0;
gz_temp=mr.splitGradientAt(gz_2,gz_2.riseTime);
gz_temp(2).delay=0;
gz_parts=mr.splitGradientAt(gz_temp(2),rf_180_4.delay+rf_180_duration);
gz_parts(1).delay=mr.calcDuration(gzSpoil_semiFinal);
rf_180_4.delay=mr.calcDuration(gzSpoil_semiFinal);

gzSpoil_Final=mr.makeExtendedTrapezoidArea('z',gz.amplitude,0,spoilMoment,system);
gxSpoil_Final=mr.makeTrapezoid('x','area',spoilMoment,'system',system);
gySpoil_Final=mr.makeTrapezoid('y','area',spoilMoment,'system',system);

gzSpoil_Final.delay=rf_180_4.delay+rf_180_duration;
gz_comb=mr.addGradients({gzSpoil_semiFinal,gz_parts(1),gzSpoil_Final},'system',system);

gxSpoil_Final.delay=rf_180_4.delay+rf_180_duration;
gySpoil_Final.delay=rf_180_4.delay+rf_180_duration;

gxSpoil_combi=mr.addGradients({gxSpoil_semiFinal,gxSpoil_Final},'system',system);
gySpoil_combi=mr.addGradients({gySpoil_semiFinal,gySpoil_Final},'system',system);

%timing calculation
lTime1=(rf_90_duration+rf_180_duration)/2+system.rfRingdownTime+mr.calcDuration(gzSpoil);
lTime2=(rf_180_duration+rf_180_duration)/2+system.rfRingdownTime+mr.calcDuration(gxSpoil);
lTime3=(rf_180_duration+rf_180_duration)/2+system.rfRingdownTime+mr.calcDuration(gxSpoil_2);

lTime4=TE/2-lTime2;
lTime5=TE/2-lTime1-lTime3;

%Define ADC events
adc = mr.makeAdc(Nx, 'Dwell', 2e-4, 'system', system);
delayTE1=lTime4-(rf_180_duration+gz.fallTime+mr.calcDuration(gySpoil_semiFinal));
delayTE2=lTime5-(rf_180_duration/2+mr.calcDuration(gySpoil_Final)-gySpoil_Final.delay)-adc.dwell/2;
adc.delay=delayTE2;

% Loop over repetitions and define sequence blocks
for i=1:Nrep
    seq.addBlock(rf_90, gx_p1);
    seq.addBlock(gzSpoil,gx_parts(2),rf_180_1,gy_p1);
    seq.addBlock(gxSpoil,gzSpoil_2,rf_180_2, gy_comb_parts(1));
    seq.addBlock(gxSpoil_2,gy_comb_parts(2),rf_180_3, gz);
    seq.addBlock(mr.makeDelay(delayTE1));
    seq.addBlock(gxSpoil_combi,gySpoil_combi,gz_comb,rf_180_4);
    seq.addBlock(adc,mr.makeDelay(mr.calcDuration(adc)+system.adcDeadTime));

    % this is realy a lazy way of defining the TR delay
    if i==1
        delayTR = TR- seq.duration();
        assert(delayTR>0);
    end

    seq.addBlock(mr.makeDelay(delayTR))
end

%% check whether the timing of the sequence is correct
[ok, error_report]=seq.checkTiming;

if (ok)
    fprintf('Timing check passed successfully\n');
else
    fprintf('Timing check failed! Error listing follows:\n');
    fprintf([error_report{:}]);
    fprintf('\n');
end
seq.setDefinition('FOV', [sliceThickness sliceThickness sliceThickness]);
seq.setDefinition('Name', 'semiLaser');

seq.write('semiLASER.seq');       % Write to pulseq file

%% do some visualizations

seq.plot('timeDisp','us','showBlocks',1,'timeRange',[0 TE*1.2]);
%seq.plot();             % Plot sequence waveforms

%% trajectory calculation
% [ktraj_adc, ktraj, t_excitation, t_refocusing, t_adc] = seq.calculateKspace();
% t_ktraj=(1:(size(ktraj,2)))*system.gradRasterTime;
[ktraj_adc, t_adc, ktraj, t_ktraj, t_excitation, t_refocusing] = seq.calculateKspacePP();
[ktraj_adc_ofs, t_adc_ofs, ktraj_ofs, t_ktraj_ofs] = seq.calculateKspacePP('gradient_offset',[1.0 .5 -1.0]*1e3); % this will help us to verify the correct echo positions

% plot "k-spaces" (gradient moments)
figure; plot(t_ktraj, ktraj'); % plot the entire k-space trajectory
hold on; plot(t_adc,ktraj_adc(1,:),'.'); % and sampling points on the kx-axis
plot(t_ktraj_ofs, ktraj_ofs');
axis([0 TE*2 -250 250]);
title('gradient moments without and with background gradients');

%% additional timing checks
% Karl Landheer et al define tau(1:5) and require tau_1+tau_2+tau_3 = tau_2+tau4
tau=diff([t_excitation t_refocusing t_adc(1)]);

if abs(sum(tau)-TE)>5e-5 % we tolerate an error of 1/2 grad rasters
    warning('TE calculation seems to be wrong, check timing!');
end

if abs(sum(tau([1 3 5]))-sum(tau([2 4])))>5e-5 % we tolerate an error of 1/2 grad rasters
    warning('spin echo condition is not fulfilled, check timing!');
end

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function [output,frequency_parameters,notes,spike_triggers,varargout] = read_Intan(fn,varargin)
% Adapted function to read intan files. Adapted from
% 'read_Intan_RHD2000file.m' as can be downloaded from the manufacuter's
% website: http://www.intantech.com/downloads.html

% This function requires user inputs, opens no dialogues and generates no
% command line output (except errors).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUTS:
% fn : Filename (including path if necessary)
% 'range' : declares that the following input is a range of data (in sec)
%           to load.
% 'ampch' : Lists the amplifier channnels to load
% 'auxch' : Lists the auxilary channels to load
% 'volt'  : Declares whether the voltage channel should be loaded or no
%           ('yes' or 'no')
% 'ADCch' : Lists the onboard ADC channels to load
% 'DIGch' : Lists the onboard digital channels to load
% 'tempch': Declares whether the temperature data should be loaded ('yes'
%           or 'no')
% 'notch' : Assign '50Hz' or '60Hz' to enable notch filtering. If left
%           undefined then no filter is applied.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OUTPUTS
% The first output is always the overview of all outputs.
% The second output is always a structure with frequency settings
% The third output is always the notes structure
% The fourth output is always the spike_triggers structure
% Then the function will provide the requested output as data (in
% timeseries format) and channel information pairs. This is always in a fixed
% order:
% 1. amplifier channels
% 2. auxiliary channels
% 3. ADC channels
% 4. Digital channels
% 5. volts
% 6. temperature
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EXAMPLES:
% [output,freqs,notes,spike_triggers,ampl_data,ampl_ch,aux_data,aux_ch]...
%    = read_Intan(fn,'auxch',[1 2 3],'ampch',[1]);
% This loads the three auxiliary streams, and channel 1 from the general
% amplifier. No range is specified, so the whole stream is loaded.
%
% [output

% [output,freqs,notes,ampl_data,ampl_channels,ADC_data,...
%   ADC_channels,temp_data,temp_ch] = read_Intan_new(fn,'ampch',...
%   [1 2 4 10:20],'tempch','yes','ADCch',[1 2],'range',[1 100.1]);
%
%

% First we parse the variable inputs to check whether the list makes sense:
    p = inputParser;
    p.addRequired('Filename',@(fn) ischar(fn));
    p.addParamValue('notch','notchfreq',@(notchfreq)...
        strcmp(notchfreq,'50Hz') || strcmp(notchfreq,'60Hz'));
    p.addParamValue('range','rng',@(rng) length(rng)==2);
    p.addParamValue('ampch','chlist',@(chlist) ~isempty(chlist));
    p.addParamValue('auxch','aux',@(aux) ~isempty(aux));
    p.addParamValue('volt','Vch', @(Vch) ...
        strcmpi(Vch,'yes') || strcmpi(Vch,'no'));
    p.addParamValue('ADCch','ADC',@(ADC) ~isempty(ADC));
    p.addParamValue('DIGch','DIG',@(DIG) ~isempty(DIG)>0);
    p.addParamValue('tempch','temp',@(temp)...
        strcmpi(temp,'yes') || strcmpi(temp,'yes'));

    p.parse(fn,varargin{:,:})

% Now we assign the values defined in the inputs to the respective
% variables and set some defaults.
nrstreamsrequested = 0;
notchfreq = 0;
channellist = 0;
range = 0;
aux = 0;
supplyvolt = 0;
ADCchannels = 0;
digitalchannels = 0;
temperature = 0;
for arg=1:2:length(varargin)
    switch (varargin{arg})
        case 'notch'
            notchfreq = varargin{arg+1};
        case 'ampch'
            channellist = varargin{arg+1};
            nrstreamsrequested=nrstreamsrequested+1;
        case 'range'
            range = varargin{arg+1};
        case 'auxch'
            aux = varargin{arg+1};
            nrstreamsrequested=nrstreamsrequested+1;
        case 'volt'
            supplyvolt = 1;
            nrstreamsrequested=nrstreamsrequested+1;
        case 'ADCch'
            ADCchannels = varargin{arg+1};
            nrstreamsrequested=nrstreamsrequested+1;
        case 'DIGch'
            digitalchannels = varargin{arg+1};
            nrstreamsrequested=nrstreamsrequested+1;
        case 'tempch'
            temperature = 1;
            nrstreamsrequested=nrstreamsrequested+1;
    end
end

% The actual start of the file opening. Here we first open the file to gain
% access to the parts in it.
    fid = fopen(fn, 'r');
    s = dir(fn);
    filesize = s.bytes;

    % Check 'magic number' at beginning of file to make sure this is an Intan
    % Technologies RHD2000 data file.
    magic_number = fread(fid, 1, 'uint32');
    if magic_number ~= hex2dec('c6912702')
        error('Unrecognized file type.');
    end

    % Read version number.
    data_file_main_version_number = fread(fid, 1, 'int16');
    data_file_secondary_version_number = fread(fid, 1, 'int16');
    % Read information of sampling rate and amplifier frequency settings.
    sample_rate = fread(fid, 1, 'single');
    dsp_enabled = fread(fid, 1, 'int16');
    actual_dsp_cutoff_frequency = fread(fid, 1, 'single');
    actual_lower_bandwidth = fread(fid, 1, 'single');
    actual_upper_bandwidth = fread(fid, 1, 'single');

    desired_dsp_cutoff_frequency = fread(fid, 1, 'single');
    desired_lower_bandwidth = fread(fid, 1, 'single');
    desired_upper_bandwidth = fread(fid, 1, 'single');

    % Enables notch filtering if requested
    notch_filter_mode = fread(fid, 1, 'int16');
    notch_filter_frequency = 0;
    if strcmp(notchfreq,'50Hz')
        notch_filter_frequency = 50;
    elseif strcmp(notchfreq,'60Hz')
        notch_filter_frequency = 60;
    end

    desired_impedance_test_frequency = fread(fid, 1, 'single');
    actual_impedance_test_frequency = fread(fid, 1, 'single');

    % Place notes in data strucure
    notes = struct( ...
    'note1', fread_QString(fid), ...
    'note2', fread_QString(fid), ...
    'note3', fread_QString(fid) );

    % If data file is from GUI v1.1 or later, see if temperature sensor data
    % was saved.
    num_temp_sensor_channels = 0;
    if ((data_file_main_version_number == 1 && data_file_secondary_version_number >= 1) ...
        || (data_file_main_version_number > 1))
        num_temp_sensor_channels = fread(fid, 1, 'int16');
    end

    % Place frequency-related information in data structure.
    frequency_parameters = struct( ...
        'amplifier_sample_rate', sample_rate, ...
        'aux_input_sample_rate', sample_rate / 4, ...
        'supply_voltage_sample_rate', sample_rate / 60, ...
        'board_adc_sample_rate', sample_rate, ...
        'board_dig_in_sample_rate', sample_rate, ...
        'desired_dsp_cutoff_frequency', desired_dsp_cutoff_frequency, ...
        'actual_dsp_cutoff_frequency', actual_dsp_cutoff_frequency, ...
        'dsp_enabled', dsp_enabled, ...
        'desired_lower_bandwidth', desired_lower_bandwidth, ...
        'actual_lower_bandwidth', actual_lower_bandwidth, ...
        'desired_upper_bandwidth', desired_upper_bandwidth, ...
        'actual_upper_bandwidth', actual_upper_bandwidth, ...
        'notch_filter_frequency', notch_filter_frequency, ...
        'desired_impedance_test_frequency', desired_impedance_test_frequency, ...
        'actual_impedance_test_frequency', actual_impedance_test_frequency );
    % Define data structure for spike trigger settings.
    spike_trigger_struct = struct( ...
        'voltage_trigger_mode', {}, ...
        'voltage_threshold', {}, ...
        'digital_trigger_channel', {}, ...
        'digital_edge_polarity', {} );

    new_trigger_channel = struct(spike_trigger_struct);
    spike_triggers = struct(spike_trigger_struct);

    % Define data structure for data channels.
    channel_struct = struct( ...
        'custom_channel_name', {}, ...
        'native_channel_name', {}, ...
        'native_order', {}, ...
        'custom_order', {}, ...
        'board_stream', {}, ...
        'chip_channel', {}, ...
        'port_name', {}, ...
        'port_prefix', {}, ...
        'port_number', {}, ...
        'electrode_impedance_magnitude', {}, ...
        'electrode_impedance_phase', {} );

    new_channel = struct(channel_struct);

    % Create structure arrays for each type of data channel.
    amplifier_channels = struct(channel_struct);
    aux_input_channels = struct(channel_struct);
    supply_voltage_channels = struct(channel_struct);
    board_adc_channels = struct(channel_struct);
    board_dig_in_channels = struct(channel_struct);

    amplifier_index = 1;
    aux_input_index = 1;
    supply_voltage_index = 1;
    board_adc_index = 1;
    board_dig_in_index = 1;

    % Read signal summary from data file header.

    number_of_signal_groups = fread(fid, 1, 'int16');

    for signal_group = 1:number_of_signal_groups
        signal_group_name = fread_QString(fid);
        signal_group_prefix = fread_QString(fid);
        signal_group_enabled = fread(fid, 1, 'int16');
        signal_group_num_channels = fread(fid, 1, 'int16');
        signal_group_num_amp_channels = fread(fid, 1, 'int16');

        if (signal_group_num_channels > 0 && signal_group_enabled > 0)
            new_channel(1).port_name = signal_group_name;
            new_channel(1).port_prefix = signal_group_prefix;
            new_channel(1).port_number = signal_group;
            for signal_channel = 1:signal_group_num_channels
                new_channel(1).custom_channel_name = fread_QString(fid);
                new_channel(1).native_channel_name = fread_QString(fid);
                new_channel(1).native_order = fread(fid, 1, 'int16');
                new_channel(1).custom_order = fread(fid, 1, 'int16');
                signal_type = fread(fid, 1, 'int16');
                channel_enabled = fread(fid, 1, 'int16');
                new_channel(1).chip_channel = fread(fid, 1, 'int16');
                new_channel(1).board_stream = fread(fid, 1, 'int16');
                new_trigger_channel(1).voltage_trigger_mode = fread(fid, 1, 'int16');
                new_trigger_channel(1).voltage_threshold = fread(fid, 1, 'int16');
                new_trigger_channel(1).digital_trigger_channel = fread(fid, 1, 'int16');
                new_trigger_channel(1).digital_edge_polarity = fread(fid, 1, 'int16');
                new_channel(1).electrode_impedance_magnitude = fread(fid, 1, 'single');
                new_channel(1).electrode_impedance_phase = fread(fid, 1, 'single');

                if (channel_enabled)
                    switch (signal_type)
                        case 0
                            amplifier_channels(amplifier_index) = new_channel;
                            spike_triggers(amplifier_index) = new_trigger_channel;
                            amplifier_index = amplifier_index + 1;
                        case 1
                            aux_input_channels(aux_input_index) = new_channel;
                            aux_input_index = aux_input_index + 1;
                        case 2
                            supply_voltage_channels(supply_voltage_index) = new_channel;
                            supply_voltage_index = supply_voltage_index + 1;
                        case 3
                            board_adc_channels(board_adc_index) = new_channel;
                            board_adc_index = board_adc_index + 1;
                        case 4
                            board_dig_in_channels(board_dig_in_index) = new_channel;
                            board_dig_in_index = board_dig_in_index + 1;
                        otherwise
                            error('Unknown channel type');
                    end
                end

            end
        end
    end

    % Summarize contents of data file.
    num_amplifier_channels = amplifier_index - 1;
    num_aux_input_channels = aux_input_index - 1;
    num_supply_voltage_channels = supply_voltage_index - 1;
    num_board_adc_channels = board_adc_index - 1;
    num_board_dig_in_channels = board_dig_in_index - 1;

    % Determine how many samples the data file contains.

    % Each data block contains 60 amplifier samples.
    bytes_per_block = 60 * 4;  % timestamp data
    bytes_per_block = bytes_per_block + 60 * 2 * num_amplifier_channels;
    % Auxiliary inputs are sampled 4x slower than amplifiers
    bytes_per_block = bytes_per_block + 15 * 2 * num_aux_input_channels;
    % Supply voltage is sampled 60x slower than amplifiers
    bytes_per_block = bytes_per_block + 1 * 2 * num_supply_voltage_channels;
    % Board analog inputs are sampled at same rate as amplifiers
    bytes_per_block = bytes_per_block + 60 * 2 * num_board_adc_channels;
    % Board digital inputs are sampled at same rate as amplifiers
    if (num_board_dig_in_channels > 0)
        bytes_per_block = bytes_per_block + 60 * 2;
    end
    % Temp sensor is sampled 60x slower than amplifiers
    if (num_temp_sensor_channels > 0)
       bytes_per_block = bytes_per_block + 1 * 2 * num_temp_sensor_channels;
    end

    % How many data blocks remain in this file?
    bytes_remaining = filesize - ftell(fid);
    num_data_blocks = bytes_remaining / bytes_per_block;

    num_amplifier_samples = 60 * num_data_blocks;
    num_aux_input_samples = 15 * num_data_blocks;
    num_supply_voltage_samples = 1 * num_data_blocks;
    num_board_adc_samples = 60 * num_data_blocks;
    num_board_dig_in_samples = 60 * num_data_blocks;

    record_time = num_amplifier_samples / sample_rate;

    % Pre-allocate memory for data.
    t_amplifier = zeros(1, num_amplifier_samples);

    amplifier_data = zeros(num_amplifier_channels, num_amplifier_samples);
    aux_input_data = zeros(num_aux_input_channels, num_aux_input_samples);
    supply_voltage_data = zeros(num_supply_voltage_channels, num_supply_voltage_samples);
    temp_sensor_data = zeros(num_temp_sensor_channels, num_supply_voltage_samples);
    board_adc_data = zeros(num_board_adc_channels, num_board_adc_samples);
    board_dig_in_data = zeros(num_board_dig_in_channels, num_board_dig_in_samples);
    board_dig_in_raw = zeros(1, num_board_dig_in_samples);
    % Read sampled data from file.

    amplifier_index = 1;
    aux_input_index = 1;
    supply_voltage_index = 1;
    board_adc_index = 1;
    board_dig_in_index = 1;

    for i=1:num_data_blocks
        t_amplifier(amplifier_index:(amplifier_index+59)) = fread(fid, 60, 'uint32');
        if (num_amplifier_channels > 0)
            amplifier_data(:, amplifier_index:(amplifier_index+59)) = fread(fid, [60, num_amplifier_channels], 'uint16')';
        end
        if (num_aux_input_channels > 0)
            aux_input_data(:, aux_input_index:(aux_input_index+14)) = fread(fid, [15, num_aux_input_channels], 'uint16')';
        end
        if (num_supply_voltage_channels > 0)
            supply_voltage_data(:, supply_voltage_index) = fread(fid, [1, num_supply_voltage_channels], 'uint16')';
        end
        if (num_temp_sensor_channels > 0)
            temp_sensor_data(:, supply_voltage_index) = fread(fid, [1, num_temp_sensor_channels], 'int16')';
        end
        if (num_board_adc_channels > 0)
            board_adc_data(:, board_adc_index:(board_adc_index+59)) = fread(fid, [60, num_board_adc_channels], 'uint16')';
        end
        if (num_board_dig_in_channels > 0)
            board_dig_in_raw(board_dig_in_index:(board_dig_in_index+59)) = fread(fid, 60, 'uint16');
        end

        amplifier_index = amplifier_index + 60;
        aux_input_index = aux_input_index + 15;
        supply_voltage_index = supply_voltage_index + 1;
        board_adc_index = board_adc_index + 60;
        board_dig_in_index = board_dig_in_index + 60;

    end

    % Make sure we have read exactly the right amount of data.
    bytes_remaining = filesize - ftell(fid);
    if (bytes_remaining ~= 0)
        %error('Error: End of file not reached.');
    end

    % Close data file.
    fclose(fid);

    % Extract digital input channels to separate variables.
    for i=1:num_board_dig_in_channels
       mask = 2^(board_dig_in_channels(i).native_order) * ones(size(board_dig_in_raw));
       board_dig_in_data(i, :) = (bitand(board_dig_in_raw, mask) > 0);
    end

    % Scale voltage levels appropriately.
    amplifier_data = 0.195 * (amplifier_data - 32768); % units = microvolts
    aux_input_data = 37.4e-6 * aux_input_data; % units = volts
    supply_voltage_data = 74.8e-6 * supply_voltage_data; % units = volts
    board_adc_data = 50.354e-6 * board_adc_data; % units = volts
    temp_sensor_data = temp_sensor_data / 100; % units = deg C

    % Check for gaps in timestamps.
    num_gaps = sum(diff(t_amplifier) ~= 1);
    if (num_gaps == 0)
    else
        fprintf(1, 'Warning: %d gaps in timestamp data found.  Time scale will not be uniform!\n', ...
            num_gaps);
    end

    % Scale time steps (units = seconds).
    t_amplifier = t_amplifier / sample_rate;
    t_aux_input = t_amplifier(1:4:end);
    t_supply_voltage = t_amplifier(1:60:end);
    t_board_adc = t_amplifier;
    t_dig_in = t_amplifier;
    t_temp_sensor = t_supply_voltage;

    % If notch filtering was requested, apply it here:
    if (notch_filter_frequency > 0)
        %disp(['Applying Notch Filter of ',notchfreq]);
        for i=1:num_amplifier_channels
            amplifier_data(i,:) = ...
                notch_filter(amplifier_data(i,:), sample_rate, notch_filter_frequency, 10);
        end
    end

    % Now that we have the whole file we select the parts we need.
    % First we select the datastreams requested, then select the data range
    % requested. With every step we check whether we do have the data in
    % the file actually. If something is requested, but not included in the
    % file, a warning is displayed in the commandline.

    if nrstreamsrequested==0    %Nothing specified, all output requested
        % Output order:
        % amplifier channels => auxiliary channels => ADC channels => Digital
        % channels => volts => temperature
        output = {...
            'ampch' [1:size(amplifier_data,1)] ...
            'auxch' [1:size(aux_input_data,1)] ...
            'ADCch' [1:size(board_adc_data,1)] ...
            'DIGch' [1:size(board_dig_in_data,1)] ...
            'volt' 'yes' ...
            'tempch' 'yes'...
            };
    else
        output = cell(1,1);
        i=1;
        if sum(cellfun(@isequal,varargin,repmat({'ampch'},size(varargin))))
            output{1,i} = 'ampch';
            output{1,i+1} = channellist;
            i=i+2;
        end
        if sum(cellfun(@isequal,varargin,repmat({'auxch'},size(varargin))))
            output{1,i} = 'auxch';
            output{1,i+1} = aux;
            i=i+2;
        end
        if sum(cellfun(@isequal,varargin,repmat({'ADCch'},size(varargin))))
            output{1,i} = 'ADCch';
            output{1,i+1} = ADCchannels;
            i=i+2;
        end
        if sum(cellfun(@isequal,varargin,repmat({'DIGch'},size(varargin))))
            output{1,i} = 'DIGch';
            output{1,i+1} = digitalchannels;
            i=i+2;
        end
        if sum(cellfun(@isequal,varargin,repmat({'volt'},size(varargin))))
            output{1,i} = 'volt';
            output{1,i+1} = supplyvolt;
            i=i+2;
        end
        if sum(cellfun(@isequal,varargin,repmat({'tempch'},size(varargin))))
            output{1,i} = 'tempch';
            output{1,i+1} = temperature;
            i=i+2;
        end
    end
    for arg=1:2:length(output)
        switch (output{arg})
            case 'ampch'
                if size(amplifier_data,1)>0
                    amplifier_channels = amplifier_channels(1,output{arg+1});
                    amplifier_data = amplifier_data(output{arg+1},:);
                    ts_amp = timeseries;
                    ts_amp.time = t_amplifier;
                    ts_amp.data = amplifier_data';
                    if range ~=0
                        ts_amp = selecttime(ts_amp,range);
                    end
                    ts_amp.DataInfo.Units = 'uV';
                    ts_amp.TimeInfo.Units = 'seconds';
                    varargout{arg} = ts_amp;
                    varargout{arg+1} = amplifier_channels;
                else
                    disp('WARNING: No amplifier channels found in file!')
                end
            case 'auxch'
                if size(aux_input_data,1)>0
                    aux_input_channels = aux_input_channels(1,output{arg+1});
                    aux_input_data = aux_input_data(output{arg+1},:);
                    ts_aux = timeseries;
                    ts_aux.time = t_aux_input;
                    ts_aux.data = aux_input_data';
                    if range~=0
                        ts_aux = selecttime(ts_aux,range);
                    end
                    ts_aux.DataInfo.Unit = 'uV';
                    ts_aux.TimeInfo.Units = 'seconds';
                    varargout{arg} = ts_aux;
                    varargout{arg+1} = aux_input_channels;
                else
                    disp('WARNING: No auxiliary inputs found in file!')
                end
            case 'ADCch'
                if size(board_adc_data,1)>0
                    board_adc_channels = board_adc_channels(1,output{arg+1});
                    board_adc_data = board_adc_data(output{arg+1},:);
                    ts_board_adc = timeseries;
                    ts_board_adc.time = t_board_adc;
                    ts_board_adc.data = board_adc_data';
                    if range ~=0
                        ts_board_adc = selecttime(ts_board_adc,range);
                    end
                    ts_board_adc.DataInfo.Unit = 'uV';
                    ts_board_adc.TimeInfo.Units = 'seconds';
                    varargout{arg} = ts_board_adc;
                    varargout{arg+1} = board_adc_channels;
                else
                    disp('WARNING: No board ADC channels found in file!');
                end
            case 'DIGch'
                if size(board_dig_in_data,1)>0
                    board_dig_in_channels = board_dig_in_channels(1,output{arg+1});
                    board_dig_in_data = board_dig_in_data(output{arg+1},:);
                    ts_dig_in = timeseries;
                    ts_dig_in.time = t_amplifier;
                    ts_dig_in.data = board_dig_in_data';
                    if range ~=0
                        ts_dig_in = selecttime(ts_dig_in,range);
                    end
                    ts_dig_in.TimeInfo.Units = 'seconds';
                    varargout{arg} = ts_dig_in;
                    varargout{arg+1} = board_dig_in_channels;
                else
                    disp('WARNING: No board digital in channels found in file!');
                end
            case 'volt'
                if size(supply_voltage_data,1)>0
                    ts_volt = timeseries;
                    ts_volt.time = t_supply_voltage;
                    ts_volt.data = supply_voltage_data';
                    if range ~=0
                        ts_volt = selecttime(ts_volt,range);
                    end
                    ts_volt.TimeInfo.Units = 'seconds';
                    varargout{arg} = ts_volt;
                    varargout{arg+1} = supply_voltage_channels;
                else
                    disp('WARNING: No supply voltage data found in file!');
                end
            case 'tempch'
                if size(temp_sensor_data,1)>0
                    ts_temp = timeseries;
                    ts_temp.time = t_temp_sensor;
                    ts_temp.data = temp_sensor_data';
                    if range ~= 0
                        ts_temp = selecttime(ts_temp,range);
                    end
                    ts_temp.TimeInfo.Units = 'seconds';
                    varargout{arg} = ts_temp;
                    varargout{arg+1} = '1';
                else
                    disp('WARNING: No temperature data found in file!');
                end
        end
    end

return

function new = selecttime(old,range)
% Enter a timeseries and a range.
% Returns a new timeseries according to range.
% Range should be a 1x2 matrix with [start end]
new = timeseries;
data = old.data(old.time>range(1,1),:);
time = old.time(old.time>range(1,1));
data = data(time<range(1,2),:);
time = time(time<range(1,2));
new.data = data;
new.time = time;

return

function a = fread_QString(fid)

% a = read_QString(fid)
%
% Read Qt style QString.  The first 32-bit unsigned number indicates
% the length of the string (in 16-bit unicode characters).  If this
% number equals 0xFFFFFFFF, the string is null.

a = '';
length = fread(fid, 1, 'uint32');
if length == hex2num('ffffffff')
    return;
end
% convert length from bytes to 16-bit unicode words
length = length / 2;

for i=1:length
    a(i) = fread(fid, 1, 'uint16');
end

return

function s = plural(n)

    % s = plural(n)
    %
    % Utility function to optionally plurailze words based on the value
    % of n.

    if (n == 1)
        s = '';
    else
        s = 's';
    end

return

function out = notch_filter(in, fSample, fNotch, Bandwidth)

    % out = notch_filter(in, fSample, fNotch, Bandwidth)
    %
    % Implements a notch filter (e.g., for 50 or 60 Hz) on vector 'in'.
    % fSample = sample rate of data (in Hz or Samples/sec)
    % fNotch = filter notch frequency (in Hz)
    % Bandwidth = notch 3-dB bandwidth (in Hz).  A bandwidth of 10 Hz is
    %   recommended for 50 or 60 Hz notch filters; narrower bandwidths lead to
    %   poor time-domain properties with an extended ringing response to
    %   transient disturbances.
    %
    % Example:  If neural data was sampled at 30 kSamples/sec
    % and you wish to implement a 60 Hz notch filter:
    %
    % out = notch_filter(in, 30000, 60, 10);

    tstep = 1/fSample;
    Fc = fNotch*tstep;

    L = length(in);

    % Calculate IIR filter parameters
    d = exp(-2*pi*(Bandwidth/2)*tstep);
    b = (1 + d*d)*cos(2*pi*Fc);
    a0 = 1;
    a1 = -b;
    a2 = d*d;
    a = (1 + d*d)/2;
    b0 = 1;
    b1 = -2*cos(2*pi*Fc);
    b2 = 1;
    out = zeros(size(in));
    out(1) = in(1);
    out(2) = in(2);
    % (If filtering a continuous data stream, change out(1) and out(2) to the
    %  previous final two values of out.)

    % Run filter
    for i=3:L
        out(i) = (a*b2*in(i-2) + a*b1*in(i-1) + a*b0*in(i) - a2*out(i-2) - a1*out(i-1))/a0;
    end
return

author: Sable

repo: mcbench-benchmarks


text:
function [to_fine, from_fine] = twoDquad_spline(x_fine,y_fine,x_knots,y_knots)
% Create 2D quadratic polynomial based spline approximation
%
% by SeHyoun Ahn, July 2016
%
% REFERENCE: To be Written
%
% PARAMTERS:
%    x_fine,y_fine = fine grid points
%    x_knots,y_knots = coarser grid to reduce to
%
% OUTPUTS:
%    to_fine = change of basis from spline basis to finer grid points
%    from_fine = change of basis from finer grid points to spline basis
%
% EXAMPLES:
%   x = linspace(-1,1,100)';
%   y = linspace(-1,1,100)';
%
%   % z = e^(-(x^2+y^2));
%   z = exp(-bsxfun(@plus,x.^2,y'.^2));
%
%   knot_x = linspace(-1,1,4)';
%   knot_y = linspace(-1,1,4)';
%
%   [to_fine, from_fine] = twoDquad_spline(x,y,knot_x,knot_y);
%
%   z_approxed = to_fine*from_fine*z(:);
%
%   surf(z);
%   hold on;
%   disp('press any key');
%   pause();
%   surf(reshape(z_approxed,100,100));
%
% SYNTAX:
% [to_fine, from_fine] = twoDquad_spline(x_fine,y_fine,x_knots,y_knots)

n_x_fine = length(x_fine);
n_y_fine = length(y_fine);

n_x_knots = length(x_knots);
n_y_knots = length(y_knots);
h_x = diff(x_knots);
h_y = diff(y_knots);

% Find the Location of Rectangle
x_locs = sum(bsxfun(@ge,x_fine,x_knots'),2);
x_locs = min(x_locs,n_x_knots-1);
x_adjusted = (x_fine - x_knots(x_locs))./h_x(x_locs);

y_locs = sum(bsxfun(@ge,y_fine,y_knots'),2);
y_locs = min(y_locs,n_y_knots-1);
y_adjusted = (y_fine - y_knots(y_locs))./h_y(y_locs);

xx = bsxfun(@times,x_adjusted,ones(1,n_y_fine));
yy = bsxfun(@times,ones(n_x_fine,1),y_adjusted');

xx_locs = bsxfun(@times,x_locs,ones(1,n_y_fine));
yy_locs = bsxfun(@times,ones(n_x_fine,1),y_locs');

% Stacked value of which rectangle the grid point is in
position = xx_locs(:) + (yy_locs(:)-1)*(n_x_knots-1);

% Evaluation of Polynomial Basis
polynomials_stacked = [ones(n_x_fine*n_y_fine,1),  xx(:),  yy(:),  ...
                        xx(:).^2,  xx(:).*yy(:),  yy(:).^2,  ...
                        xx(:).^2.*yy(:),  xx(:).*yy(:).^2];

evaluation = zeros(n_x_fine*n_y_fine,8*(n_x_knots-1)*(n_y_knots-1));

for i=(n_x_knots-1)*(n_y_knots-1):-1:1
    locs = find(position==i);
    evaluation(locs,(i-1)*8+1:8*i) = polynomials_stacked(locs,:);
end

evaluation = sparse(evaluation);

coeff_change = [1, 0, 0, 0, 0, 0, 0, 0;         % f(0,0)
                1, 1, 0, 1, 0, 0, 0, 0;         % f(1,0)
                1, 0, 1, 0, 0, 1, 0, 0;         % f(0,1)
                1, 1, 1, 1, 1, 1, 1, 1;         % f(1,1)
                0, 1, 0, 0, 0, 0, 0, 0;         % dx(0,0)
                0, 0, 1, 0, 0, 0, 0, 0;         % dy(0,0)
                0, 0, 1, 0, 1, 0, 1, 0;         % dy(1,0)
                0, 1, 0, 0, 1, 0, 0, 1];        % dx(0,1)

revert_change = sparse(inv(coeff_change));
aux = kron(speye(n_y_knots-1),spdiags(ones(n_x_knots-1,1),0,n_x_knots-1,n_x_knots));

%% Unpack Values-Derivatives to Proper Location
f00 = [aux,sparse((n_x_knots-1)*(n_y_knots-1),2*n_x_knots*n_y_knots-n_y_knots)]';
f10 = [sparse((n_x_knots-1)*(n_y_knots-1),1),aux,sparse((n_x_knots-1)*(n_y_knots-1),2*n_x_knots*n_y_knots-n_y_knots-1)]';
f01 = [sparse((n_x_knots-1)*(n_y_knots-1),n_x_knots),aux,sparse((n_x_knots-1)*(n_y_knots-1),2*n_x_knots*n_y_knots-n_x_knots-n_y_knots)]';
f11 = [sparse((n_x_knots-1)*(n_y_knots-1),n_x_knots+1),aux,sparse((n_x_knots-1)*(n_y_knots-1),2*n_x_knots*n_y_knots-n_x_knots-n_y_knots-1)]';
dx00 = [sparse((n_x_knots-1)*(n_y_knots-1),n_x_knots*n_y_knots),spdiags(repmat(h_x,n_y_knots-1,1),0,(n_x_knots-1)*(n_y_knots-1),(n_x_knots-1)*(n_y_knots-1)),sparse((n_x_knots-1)*(n_y_knots-1),n_x_knots*n_y_knots-1)]';
dx01 = [sparse((n_x_knots-1)*(n_y_knots-1),n_x_knots*n_y_knots+n_x_knots-1),spdiags(repmat(h_x,n_y_knots-1,1),0,(n_x_knots-1)*(n_y_knots-1),(n_x_knots-1)*(n_y_knots-1)),sparse((n_x_knots-1)*(n_y_knots-1),n_x_knots*n_y_knots-n_x_knots)]';
aux = kron(spdiags(h_y,0,n_y_knots-1,n_y_knots-1),spdiags(ones(n_x_knots-1,1),0,n_x_knots-1,n_x_knots));
dy00 = [sparse((n_x_knots-1)*(n_y_knots-1),2*n_x_knots*n_y_knots-n_y_knots),aux]';
dy10 = [sparse((n_x_knots-1)*(n_y_knots-1),2*n_x_knots*n_y_knots-n_y_knots+1),aux(:,1:end-1)]';

% Warning: Overwriting aux matrix here
aux = [f00(:),f10(:),f01(:),f11(:),dx00(:),dy00(:),dy10(:),dx01(:)]';

% Change Value/Derivatives to Coefficient Basis
aux = revert_change*aux;

%% A nightmare set of transformations. Good luck!
% These are supposed to reorder things to match location of the basis
% expansion
aux = reshape(aux,8*(3*n_x_knots*n_y_knots-n_x_knots-n_y_knots),(n_x_knots-1)*(n_y_knots-1));
permute = reshape(1:8*(3*n_x_knots*n_y_knots-n_x_knots-n_y_knots),8,3*n_x_knots*n_y_knots-n_x_knots-n_y_knots)';
aux = reshape(aux(permute(:),:),3*n_x_knots*n_y_knots-n_x_knots-n_y_knots,8*(n_x_knots-1)*(n_y_knots-1))';

%% Finally Basis Change
to_fine = evaluation*aux;

% Sometimes less than full basis are required up to machine accuarcy, so reduce with SVD
[u,d,~] = svd(full(to_fine),'econ');
singles = find(diag(d)>d(1,1)*eps*10);

% Basis Change Functions
to_fine = u(:,singles)*d(singles,singles);
% Since we took an SVD, so we can just invert by multiplication
from_fine = (d(singles,singles)\u(:,singles)');

author: gregkaplan

repo: phact

save_path: github-repos/MATLAB/gregkaplan-phact


text:
function r = minus(p,q)
%MINUS: this function subtracts the x, y, and z values of the structures that I use in
% the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

% First implementation:
% m = size(p.xvalue);
% d = ndims(p.xvalue);
% if (min(m) == 1) % 1D: only x value
%     r.xvalue = p.xvalue-q.xvalue;
% elseif (d == 2) % 2D: x and y values
%     r.xvalue = p.xvalue-q.xvalue;
%     r.yvalue = p.yvalue-q.yvalue;
% else % 3D:x, y, and z values
%     r.xvalue = p.xvalue-q.xvalue;
%     r.yvalue = p.yvalue-q.yvalue;
%     r.zvalue = p.zvalue-q.zvalue;
% end
r = p+(-q);

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function [sr_img] = superresolve_ebsr(slidingWindows, magFactor)
%
% LMSSR_GENERATESRUSINGEBSR Performs a Single Image Super-Resolution (SISR) approach according to "Example-based Learning for Single-Image
%                      Super-resolution". Supports GS/Y and RGB.
%    [sr_img] = LMSSR_GENERATESRUSINGEBSR(lr_img, upscaling)
%
% Parameters: lr_img          -   LR RGB or GS/Y image to be superresolved (height,width,dye)
%             upscaling       -   Scalar value containing the upscaling factor for both Y- and X-coordinates
%
% Important note: YUV images are also supported but yield bad results!
%
% Author: Michel Bätz (LMS)
%
% See also: lmsSR_framework
%

lr_img    = slidingWindows.referenceFrame;
upscaling = magFactor;

%Low = lr_ref_img;
%Low(:,:,2) = Low(:,:,1);
%Low(:,:,3) = Low(:,:,1);
%Low = im2uint8(Low);

if size(lr_img,3) == 3, % Color image case
    lr_img = lmsSR_convertYUV2RGB(lr_img);
end

sr_img = SuperresCode(lr_img,upscaling); % EXTERNAL CODE
sr_img = uint8(sr_img); % Required as output contains uint8 values

sr_img = im2double(sr_img);

if size(lr_img,3) == 3, % Color image case
    sr_img = lmsSR_convertRGB2YUV(sr_img);
end

author: thomas-koehler

repo: SupER


text:
function [grad, grad_W, grad_b] = B_tconv(prev_layer, curr_layer, future_layers, skip_grad)
input = prev_layer{1}.a;
[D,T,N] = size(input);
% Each row of W is a DxP filter to be applied along the time axis. We have
% totally H such filters.
W = curr_layer.W;
b = curr_layer.b;
[H,DP] = size(W);
P = DP/D;
halfP = (P-1)/2;

if 0
    % Instead of using conv2, we splice input and use matrix transform
    % We first put the N samples into a big sample, separated by zeros
    if strcmpi(class(input), 'gpuArray')
        input2 = gpuArray.zeros(D,T+P,N);
    else
        input2 = zeros(D,T+P,N);
    end
    input2(:,halfP+1:halfP+T,:) = input;
    X = reshape(input2, D, (T+P)*N);
    context = -halfP : halfP;
    X2 = ExpandContext_v2(X, context);
else
    X2 = curr_layer.X2;
end

X3 = reshape(X2, size(X2,1), size(X2,2)/N,N);

if isfield(curr_layer, 'range') && strcmpi(curr_layer.range, 'valid') == 0
    X3 = X3(:,halfP+1:halfP+T,:);
else
    X3 = X3(:,halfP*2+1:T,:);
end
X4 = reshape(X3, size(X3,1), size(X3,2)*N);

future_grad = GetFutureGrad(future_layers, curr_layer);
if strcmpi(class(future_grad), 'gpuArray')
    grad = gpuArray.zeros(D,T,N);
else
    grad = zeros(D,T,N);
end

fg = reshape(future_grad, H, size(future_grad,2)*N);

% Y2 = F_affine_transform(X2, W, b);
grad_W = fg * X4';
grad_b = sum(fg,2);

if skip_grad==0
    gradX4 = W' * fg;
    gradX4 = reshape(gradX4, size(gradX4,1), size(gradX4,2)/N, N);

    % now we need to distribute gradX4 to the input
    if strcmpi(class(input), 'gpuArray')
        gradX2 = gpuArray.zeros(size(X2,1), size(X2,2)/N,N);
    else
        gradX2 = zeros(size(X2,1), size(X2,2)/N,N);
    end
    if isfield(curr_layer, 'range') && strcmpi(curr_layer.range, 'valid') == 0
        gradX2(:,halfP+1:halfP+T,:) = gradX4;
    else
        gradX2(:,halfP*2+1:T,:) = gradX4;
    end
    gradX2 = reshape(gradX2, size(X2,1), size(X2,2));
    fl{1}.grad = gradX2;
    fl{1}.name = 'dummy';
    gradX = B_splice(fl, P);
    grad = reshape(gradX, D, T+P, N);
    grad = grad(:, halfP+1:halfP+T,:);
end

end

author: singaxiong

repo: SignalGraph


text:
function x = blend_rst_0dn ( r, s, t, n, bound_rst )

%*****************************************************************************80
%
%% BLEND_RST_0DN extends vector data at corners into a cube.
%
%  Diagram:
%
%    010-----r10-----110        011-----r11-----111
%      |       .       |          |       .       |
%      |       .       |          |       .       |
%    0s0.....rs0.....1s0        0s1.....rs1.....1s1     S
%      |       .       |          |       .       |     |
%      |       .       |          |       .       |     |
%    000-----r00-----100        001-----r01-----101     +----R
%           BOTTOM                      TOP
%
%    011-----0s1-----001        111-----1s1-----101
%      |       .       |          |       .       |
%      |       .       |          |       .       |
%    01t.....0st.....00t        11t.....1st.....10t          T
%      |       .       |          |       .       |          |
%      |       .       |          |       .       |          |
%    010-----0s0-----000        110-----1s0-----100     S----+
%           LEFT                       RIGHT
%
%    001-----r01-----101        011-----r11-----111
%      |       .       |          |       .       |
%      |       .       |          |       .       |
%    00t.....r0t.....100        01t.....r1t.....11t     T
%      |       .       |          |       .       |     |
%      |       .       |          |       .       |     |
%    000-----r00-----100        010-----r10-----110     +----R
%           FRONT                       BACK
%
%  Discussion:
%
%    BLEND_RST_0DN is equivalent to a trilinear finite element method.
%    Data along the edges, faces, and interior of the cube is
%    interpolated from the data at the corners.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    William Gordon,
%    Blending-Function Methods of Bivariate and Multivariate Interpolation
%    and Approximation,
%    SIAM Journal on Numerical Analysis,
%    Volume 8, Number 1, March 1971, pages 158-177.
%
%    William Gordon and Charles Hall,
%    Transfinite Element Methods: Blending-Function Interpolation over
%    Arbitrary Curved Element Domains,
%    Numerische Mathematik,
%    Volume 21, Number 1, 1973, pages 109-129.
%
%    William Gordon and Charles Hall,
%    Construction of Curvilinear Coordinate Systems and Application to
%    Mesh Generation,
%    International Journal of Numerical Methods in Engineering,
%    Volume 7, 1973, pages 461-477.
%
%    Joe Thompson, Bharat Soni, Nigel Weatherill,
%    Handbook of Grid Generation,
%    CRC Press, 1999.
%
%  Parameters:
%
%    Input, real R, S, T, the (R,S,T) coordinates of the
%    point to be evaluated.
%
%    Input, integer N, the dimension of the vector space.
%
%    External, BOUND_RST, is a function which is given (R,S,T)
%    coordinates and an component value I, and returns XI, the value
%    of the I-th component of the N-vector at that point.  BOUND_RST
%    will only be called for "corners", that is, for values (R,S,T)
%    where R, S and T are either 0.0 or 1.0.  BOUND_RST has the form:
%      function xi = bound_rst ( r, s, t, i )
%
%    Output, real X(N), the interpolated value at the
%    point (R,S,T).
%
  for i = 1 : n
%
%  Get the I-th coordinate component at the corners.
%
    x000 = bound_rst ( 0.0, 0.0, 0.0, i );
    x001 = bound_rst ( 0.0, 0.0, 1.0, i );
    x010 = bound_rst ( 0.0, 1.0, 0.0, i );
    x011 = bound_rst ( 0.0, 1.0, 1.0, i );
    x100 = bound_rst ( 1.0, 0.0, 0.0, i );
    x101 = bound_rst ( 1.0, 0.0, 1.0, i );
    x110 = bound_rst ( 1.0, 1.0, 0.0, i );
    x111 = bound_rst ( 1.0, 1.0, 1.0, i );
%
%  Interpolate the I-th coordinate component at the edges.
%
    xr00 = blend_101 ( r, x000, x100 );
    xr01 = blend_101 ( r, x001, x101 );
    xr10 = blend_101 ( r, x010, x110 );
    xr11 = blend_101 ( r, x011, x111 );

    x0s0 = blend_101 ( s, x000, x010 );
    x0s1 = blend_101 ( s, x001, x011 );
    x1s0 = blend_101 ( s, x100, x110 );
    x1s1 = blend_101 ( s, x101, x111 );

    x00t = blend_101 ( t, x000, x001 );
    x01t = blend_101 ( t, x010, x011 );
    x10t = blend_101 ( t, x100, x101 );
    x11t = blend_101 ( t, x110, x111 );
%
%  Interpolate the I-th component on the faces.
%
    x0st = blend_112 ( s, t, x000, x001, x010, x011, x0s0, x0s1, x00t, x01t );

    x1st = blend_112 ( s, t, x100, x101, x110, x111, x1s0, x1s1, x10t, x11t );

    xr0t = blend_112 ( r, t, x000, x001, x100, x101, xr00, xr01, x00t, x10t );

    xr1t = blend_112 ( r, t, x010, x011, x110, x111, xr10, xr11, x01t, x11t );

    xrs0 = blend_112 ( r, s, x000, x010, x100, x110, xr00, xr10, x0s0, x1s0 );

    xrs1 = blend_112 ( r, s, x001, x011, x101, x111, xr01, xr11, x0s1, x1s1 );
%
%  Interpolate the I-th coordinate component of the interior point.
%
    x(i) = blend_123 ( r, s, t, x000, x001, x010, x011, x100, x101, x110, x111, ...
      xr00, xr01, xr10, xr11, x0s0, x0s1, x1s0, x1s1, x00t, x01t, x10t, x11t, ...
      x0st, x1st, xr0t, xr1t, xrs0, xrs1 );

  end

  return
end


text:
function img = dtiSmooth3(img, sigma)
% img = dtiSmooth3(img, sigma)
%
% Gaussian smoothing (faster than matlab's smooth3)
%
% NOTE: We should add some zero-padding to avoid wrap-around issues.

persistent gauss;
persistent sigmaCache;
persistent dimsCache;

if(numel(sigma==1))
    sigma = [sigma sigma sigma];
end
if(all(sigma<=0))
    return;
end
dims = size(img);
if ~isa(img, 'double'), type = class(img); img = double(img); end
% there should be a more elegant way of preserving the image intensity
% range!
scale = max(abs(img(:)));
if(scale==0), return; end

if(isempty(gauss) || numel(sigma)~=numel(sigmaCache) || any(sigma~=sigmaCache) || numel(dims)~=numel(dimsCache) || any(dims~=dimsCache))
    gauss = newGauss(dims, sigma);
    sigmaCache = sigma;
    dimsCache = dims;
end

% Allow for 4d timeseries or 5d tensor data.
for(ii=1:size(img,4))
    for(jj=1:size(img,5))
        ft = fftshift(fftn(img(:,:,:,ii,jj)));
        ft = gauss.*ft;
        img(:,:,:,ii,jj) = real(ifftn(ifftshift(ft)));
    end
end
img = img./max(abs(img(:))).*scale;

if exist('type', 'var')
	% we had a non-double matrix, convert back
	img = feval(type, img);
end
return;

function gauss = newGauss(dimSrc,variance)
gaussx = myGausswin(dimSrc(1),variance(1));
gaussy = myGausswin(dimSrc(2),variance(2));
gaussz = myGausswin(dimSrc(3),variance(3));
gaussxy = gaussx*gaussy';
gauss = repmat(gaussxy,[1 1 dimSrc(3)]);
for i = 1:dimSrc(1)
    for j = 1:dimSrc(2)
        gauss(i,j,:) = gaussxy(i,j)*gaussz;
    end
end
return

function w = myGausswin(n, alpha)
% alpha is 1/stdev (width of the gaussian in Fourier space).
k = [-(n-1)/2:(n-1)/2];
w = exp((-1/2)*(alpha * k/(n/2)).^2)';
% Make sure the peak is at 1 so that we'll preserve to DC level
w(w==max(w)) = 1;
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% S = SKEW2(MTX,MEAN,VAR)
%
% Sample skew (third moment divided by variance^3/2) of a matrix.
%  MEAN (optional) and VAR (optional) make the computation faster.

function res = skew2(mtx, mn, v)

if (exist('mn') ~= 1)
  mn =  mean2(mtx);
end

if (exist('v') ~= 1)
  v =  var2(mtx,mn);
end

if (isreal(mtx))
  res = mean(mean((mtx-mn).^3)) / (v^(3/2));
else
  res = mean(mean(real(mtx-mn).^3)) / (real(v)^(3/2)) + ...
      i * mean(mean(imag(mtx-mn).^3)) / (imag(v)^(3/2));
end

repo: SelfExSR