text:
function test_ft_timelockanalysis_new(datainfo,writeflag)

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_timelockanalysis_new ft_timelockanalysis ref_datasets

% the optional writeflag determines whether the output should be saved to
% disk

%%
% this is a function for testing ft_timelockanalysis_new, which is not official yet
%
% This function is testing a new ft_timelockanalysis_new which is something
% Johanna is working on and not in SVN yet.

return

%%

if nargin<2
  writeflag = 0;
end
if nargin<1
  datainfo = ref_datasets;
end

% for k = 1:numel(datainfo)
for k = 1:10
  datanew = timelockanalysis10trials(datainfo(k), writeflag);

  fname = fullfile(datainfo(k).origdir,'latest/timelock',datainfo(k).type,'timelock_',datainfo(k).datatype);
  tmp = load(fname);
  if isfield(tmp, 'data')
    data = tmp.data;
  elseif isfield(tmp, 'datanew')
    data = tmp.datanew;
  else isfield(tmp, 'timelock')
    data = tmp.timelock;
  end

  datanew = rmfield(datanew, 'cfg'); % these are per construction different if writeflag = 0;
  data    = rmfield(data,    'cfg');
  assert(isequaln(data, datanew));
end

test_cfg_options;

function [tlck] = timelockanalysis10trials(dataset, writeflag)

cfg        = [];
cfg.inputfile  = fullfile(dataset.origdir,'latest/raw',dataset.type,['preproc_',dataset.datatype]);
if writeflag
  cfg.outputfile = fullfile(dataset.origdir,'latest/timelock',dataset.type,'timelock_',dataset.datatype);
end
tlck = ft_timelockanalysis(cfg);
tlck1 = ft_timelockanalysis_new(cfg);
return

function test_cfg_options
load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/eventrelatedaveraging/dataFC_LP.mat'));
data = dataFC_LP;
clear dataFC_LP

data.time{2}=data.time{2}+.5; % purposely add some jitter to time window
data.time{3}=data.time{3}-.5;

cfg=[];
try
    tlock=ft_timelockanalysis_new(cfg,data);
catch me
%     if ~strcmp(me.message,'the option "output" was not specified or was empty');
        error(me.message)
%     end
end
cfg=[];
cfg.output='rubbish';
try
    tlock=ft_timelockanalysis_new(cfg,data);
catch me
    if ~strcmp(me.message,'the value of cfg.output is not set correctly');
        error(me.message)
    end
end

% no latency or covlatency given, use defaults
cfg=[];
cfg.output='avg';
cfg.feedback='none';
cfg.preproc.feedback='textbar';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
% cfg.vartrllength=1;
% tlocko=ft_timelockanalysis(cfg,data);
cfg=[];
cfg.output='cov';
cfg.feedback='none';
cfg.preproc.feedback='textbar';
cfg.keeptrials='no';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end
cfg=[];
cfg.output='avgandcov';
cfg.feedback='none';
cfg.preproc.feedback='textbar';
cfg.keeptrials='no';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end
% cfg.covariance='yes';
% cfg.vartrllength=1;
% tlocko=ft_timelockanalysis(cfg,data);
cfg=[];
cfg.output='avg';
cfg.feedback='none';
cfg.preproc.feedback='textbar';
cfg.keeptrials='yes';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end
% cfg.vartrllength=1;
% tlocko=ft_timelockanalysis(cfg,data);
cfg=[];
cfg.output='cov';
cfg.feedback='none';
cfg.preproc.feedback='textbar';
cfg.keeptrials='yes';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end
cfg=[];
cfg.output='avgandcov';
cfg.feedback='none';
cfg.preproc.feedback='textbar';
cfg.keeptrials='yes';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end
% cfg.covariance='yes';
% cfg.vartrllength=1;
% tlocko=ft_timelockanalysis(cfg,data);

% options with .latency and .covlatency specified
cfg=[];
cfg.output='avg';
cfg.feedback='none';
cfg.preproc.feedback='textbar';
cfg.keeptrials='no';
cfg.latency=[min(data.time{1}) max(data.time{1})];
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end
% cfg.vartrllength=1;
% tlocko=ft_timelockanalysis(cfg,data);
cfg.latency='maxperlength';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end
cfg.latency='minperlength';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end
cfg.latency='prestim';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end
cfg.latency='poststim';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end

cfg=[];
cfg.feedback='none';
cfg.preproc.feedback='textbar';
cfg.output='avg';
cfg.keeptrials='yes';
cfg.latency=[min(data.time{1}) max(data.time{1})];
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end
cfg.latency='maxperlength';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end
cfg.latency='minperlength';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end
cfg.latency='prestim';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end
cfg.latency='poststim';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end

cfg=[];
cfg.feedback='none';
cfg.preproc.feedback='textbar';
cfg.output='avgandcov';
cfg.keeptrials='no';
cfg.latency=[min(data.time{1}) max(data.time{1})];
cfg.covlatency=[min(data.time{1}) max(data.time{1})];
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end
% cfg.covariance='yes';
% cfg.vartrllength=1;
% tlocko=ft_timelockanalysis(cfg,data);
cfg.latency='maxperlength';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end

cfg=[];
cfg.feedback='none';
cfg.preproc.feedback='textbar';
cfg.output='avgandcov';
cfg.keeptrials='yes';
cfg.latency=[min(data.time{1}) max(data.time{1})];
cfg.covlatency=[min(data.time{1}) max(data.time{1})];
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end
% cfg.covariance='yes';
% cfg.vartrllength=1;
% tlocko=ft_timelockanalysis(cfg,data);

cfg=[];
cfg.feedback='none';
cfg.preproc.feedback='textbar';
cfg.output='cov';
cfg.keeptrials='no';
cfg.latency=[min(data.time{1}) max(data.time{1})];
cfg.covlatency=[min(data.time{1}) max(data.time{1})];
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end

cfg=[];
cfg.feedback='none';
cfg.preproc.feedback='textbar';
cfg.output='cov';
cfg.keeptrials='yes';
cfg.latency=[min(data.time{1}) max(data.time{1})];
cfg.covlatency=[min(data.time{1}) max(data.time{1})];
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end

% check toi options
cfg=[];
cfg.feedback='none';
cfg.preproc.feedback='textbar';
cfg.output='cov';
cfg.latency=[min(data.time{1}) max(data.time{1})];
cfg.covlatency='minperlength';
cfg.toi=[-0.5 0.7];
cfg.timwin=1;
cfg.equatenumtrials='yes';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end

cfg=[];
cfg.feedback='none';
cfg.preproc.feedback='textbar';
cfg.output='cov';
cfg.covlatency='minperlength';
cfg.toi=[-.8:.3:0.1];
cfg.timwin=1;
cfg.equatenumtrials='no';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end

cfg=[];
cfg.feedback='none';
cfg.preproc.feedback='textbar';
cfg.output='cov';
cfg.covlatency='minperlength';
cfg.toi=[-.8:.3:0.1];
cfg.timwin=1;
cfg.equatenumtrials='no';
cfg.keeptrials='yes';
try
    tlock=ft_timelockanalysis_new(cfg,data);
    if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end
catch me
    if ~strcmp(me.message,'sorry, if keeping trials and computing cov, cfg.equatenumtrials should be yes')
        error(me.message)
    end
end

cfg=[];
cfg.feedback='none';
cfg.preproc.feedback='textbar';
cfg.output='cov';
cfg.covlatency='minperlength';
cfg.toi=[-.8:.3:0.1];
cfg.timwin=1;
cfg.equatenumtrials='yes';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end

cfg=[];
cfg.feedback='none';
cfg.preproc.feedback='textbar';
cfg.output='cov';
cfg.covlatency='minperlength';
cfg.toi=[-.8:.3:0.1];
cfg.timwin=1;
cfg.equatenumtrials='yes';
cfg.keeptrials='yes';
tlock=ft_timelockanalysis_new(cfg,data);
if ~strmatch(ft_datatype(tlock),'timelock'),error('datatype');end;
if ~strcmp(ft_datatype(ft_checkdata(tlock,'datatype','raw')),'raw') || ~strcmp(ft_datatype(ft_checkdata(tlock)),'timelock'),error('checkdata');end

%% ignore for svn testing, but used to test output of new function in further functions

% test ft_sourceanalysis versus MNE event related tutorial
if 0

    cfg=[];
    cfg.output='avg';
    cfg.feedback='none';
    cfg.preproc.feedback='textbar';
    tlock=ft_timelockanalysis_new(cfg,data);
    cfg=[];
    cfg.method='lcmv';
    cfg.hdmfile=dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.hdm');
    cfg.grad=data.grad;
    source=ft_sourceanalysis(cfg,tlock);

    cfg = [];
    cfg.covariance = 'yes';
    cfg.vartrllength=1;
    cfg.covariancewindow = [-inf 0]; %it will calculate the covariance matrix
    % on the timepoints that are
    % before the zero-time point in the trials
    tlckFC = ft_timelockanalysis(cfg, data);
%     tlckFIC = ft_timelockanalysis(cfg, dataFIC_LP);
    save tlck tlckFC tlckFIC;
    cfg=[];
    cfg.method='lcmv';
    cfg.hdmfile=dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.hdm');
    cfg.grad=data.grad;
    sourceFC=ft_sourceanalysis(cfg,tlckFC);
end

% spin off of beamformer tutorial but in time-domain, results won't match exactly
if 0
    cfg=[];
    cfg.output='cov';
    cfg.feedback='none';
    cfg.preproc.feedback='textbar';
    cfg.covlatency=[0.8 1.3];
    cfg.preproc.bpfilter='yes';
    cfg.preproc.bpfreq=[16 20];
    tlock=ft_timelockanalysis_new(cfg,data);
    load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/beamformer/segmentedmri.mat'))
    cfg=[];
    vol=ft_prepare_singleshell(cfg,segmentedmri);
    cfg=[];
    cfg.headmodel=vol;
    cfg.reducerank = 2;
    cfg.grad=tlock.grad;
    cfg.sourcemodel.resolution=1;
    cfg.channel={'MEG','-MLP31','-MLO12'};
    grid=ft_prepare_leadfield(cfg);
    cfg=[];
    cfg.method='lcmv';
    cfg.projectnoise='yes';
    cfg.sourcemodel=grid;
    cfg.headmodel=vol;
    source=ft_sourceanalysis(cfg,tlock);

    mri=ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.mri'));
    sourcediff=source;
    % sourcediff.avg.pow=(source.avg.pow-source.avg.noise)./source.avg.noise;
    sourcediff.avg.pow=(source.avg.pow)./source.avg.noise;
    cfg=[];
    cfg.downsample=2;
    sourcediffint=ft_sourceinterpolate(cfg,sourcediff,mri);
    cfg=[];
    cfg.method='slice';
    cfg.funparameter='avg.pow';
    cfg.maskparameter=cfg.funparameter;
    cfg.funcolorlim=[5 6.2];
    cfg.opacitylim=[5 6.2];
    cfg.opacitymap='rampup';
    figure;ft_sourceplot(cfg,sourcediffint); %ok

    cfg=[];
    cfg.downsample=2;
    sourceint=ft_sourceinterpolate(cfg,source,mri);
    cfg=[];
    cfg.method='slice';
    cfg.funparameter='avg.pow';
    % cfg.maskparameter=cfg.funparameter;
    % cfg.funcolorlim=[5 6.2];
    % cfg.opacitylim=[5 6.2];
    % cfg.opacitymap='rampup';
    figure;ft_sourceplot(cfg,sourceint); %ok

end

% test ft_timelockstatistics
if 0

cfg=[];
cfg.output='cov';
cfg.feedback='none';
cfg.keeptrials='yes';
cfg.preproc.feedback='textbar';
tlock=ft_timelockanalysis_new(cfg,data);

tlock1 = ft_selectdata(tlock, 'rpt',1:36)
tlock2 = ft_selectdata(tlock, 'rpt',37:72)

cfg=[];
cfg.method='analytic';
   cfg.statistic='ft_statfun_indepsamplesT' ;
stat=ft_timelockstatistics(cfg,tlock1,tlock2);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function SPGP
clc;
clear;
close all;
%% the example fuction
datasize=200;
var_noise=0.1;
x=linspace(-10,10,datasize)';
    f=sin(x)./x;
noise=var_noise*randn(datasize,1);
y=f+noise;
x_raw=x;
plot(x_raw,f,'r-','LineWidth',2),hold on;
%%
Inset=30;
dataindex=randsample(datasize,Inset);
xI=x(dataindex);
yI=y(dataindex);
%%
%%test
testnum=500;
x_test=linspace(-15,15,testnum)';
K_nn=Kernel_func(x,x,0.7,1.2);
K_mn=Kernel_func(xI,x,0.7,1.2);
K_mm=Kernel_func(xI,xI,0.7,1.2);
lambda=diag(diag(K_nn)-diag(K_mn'*inv(K_mm)*K_mn));
K_mat=K_mn*inv(var_noise*eye(datasize)+lambda)*K_mn'+K_mm;
K_sm=Kernel_func(x_test,xI,0.7,1.2);
L=chol(K_mat,'lower');
b=K_mn*inv(var_noise*eye(datasize)+lambda)*y;
f_mean=K_sm*(L'\(L\b));

for i=1:testnum
K_sm=Kernel_func(x_test(i),xI,0.7,1.2);
K_ss=Kernel_func(x_test(i),x_test(i),0.7,1.2);
v=L\(K_sm');
  f_var(i)=K_ss-K_sm*inv(K_mm)*K_sm'+v'*v;
end
f_var=f_var';

plot(x_test,f_mean,'b-','LineWidth',2),hold on;

x_area=[x_test;flipdim(x_test,1)];
ybound_std=[(f_mean+1.96*sqrt(f_var));flipdim(f_mean-1.96*sqrt(f_var),1)];
fill(x_area,ybound_std,[187 255 0]/255,'EdgeColor','none');
plot(x_raw,f,'r-','LineWidth',2),hold on;
plot(xI,yI,'b+'),hold on;
plot(x_test,f_mean,'b-','LineWidth',2),hold on;
%%
legend1=legend('target function','training points','preditive distribution','95% confidence interval');
 set(legend1,'Box','off','Color','none',...
     'Location','NorthEast');
xlabel('(g)SPGP with 30 points ')
xlim([-15,15]);
ylim([-2,2]);
set(gca,'XTick',[-15:3:15])
matlab2tikz( 'SPGP30.tex' )

author: epfl-lasa


text:
%  Inspired by a question on the dynare forum at
%  http://www.dynare.org/phpBB3/viewtopic.php?f=1&t=3206
%

author: jmaih


text:
%IMM_UPDATE  UKF based Interacting Multiple Model (IMM) Filter update step
%
% Syntax:
%   [X_i,P_i,MU,X,P] = IMM_UPDATE(X_p,P_p,c_j,ind,dims,Y,H,R)
%
% In:
%   X_p  - Cell array containing N^j x 1 mean state estimate vector for
%          each model j after prediction step
%   P_p  - Cell array containing N^j x N^j state covariance matrix for
%          each model j after prediction step
%   c_j  - Normalizing factors for mixing probabilities
%   ind  - Indices of state components for each model as a cell array
%   dims - Total number of different state components in the combined system
%   Y    - Dx1 measurement vector.
%   H    - Measurement matrices for each model as a cell array.
%   h    - Measurement mean
%   param - parameters
%   R    - Measurement noise covariances for each model as a cell array.
%
% Out:
%   X_i  - Updated state mean estimate for each model as a cell array
%   P_i  - Updated state covariance estimate for each model as a cell array
%   MU   - Probabilities of each model
%   X    - Combined state mean estimate
%   P    - Combined state covariance estimate
%
% Description:
%   IMM-UKF filter measurement update step. If some of the models have linear
%   measurements standard Kalman filter update step is used for those.
%
% See also:
%   IMM_PREDICT, IMM_SMOOTH, IMM_FILTER

% History:
%   01.11.2007 JH The first official version.
%
% Copyright (C) 2007 Jouni Hartikainen
%
% $Id: imm_update.m 111 2007-11-01 12:09:23Z jmjharti $
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.

function [X_i,P_i,MU,X,P] = uimm_update(X_p,P_p,c_j,ind,dims,Y,H,h,R,param)
    % Number of models
    m = length(X_p);

    % Space for update state mean, covariance and likelihood of measurements
    X_i = cell(1,m);
    P_i = cell(1,m);
    lambda = zeros(1,m);

    % Update for each model
    for i = 1:m
        % Update the state estimates
        if isempty(h) | isempty(h{i})
            [X_i{i}, P_i{i}, K, IM, IS, lambda(i)] = kf_update(X_p{i},P_p{i},Y,H{i},R{i});
        else
            [X_i{i}, P_i{i}, K, IM, IS, lambda(i)] = ukf_update1(X_p{i},P_p{i},Y,h{i},R{i},param{i});
        end
    end

    % Calculate the model probabilities
    MU = zeros(1,m);
    c = sum(lambda.*c_j);
    MU = c_j.*lambda/c;

    % In case lambda's happen to be zero
    if c == 0
        MU = c_j;
    end

    % Output the combined updated state mean and covariance, if wanted.
    if nargout > 3
        % Space for estimates
        X = zeros(dims,1);
        P = zeros(dims,dims);
        % Updated state mean
        for i = 1:m
            X(ind{i}) = X(ind{i}) + MU(i)*X_i{i};
        end
        % Updated state covariance
        for i = 1:m
            P(ind{i},ind{i}) = P(ind{i},ind{i}) + MU(i)*(P_i{i} + (X_i{i}-X(ind{i}))*(X_i{i}-X(ind{i}))');
        end
    end

author: EEA-sensors

repo: ekfukf

save_path: github-repos/MATLAB/EEA-sensors-ekfukf


text:
% Example for the adjacency matrix of
%      0     0     1     1     1    -1     1     0
%      0     0     1    -1     0     1     1    -1
%      1     1     0     0     0     1     0     0
%      1    -1     0     0     0     0    -1     0
%      1     0     0     0     0     1    -1     1
%     -1     1     1     0     1     0     1    -1
%      1     1     0    -1    -1     1     0     0
%      0    -1     0     0     1    -1     0     0
% and 100 iterations of simulations

load adj2
heider(adj2,100);

author: Sable

repo: mcbench-benchmarks


text:
2     % problem
5     % outflow boundary
5     % grid parameter
1     % discretisation
0.02  % viscosity parameter
1     % Picard/Newton/hybrid linearization
9     % number of Picard iterations
1.d-5 % nonlinear tolerance

%% Data file for backward facing step

author: johannesgerer

repo: jburkardt-m


text:
function [spos_n, svel_b, sCsn, simu]=la_imu(acc, gyro, gyro_der, pos_n, vel_b, Cbn, imu_la, imu_or, typ, gyro_vib, imu_vib)

nimu=size(imu_la,2);
spos_n=zeros(3,nimu);
svel_b=zeros(3,nimu);
sCsn=zeros(3,3,nimu);
simu=zeros(6,nimu);

%vibration information is not compulsory
if (isempty(gyro_vib))
    gyro_vib=zeros(3,nimu);
end

if (isempty(imu_vib))
    imu_vib=zeros(3,3,nimu);
    for in=1:nimu
        imu_vib(:,:,in)=eye(3);
    end
end

if (typ==0)
    [Rn, Re, g, sL, cL, WIE]=geo_param(pos_n);
end
for in=1:nimu
    %Attitude
    Crig=imu_or(:,:,in);
    Cvib=imu_vib(:,:,in);
    Csm=Crig*Cvib;
    Csn=Cbn*Csm;
    sCsn(:,:,in)=Csn;

    %position
    la_m=imu_la(:,in);
    if (typ==0)
        spos_n(:,in)=pos_n+[1/(Rn+pos_n(3)) 0 0;0 1/cL/(Re+pos_n(3)) 0; 0 0 -1]*Cbn*la_m;
    elseif (typ==1)
        spos_n(:,in)=pos_n+Cbn*la_m;
    end

    %Velocity
    if (typ==0)
        wie_m=Cbn'*[WIE*cL; 0; -WIE*sL];
    else
        wie_m=zeros(3,1);
    end
    svel_b(:,in)=Csm'*(vel_b+cross(gyro-wie_m,la_m));

    %Acceleration
    simu(1:3,in)=Csm'*(acc+cross(gyro,cross(gyro,la_m))+cross(gyro_der,la_m));

    %Rotation rate
    simu(4:6,in)=Csm'*(gyro+gyro_vib(:,in));
end

repo: ekfmonoslam


text:
clear;
beta1=[2 0.9];
beta0=[3 2   ];
% fitfunc=@(x,betav)(betav(1)*x.^2+betav(2)*x+betav(3));
fitfunc=@(x,betav) betav(1)*x.^3.*1./(exp(x./betav(2))-1);
% fitfunc=@(x,betav) (sin(betav(1).*x)+0.1*cos(betav(2).*x)+0.15*sin(betav(3).*x));
% fitfunc=@(x,betav) betav(1).*x+betav(2)+betav(3)*x.^2+betav(4).*x.^3;
% fitfunc=@(x,betav) betav(3)*sqrt(2*pi)*abs(betav(2)) *normpdf(x,betav(1),abs(betav(2)));
% fitfunc=@(x,betav) betav(1).*x+betav(2)+betav(3)*x.^2+betav(4).*x.^3+...
% betav(5).*x.^4+betav(6).*x.^5+betav(7).*x.^6+betav(8).*x.^7;

error=0.01 ;
x=[10*eps:0.5:10];
y=fitfunc(x,beta1);
y=y+(randn(size(y)))*error;
yerr=error.*ones(size(y));

errorbarm([x; y; yerr; ])
% pause
textposition=[];
plotbool=true;
axisin=[-10 15 0 12];
axisin=[];
headercell={'Nonlinear Fit'  '$f(x)=\frac{A_0}{\sqrt{2\pi\sigma^2}}\exp{\frac{-(x-\mu)^2}{2\sigma}}$' ''};
mylabel={'xaxis','yaxis [cm]','$\mu$','$\sigma$' '$A_0$'};

[beta betaerr chi prob chiminvec]=wnonlinfit(x,y,yerr,fitfunc,beta0...
    ,'chitol',5,'label',mylabel,'position',textposition,...
    'header',headercell,...
    'printchi','off','axis',axisin,'errprec',2);

 print -dpdf -cmyk test.pdf

% figure(2)
% print -dpdf -cmyk testres.pdf


text:
function [ output_meta ] = meta2sicd_rs_xml( xml_domnode, beta_domnode, noise_domnode )
%META2SICD_RS_XML Converts Radarsat product.xml description into a SICD-style metadata structure
%
% Takes as input a Document Object Model (DOM) node from the RS2
% product.xml descriptor file.
%
% This function does NOT handle ScanSAR datasets.
%
% There is an outstanding question with regard to the meaning of the
% pulseRepetitionFrequency as provided.  See comments in Timeline section
% below.
%
% Written by: Wade Schwartzkopf, NGA/Research
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Setup
SECONDS_IN_A_DAY = 24*60*60;
xp=javax.xml.xpath.XPathFactory.newInstance.newXPath();

%% CollectionInfo
output_meta.CollectionInfo.CollectorName=char(xp.evaluate(...
    xpath_str({'product','sourceAttributes','satellite'}),xml_domnode));
if strcmpi(output_meta.CollectionInfo.CollectorName, 'RADARSAT-2')
    gen = 'RS2';
elseif strncmpi(output_meta.CollectionInfo.CollectorName, 'RCM', 3)
    gen = 'RCM';
end
[rawDataStartTime, rawDataStartTimeFrac] = datenum_w_frac(char(xp.evaluate(...
    xpath_str({'product','sourceAttributes','rawDataStartTime'}),...
    xml_domnode)));
if strcmp(gen,'RS2')
    output_meta.CollectionInfo.CoreName = [... % Start with NGA-like prefix
        upper(datestr(rawDataStartTime,'ddmmmyy')) 'RS02' ...
        char(xp.evaluate(xpath_str({'product','sourceAttributes','imageId'}),xml_domnode))];
elseif strcmp(gen, 'RCM')
    output_meta.CollectionInfo.CoreName = [... % Start with NGA-like prefix
        upper(datestr(rawDataStartTime,'ddmmmyy')) ...
        'RCM' output_meta.CollectionInfo.CollectorName(end) ...
        datestr(rawDataStartTime,'HHMMSS')]; % Make time of day unique identier within day
    % ScanSAR might need multiple CoreNames, one for each burst?
end
output_meta.CollectionInfo.CollectType='MONOSTATIC';
output_meta.CollectionInfo.RadarMode.ModeID=char(xp.evaluate(...
    xpath_str({'product','sourceAttributes','beamModeMnemonic'}),xml_domnode));
beamMode = char(xp.evaluate(xpath_str({...
    'product','sourceAttributes','beamMode'}),xml_domnode));
acqType = char(xp.evaluate(xpath_str({...
    'product','sourceAttributes','radarParameters','acquisitionType'}),xml_domnode));
if ((~isempty(beamMode) && strncmpi(beamMode, 'SPOTLIGHT', 9)) || ...
   (~isempty(acqType) && strncmpi(acqType, 'SPOTLIGHT', 9)) || ...
   ~isempty(strfind(output_meta.CollectionInfo.RadarMode.ModeID,'SL')))
    output_meta.CollectionInfo.RadarMode.ModeType = 'SPOTLIGHT';
elseif strcmpi(output_meta.CollectionInfo.RadarMode.ModeID(1:2), 'SC')
    error('META2SICD_RS_XML:RS_SCANSAR', 'ScanSAR mode data is not currently handled.');
else % Finally assume it's stripmap
    output_meta.CollectionInfo.RadarMode.ModeType = 'STRIPMAP';
end
if strcmp(gen, 'RS2') % All RS2 data is unclassified
    output_meta.CollectionInfo.Classification='UNCLASSIFIED';
elseif strcmp(gen, 'RCM') % RCM has this as a specific field
    output_meta.CollectionInfo.Classification=upper(char(xp.evaluate(...
        xpath_str({'product', 'securityAttributes', 'securityClassification'}), xml_domnode)));
end

%% ImageCreation
output_meta.ImageCreation.Application=char(xp.evaluate(...
    xpath_str({'product','imageGenerationParameters','generalProcessingInformation','softwareVersion'}),...
    xml_domnode));
output_meta.ImageCreation.DateTime=datenum(char(xp.evaluate(...
    xpath_str({'product','imageGenerationParameters','generalProcessingInformation','processingTime'}),...
    xml_domnode)),'yyyy-mm-ddTHH:MM:SS.FFF');
output_meta.ImageCreation.Site=char(xp.evaluate(...
    xpath_str({'product','imageGenerationParameters','generalProcessingInformation','processingFacility'}),...
    xml_domnode));
output_meta.ImageCreation.Profile='Prototype';

%% ImageData
if strcmp(gen, 'RS2')
    output_meta.ImageData.NumRows=uint32(str2double(xp.evaluate(...
       xpath_str({'product','imageAttributes','rasterAttributes','numberOfSamplesPerLine'}),...
       xml_domnode)));
    output_meta.ImageData.NumCols=uint32(str2double(xp.evaluate(...
       xpath_str({'product','imageAttributes','rasterAttributes','numberOfLines'}),...
        xml_domnode)));
elseif strcmp(gen, 'RCM')
    output_meta.ImageData.NumRows=uint32(str2double(xp.evaluate(...
       xpath_str({'product','sceneAttributes','imageAttributes','samplesPerLine'}),...
       xml_domnode)));
    output_meta.ImageData.NumCols=uint32(str2double(xp.evaluate(...
       xpath_str({'product','sceneAttributes','imageAttributes','numLines'}),...
        xml_domnode)));
end
output_meta.ImageData.FullImage=output_meta.ImageData;
output_meta.ImageData.FirstRow=uint32(0); output_meta.ImageData.FirstCol=uint32(0);
output_meta.ImageData.PixelType='RE16I_IM16I';  % RS2 always 16-bit
if strcmp(gen, 'RCM') && str2double(xp.evaluate(...  % RCM can be 16 or 32
       xpath_str({'product','imageReferenceAttributes','rasterAttributes','bitsPerSample'}),...
        xml_domnode)) == 32
    output_meta.ImageData.PixelType='RE32F_IM32F';
end
% Seems that all pixels are always valid
output_meta.ImageData.ValidData.Vertex(1).Row = uint32(0);
output_meta.ImageData.ValidData.Vertex(1).Col = uint32(0);
output_meta.ImageData.ValidData.Vertex(2).Row = uint32(0);
output_meta.ImageData.ValidData.Vertex(2).Col = output_meta.ImageData.NumCols-1;
output_meta.ImageData.ValidData.Vertex(3).Row = output_meta.ImageData.NumRows-1;
output_meta.ImageData.ValidData.Vertex(3).Col = output_meta.ImageData.NumCols-1;
output_meta.ImageData.ValidData.Vertex(4).Row = output_meta.ImageData.NumRows-1;
output_meta.ImageData.ValidData.Vertex(4).Col = uint32(0);

%% SCP
if strcmp(gen, 'RS2')
    im_at_str = 'imageAttributes';
elseif strcmp(gen, 'RCM')
    im_at_str = 'imageReferenceAttributes';
end
% There are many different equally valid options for picking the SCP point.
% One way is to chose the tie point that is closest to the image center.
num_tie_points=str2double(xp.evaluate(...
    ['count(' xpath_str({'product',im_at_str,'geographicInformation','geolocationGrid','imageTiePoint'}) ')'],...
    xml_domnode));
tiePointPixels = zeros(2,num_tie_points);
tiePointGeo    = zeros(3,num_tie_points);
for i=1:num_tie_points
    tiePointPixels(1,i) = str2double(xp.evaluate(...
            [xpath_str({'product',im_at_str,'geographicInformation','geolocationGrid','imageTiePoint'})...
            '[' num2str(i) ']' xpath_str({'imageCoordinate','pixel'})],...
            xml_domnode));
    tiePointPixels(2,i) = str2double(xp.evaluate(...
            [xpath_str({'product',im_at_str,'geographicInformation','geolocationGrid','imageTiePoint'})...
            '[' num2str(i) ']' xpath_str({'imageCoordinate','line'})],...
            xml_domnode));
    tiePointGeo(1,i) = str2double(xp.evaluate(...
            [xpath_str({'product',im_at_str,'geographicInformation','geolocationGrid','imageTiePoint'})...
            '[' num2str(i) ']' xpath_str({'geodeticCoordinate','latitude'})],...
            xml_domnode));
    tiePointGeo(2,i) = str2double(xp.evaluate(...
            [xpath_str({'product',im_at_str,'geographicInformation','geolocationGrid','imageTiePoint'})...
            '[' num2str(i) ']' xpath_str({'geodeticCoordinate','longitude'})],...
            xml_domnode));
    tiePointGeo(3,i) = str2double(xp.evaluate(...
            [xpath_str({'product',im_at_str,'geographicInformation','geolocationGrid','imageTiePoint'})...
            '[' num2str(i) ']' xpath_str({'geodeticCoordinate','height'})],...
            xml_domnode));
end
% Pick tie point closest to center for SCP
centerPoint = [double(output_meta.ImageData.NumRows-1)/2.0; ...
               double(output_meta.ImageData.NumCols-1)/2.0];
D = (tiePointPixels - repmat(centerPoint,1,num_tie_points));
[C,scp_index] = min( sqrt( sum(D.^2) ) );
output_meta.ImageData.SCPPixel.Row = uint32(tiePointPixels(1,scp_index));
output_meta.ImageData.SCPPixel.Col = uint32(tiePointPixels(2,scp_index));

% Sometimes lines up with SCP in Lockheed SICDs (and sometimes not):
% output_meta.ImageData.SCPPixel.Col = ceil(output_meta.ImageData.NumCols/2)-1;
% output_meta.ImageData.SCPPixel.Row = floor(output_meta.ImageData.NumRows/2);

%% GeoData
% All RS2 and RCM data we know use the WGS84 model, although it is stated
% in slightly different XML fields in the RS2 and RCM product.xml
output_meta.GeoData.EarthModel='WGS_84';
% Initially, we just seed this with a rough value.  Later we will put in
% something more precise.
output_meta.GeoData.SCP.LLH.Lat = tiePointGeo(1,scp_index);
output_meta.GeoData.SCP.LLH.Lon = tiePointGeo(2,scp_index);
output_meta.GeoData.SCP.LLH.HAE = tiePointGeo(3,scp_index);
pos_ecf = geodetic_to_ecf(tiePointGeo(:,scp_index));
output_meta.GeoData.SCP.ECF.X = pos_ecf(1);
output_meta.GeoData.SCP.ECF.Y = pos_ecf(2);
output_meta.GeoData.SCP.ECF.Z = pos_ecf(3);
% Corner coordinates will be computed later by derived_sicd_fields.
% Corners
% min_row=min(tiePointPixels(1,:));
% min_col=min(tiePointPixels(2,:));
% max_row=max(tiePointPixels(1,:));
% max_col=max(tiePointPixels(2,:));
% ll_index=find((tiePointPixels(1,:)==min_row)&(tiePointPixels(2,:)==min_col), 1);
% ul_index=find((tiePointPixels(1,:)==max_row)&(tiePointPixels(2,:)==min_col), 1);
% ur_index=find((tiePointPixels(1,:)==max_row)&(tiePointPixels(2,:)==max_col), 1);
% lr_index=find((tiePointPixels(1,:)==min_row)&(tiePointPixels(2,:)==max_col), 1);
% output_meta.GeoData.ImageCorners.ICP.FRFC.Lat=tiePointGeo(1,ul_index);
% output_meta.GeoData.ImageCorners.ICP.FRFC.Lon=tiePointGeo(2,ul_index);
% output_meta.GeoData.ImageCorners.ICP.FRLC.Lat=tiePointGeo(1,ur_index);
% output_meta.GeoData.ImageCorners.ICP.FRLC.Lon=tiePointGeo(2,ur_index);
% output_meta.GeoData.ImageCorners.ICP.LRLC.Lat=tiePointGeo(1,lr_index);
% output_meta.GeoData.ImageCorners.ICP.LRLC.Lon=tiePointGeo(2,lr_index);
% output_meta.GeoData.ImageCorners.ICP.LRFC.Lat=tiePointGeo(1,ll_index);
% output_meta.GeoData.ImageCorners.ICP.LRFC.Lon=tiePointGeo(2,ll_index);

%% State vectors
num_state_vectors=str2double(xp.evaluate(...
    ['count(' xpath_str({'product','sourceAttributes','orbitAndAttitude','orbitInformation','stateVector'}) ')'],...
    xml_domnode));
state_vector_T  = zeros(1,num_state_vectors);
state_vector_T_frac  = zeros(1,num_state_vectors);
state_vector_X  = zeros(1,num_state_vectors);
state_vector_Y  = zeros(1,num_state_vectors);
state_vector_Z  = zeros(1,num_state_vectors);
% state_vector_VX = zeros(1,num_state_vectors);
% state_vector_VY = zeros(1,num_state_vectors);
% state_vector_VZ = zeros(1,num_state_vectors);
for i=1:num_state_vectors
    timeStamp = char(xp.evaluate(...
        [xpath_str({'product','sourceAttributes','orbitAndAttitude','orbitInformation','stateVector'})...
        '[' num2str(i) ']' xpath_str({'timeStamp'})],...
        xml_domnode));
    [state_vector_T(i), state_vector_T_frac(i)] = datenum_w_frac(timeStamp);

    state_vector_X(i) = str2double(xp.evaluate(...
        [xpath_str({'product','sourceAttributes','orbitAndAttitude','orbitInformation','stateVector'})...
        '[' num2str(i) ']' xpath_str({'xPosition'})],...
        xml_domnode));
    state_vector_Y(i) = str2double(xp.evaluate(...
        [xpath_str({'product','sourceAttributes','orbitAndAttitude','orbitInformation','stateVector'})...
        '[' num2str(i) ']' xpath_str({'yPosition'})],...
        xml_domnode));
    state_vector_Z(i) = str2double(xp.evaluate(...
        [xpath_str({'product','sourceAttributes','orbitAndAttitude','orbitInformation','stateVector'})...
        '[' num2str(i) ']' xpath_str({'zPosition'})],...
        xml_domnode));
%     state_vector_VX(i) = str2double(xp.evaluate(...
%         [xpath_str({'product','sourceAttributes','orbitAndAttitude','orbitInformation','stateVector'})...
%         '[' num2str(i) ']' xpath_str({'xVelocity'})],...
%         xml_domnode));
%     state_vector_VY(i) = str2double(xp.evaluate(...
%         [xpath_str({'product','sourceAttributes','orbitAndAttitude','orbitInformation','stateVector'})...
%         '[' num2str(i) ']' xpath_str({'yVelocity'})],...
%         xml_domnode));
%     state_vector_VZ(i) = str2double(xp.evaluate(...
%         [xpath_str({'product','sourceAttributes','orbitAndAttitude','orbitInformation','stateVector'})...
%         '[' num2str(i) ']' xpath_str({'zVelocity'})],...
%         xml_domnode));
end
state_vector_T = round((state_vector_T-rawDataStartTime)*SECONDS_IN_A_DAY) + ... % Convert from days to secs
        (state_vector_T_frac-rawDataStartTimeFrac); % Handle fractional seconds
% sv2poly.m shows ways to determine best polynomial order, but 5th is almost always best
polyorder = min(5, numel(state_vector_T) - 1);
P_x = polyfit(state_vector_T, state_vector_X, polyorder);
P_y = polyfit(state_vector_T, state_vector_Y, polyorder);
P_z = polyfit(state_vector_T, state_vector_Z, polyorder);
output_meta.Position.ARPPoly.X = P_x(end:-1:1).';
output_meta.Position.ARPPoly.Y = P_y(end:-1:1).';
output_meta.Position.ARPPoly.Z = P_z(end:-1:1).';

%% Grid
if strcmp(char(xp.evaluate(...
    xpath_str({'product','imageGenerationParameters','generalProcessingInformation','productType'}),...
    xml_domnode)),'SLC')
    output_meta.Grid.ImagePlane = 'SLANT';
    output_meta.Grid.Type = 'RGZERO';
else
    output_meta.Grid.ImagePlane = 'GROUND';
end
output_meta.Grid.Row.SS = str2double(xp.evaluate(...
    xpath_str({'product',im_at_str,'rasterAttributes','sampledPixelSpacing'}),...
    xml_domnode));
% Col.SS is derived after DRateSFPoly below, rather than used from this
% given field, so that SICD metadata can be internally consistent.
% output_meta.Grid.Col.SS = str2double(xp.evaluate(...
%     xpath_str({'product',im_at_str,'rasterAttributes','sampledLineSpacing'}),...
%     xml_domnode));
output_meta.Grid.Row.Sgn = -1; % Always true for RS2
output_meta.Grid.Col.Sgn = -1; % Always true for RS2
fc = str2double(xp.evaluate(...
    xpath_str({'product','sourceAttributes','radarParameters','radarCenterFrequency'}),...
    xml_domnode)); % Center frequency
output_meta.Grid.Row.ImpRespBW = 2*str2double(xp.evaluate(...
    xpath_str({'product','imageGenerationParameters','sarProcessingInformation','totalProcessedRangeBandwidth'}),...
    xml_domnode))/SPEED_OF_LIGHT;
dop_bw = str2double(xp.evaluate(...
    xpath_str({'product','imageGenerationParameters','sarProcessingInformation','totalProcessedAzimuthBandwidth'}),...
    xml_domnode)); % Doppler bandwidth
[zd_last, zd_last_frac] = datenum_w_frac(char(xp.evaluate(...
    xpath_str({'product','imageGenerationParameters','sarProcessingInformation','zeroDopplerTimeLastLine'}),...
    xml_domnode)));
[zd_first, zd_first_frac] = datenum_w_frac(char(xp.evaluate(...
    xpath_str({'product','imageGenerationParameters','sarProcessingInformation','zeroDopplerTimeFirstLine'}),...
    xml_domnode)));
ss_zd_s = abs(round((zd_last-zd_first)*SECONDS_IN_A_DAY) + ... % Convert from days to secs
        (zd_last_frac-zd_first_frac))/... % Handle fractional seconds
        double(output_meta.ImageData.NumCols-1); % Image column spacing in zero doppler time (seconds)
output_meta.Grid.Row.KCtr = 2*fc/SPEED_OF_LIGHT;
output_meta.Grid.Col.KCtr = 0;
output_meta.Grid.Row.DeltaK1 = -output_meta.Grid.Row.ImpRespBW/2;
output_meta.Grid.Row.DeltaK2 = -output_meta.Grid.Row.DeltaK1;
output_meta.Grid.Row.DeltaKCOAPoly = 0;
% Constants used to compute weighting parameters
NUM_SAMPLES = 512;
OVERSAMPLE = 1024;
output_meta.Grid.Row.WgtType.WindowName = upper(char(xp.evaluate(...
    xpath_str({'product','imageGenerationParameters','sarProcessingInformation','rangeWindow','windowName'}),...
    xml_domnode)));
if strcmpi(output_meta.Grid.Row.WgtType.WindowName,'KAISER') % The usual RS2 weigting
    output_meta.Grid.Row.WgtType.Parameter.name = 'BETA';
    output_meta.Grid.Row.WgtType.Parameter.value = char(xp.evaluate(...
        xpath_str({'product','imageGenerationParameters','sarProcessingInformation','rangeWindow','windowCoefficient'}),...
        xml_domnode));
    beta_row = str2double(output_meta.Grid.Row.WgtType.Parameter.value);
    % We don't use the Mathworks Kaiser function, so we won't be dependent on the Signal Processing Toolbox
    output_meta.Grid.Row.WgtFunct = kaiser_nosptb(NUM_SAMPLES,beta_row);
    imp_resp = abs(fft(output_meta.Grid.Row.WgtFunct, round(NUM_SAMPLES*OVERSAMPLE))); % Oversampled response function
    imp_resp = imp_resp/sum(output_meta.Grid.Row.WgtFunct); % Normalize to unit peak
    ind = find(imp_resp<1/sqrt(2),1,'first')+[-1 -0]; % Samples surrounding half-power point
    ind = interp1(imp_resp(ind), ind, 1/sqrt(2)); % Linear interpolation to solve for half-power point
    row_broadening_factor = 2*(ind - 1)/OVERSAMPLE;
    output_meta.Grid.Row.ImpRespWid = row_broadening_factor/output_meta.Grid.Row.ImpRespBW;
end
output_meta.Grid.Col.WgtType.WindowName = upper(char(xp.evaluate(...
    xpath_str({'product','imageGenerationParameters','sarProcessingInformation','azimuthWindow','windowName'}),...
    xml_domnode)));

%% Radar Collection
% Ultrafine and spotlight modes have "lower" and "upper" parts to the
% pulse.
% output_meta.RadarCollection.RefFreqIndex=uint32(0); % Absence of this field means all frequencies are true values
num_pulse_parts = str2double(xp.evaluate(...
    ['count(' xpath_str({'product','sourceAttributes','radarParameters','pulseLength'}) ')'],...
    xml_domnode));
for i=1:num_pulse_parts
    output_meta.RadarCollection.Waveform.WFParameters(i).TxRFBandwidth = str2double(xp.evaluate(...
        [xpath_str({'product','sourceAttributes','radarParameters','pulseBandwidth'}) '[' num2str(i) ']'],...
        xml_domnode)); % Bandwidth
    output_meta.RadarCollection.Waveform.WFParameters(i).TxPulseLength = str2double(xp.evaluate(...
        [xpath_str({'product','sourceAttributes','radarParameters','pulseLength'}) '[' num2str(i) ']'],...
        xml_domnode));
    output_meta.RadarCollection.Waveform.WFParameters(i).RcvDemodType='CHIRP';
    sample_rate = str2double(xp.evaluate(...
        [xpath_str({'product','sourceAttributes','radarParameters','adcSamplingRate'}) '[' num2str(i) ']'],...
        xml_domnode));
    output_meta.RadarCollection.Waveform.WFParameters(i).RcvWindowLength = str2double(xp.evaluate(...
        xpath_str({'product','sourceAttributes','radarParameters','samplesPerEchoLine'}),...
        xml_domnode))/sample_rate;
    output_meta.RadarCollection.Waveform.WFParameters(i).ADCSampleRate = sample_rate;
    output_meta.RadarCollection.Waveform.WFParameters(i).RcvFMRate = 0; % True for RcvDemodType='CHIRP'
end
bw = sum([output_meta.RadarCollection.Waveform.WFParameters.TxRFBandwidth]);
output_meta.RadarCollection.TxFrequency.Min = fc-(bw/2); % fc calculated in Grid section
output_meta.RadarCollection.TxFrequency.Max = fc+(bw/2);
% Assumes pulse parts are exactly adjacent in bandwidth
output_meta.RadarCollection.Waveform.WFParameters(1).TxFreqStart = ...
    output_meta.RadarCollection.TxFrequency.Min;
for i=2:num_pulse_parts
    output_meta.RadarCollection.Waveform.WFParameters(i).TxFreqStart = ...
        output_meta.RadarCollection.Waveform.WFParameters(i-1).TxFreqStart + ...
        output_meta.RadarCollection.Waveform.WFParameters(i-1).TxRFBandwidth;
end
% Polarization
pols = textscan(char(xp.evaluate(...
    xpath_str({'product','sourceAttributes','radarParameters','polarizations'}),...
    xml_domnode)),'%s');
pols = pols{1};
M = struct('H','H','V','V','C','RHC');
tx_pols = unique(cellfun(@(x) x(1), pols));
for i=1:numel(pols)
    output_meta.RadarCollection.RcvChannels.ChanParameters(i).TxRcvPolarization = ...
        [M.(pols{i}(1)) ':' M.(pols{i}(2))];
end
if isscalar(tx_pols) % Only one transmit polarization
    output_meta.RadarCollection.TxPolarization = M.(tx_pols);
else % Multiple transmit polarizations
    output_meta.RadarCollection.TxPolarization = 'SEQUENCE';
    for i = 1:numel(tx_pols)
        output_meta.RadarCollection.TxSequence.TxStep(i).TxPolarization = M.(tx_pols(i));
    end
end
% Another way to get polarimetric channels:
% num_pol_bands = str2double(xp.evaluate(...
%     ['count(' xpath_str({'product','imageAttributes','fullResolutionImageData'}) ')'],...
%     xml_domnode));
% for i=1:num_pol_bands
%     pol = char(xp.evaluate(...
%         [xpath_str({'product','imageAttributes','fullResolutionImageData'}) '[' num2str(i) ']/@pole'],...
%         xml_domnode));
% end

%% Timeline
output_meta.Timeline.CollectStart = rawDataStartTime + (rawDataStartTimeFrac/SECONDS_IN_A_DAY);
if strcmp(gen, 'RS2')
    prf_xp_str = {'pulseRepetitionFrequency'};
elseif strcmp(gen, 'RCM')
    prf_xp_str = {'prfInformation', 'pulseRepetitionFrequency'};
end
prf = str2double(xp.evaluate(...
    xpath_str([{'product','sourceAttributes', 'radarParameters'} prf_xp_str]),...
    xml_domnode));
num_lines_entries = str2double(xp.evaluate(...
    ['count(' xpath_str({'product','imageGenerationParameters','sarProcessingInformation','numberOfLinesProcessed'}) ')'],...
    xml_domnode));
num_lines_processed = zeros(num_lines_entries,1);
for i = 1:num_lines_entries
    num_lines_processed(i) = str2double(xp.evaluate(...
        [xpath_str({'product','imageGenerationParameters','sarProcessingInformation','numberOfLinesProcessed'}) '[@pole="' pols{i} '"]'],...
        xml_domnode));
end
num_lines_processed = num_lines_processed(1) * numel(tx_pols);
if num_lines_entries ~= numel(pols) || ~all(num_lines_processed == num_lines_processed(1))
    % This should never happen, but we'll throw an error if it does.
    warning('META2SICD_RS_XML:UnableToComputeCollectDuration', 'Unhandled data condition.');
end
prf = prf * num_pulse_parts;
if num_pulse_parts==2 && ...
        strcmp(output_meta.CollectionInfo.RadarMode.ModeType,'STRIPMAP')
    % Why????
    % This seems to be necessary to make CollectDuration match CA ranges
    % and to make 1/prf roughly equal to ss_zd_s (which is generally true
    % for STRIPMAP).  But we already doubled the prf above to account for
    % the pulse parts (to make real vs effective prf), so why do we have to
    % do it again? And why don't we have to do it for SPOTLIGHT?
    prf = 2*prf;
end
output_meta.Timeline.CollectDuration = num_lines_processed/prf;
output_meta.Timeline.IPP.Set.TStart = 0;
output_meta.Timeline.IPP.Set.TEnd = 0; % Apply real value later.  Just a placeholder.
output_meta.Timeline.IPP.Set.IPPStart = uint32(0);
output_meta.Timeline.IPP.Set.IPPEnd = uint32(num_lines_processed);
output_meta.Timeline.IPP.Set.IPPPoly = [0; prf];
output_meta.Timeline.IPP.Set.TEnd = output_meta.Timeline.CollectDuration;

%% Image Formation
output_meta.ImageFormation.RcvChanProc = ...
    struct('NumChanProc', uint32(1), ... % Assumes not a MODEX collect
    'PRFScaleFactor', 1/max(num_pulse_parts, numel(tx_pols))); % Either polarimetric or multi-step, but not both.
output_meta.ImageFormation.ImageFormAlgo = 'RMA';
output_meta.ImageFormation.TStartProc = 0;
output_meta.ImageFormation.TEndProc = output_meta.Timeline.CollectDuration;
output_meta.ImageFormation.TxFrequencyProc.MinProc = ...
    output_meta.RadarCollection.TxFrequency.Min;
output_meta.ImageFormation.TxFrequencyProc.MaxProc = ...
    output_meta.RadarCollection.TxFrequency.Max;
output_meta.ImageFormation.STBeamComp = 'GLOBAL';
output_meta.ImageFormation.ImageBeamComp = 'SV';
output_meta.ImageFormation.AzAutofocus = 'NO';
output_meta.ImageFormation.RgAutofocus = 'NO';

%% RMA.INCA
output_meta.RMA.RMAlgoType = 'OMEGA_K';
output_meta.RMA.ImageType = 'INCA';
output_meta.SCPCOA.SideOfTrack = char(xp.evaluate(...
    xpath_str({'product','sourceAttributes','radarParameters','antennaPointing'}),...
    xml_domnode));  % Should always be right looking for RCM
output_meta.SCPCOA.SideOfTrack = upper(output_meta.SCPCOA.SideOfTrack(1));
if output_meta.SCPCOA.SideOfTrack=='L'
    ss_zd_s = -ss_zd_s;
    % In addition to left/right, RS2 data can independently be in
    % increasing/decreasing line order.
    if (round((zd_first-zd_last)*SECONDS_IN_A_DAY) + ...
        (zd_first_frac-zd_last_frac)) < 0 % zd_last occurred after zd_first
        zd_first = zd_last;
        zd_first_frac = zd_last_frac;
    end
    look = 1;
else
    if (round((zd_first-zd_last)*SECONDS_IN_A_DAY) + ...
        (zd_first_frac-zd_last_frac)) > 0 % zd_last occurred before zd_first
        zd_first = zd_last;
        zd_first_frac = zd_last_frac;
    end
    look = -1;
end
% Zero doppler time of SCP relative to collect start
zd_t_scp = round((zd_first-rawDataStartTime)*SECONDS_IN_A_DAY) + ... % Convert days to seconds
    (zd_first_frac-rawDataStartTimeFrac) + ... % Handle fractional seconds
    (double(output_meta.ImageData.SCPPixel.Col) * ss_zd_s);
if strcmp(gen, 'RS2')
    near_range = str2double(xp.evaluate(...
        xpath_str({'product','imageGenerationParameters','sarProcessingInformation','slantRangeNearEdge'}),...
        xml_domnode)); % in meters
elseif strcmp(gen, 'RCM')
    near_range = str2double(xp.evaluate(...
        xpath_str({'product','sceneAttributes','imageAttributes','slantRangeNearEdge'}),...
        xml_domnode)); % in meters
end
output_meta.RMA.INCA.R_CA_SCP = near_range + ...
    (double(output_meta.ImageData.SCPPixel.Row)*output_meta.Grid.Row.SS);
output_meta.RMA.INCA.FreqZero = fc;
% Doppler Rate (We do this first since some other things are dependent on
% it.)
% For the purposes of the DRateSFPoly computation, we ignore any
% changes in velocity over the azimuth dimension.
pos_coefs = [P_x(:) P_y(:) P_z(:)];
% Velocity is derivate of position.
vel_coefs=pos_coefs(1:end-1,:).*repmat(((size(pos_coefs,1)-1):-1:1)',[1 3]);
vel_x = polyval(vel_coefs(:,1), zd_t_scp);
vel_y = polyval(vel_coefs(:,2), zd_t_scp);
vel_z = polyval(vel_coefs(:,3), zd_t_scp);
vm_ca_sq = vel_x.^2 + vel_y.^2 + vel_z.^2; % Magnitude of the velocity squared
r_ca = [output_meta.RMA.INCA.R_CA_SCP; 1]; % Polynomial representing range as a function of range distance from SCP
if strcmp(gen, 'RS2')
    drc_xp_str = {'product','imageGenerationParameters','dopplerRateValues', 'dopplerRateValuesCoefficients'};
elseif strcmp(gen, 'RCM')
    drc_xp_str = {'product','dopplerRate','dopplerRateEstimate', 'dopplerRateCoefficients'};
end
drr_xp_str = [drc_xp_str{1:(end-1)} {'dopplerRateReferenceTime'}];
dop_rate_coefs = str2num(xp.evaluate(xpath_str(drc_xp_str),... % Multiple numbers.  We need str2num instead of str2double
    xml_domnode)); % Shifted (origin at dop_rate_ref_t, not SCP) and scaled (sec, not m) version of SICD DopCentroidPoly
dop_rate_ref_t = str2double(xp.evaluate(xpath_str(drr_xp_str),...
    xml_domnode)); % Reference time of Doppler rate polynomial
dop_rate_coefs_shifted = polyshift(dop_rate_coefs, ... % Shift so SCP is reference
    (output_meta.RMA.INCA.R_CA_SCP*2/SPEED_OF_LIGHT) - ... % SICD reference time (SCP)
    dop_rate_ref_t); % Reference time of native Doppler Centroid polynomial
dop_rate_coefs_scaled = dop_rate_coefs_shifted .*  ... % Scale from seconds to meters
    (2/SPEED_OF_LIGHT) .^ (0:(length(dop_rate_coefs)-1));
output_meta.RMA.INCA.DRateSFPoly = - conv(dop_rate_coefs_scaled.',r_ca) * ... % Multiplication of two polynomials is just a convolution of their coefficients
    SPEED_OF_LIGHT / (2 * fc * vm_ca_sq(1)); % Assumes a SGN of -1

%% Fields dependent on Doppler rate
% This computation of SS is actually better than the claimed SS
% (sampledLineSpacing) in many ways, because this makes all of the metadata
% internally consistent.  This must be the sample spacing exactly at SCP
% (which is the definition for SS in SICD), if the other metadata from
% which is it computed is correct and consistent. Since column SS can vary
% slightly over a RGZERO image, we don't know if the claimed sample spacing
% in the native metadata is at our chosen SCP, or another point, or an
% average across image or something else.
output_meta.Grid.Col.SS = sqrt(vm_ca_sq(1)) * abs(ss_zd_s) * ...
    output_meta.RMA.INCA.DRateSFPoly(1,1);
output_meta.Grid.Col.ImpRespBW = dop_bw*abs(ss_zd_s)/output_meta.Grid.Col.SS; % Convert to azimuth spatial bandwidth (cycles per meter)
output_meta.RMA.INCA.TimeCAPoly = [zd_t_scp; ss_zd_s/output_meta.Grid.Col.SS];
if strcmpi(output_meta.Grid.Col.WgtType.WindowName,'KAISER') % The usual RS2 weigting
    output_meta.Grid.Col.WgtType.Parameter.name = 'BETA';
    output_meta.Grid.Col.WgtType.Parameter.value = char(xp.evaluate(...
        xpath_str({'product','imageGenerationParameters','sarProcessingInformation','azimuthWindow','windowCoefficient'}),...
        xml_domnode));
    beta_col = str2double(output_meta.Grid.Col.WgtType.Parameter.value);
    % We don't use the Mathworks Kaiser function, so we won't be dependent on the Signal Processing Toolbox
    output_meta.Grid.Col.WgtFunct = kaiser_nosptb(NUM_SAMPLES,beta_col);
    imp_resp = abs(fft(output_meta.Grid.Col.WgtFunct, round(NUM_SAMPLES*OVERSAMPLE))); % Oversampled response function
    imp_resp = imp_resp/sum(output_meta.Grid.Col.WgtFunct); % Normalize to unit peak
    ind = find(imp_resp<1/sqrt(2),1,'first')+[-1 -0]; % Samples surrounding half-power point
    ind = interp1(imp_resp(ind), ind, 1/sqrt(2)); % Linear interpolation to solve for half-power point
    col_broadening_factor = 2*(ind - 1)/OVERSAMPLE;
    output_meta.Grid.Col.ImpRespWid = col_broadening_factor/output_meta.Grid.Col.ImpRespBW;
end

%% Doppler Centroid
if strcmp(gen, 'RS2')
    dc_xp_str = {'product','imageGenerationParameters','dopplerCentroid'};
elseif strcmp(gen, 'RCM')
    dc_xp_str = {'product','dopplerCentroid','dopplerCentroidEstimate'};
end
dop_cent_coefs = str2num(xp.evaluate(... % Multiple numbers.  We need str2num instead of str2double
    xpath_str([dc_xp_str,{'dopplerCentroidCoefficients'}]),...
    xml_domnode)); % Shifted (origin at dop_cent_ref_t, not SCP) and scaled (sec, not m) version of SICD DopCentroidPoly
dop_cent_ref_t = str2double(xp.evaluate(...
    xpath_str([dc_xp_str,{'dopplerCentroidReferenceTime'}]),...
    xml_domnode)); % Reference time of Doppler Centroid polynomial
dop_cent_coefs_shifted = polyshift(dop_cent_coefs, ... % Shift so SCP is reference
    (output_meta.RMA.INCA.R_CA_SCP*2/SPEED_OF_LIGHT) - ... % SICD reference time (SCP)
    dop_cent_ref_t); % Reference time of native Doppler Centroid polynomial
dop_cent_coefs_scaled = dop_cent_coefs_shifted .*  ... % Scale from seconds to meters
    (2/SPEED_OF_LIGHT) .^ (0:(length(dop_cent_coefs)-1));
output_meta.RMA.INCA.DopCentroidPoly=dop_cent_coefs_scaled.';
% Adjust Doppler Centroid for spotlight
if strcmp(output_meta.CollectionInfo.RadarMode.ModeType,'SPOTLIGHT')
    [dop_est, dop_est_frac] = datenum_w_frac(char(xp.evaluate(...  % Doppler estimate time
        xpath_str([dc_xp_str 'timeOfDopplerCentroidEstimate']), xml_domnode)));
    dop_est_t = abs(round((dop_est - rawDataStartTime)*SECONDS_IN_A_DAY) + ... % Convert from days to secs
        (dop_est_frac - rawDataStartTimeFrac)); % Handle fractional seconds
    dop_est_col = (dop_est_t - zd_t_scp)/ss_zd_s; % This is the column where the doppler centroid was computed.
    % Column-dependent variation in DopCentroidPoly due to spotlight
    output_meta.RMA.INCA.DopCentroidPoly(1,2) = ...
        -look * fc * (2 / SPEED_OF_LIGHT) * sqrt(vm_ca_sq(1)) / ...
        output_meta.RMA.INCA.R_CA_SCP;
    % dopplerCentroid in native metadata was defined at specific column,
    % which might not be our SCP column.  Adjust so that SCP column is
    % correct.
    output_meta.RMA.INCA.DopCentroidPoly(1,1) = ...
        output_meta.RMA.INCA.DopCentroidPoly(1,1) - ...
        (output_meta.RMA.INCA.DopCentroidPoly(1,2) * ...
        dop_est_col * output_meta.Grid.Col.SS);
end
output_meta.Grid.Col.DeltaKCOAPoly = ...
    output_meta.RMA.INCA.DopCentroidPoly * ss_zd_s / output_meta.Grid.Col.SS;
% Compute Col.DeltaK1/K2 from DeltaKCOAPoly
% This is not always straightforward to do generically for any possible
% DeltaKCOAPoly 2D polynomial, since you would have to compute all 2D roots
% and edge cases.  However, for the RS case, this can be solved exactly,
% since its usually a 1D polynomial.  Even the spotlight case only brings
% in a linear variation in the second dimensions, so its still easily
% solved.
% Min/max in row/range must exist at edges or internal local min/max
minmax = roots(polyder(output_meta.Grid.Col.DeltaKCOAPoly(end:-1:1,1)));
rg_bounds_m = (double([0 (output_meta.ImageData.NumRows-1)]) - ...
     double(output_meta.ImageData.SCPPixel.Row)) * output_meta.Grid.Row.SS;
possible_bounds_rg = [rg_bounds_m minmax(minmax>min(rg_bounds_m) & minmax<max(rg_bounds_m))];
% Constant or (linearly increasing\decreasing for spotlight) in column, so
% edges must contain max/min.
possible_bounds_az = (double([0 (output_meta.ImageData.NumCols-1)]) - ...
     double(output_meta.ImageData.SCPPixel.Col)) * output_meta.Grid.Col.SS;
possible_bounds_deltak = sicd_polyval2d(output_meta.Grid.Col.DeltaKCOAPoly, ...
    possible_bounds_az, possible_bounds_rg);
output_meta.Grid.Col.DeltaK1 = min(possible_bounds_deltak(:)) - ...
    (output_meta.Grid.Col.ImpRespBW/2);
output_meta.Grid.Col.DeltaK2 = max(possible_bounds_deltak(:)) + ...
    (output_meta.Grid.Col.ImpRespBW/2);
% Wrapped spectrum
if (output_meta.Grid.Col.DeltaK1 < -(1/output_meta.Grid.Col.SS)/2) || ...
        (output_meta.Grid.Col.DeltaK2 > (1/output_meta.Grid.Col.SS)/2)
    output_meta.Grid.Col.DeltaK1 = -(1/output_meta.Grid.Col.SS)/2;
    output_meta.Grid.Col.DeltaK2 = -output_meta.Grid.Col.DeltaK1;
end
%% TimeCOAPoly
% TimeCOAPoly=TimeCA+(DopCentroid/dop_rate)
% Since we can't evaluate this equation analytically, we will evaluate
% samples of it across our image and fit a 2D polynomial to it.
POLY_ORDER = 2; % Order of polynomial which we want to compute
grid_samples = POLY_ORDER + 1; % in each dimension
coords_az_m = linspace(-double(output_meta.ImageData.SCPPixel.Col),...
    double(output_meta.ImageData.NumCols-output_meta.ImageData.SCPPixel.Col-1), grid_samples) * ...
    output_meta.Grid.Col.SS;
coords_rg_m = linspace(-double(output_meta.ImageData.SCPPixel.Row),...
    double(output_meta.ImageData.NumRows-output_meta.ImageData.SCPPixel.Row-1), grid_samples) * ...
    output_meta.Grid.Row.SS;
timeca_sampled = sicd_polyval2d(output_meta.RMA.INCA.TimeCAPoly(:).',coords_az_m,coords_rg_m);
dopcentroid_sampled = sicd_polyval2d(output_meta.RMA.INCA.DopCentroidPoly,coords_az_m,coords_rg_m);
doprate_sampled = sicd_polyval2d(dop_rate_coefs_scaled.',coords_az_m,coords_rg_m);
timecoapoly_sampled = timeca_sampled+(dopcentroid_sampled./doprate_sampled);
% Least squares fit for 2D polynomial
% A*x = b
[coords_az_m, coords_rg_m] = ndgrid(coords_az_m, coords_rg_m);
a = zeros(grid_samples^2, (POLY_ORDER+1)^2);
for i = 0:POLY_ORDER
    for j = 0:POLY_ORDER
        a(:,i*(POLY_ORDER+1)+j+1) = (coords_rg_m(:).^j).*(coords_az_m(:).^i);
    end
end
b_coa = zeros((POLY_ORDER+1)^2,1);
for i=1:((POLY_ORDER+1)^2)
   b_coa(i)=sum(timecoapoly_sampled(:).*a(:,i)); % center of aperture
end
A=zeros((POLY_ORDER+1)^2);
for i=1:((POLY_ORDER+1)^2)
    for j=1:((POLY_ORDER+1)^2)
        A(i,j)=sum(a(:,i).*a(:,j));
    end
end
old_warning_state=warning('off','MATLAB:nearlySingularMatrix');
x=A\b_coa; % MATLAB often flags this as badly scaled, but results still appear valid
warning(old_warning_state);
output_meta.Grid.TimeCOAPoly=reshape(x, POLY_ORDER+1, POLY_ORDER+1);
if strcmp(output_meta.CollectionInfo.RadarMode.ModeType,'SPOTLIGHT')
    output_meta.Grid.TimeCOAPoly = output_meta.Grid.TimeCOAPoly(1);
    % This field required to compute TimeCOAPoly, but not allowed for
    % spotlight in SICD.
    output_meta.RMA.INCA = rmfield(output_meta.RMA.INCA, 'DopCentroidPoly');
else % This field also not allowed for spotlight in SICD.
    output_meta.RMA.INCA.DopCentroidCOA=true;
end

%% GeoData
% Now that sensor model fields have been populated, we can populate
% GeoData.SCP more precisely.
ecf = point_image_to_ground([output_meta.ImageData.SCPPixel.Row;output_meta.ImageData.SCPPixel.Col],output_meta);
output_meta.GeoData.SCP.ECF.X=ecf(1);
output_meta.GeoData.SCP.ECF.Y=ecf(2);
output_meta.GeoData.SCP.ECF.Z=ecf(3);
llh=ecf_to_geodetic([output_meta.GeoData.SCP.ECF.X output_meta.GeoData.SCP.ECF.Y output_meta.GeoData.SCP.ECF.Z]);
output_meta.GeoData.SCP.LLH.Lat=llh(1);
output_meta.GeoData.SCP.LLH.Lon=llh(2);
output_meta.GeoData.SCP.LLH.HAE=llh(3);

%% Radiometric
if exist('beta_domnode','var')
    % Offset always zero for SLC, so we only need gains
    betas = str2num(xp.evaluate(xpath_str({'lut','gains'}),beta_domnode)); %#ok<ST2NM>
    % Of the provided LUTs, we really only work with beta here, since it is
    % the radiometric term most commonly kept constant, and the others
    % (sigma/gamma) can always be derived from it.
    % if any(strcmp({'Constant-beta', 'Point Target', 'Point Target-1', ...
    %         'Calibration-1', 'Calibration-2', 'Ship-1', 'Ship-2', ...
    %         'Ship-3', 'Unity'}, ... % Known modes with constant beta
    %         xp.evaluate(xpath_str({'product','imageGenerationParameters', ...
    %         'sarProcessingInformation','lutApplied'}), xml_domnode)))
    if all(betas==betas(1)) % In case we missed some modes, this condition may be more reliable
        output_meta.Radiometric.BetaZeroSFPoly = 1/betas(1)^2;
    else  % Otherwise fit a 1D polynomial in range
        % RS2 has value for every row
        coords_rg_m = (double(0:(output_meta.ImageData.NumRows-1)) - ...
            double(output_meta.ImageData.SCPPixel.Row)) * output_meta.Grid.Row.SS;
        if strcmp(gen, 'RCM') % RCM subsamples the rows
            rng_indices = ((str2double(xp.evaluate(xpath_str({'lut','pixelFirstAnglesValue'}), beta_domnode)):...
                ... % Should be this, but simulated datasets are inconsistent with spec document
                ... % crm_indices = (xp.evaluate(xpath_str({'lut','pixelFirstLutValue'}), beta_domnode):...
                (str2double(xp.evaluate(xpath_str({'lut','numberOfValues'}), beta_domnode))-1)) * ... % First row is zero
                str2double(xp.evaluate(xpath_str({'lut','stepSize'}), beta_domnode))) + 1;
            coords_rg_m = coords_rg_m(rng_indices);
        end
        % For the datasets we have seen, this function is very close to
        % linear.  For the "Mixed" LUT, there is a tiny linear piecewise
        % deviation from a single overall linear.
        betapoly = polyfit(coords_rg_m, 1./betas.^2, 2);
        output_meta.Radiometric.BetaZeroSFPoly = betapoly(end:-1:1).';
    end
    % RCS, Sigma, and Gamma will be computed below in derived_sicd_fields
    % derived_sicd_fields.
    output_meta.Radiometric.NoiseLevel.NoiseLevelType = 'ABSOLUTE';
    if strcmp(gen, 'RS2')
        beta0_str = [xpath_str({'product', 'sourceAttributes', ...
            'radarParameters','referenceNoiseLevel'}) ...
            '[@incidenceAngleCorrection="Beta Nought"]'];
        noise_domnode = xml_domnode;  % RS2 noise is in main product.xml
    elseif strcmp(gen, 'RCM') % RCM noise is in separate file
        beta0_str = [xpath_str({'noiseLevels', 'referenceNoiseLevel'}) ...
            '/*[text()=''Beta Nought'']/..'];
    end
    if exist('noise_domnode','var')
        pfv = str2double(xp.evaluate([beta0_str ...
            '/*[local-name()=''pixelFirstNoiseValue'']'], noise_domnode)); % Index of first row defined to be zero
        step = str2double(xp.evaluate([beta0_str ...
            '/*[local-name()=''stepSize'']'], noise_domnode));
        beta0s = str2num(xp.evaluate([beta0_str ...
            '/*[local-name()=''noiseLevelValues'']'], noise_domnode)); %#ok<ST2NM>
        range_coords = output_meta.Grid.Row.SS * ...
            ((((1:numel(beta0s))-1) * step) + pfv - double(output_meta.ImageData.SCPPixel.Row));
        noisepoly = polyfit(range_coords, beta0s - (10*log10(polyval(...
            output_meta.Radiometric.BetaZeroSFPoly(end:-1:1), range_coords))), 2);
        output_meta.Radiometric.NoiseLevel.NoisePoly = noisepoly(end:-1:1).';
    end
end

%% SCPCOA
% All of these fields are derivable for more fundamental fields.
output_meta = derived_sicd_fields(output_meta);

%% Process fields specific to each polarimetric band
band_independent_meta = output_meta; % Values that are consistent across all bands
grouped_meta = cell(numel(pols),1);
for i=1:numel(pols)
    output_meta = band_independent_meta;

    output_meta.ImageFormation.RcvChanProc.ChanIndex = i;
    output_meta.ImageFormation.TxRcvPolarizationProc = ...
        output_meta.RadarCollection.RcvChannels.ChanParameters(i).TxRcvPolarization;

    grouped_meta{i} = output_meta;
end
output_meta = grouped_meta; % Cell array with metadata struct for each band

end

% Creates an XPath query that is a namespace insensitive search for an XML
% heirachy.  Input is a cell array of strings specifying the XML fields,
% starting with the root node and working its way down the XML tree to the
% leaf nodes.  This is required since RS2 XML files specify a default
% namespace, so XPath queries with no namespace speficied will not work.
% (That only searches for elements associated with no namespace at all.)
% Note: This seems to be an issue only in MATLAB 2010a and above, as
% previous versions used another XPATH library that was more tolerant.
function out_str = xpath_str(in_cell_array)
    out_str='';
    for j=1:length(in_cell_array)
        out_str=[out_str '/*[local-name()=''' in_cell_array{j} ''']'];
    end
end

% MATLAB's datenum function won't handle precise times down under a
% millisecond, because 1) It won't accept a format with more than 3 .FFF in
% the string description of the date format and 2) the resulting serial
% date number is stored in days from 00-JAN-0000 and just doesn't have
% enough bits to handle fractional seconds to the level we want.  Here we
% handle the fractional seconds separately so we can read date with the
% precision we need.
function [datenum_s, datenum_frac] = datenum_w_frac(datestring)
    datenum_s = datenum(datestring,'yyyy-mm-ddTHH:MM:SS');
    datenum_frac = str2double(regexp(datestring,'\.\d*','match'));
    if isnan(datenum_frac), datenum_frac = 0; end;
end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
% time varying Rayleigh flat fading channel with Clark's model
% using filter method
% Ns number of channel samples
% fDTs: fD*Ts doppler spread fD*Ts
% fdTs= .001 slow, .01 medioum, .03 fast

function out=flat_filter(Ns,fD,fDTs)

% parameters
Ts=fDTs/fD;
% filter window length at least three side lobs
Ng=ceil(2/fDTs);

% Complex Gaussian input
x=cxn(Ns,1);

% compute g(t) and g_hat(t)
t=(-Ng:Ng)*Ts;
g=besselj(1/4,2*pi*fD*abs(t))./((abs(t)).^(1/4));
g(Ng+1)=((pi*fD)^(1/4))/gamma(5/4);
%verify filter
%plot(abs(g));

% compute K such that ...
g_hat=g./sqrt(sum(g.^2));
% verify K
%sum_g_hat=sum(g_hat.^2);

% channel impulse response
cn=conv(x,g_hat)*Ts;

% reject transient
cn=cn(2*Ng+1:length(cn)-2*Ng);
out=cn/sqrt(var(cn));

end


text:
clear all; close all; clc;

M = fixedTTrankfactory([5 5 5], [1 5 5 1]);

checkmanifold(M);

author: NicolasBoumal

repo: manopt


text:
function varargout = grad(varargin)
% VL_GRAD Compute the gradient of an image
%   [IX,IY] = VL_GRAD(I) returns the gradient components IX,IY of the
%   2-D discrete function I. I must be a two-dimensional
%   matrix. VL_GRAD() computes the gradient by using finite
%   differences; specifically, it uses central differences for all but
%   the boundary pixels, for which it uses forward/backward
%   differences as appropriate.
%
%   Remark::
%     VL_GRAD() is similar to the MATLAB built-in GRADIENT() function,
%     excepts that it supports different gradient approximations.
%
%   VL_GRAD() accepts the following options:
%
%   Type:: central
%     Specify which type of finite differences to use for all but the
%     boundary samples. TYPE can be one of 'central', 'forward', or
%     'backward'.
%
%   See also: GRADIENT(), VL_HELP().
[varargout{1:nargout}] = vl_grad(varargin{:});

author: yihui-he

repo: panorama

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


text:
function indexFoundPattern = tapas_physio_find_string(stringArray, searchPattern, isExact)
% Finds a string pattern in a cell of strings and returns found cell indices
%
% NOTE: this function behaves the same as find should behave
%
%   indexFoundPattern = find_string(stringArray, searchPattern, isExact)
%
% IN
%       stringArray     cell(nStrings,1) of cells to be checked for search
%                       pattern
%       searchPattern   string to be searched for (regular expressions
%                       possible
%                       OR cell array of search patterns
%       isExact         default:false; if true, only exact matched of
%                       searchPattern to strings are returned
% OUT
%       indexFoundPattern vector of indices in stringArray that match
%                       searchPattern
%                       OR, if searchPattern was cell array:
%                       cell(nPatterns,1) of index vectors, one cell
%                       element for each pattern
% EXAMPLE
%   find_string
%
%   See also

% Author: Lars Kasper
% Created: 2014-11-05
% Copyright (C) 2014 Institute for Biomedical Engineering, ETH/Uni Zurich.

if ~iscell(stringArray)
    stringArray = {stringArray};
end

if nargin < 3
    isExact = false;
end

if ~iscell(searchPattern)

    if isExact
        % enforce start and end of string match that in sought-through
        searchPattern = ['^' searchPattern '$'];
    end

    indexFoundPattern = find(~cell2mat(cellfun(@isempty, regexp(stringArray, ...
        searchPattern), ...
        'UniformOutput', false)));
else

    nPatterns = numel(searchPattern);
    indexFoundPattern = cell(nPatterns,1);
    for iPattern = 1:nPatterns

        if isExact
            % enforce start and end of string match that in sought-through
            searchPattern{iPattern} = ['^' searchPattern{iPattern} '$'];
        end

        indexFoundPattern{iPattern} = find(~cell2mat(cellfun(@isempty, regexp(stringArray, ...
            searchPattern{iPattern}), ...
            'UniformOutput', false)));
    end
end

end

author: translationalneuromodeling

repo: tapas


text:
% Synchronoise the streams such that they will work with the network.
% Things did in this function:
%   1. upsample streams to match other streams of higher frame rate. This
%   is only done when streams are specified by DataSyncSet to have the same
%   frame rate. Also discarded last frames longer streams to make the
%   streams have the same number of frames.
%   2. apply VAD to some streams. Here VAD is stored in one stream
%   (normally the last stream). We first upsample/downsample VAD framerate
%   to match the stream to be VADed. Then remove nonspeech segments. VAD
%   should be stored as 1-D feature stream.
%
% Author: Xiong Xiao, Temasek labs, NTU, Singapore.
% Date Created: 10 Oct 2015
% Last Modified: 21 Apr 2016
%
function [data_sync, isVAD] = SynchronizeDataStreams3(alldata, para)
[nStream, nUtt] = size(alldata);

isTensor = para.IO.isTensor;
isVAD = zeros(nStream,1);
data_sync = alldata;

ApplyVADSet = para.IO.ApplyVADSet;
DataSyncSet = para.IO.DataSyncSet;

for vi=1:length(ApplyVADSet)
    if ~isempty(ApplyVADSet{vi})
        isVAD(ApplyVADSet{vi}(end))=1;      % set whether a stream is VAD stream or not
    end
end

if isempty(DataSyncSet) && isempty(ApplyVADSet)
    return;
end

for utt_i = 1:nUtt
    % Apply VAD and segmentation first, if required
    for vi=1:length(ApplyVADSet)
        curr_vad_set = ApplyVADSet{vi};
        if isempty(curr_vad_set); continue; end
        data_sync{curr_vad_set(1), utt_i} = SynchronizeStreamsVAD( data_sync(curr_vad_set, utt_i), para.IO.vadAction{vi} );
    end

    % Apply Frame rate synchronization to all Sync Sets
    for di=1:length(DataSyncSet)
        curr_sync_set = DataSyncSet{di};
        if isempty(curr_sync_set); continue; end
        data_sync(curr_sync_set, utt_i) = SynchronizeStreamsFrameRate2(  data_sync(curr_sync_set, utt_i) , para.IO.frame_rate(curr_sync_set) , isTensor(curr_sync_set)  );
    end
end
end

author: singaxiong

repo: SignalGraph


text:
function jed = epoch_to_jed_bahai ( )

%*****************************************************************************80
%
%% EPOCH_TO_JED_BAHAI returns the epoch of the Bahai calendar as a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real JED, the Julian Ephemeris Date of the epoch.
%
  jed = 2394646.5;

  return
end


text:
function  [volselpts] = mrSelLineOblVol(obPts,obSlice,obSize,obMin,obMax,sagX,sagY,volselpts,sagSize)
%
% MRSELLINEOBLVOL
%
%	 [volselpts] = mrSelLineOblVol(obSlice,obSize,obMin,obMax,sagX,sagY,volselpts,sagSize)
%
%	Select ROI points from the current oblique image.
%

[thexs,theys] = mrSelLineVol(obSlice,obSize,obMin,obMax,sagX,sagY);
obcoords = round(theys +round(thexs-1)*obSize(1));
obsel = obPts(obcoords,:);
sagcoords = round(obsel(:,2) +round(obsel(:,1)-1)*sagSize(1));
plcoords = round(obsel(:,3));
newsel = sagcoords+(plcoords-1)*prod(sagSize);
volselpts = [volselpts,newsel];

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function cheby_t_poly_test ( )

%*****************************************************************************80
%
%% CHEBY_T_POLY_TEST tests CHEBY_T_POLY.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 March 2012
%
%  Author:
%
%    John Burkardt
%
  n_max = 12;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHEBY_T_POLY_TEST:\n' );
  fprintf ( 1, '  CHEBY_T_POLY evaluates the Chebyshev T polynomial.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N      X        Exact F       T(N)(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, x, fx ] = cheby_t_poly_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = cheby_t_poly ( 1, n, x );

    fprintf ( 1, '  %6d  %8f  %12f  %12f\n', n, x, fx, fx2(n+1) );

  end

  return
end


text:
function f = plus(f, g)
%+   Addition of two TRIGTECH objects.
%   F + G adds F and G, where F and G may be TRIGTECH objects or scalars.
%
%   If F is an array-valued TRIGTECH, then F + C is supported if C is a row
%   vector of doubles with the same number of columns as F.
%
% See also MINUS, UPLUS.

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

if ( isempty(f) || isempty(g) ) % TRIGTECH + [] = [].

    f = [];

elseif ( isa(g, 'double') ) % TRIGTECH + double.

    % Update values (use bsxfun() to handle the case in which g is a vector
    % and f is an array-valued TRIGTECH):
    f.values = bsxfun(@plus, f.values, g);

    % Update coeffs:
    if ( (size(g, 2) > 1) && (size(f.coeffs, 2) == 1) )
        % Perform singleton expansion of f:
        f.coeffs = repmat(f.coeffs, 1, size(g, 2));
    end
    N = size(f.coeffs,1);

    % Determine the index in the coefficient vector where the constant term
    % is stored. The way we have arranged the coefficients means it should
    % be in the middle of the array, but that depends on whether the number
    % of coefficients is even or odd.
    if mod(N,2) == 1
        const_index = (N+1)/2;
    else
        const_index = N/2+1;
    end
    f.coeffs(const_index,:) = f.coeffs(const_index,:) + g;

    % Update isReal:
    f.isReal = f.isReal & isreal(g);

elseif ( isa(f, 'double') ) % double + TRIGTECH.

    % Switch argument order and call TRIGTECH/PLUS again:
    f = plus(g, f);

elseif ( isa(f, 'trigtech') && isa(g, 'trigtech') )  % TRIGTECH + TRIGTECH.

    % We will simply add the values together then compute the coefficients
    % of the result. This is probably not the most efficient means of
    % determing the sum.
    nf = size(f.values, 1);
    ng = size(g.values, 1);
    if ( nf > ng )
        % Increase the length of g (via PROLONG):
        g = prolong(g, nf);
    elseif ( nf < ng )
        % Increase the length of f (via PROLONG):
        f = prolong(f, ng);
    end

    % Update values and coefficients:
    f.values = f.values + g.values;
    f.coeffs = f.vals2coeffs(f.values);

    % Update isReal:
    f.isReal = f.isReal & g.isReal;

    % Force the values to be real where f is real.
    f.values(:,f.isReal) = real(f.values(:,f.isReal));

    % Look for a zero output:
    if ( ~any(f.values(:)) || ~any(f.coeffs(:)) )
        % Create a zero TRIGTECH:
        ishappy = f.ishappy && g.ishappy;
        z = zeros(1, size(f.values, 2));
        data.vscale = z;
        f = f.make(z, data);
        f.ishappy = ishappy;
    else
        f.ishappy = f.ishappy && g.ishappy;
    end

else    % Don't know how to do the addition of the objects.

    error('CHEBFUN:TRIGTECH:plus:typeMismatch',['Incompatible operation between objects.\n', ...
    'Make sure functions are of the same type.']);

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ptheta] = tapas_sem_prosa_gaussian_priors()
%% Generates standard priors for a gaussian distribution.
%
%   Input
%   ptheta      -- Structure with the priors.
%

% aponteeduardo@gmail.com
% copyright (C) 2016
%

DIM_THETA = tapas_sem_prosa_ndims();

ptheta = struct();

[mmu, mvu, vmu, vvu, me, ve, ml, vl, p0m, p0v] = ...
    tapas_sem_unified_gaussian_priors();

mu = repmat([mmu, vmu], 1, 3);
ptheta.mu = [mu, me, ml, p0m]';

pm = repmat([mvu, vvu], 1, 3);
ptheta.pm = [1./[pm, ve, vl] p0v]';

ptheta.p0 = ptheta.mu;
ptheta.p0(9) = tan(pi * (-0.4));
% Eta is beta distributed
ptheta.bdist = [9];

end


text:
function pass = test_size()
% Test SIZE

% Check two components:
F = chebfun3v(@(x,y,z) cos(x), @(x,y,z) sin(x+y+z));
pass(1) = all(size(F) == [2 inf inf inf]);

% Check three components:
F = chebfun3v(@(x,y,z) cos(x), @(x,y,z) sin(y), @(x,y,z) cos(z));
pass(2) = all(size(F) == [3 inf inf inf]);

end


text:
function [E,v] = km_kpca(X,m,ktype,kpar)
% KM_KPCA performs kernel principal component analysis (KPCA) on a data set
% X.
% Input:	- X: data matrix in column format (each data point is a row)
%			- m: the number of principal components to return. If m is
%			smaller than 1, it is interpreted as the fraction of the signal
%			energy that is to be contained within the returned principal
%			components.
%			- ktype: string representing kernel type.
%			- kpar: vector containing the kernel parameters.
% Output:	- E: matrix containing the principal components.
%			- v: array containing the eigenvalues.
%			- Xep: projections of Xe on principal directions
% USAGE: [E,v] = km_kpca(X,m,ktype,kpar)
%
% Author: Steven Van Vaerenbergh (steven *at* gtas.dicom.unican.es), 2010.
%
% This file is part of the Kernel Methods Toolbox for MATLAB.
% https://github.com/steven2358/kmbox

n = size(X,1);

K = km_kernel(X,X,ktype,kpar);
[E,V] = eig(K);
v = diag(V);	% eigenvalues
[v,ind] = sort(v,'descend');
v = v(1:m);
E = E(:,ind(1:m));	% principal components
for i=1:m
	E(:,i) = E(:,i)/sqrt(n*v(i));	% normalization
end


function [d,fsign] = checkrestrictions2(restriction,y)

% Check the restrictions
d     = 0;
count = 0;
fsign = 1;

for ii = 1 : size(restriction,2)
    tmp = eval(restriction{ii});
    count = count + min(tmp);
end
if isempty(count)==1
    error('There is a nan in the narrative restrictions.');
end
% check if -1*IRF are verified
y        = -y;
count1   = 0;
for ii = 1 : size(restriction,2)
    tmp = eval(restriction{ii});
    count1 = count1 + min(tmp);
end
if count == size(restriction,2) || count1 == size(restriction,2) % if all signs are verified stop
    d=1;
end
if count1 == size(restriction,2)
    fsign = -1;
end


text:
% Read single/multi channel waveforms
% Inputs:
%   files: a cell array of files to be read. Each cell represents one
%     sentence could be a cell and could be array of files itself if it is
%     multi-channel recording. A filename may also contains the starting
%     and ending time to be read, e.g. "abc.wav 0.2 0.5" will read the
%     segment from 0.2s to 0.5s of abc.wav.
%   reader: a structure that may contains following fields
%       fs: sampling frequency, if known in advance.
%       array: set to 1 if the sentence has multiple channels
%       multiArrayFiles: set to 1 if each channel is stored in one file.
%       useChannel: an array of integers specifying which channel(s) to be
%           returned.
%
%   Author: Xiong Xiao, Nanyang Technological University, Singapore
%   Last modified: 27 Jul 2016
%
function [all_wav] = Reader_binary(files, reader)

for i=1:length(files)
    [wav] = Reader_binary_core(files{i}, reader);
    if isfield(reader, 'precision')
        switch lower(reader.precision)
            case 'int16'
                wav = StoreWavInt16(wav);
            case 'single'
                wav = single(wav);
            case 'double'
                wav = double(wav);
        end
    end
    all_wav{i} = wav;
end
end

%%
function [data] = Reader_binary_core(file, reader)

words = ExtractWordsFromString_v2(file);

selectChannel = 0;
selectTime = 0;

filename = words{1};
if length(words)==4 || length(words)==2
    channelID = str2num(words{end});
    selectChannel = 1;
end
if length(words)>=3
    time1 = str2num(words{2});
    time2 = str2num(words{3});
    selectTime = 1;
end

big_endian = ReturnFieldWithDefaultValue(reader, 'big_endian', 0);

switch lower(reader.name)
    case 'raw'
        [data] = read08(filename, big_endian);   % you have to set fs correctly in the reader.
    otherwise
        if selectTime
            [data] = reader.reader.read(filename, [time1 time2]);
        else
            [data] = reader.reader.read(filename);
        end
end

if selectChannel
    data = data(channelID,:);
end

end


text:
% very basic generic non-Cartesian recon using BART
% assumes raw data *.mat (or *.dat) fiels are stored together with the
% *.seq Pulseq files withe the identical(!) names
%
% optionally needs mapVBVD in the path (if raw data have not yet been converted to *.mat)
% requires BART, modify or remove the line below depending on your path settings
bart_path='~/pulseq_home/bart/bart-0.6.00';
addpath([bart_path '/matlab']);
cur_dir=pwd;
cd(bart_path);
setenv('TOOLBOX_PATH', pwd);
cd(cur_dir);
%% Load the latest file from a dir
path='../IceNIH_RawSend/'; % directory to be scanned for data files
%path='/data/zte_petra/';
%path='~/Dropbox/shared/data/siemens/';
%path='~/20211025-AMR/data';

pattern='*.seq';

if path(end)~=filesep, path=[path filesep]; end

D=dir([path pattern]);
[~,I]=sort([D(:).datenum]);
seq_file_path = [path D(I(end-0)).name]; % use end-1 to reconstruct the second-last data set, etc.

%% alternatively just provide the path to the .seq file

%seq_file_path = '../interpreters/siemens/data_example/gre_example.seq'
%seq_file_path = [path '2020-11-10-073628.seq'];

%% Load sequence from the file with the same name as the raw data
fprintf(['loading `' seq_file_path '´ ...\n']);
seq = mr.Sequence();              % Create a new sequence object
seq.read(seq_file_path,'detectRFuse');

%% keep basic filename without the extension
[p,n,e] = fileparts(seq_file_path);
basic_file_path=fullfile(p,n);

%% load the raw data file
data_file_path= [basic_file_path '.mat']; % try to load a matlab file with raw data...
try
    fprintf(['loading `' data_file_path '´ ...\n']);
    data_unsorted = load(data_file_path);
    if isstruct(data_unsorted)
        fn=fieldnames(data_unsorted);
        assert(length(fn)==1); % we only expect a single variable
        data_unsorted=data_unsorted.(fn{1});
    end
catch
    data_file_path= [basic_file_path '.dat']; % now try to load a raw data file...
    fprintf(['falling back to `' data_file_path '´ ...\n']);
    twix_obj = mapVBVD(data_file_path);
    if iscell(twix_obj)
        data_unsorted = twix_obj{end}.image.unsorted();
    else
        data_unsorted = twix_obj.image.unsorted();
    end
    seqHash_twix=twix_obj.hdr.Dicom.tSequenceVariant;
    if length(seqHash_twix)==32
        fprintf(['raw data contain pulseq-file signature ' seqHash_twix '\n']);
    end
    clear twix_obj
end

%% calculate k-space trajectory
traj_recon_delay=0*1e-6;%1.75e-6; % adjust this parameter to potentially improve resolution & geometric accuracy. It can be calibrated by inverting the spiral revolution dimension and making two images match. for our Prisma and a particular trajectory we found 1.75e-6
samples_to_mask=0; % each ADC may contain damaged samples in the begining
fprintf('calculating k-space trajectory ...');

%[ktraj_adc, t_adc, ktraj, t, t_excitation, t_refocusing]=seq.calculateKspacePP('trajectory_delay', traj_recon_delay);
%[ktraj_adc, ktraj, t_excitation, t_refocusing, t_adc] = seq.calculateKspace('trajectory_delay', traj_recon_delay);
[ktraj_adc, t_adc, ktraj, t_ktraj, t_excitation] = seq.calculateKspacePP('trajectory_delay', traj_recon_delay);

% optionally transform or patch the trajectory
ktraj_adc(2,:)=-ktraj_adc(2,:); % seems to be needed for the correct orientaton
%ktraj_adc(3,:)=0;
fprintf(' done\n');

figure; plot3(ktraj_adc(1,:),ktraj_adc(2,:),ktraj_adc(3,:),'r.');
axis('equal');
drawnow;

fov=seq.getDefinition('FOV');
if isempty(fov)
    fov=[256 256 256]; % default FOV
    fprintf('WARNING: no FOV defined in the pulseq file, assuming FOV of %g x %g x %g m, matrix size will be wrong if these are incorrect\n', fov);
end

%fov=fov*1.5;
%% raw data preparation

adc_len=size(data_unsorted,1);
channels=size(data_unsorted,2);
readouts=size(data_unsorted,3);

%% average repetitions
na=numel(data_unsorted)/channels/numel(t_adc);
if (na>1)
    fprintf('averaging over %d acquisitions ...\n',na);
    data_unsorted=reshape(sum(reshape(data_unsorted,[adc_len,channels,readouts/na,na]),4),[adc_len,channels,readouts/na]);
    %data_unsorted=reshape(data_unsorted,[adc_len,channels,readouts/na,na]);
    %data_unsorted=reshape(data_unsorted(:,:,:,1),[adc_len,channels,readouts/na]);
    readouts=readouts/na;
end

%% mask data if requested
if samples_to_mask>0
    ktraj_adc=reshape(ktraj_adc,[3,adc_len,readouts]);
    ktraj_adc(:,1:samples_to_mask,:)=[];
    t_adc=reshape(t_adc,[adc_len,readouts]);
    t_adc(1:samples_to_mask,:)=[];
    data_unsorted(1:samples_to_mask,:,:)=[];
    adc_len=adc_len-samples_to_mask;
    ktraj_adc=reshape(ktraj_adc,[3,adc_len*readouts]);
    t_adc=reshape(t_adc,[1,adc_len*readouts]);
end

%% data preparation and petra/zte specific stuff
rawdata = permute(data_unsorted, [1,3,2]);

if strcmp('petra',seq.getDefinition('Name'))
    fprintf('performing special data handling for petra ...\n');
    SamplesPerShell=seq.getDefinition('SamplesPerShell');
    waspi_samples=1; % 1 16 32 tested; BART (nlinv) does not seem to like waspi (background noise goes up)
    if ~isempty(SamplesPerShell)
        SamplesSPI=sum(SamplesPerShell(2:end));
        data_mask=boolean([ones(SamplesPerShell(1),adc_len); ...
                          ones(SamplesSPI,waspi_samples), zeros(SamplesSPI,adc_len-waspi_samples)])';
    else
        data_mask=boolean(ones(adc_len,readouts));
    end
    ktraj_adc=ktraj_adc(:,data_mask(:));

    rdtmp=reshape(rawdata, [adc_len*readouts,channels]);
    rawdata=rdtmp(data_mask(:),:);
    clear rdtmp;
    % reshape to 4D (1 samples 1 channels) as BART needs it
    rawdata=reshape(rawdata,[1 size(rawdata,1) 1 channels]);

    kyz_plane=abs(ktraj_adc(1,:))<0.9; % *(1/fov(1)) -- for BART-like trajectory not needed
    figure; plot(ktraj_adc(2,kyz_plane),ktraj_adc(3,kyz_plane),'r.','MarkerSize',0.5); axis equal;title('Kyz plane');
    %figure; plot3(kspace(kxz_plane,1),kspace(kxz_plane,2),kspace(kxz_plane,3),'r.'); axis equal;
    drawnow;
elseif strcmp('ute_rs',seq.getDefinition('Name'))
    fprintf('performing special data handling for ute ...\n');
    % average every 2nd spoke because of the half-rf excitation
    ktraj_adc=reshape(ktraj_adc,[3,adc_len,readouts]);
    ktraj_adc=ktraj_adc(:,:,1:2:end-1);
    %rawdata=rawdata(:,1:2:end-1,:);
    rawdata=rawdata(:,1:2:end-1,:)+rawdata(:,2:2:end,:);
    readouts=readouts/2;
    ktraj_adc=reshape(ktraj_adc,[3,adc_len*readouts]);
    % make it 2D
    ktraj_adc(3,:)=0;
    % reshape raw data to 4D (1 samples 1 channels) as BART needs it
    rawdata=reshape(rawdata,[1 adc_len*readouts 1 channels]);
else
    % just a general reshape to 4D (1 samples 1 channels) as BART needs it
    rawdata = reshape(rawdata, [1 adc_len*readouts 1 channels]);
end

%% BART expects the trajectory in units of 1/FOV, Pulseq's trajectory is in inverse meters, need to rescale
for i=1:3, ktraj_adc(i,:)=ktraj_adc(i,:)*fov(i); end

% left-right swap (tested on TRIO)
ktraj_adc(1,:)=-ktraj_adc(1,:);

%% detect and clean-up 2D trajectory
kmax=max(abs(ktraj_adc'));
klim_almost_zero=1e-3; % why?
if any(kmax<klim_almost_zero)
    fprintf('2D trajectory detected\n');
    ktraj_adc(kmax<klim_almost_zero,:)=0;
end

%% calculate signal power in respective channels
datapwr=squeeze(sum(rawdata.*conj(rawdata),2));
[~,pwrsort]=sort(-datapwr);

%% call bart for the actual recon -- first try inverse gridding
igrid = bart('nufft -i -t -lowmem1', ktraj_adc, rawdata); % inverse interative gridding (-i) % -lowmem1 may be of help with large data
save([basic_file_path '_recon_igrid'],'igrid');
% display (needs imab for displaying multiple slices/partitions)
figure; imab(sum(abs(igrid).^2,4).^0.5); colormap gray;
title('BART recon, inverse gridding'); drawnow;

%% call bart for the actual recon -- now try iterative non-linear inversion which also uses multiple coils in a parallel imaging way
nlinv = bart('nlinv -n -d4 -m1 -i16 -t', ktraj_adc, rawdata); % there is a -R<n> key, but I did not notice any effect
%nlinv = bart('nlinv -n -d4 -m1 -i16 -t', ktraj_adc, rawdata(:,:,:,pwrsort(1:3))); % only include N most important channels (to save time and/or RAM)
save([basic_file_path '_recon_nlinv'],'nlinv');
figure; imab(abs(nlinv)/max(abs(nlinv(:)))); colormap gray;
title('BART recon, nonlinear inversion'); drawnow;

%% call bart for the pics algorithm with compressed sensing (with dummy coil sensitivities)
pics = bart('pics -S -l2 -r5 -t', ktraj_adc, rawdata, ones(size(igrid))); % inverse interative gridding (-i) % -lowmem1 may be of help with large data
%pics_l1 = bart('pics -S -l1 -i100 -s0.0002 -r0.003 -t', ktraj_adc, rawdata, ones(size(igrid))); % inverse interative gridding (-i) % -lowmem1 may be of help with large data
% tested r 0.002 0.004 and 0.01 (l1, l1s or ++, l1ss or +++)
save([basic_file_path '_recon_pics'],'pics');
% display (needs imab for displaying multiple slices/partitions)
figure; imab(sum(abs(pics).^2,4).^0.5); colormap gray;
title('BART recon, pics-l2'); drawnow;

% %%
% pics_l2 = bart('pics -S -l2 -r5 -t', ktraj_adc, rawdata, ones(size(igrid))); % inverse interative gridding (-i) % -lowmem1 may be of help with large data
% save([basic_file_path '_recon_picsL2'],'pics_l2');
% % display (needs imab for displaying multiple slices/partitions)
% figure; imab(sum(abs(pics_l2).^2,4).^0.5); colormap gray; drawnow;

%% 3D viewer
% img3d=abs(pics_l1(:,:,end:-1:1));
% %img3d=imag(pics);
% %img3d(img3d<0)=0;
% % img3d=sum(abs(igrid).^2,4).^0.5;
%
% s=size(img3d);
% [xx,yy,zz]=ndgrid(-s(1)/2:s(1)/2-1,-s(2)/2:s(2)/2-1,-s(3)/2:s(3)/2-1);
% img3d(vecnorm([xx(:),yy(:),zz(:)]')>min(s)/2)=0;
% figure; imab(permute(img3d(116:140,:,:),[2,3,1])); colormap gray;
% volumeViewer(img3d);

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function [S, neighbours] = findAcyclicNeighbours2(Y, k)

% FINDACYCLICNEIGHBOURS2 find the k nearest neighbours for each point in Y preventing cycles in the graph.
% FORMAT
% DESC returns the indices of the k nearest neighbours to each point in
% the given data matrix Y.
% ARG y : the data in which neighbours need to be found.
% ARG k : the number of neighbours that need to be found.
% RETURN ind : the indices of each points neighbours.
% RETURN D : the squared distance to each of the neighbours.
%
% COPYRIGHT : Neil D. Lawrence, 2010
%
% SEEALSO : lleOptimise, fmvuOptimise, isomapCreate

% MLTOOLS

  [neighboursInd, A] = findNeighbours(Y, k);
  N = size(Y, 1);
  W = spalloc(N, N, 2*size(Y, 1)*k);
  for i = 1:N
    for j = 1:k
      W(i, neighboursInd(i, j)) = -1;
      W(neighboursInd(i, j), i) = -1;
    end
  end
  L = W;
  jitter = 1e-6;
  L(1:N+1:end) = -sum(W)+jitter;
  %P = amd(L);
  %Y = Y(P, :);
  P = 1:size(L, 1);
  [UT, p, S] = chol(L, 'lower', 'vector');
  %UT = chol(L(P, P))';
  UT(1:N+1:end)=0;
  neighbours{1} = [];
  for i = 1:N-1
    neighbours{i} = find(UT(:, i));
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function fullErrMsg = bst_error(errMsg, errTitle, isStack)
% BST_ERROR: Process an error (Display error message, and send a bug report)
%
% Usage : bst_error()                          : Display lasterror
%         bst_error(errMsg, errTitle, isStack) : Report an error
%         bst_error(errMsg, errTitle)          : Report an error (stack is displayed)
%         bst_error(errMsg)                    : Report an error (errTitle is the name of the calling script)
%         fullErrMsg = bst_error(...)          : Do not display anything, just return the full error message

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2008-2022

global GlobalData;

% Display lasterror
if (nargin == 0)
    e = lasterror;
    % If no error available: create an empty error
    if isempty(e)
        try
            error('Brainstorm:NoError', 'NoError');
        catch
        end
        e = lasterror;
    end
    % Catch OUT OF MEMORY errors
    if strcmpi(e.identifier, 'MATLAB:nomem')
        errStack = e.stack;
        errMsg   = ['Out of memory. ' 10 10 'Please unload some datasets (or restart Matlab) an try again.'];
        errTitle = '';
        isStack  = 1;
    % Get lasterror stack and error message
    else
        % If error stack is empty, use current call stack
        if isempty(e.stack)
            % Get current stack trace
            errStack = dbstack();
            % Remove call to bst_error() from stack
            errStack(1) = [];
        else
            errStack = e.stack;
        end
        errMsg   = e.message;
        errTitle = '';
        isStack  = 1;
        % Remove HTML tags
        errMsg = str_striptag(errMsg);
        % Remove "Error using"
        errMsg = strrep(errMsg, 'Error using ==>', '');
        errMsg = strrep(errMsg, 'Error using', '');
    end
% Else : display bst_error line
else
    % Get current stack trace
    errStack = dbstack();
    % Remove call to bst_error() from stack
    errStack(1) = [];
    % Get error title if not defined
    if (nargin < 2)
        errTitle = '';
    end
    % If isStack not defined : set default to 1
    if (nargin < 3)
        isStack = 1;
    end
end

if isempty(errStack)
    return
end

% Set default default title
if isempty(errTitle)
    [tmp__, errFile] = fileparts(errStack(1).file);
    errTitle = ['Error in ' errFile '.m'];
end

% Print line number
strLine = [];
if isStack
    strLine = sprintf('Line %d: ', errStack(1).line);
end
% Print stack trace
strStack = [];
if isStack
    % For each stack level
    for i = 1:length(errStack)
        % Get filename
        [tmp__, errFile] = fileparts(errStack(i).file);
        errFile = [errFile, '.m']; %#ok<AGROW>
        % If s.name = s.file : do not display name
        if strcmpi(errFile, [errStack(i).name '.m'])
            strStack = sprintf('%s>%s at %d\n', strStack, errFile, errStack(i).line);
        else
            strStack = sprintf('%s>%s>%s at %d\n', strStack, errFile, errStack(i).name, errStack(i).line);
        end
    end
    strStack = sprintf('\n_______________________________________________\nCall stack:\n%s_______________________________________________\n', strStack);
end

% Full dialog error message
fullErrMsg = sprintf('%s%s\n%s', strLine, errMsg, strStack);
% Full console error message
consoleMsg = strrep(fullErrMsg(1:end-1), '_______________________________________________', '');
consoleMsg = strrep(consoleMsg, char([10 10]), char(10));
consoleMsg = strrep(consoleMsg, char(10), [10 '** ']);
consoleMsg = [10 '***************************************************************************' 10 ...
              '** Error: ' consoleMsg ...
              10 '***************************************************************************' 10];

% If no text output: display messages
if (nargout == 0)
    % Do not display java dialog in case of server mode
    if isempty(GlobalData) || ~isfield(GlobalData, 'Program') || ~isfield(GlobalData.Program, 'GuiLevel') || (GlobalData.Program.GuiLevel >= 0)
        % Hide progress bar
        bst_progress('stop');
        % If there are too many lines, truncate the message
        MAX_LINES = 30;
        if (nnz(fullErrMsg == 10) > MAX_LINES)
            iCut = find(fullErrMsg == 10, MAX_LINES);
            fullErrMsgDlg = [fullErrMsg(1:iCut(end)-1) 10 '...' 10 '[See command window for full error message]'];
        else
            fullErrMsgDlg = fullErrMsg;
        end
        % Display error message window
        java_dialog('errorhelp', fullErrMsgDlg, errTitle);
    end
    % Display error message in Matlab console
    disp(consoleMsg);
    % Close all the open files
    fclose all;
end

% %% ===== SENDING BUG REPORT =====
% % Send email only if error is not a known error (arguments are given to the function)
% if (nargin == 0)
%     % Send bug report (if configuration is ok...)
%     panel_bug('SendBugReport', consoleMsg);
% end

% Close all files
%fclose('all');

author: brainstorm-tools