text:
function varargout = spm_BIDS(varargin)
% Parse and query a directory structure formated according to the BIDS standard
% FORMAT BIDS = spm_BIDS(root)
% root   - directory formated according to BIDS [Default: pwd]
% BIDS   - structure containing the BIDS file layout
%
% FORMAT result = spm_BIDS(BIDS,query,...)
% BIDS   - BIDS directory name or BIDS structure (from spm_BIDS)
% query  - type of query: {'data', 'metadata', 'sessions', 'subjects',
%          'runs', 'tasks', 'runs', 'types', 'modalities'}
% result - outcome of query
%__________________________________________________________________________
%
% BIDS (Brain Imaging Data Structure): https://bids.neuroimaging.io/
%   The brain imaging data structure, a format for organizing and
%   describing outputs of neuroimaging experiments.
%   K. J. Gorgolewski et al, Scientific Data, 2016.
%__________________________________________________________________________
% Copyright (C) 2016-2018 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_BIDS.m 7642 2019-07-19 15:04:40Z guillaume $

%-Validate input arguments
%==========================================================================
if ~nargin
    root = pwd;
elseif nargin == 1
    if ischar(varargin{1})
        root = spm_select('CPath',varargin{1});
    else
        varargout = varargin(1);
        return;
    end
else
    BIDS = spm_BIDS(varargin{1});
    varargout{1} = BIDS_query(BIDS,varargin{2:end});
    return;
end

%-BIDS structure
%==========================================================================

BIDS = struct(...
    'dir',root, ...               % BIDS directory
    'description',struct([]), ... % content of dataset_description.json
    'sessions',{{}},...           % cellstr of sessions
    'scans',struct([]),...        % content of sub-<participant_label>_scans.tsv (should go within subjects)
    'sess',struct([]),...         % content of sub-participants_label>_sessions.tsv (should go within subjects)
    'participants',struct([]),... % content of participants.tsv
    'subjects',struct([]));       % structure array of subjects

%-Validation of BIDS root directory
%==========================================================================
if isempty(BIDS.dir)
    error('A BIDS directory has to be specified.');
elseif ~exist(BIDS.dir,'dir')
    error('BIDS directory does not exist.');
elseif ~exist(fullfile(BIDS.dir,'dataset_description.json'),'file')
    error('BIDS directory not valid: missing dataset_description.json.');
end

%-Dataset description
%==========================================================================
try
    BIDS.description = spm_jsonread(fullfile(BIDS.dir,'dataset_description.json'));
catch
    error('BIDS dataset description could not be read.');
end
if ~isfield(BIDS.description,'BIDSVersion') || ~isfield(BIDS.description,'Name')
    error('BIDS dataset description not valid.');
end
% See also optional README and CHANGES files

%-Optional directories
%==========================================================================
% [code/]
% [derivatives/]
% [stimuli/]
% [sourcedata/]
% [phenotype]

%-Scans key file
%==========================================================================

% sub-<participant_label>/[ses-<session_label>/]
%     sub-<participant_label>_scans.tsv

%-Participant key file
%==========================================================================
p = spm_select('FPList',BIDS.dir,'^participants\.tsv$');
if ~isempty(p)
    BIDS.participants = spm_load(p);
end
p = spm_select('FPList',BIDS.dir,'^participants\.json$');
if ~isempty(p)
    BIDS.participants.meta = spm_jsonread(p);
end

%-Sessions file
%==========================================================================

% sub-<participant_label>/[ses-<session_label>/]
%      sub-<participant_label>[_ses-<session_label>]_sessions.tsv

%-Tasks: JSON files are accessed through metadata
%==========================================================================
%t = spm_select('FPList',BIDS.dir,...
%    '^task-.*_(beh|bold|events|channels|physio|stim|meg)\.(json|tsv)$');

%-Subjects
%==========================================================================
sub = cellstr(spm_select('List',BIDS.dir,'dir','^sub-.*$'));
if isequal(sub,{''})
    error('No subjects found in BIDS directory.');
end

for su=1:numel(sub)
    sess = cellstr(spm_select('List',fullfile(BIDS.dir,sub{su}),'dir','^ses-.*$'));
    for se=1:numel(sess)
        if isempty(BIDS.subjects)
            BIDS.subjects = parse_subject(BIDS.dir, sub{su}, sess{se});
        else
            BIDS.subjects(end+1) = parse_subject(BIDS.dir, sub{su}, sess{se});
        end
    end
end

varargout = { BIDS };

%==========================================================================
%-Parse a subject's directory
%==========================================================================
function subject = parse_subject(p, subjname, sesname)

subject.name    = subjname;   % subject name ('sub-<participant_label>')
subject.path    = fullfile(p,subjname,sesname); % full path to subject directory
subject.session = sesname; % session name ('' or 'ses-<label>')
subject.anat    = struct([]); % anatomy imaging data
subject.func    = struct([]); % task imaging data
subject.fmap    = struct([]); % fieldmap data
subject.beh     = struct([]); % behavioral experiment data
subject.dwi     = struct([]); % diffusion imaging data
subject.eeg     = struct([]); % EEG data
subject.meg     = struct([]); % MEG data
subject.pet     = struct([]); % PET imaging data

%--------------------------------------------------------------------------
%-Anatomy imaging data
%--------------------------------------------------------------------------
pth = fullfile(subject.path,'anat');
if exist(pth,'dir')
    f = spm_select('List',pth,...
        sprintf('^%s.*_([a-zA-Z0-9]+){1}\\.nii(\\.gz)?$',subject.name));
    if isempty(f), f = {}; else f = cellstr(f); end
    for i=1:numel(f)

        %-Anatomy imaging data file
        %------------------------------------------------------------------
        p = parse_filename(f{i}, {'sub','ses','acq','ce','rec','fa','echo','inv','run'});
        subject.anat = [subject.anat p];

    end
end

%--------------------------------------------------------------------------
%-Task imaging data
%--------------------------------------------------------------------------
pth = fullfile(subject.path,'func');
if exist(pth,'dir')

    %-Task imaging data file
    %----------------------------------------------------------------------
    f = spm_select('List',pth,...
        sprintf('^%s.*_task-.*_bold\\.nii(\\.gz)?$',subject.name));
    if isempty(f), f = {}; else f = cellstr(f); end
    for i=1:numel(f)

        p = parse_filename(f{i}, {'sub','ses','task','acq','rec','fa','echo','inv','run','recording', 'meta'});
        subject.func = [subject.func p];
        subject.func(end).meta = struct([]); % ?

    end

    %-Task events file
    %----------------------------------------------------------------------
    % (!) TODO: events file can also be stored at higher levels (inheritance principle)
    f = spm_select('List',pth,...
        sprintf('^%s.*_task-.*_events\\.tsv$',subject.name));
    if isempty(f), f = {}; else f = cellstr(f); end
    for i=1:numel(f)

        p = parse_filename(f{i}, {'sub','ses','task','acq','rec','fa','echo','inv','run','recording', 'meta'});
        subject.func = [subject.func p];
        subject.func(end).meta = spm_load(fullfile(pth,f{i})); % ?

    end

    %-Physiological and other continuous recordings file
    %----------------------------------------------------------------------
    % (!) TODO: stim file can also be stored at higher levels (inheritance principle)
    f = spm_select('List',pth,...
        sprintf('^%s.*_task-.*_(physio|stim)\\.tsv\\.gz$',subject.name));
    % see also [_recording-<label>]
    if isempty(f), f = {}; else f = cellstr(f); end
    for i=1:numel(f)

        p = parse_filename(f{i}, {'sub','ses','task','acq','rec','fa','echo','inv','run','recording', 'meta'});
        subject.func = [subject.func p];
        subject.func(end).meta = struct([]); % ?

    end
end

%--------------------------------------------------------------------------
%-Fieldmap data
%--------------------------------------------------------------------------
pth = fullfile(subject.path,'fmap');
if exist(pth,'dir')
    f = spm_select('List',pth,...
        sprintf('^%s.*\\.nii(\\.gz)?$',subject.name));
    if isempty(f), f = {}; else f = cellstr(f); end
    j = 1;

    %-Phase difference image and at least one magnitude image
    %----------------------------------------------------------------------
    labels = regexp(f,[...
        '^sub-[a-zA-Z0-9]+' ...              % sub-<participant_label>
        '(?<ses>_ses-[a-zA-Z0-9]+)?' ...     % ses-<label>
        '(?<acq>_acq-[a-zA-Z0-9]+)?' ...     % acq-<label>
        '(?<run>_run-[a-zA-Z0-9]+)?' ...     % run-<index>
        '_phasediff\.nii(\.gz)?$'],'names'); % NIfTI file extension
    if any(~cellfun(@isempty,labels))
        idx = find(~cellfun(@isempty,labels));
        for i=1:numel(idx)
            fb = spm_file(spm_file(f{idx(i)},'basename'),'basename');
            metafile = fullfile(pth,spm_file(fb,'ext','json'));
            subject.fmap(j).type = 'phasediff';
            subject.fmap(j).filename = f{idx(i)};
            subject.fmap(j).magnitude = {...
                strrep(f{idx(i)},'_phasediff.nii','_magnitude1.nii'),...
                strrep(f{idx(i)},'_phasediff.nii','_magnitude2.nii')}; % optional
            subject.fmap(j).ses = regexprep(labels{idx(i)}.ses,'^_[a-zA-Z0-9]+-','');
            subject.fmap(j).acq = regexprep(labels{idx(i)}.acq,'^_[a-zA-Z0-9]+-','');
            subject.fmap(j).run = regexprep(labels{idx(i)}.run,'^_[a-zA-Z0-9]+-','');
            if spm_existfile(metafile)
                subject.fmap(j).meta = spm_jsonread(metafile);
            else
                % (!) TODO: file can also be stored at higher levels (inheritance principle)
                subject.fmap(j).meta = struct([]); % ?
            end
            j = j + 1;
        end
    end

    %-Two phase images and two magnitude images
    %----------------------------------------------------------------------
    labels = regexp(f,[...
        '^sub-[a-zA-Z0-9]+' ...           % sub-<participant_label>
        '(?<ses>_ses-[a-zA-Z0-9]+)?' ...  % ses-<label>
        '(?<acq>_acq-[a-zA-Z0-9]+)?' ...  % acq-<label>
        '(?<run>_run-[a-zA-Z0-9]+)?' ...  % run-<index>
        '_phase1\.nii(\.gz)?$'],'names'); % NIfTI file extension
    if any(~cellfun(@isempty,labels))
        idx = find(~cellfun(@isempty,labels));
        for i=1:numel(idx)
            fb = spm_file(spm_file(f{idx(i)},'basename'),'basename');
            metafile = fullfile(pth,spm_file(fb,'ext','json'));
            subject.fmap(j).type = 'phase12';
            subject.fmap(j).filename = {...
                f{idx(i)},...
                strrep(f{idx(i)},'_phase1.nii','_phase2.nii')};
            subject.fmap(j).magnitude = {...
                strrep(f{idx(i)},'_phase1.nii','_magnitude1.nii'),...
                strrep(f{idx(i)},'_phase1.nii','_magnitude2.nii')};
            subject.fmap(j).ses = regexprep(labels{idx(i)}.ses,'^_[a-zA-Z0-9]+-','');
            subject.fmap(j).acq = regexprep(labels{idx(i)}.acq,'^_[a-zA-Z0-9]+-','');
            subject.fmap(j).run = regexprep(labels{idx(i)}.run,'^_[a-zA-Z0-9]+-','');
            if spm_existfile(metafile)
                subject.fmap(j).meta = {...
                    spm_jsonread(metafile),...
                    spm_jsonread(strrep(metafile,'_phase1.json','_phase2.json'))};
            else
                % (!) TODO: file can also be stored at higher levels (inheritance principle)
                subject.fmap(j).meta = struct([]); % ?
            end
            j = j + 1;
        end
    end

    %-A single, real fieldmap image
    %----------------------------------------------------------------------
    labels = regexp(f,[...
        '^sub-[a-zA-Z0-9]+' ...             % sub-<participant_label>
        '(?<ses>_ses-[a-zA-Z0-9]+)?' ...    % ses-<label>
        '(?<acq>_acq-[a-zA-Z0-9]+)?' ...    % acq-<label>
        '(?<run>_run-[a-zA-Z0-9]+)?' ...    % run-<index>
        '_fieldmap\.nii(\.gz)?$'],'names'); % NIfTI file extension
    if any(~cellfun(@isempty,labels))
        idx = find(~cellfun(@isempty,labels));
        for i=1:numel(idx)
            fb = spm_file(spm_file(f{idx(i)},'basename'),'basename');
            metafile = fullfile(pth,spm_file(fb,'ext','json'));
            subject.fmap(j).type = 'fieldmap';
            subject.fmap(j).filename = f{idx(i)};
            subject.fmap(j).magnitude = strrep(f{idx(i)},'_fieldmap.nii','_magnitude.nii');
            subject.fmap(j).ses = regexprep(labels{idx(i)}.ses,'^_[a-zA-Z0-9]+-','');
            subject.fmap(j).acq = regexprep(labels{idx(i)}.acq,'^_[a-zA-Z0-9]+-','');
            subject.fmap(j).run = regexprep(labels{idx(i)}.run,'^_[a-zA-Z0-9]+-','');
            if spm_existfile(metafile)
                subject.fmap(j).meta = spm_jsonread(metafile);
            else
                % (!) TODO: file can also be stored at higher levels (inheritance principle)
                subject.fmap(j).meta = struct([]); % ?
            end
            j = j + 1;
        end
    end

    %-Multiple phase encoded directions (topup)
    %----------------------------------------------------------------------
    labels = regexp(f,[...
        '^sub-[a-zA-Z0-9]+' ...          % sub-<participant_label>
        '(?<ses>_ses-[a-zA-Z0-9]+)?' ... % ses-<label>
        '(?<acq>_acq-[a-zA-Z0-9]+)?' ... % acq-<label>
        '_dir-(?<dir>[a-zA-Z0-9]+)?' ... % dir-<index>
        '(?<run>_run-[a-zA-Z0-9]+)?' ... % run-<index>
        '_epi\.nii(\.gz)?$'],'names');   % NIfTI file extension
    if any(~cellfun(@isempty,labels))
        idx = find(~cellfun(@isempty,labels));
        for i=1:numel(idx)
            fb = spm_file(spm_file(f{idx(i)},'basename'),'basename');
            metafile = fullfile(pth,spm_file(fb,'ext','json'));
            subject.fmap(j).type = 'epi';
            subject.fmap(j).filename = f{idx(i)};
            subject.fmap(j).ses = regexprep(labels{idx(i)}.ses,'^_[a-zA-Z0-9]+-','');
            subject.fmap(j).acq = regexprep(labels{idx(i)}.acq,'^_[a-zA-Z0-9]+-','');
            subject.fmap(j).dir = labels{idx(i)}.dir;
            subject.fmap(j).run = regexprep(labels{idx(i)}.run,'^_[a-zA-Z0-9]+-','');
            if spm_existfile(metafile)
                subject.fmap(j).meta = spm_jsonread(metafile);
            else
                % (!) TODO: file can also be stored at higher levels (inheritance principle)
                subject.fmap(j).meta = struct([]); % ?
            end
            j = j + 1;
        end
    end
end

%--------------------------------------------------------------------------
%-EEG data
%--------------------------------------------------------------------------
pth = fullfile(subject.path,'eeg');
if exist(pth,'dir')

    %-EEG data file
    %----------------------------------------------------------------------
    f = spm_select('List',pth,...
        sprintf('^%s.*_task-.*_eeg\\..*[^json]$',subject.name));
    if isempty(f), f = {}; else f = cellstr(f); end
    for i=1:numel(f)

        p = parse_filename(f{i}, {'sub','ses','task','acq','run','meta'});
        subject.eeg = [subject.eeg p];
        subject.eeg(end).meta = struct([]); % ?

    end

    %-EEG events file
    %----------------------------------------------------------------------
    f = spm_select('List',pth,...
        sprintf('^%s.*_task-.*_events\\.tsv$',subject.name));
    if isempty(f), f = {}; else f = cellstr(f); end
    for i=1:numel(f)

        p = parse_filename(f{i}, {'sub','ses','task','acq','run','meta'});
        subject.eeg = [subject.eeg p];
        subject.eeg(end).meta = spm_load(fullfile(pth,f{i})); % ?

    end

    %-Channel description table
    %----------------------------------------------------------------------
    f = spm_select('List',pth,...
        sprintf('^%s.*_task-.*_channels\\.tsv$',subject.name));
    if isempty(f), f = {}; else f = cellstr(f); end
    for i=1:numel(f)

        p = parse_filename(f{i}, {'sub','ses','task','acq','run','meta'});
        subject.eeg = [subject.eeg p];
        subject.eeg(end).meta = spm_load(fullfile(pth,f{i})); % ?

    end

    %-Session-specific file
    %----------------------------------------------------------------------
    f = spm_select('List',pth,...
        sprintf('^%s(_ses-[a-zA-Z0-9]+)?.*_(electrodes\\.tsv|photo\\.jpg|coordsystem\\.json|headshape\\..*)$',subject.name));
    if isempty(f), f = {}; else f = cellstr(f); end
    for i=1:numel(f)

        p = parse_filename(f{i}, {'sub','ses','acq','meta'});
        subject.eeg = [subject.eeg p];
        subject.eeg(end).meta = struct([]); % ?

    end

end

%--------------------------------------------------------------------------
%-MEG data
%--------------------------------------------------------------------------
pth = fullfile(subject.path,'meg');
if exist(pth,'dir')

    %-MEG data file
    %----------------------------------------------------------------------
    [f,d] = spm_select('List',pth,...
        sprintf('^%s.*_task-.*_meg\\..*[^json]$',subject.name));
    if isempty(f), f = d; end
    if isempty(f), f = {}; else f = cellstr(f); end
    for i=1:numel(f)

        p = parse_filename(f{i}, {'sub','ses','task','acq','run','proc', 'meta'});
        subject.meg = [subject.meg p];
        subject.meg(end).meta = struct([]); % ?

    end

    %-MEG events file
    %----------------------------------------------------------------------
    % (!) TODO: events file can also be stored at higher levels (inheritance principle)
    f = spm_select('List',pth,...
        sprintf('^%s.*_task-.*_events\\.tsv$',subject.name));
    if isempty(f), f = {}; else f = cellstr(f); end
    for i=1:numel(f)

        p = parse_filename(f{i}, {'sub','ses','task','acq','run','proc', 'meta'});
        subject.meg = [subject.meg p];
        subject.meg(end).meta = spm_load(fullfile(pth,f{i})); % ?

    end

    %-Channels description table
    %----------------------------------------------------------------------
    % (!) TODO: channels file can also be stored at higher levels (inheritance principle)
    f = spm_select('List',pth,...
        sprintf('^%s.*_task-.*_channels\\.tsv$',subject.name));
    if isempty(f), f = {}; else f = cellstr(f); end
    for i=1:numel(f)

        p = parse_filename(f{i}, {'sub','ses','task','acq','run','proc', 'meta'});
        subject.meg = [subject.meg p];
        subject.meg(end).meta = spm_load(fullfile(pth,f{i})); % ?

    end

    %-Session-specific file
    %----------------------------------------------------------------------
    f = spm_select('List',pth,...
        sprintf('^%s(_ses-[a-zA-Z0-9]+)?.*_(photo\\.jpg|coordsystem\\.json|headshape\\..*)$',subject.name));
    if isempty(f), f = {}; else f = cellstr(f); end
    for i=1:numel(f)

        p = parse_filename(f{i}, {'sub','ses','task','acq','run','proc', 'meta'});
        subject.meg = [subject.meg p];
        subject.meg(end).meta = struct([]); % ?

    end

end

%--------------------------------------------------------------------------
%-Behavioral experiments data
%--------------------------------------------------------------------------
pth = fullfile(subject.path,'beh');
if exist(pth,'dir')
    f = spm_select('FPList',pth,...
        sprintf('^%s.*_(events\\.tsv|beh\\.json|physio\\.tsv\\.gz|stim\\.tsv\\.gz)$',subject.name));
    if isempty(f), f = {}; else f = cellstr(f); end
    for i=1:numel(f)

        %-Event timing, metadata, physiological and other continuous
        % recordings
        %------------------------------------------------------------------
        p = parse_filename(f{i}, {'sub','ses','task'});
        subject.beh = [subject.beh p];

    end
end

%--------------------------------------------------------------------------
%-Diffusion imaging data
%--------------------------------------------------------------------------
pth = fullfile(subject.path,'dwi');
if exist(pth,'dir')
    f = spm_select('FPList',pth,...
        sprintf('^%s.*_([a-zA-Z0-9]+){1}\\.nii(\\.gz)?$',subject.name));
    if isempty(f), f = {}; else f = cellstr(f); end
    for i=1:numel(f)

        %-Diffusion imaging file
        %------------------------------------------------------------------
        p = parse_filename(f{i}, {'sub','ses','acq','run', 'bval','bvec'});
        subject.dwi = [subject.dwi p];

        %-bval file
        %------------------------------------------------------------------
        % bval file can also be stored at higher levels (inheritance principle)
        bvalfile = get_metadata(f{i},'^.*%s\\.bval$');
        if isfield(bvalfile,'filename')
            subject.dwi(end).bval = spm_load(bvalfile.filename); % ?
        end

        %-bvec file
        %------------------------------------------------------------------
        % bvec file can also be stored at higher levels (inheritance principle)
        bvecfile = get_metadata(f{i},'^.*%s\\.bvec$');
        if isfield(bvalfile,'filename')
            subject.dwi(end).bvec = spm_load(bvecfile.filename); % ?
        end

    end
end

%--------------------------------------------------------------------------
%-Positron Emission Tomography imaging data
%--------------------------------------------------------------------------
pth = fullfile(subject.path,'pet');
if exist(pth,'dir')
    f = spm_select('List',pth,...
        sprintf('^%s.*_task-.*_pet\\.nii(\\.gz)?$',subject.name));
    if isempty(f), f = {}; else f = cellstr(f); end
    for i=1:numel(f)

        %-PET imaging file
        %------------------------------------------------------------------
        p = parse_filename(f{i}, {'sub','ses','task','acq','rec','run'});
        subject.pet = [subject.pet p];

    end
end

%==========================================================================
%-Perform a BIDS query
%==========================================================================
function result = BIDS_query(BIDS,query,varargin)
opts = parse_query(varargin);
switch query
%   case 'subjects'
%       result = regexprep(unique({BIDS.subjects.name}),'^[a-zA-Z0-9]+-','');
    case 'sessions'
        result = unique({BIDS.subjects.session});
        result = regexprep(result,'^[a-zA-Z0-9]+-','');
        result = unique(result);
        result(cellfun('isempty',result)) = [];
    case 'modalities'
        hasmod = arrayfun(@(y) structfun(@(x) isstruct(x) & ~isempty(x),y),...
            BIDS.subjects,'UniformOutput',false);
        hasmod = any([hasmod{:}],2);
        mods   = fieldnames(BIDS.subjects)';
        result = mods(hasmod);
    case {'subjects', 'tasks', 'runs', 'types', 'data', 'metadata'}
        %-Initialise output variable
        result = {};
        %-Filter according to subjects
        if any(ismember(opts(:,1),'sub'))
            subs = opts{ismember(opts(:,1),'sub'),2};
            opts(ismember(opts(:,1),'sub'),:) = [];
        else
            subs = unique({BIDS.subjects.name});
            subs = regexprep(subs,'^[a-zA-Z0-9]+-','');
        end
        %-Filter subjects according to participants.tsv
        if any(ismember(opts(:,1),'participants'))
            p = find(ismember(opts(:,1),'participants'));
            for i=1:numel(p)
                q = opts{p(i),2};
                fn = fieldnames(q);
                for j=1:numel(fn)
                    if iscell(BIDS.participants.(fn{j}))
                        idx = ismember(BIDS.participants.(fn{j}),q.(fn{j}));
                    else
                        if isa(q.(fn{j}),'function_handle')
                            idx = feval(q.(fn{j}),BIDS.participants.(fn{j}));
                        else
                            idx = BIDS.participants.(fn{j}) == q.(fn{j});
                        end
                    end
                    s = regexprep(BIDS.participants.participant_id(idx),'^[a-zA-Z0-9]+-','');
                    subs = intersect(subs, s);
                end
            end
            opts(p,:) = [];
        end
        %-Filter according to modality
        if any(ismember(opts(:,1),'modality'))
            mods = opts{ismember(opts(:,1),'modality'),2};
            opts(ismember(opts(:,1),'modality'),:) = [];
        else
            mods = BIDS_query(BIDS,'modalities');
        end
        %-Get optional target option for metadata query
        if strcmp(query,'metadata') && any(ismember(opts(:,1),'target'))
            target = opts{ismember(opts(:,1),'target'),2};
            opts(ismember(opts(:,1),'target'),:) = [];
            if iscellstr(target)
                target = substruct('.',target{1});
            end
        else
            target = [];
        end
        %-Perform query
        for i=1:numel(BIDS.subjects)
            if ~ismember(BIDS.subjects(i).name(5:end),subs), continue; end
            for j=1:numel(mods)
                d = BIDS.subjects(i).(mods{j});
                for k=1:numel(d)
                    sts = true;
                    for l=1:size(opts,1)
                        if ~isfield(d(k),opts{l,1}) || ~ismember(d(k).(opts{l,1}),opts{l,2})
                            sts = false;
                        end
                    end
                    switch query
                        case 'subjects'
                            if sts
                                result{end+1} = BIDS.subjects(i).name;
                            end
                        case 'data'
                            if sts && isfield(d(k),'filename')
                                result{end+1} = fullfile(BIDS.subjects(i).path,mods{j},d(k).filename);
                            end
                        case 'metadata'
                            if sts && isfield(d(k),'filename')
                                f = fullfile(BIDS.subjects(i).path,mods{j},d(k).filename);
                                result{end+1} = get_metadata(f);
                                if ~isempty(target)
                                    try
                                        result{end} = subsref(result{end},target);
                                    catch
                                        warning('Non-existent field for metadata.');
                                        result{end} = [];
                                    end
                                end
                            end
%                             if sts && isfield(d(k),'meta')
%                                 result{end+1} = d(k).meta;
%                             end
                        case 'runs'
                            if sts && isfield(d(k),'run')
                                result{end+1} = d(k).run;
                            end
                        case 'tasks'
                            if sts && isfield(d(k),'task')
                                result{end+1} = d(k).task;
                            end
                        case 'types'
                            if sts && isfield(d(k),'type')
                                result{end+1} = d(k).type;
                            end
                    end
                end
            end
        end
        %-Postprocessing output variable
        switch query
            case 'subjects'
                result = unique(result);
                result = regexprep(result,'^[a-zA-Z0-9]+-','');
            case 'data'
                result = result';
            case 'metadata'
                if numel(result) == 1
                    result = result{1};
                end
            case {'tasks','runs','types'}
                result = unique(result);
                result(cellfun('isempty',result)) = [];
        end
    otherwise
        error('Unable to perform BIDS query.');
end

%==========================================================================
%-Parse BIDS query
%==========================================================================
function query = parse_query(query)
if numel(query) == 1 && isstruct(query{1})
    query = [fieldnames(query{1}), struct2cell(query{1})];
else
    if mod(numel(query),2)
        error('Invalid input syntax.');
    end
    query = reshape(query,2,[])';
end
for i=1:size(query,1)
    if ischar(query{i,2})
        query{i,2} = cellstr(query{i,2});
    end
    for j=1:numel(query{i,2})
        if iscellstr(query{i,2})
            query{i,2}{j} = regexprep(query{i,2}{j},sprintf('^%s-',query{i,1}),'');
        end
    end
end

%==========================================================================
%-Parse filename
%==========================================================================
function p = parse_filename(filename,fields)
filename = spm_file(filename,'filename');
[parts, dummy] = regexp(filename,'(?:_)+','split','match');
p.filename = filename;
[p.type, p.ext] = strtok(parts{end},'.');
for i=1:numel(parts)-1
    [d, dummy] = regexp(parts{i},'(?:\-)+','split','match');
    p.(d{1}) = d{2};
end
if nargin == 2
    for i=1:numel(fields)
        if ~isfield(p,fields{i})
            p.(fields{i}) = '';
        end
    end
    try
        p = orderfields(p,['filename','ext','type',fields]);
    catch
        warning('Ignoring file "%s" not matching template.',filename);
        p = struct([]);
    end
end

%==========================================================================
%-Get metadata
%==========================================================================
function meta = get_metadata(filename, pattern)
if nargin == 1, pattern = '^.*%s\\.json$'; end
pth = fileparts(filename);
p = parse_filename(filename);

meta = struct();

if isfield(p,'ses') && ~isempty(p.ses)
    N = 4; % there is a session level in the hierarchy
else
    N = 3;
end

for n=1:N
    metafile = spm_select('FPList',pth, sprintf(pattern,p.type));
    if isempty(metafile), metafile = {}; else metafile = cellstr(metafile); end
    for i=1:numel(metafile)
        p2 = parse_filename(metafile{i});
        fn = setdiff(fieldnames(p2),{'filename','ext','type'});
        ismeta = true;
        for j=1:numel(fn)
            if ~isfield(p,fn{j}) || ~strcmp(p.(fn{j}),p2.(fn{j}))
                ismeta = false;
                break;
            end
        end
        if ismeta
            if strcmp(p2.ext,'.json')
                meta = update_metadata(meta,spm_jsonread(metafile{i}));
            else
                meta.filename = metafile{i};
            end
        end
    end
    pth = fullfile(pth,'..');
end

%==========================================================================
%-Inheritance principle
%==========================================================================
function s1 = update_metadata(s1,s2)
fn = fieldnames(s2);
for i=1:numel(fn)
    if ~isfield(s1,fn{i})
        s1.(fn{i}) = s2.(fn{i});
    end
end

author: spm


text:
clear all;
close all;
clc;

display('Nth (3) order Reconstruction');
N = 3;                  % Nth order nonuniform sampling
TQ = 1;Fs = 1/TQ;                 % Nyquist Period

T = [2*TQ 3*TQ 6*TQ];      % Decimation Periods
K = 0.5*lcm(2*T(1), 2*T(2)); K = 0.5*lcm(2*K, 2*T(3))/TQ;
capT = K*TQ; M = capT./T;
capM = lcm(M(1), M(2)); capM = lcm(capM, M(3));
excess = ceil((K-1)/capM);
maxf = K/(excess*capM+1);
TQ1 = maxf*TQ;
K1 = K/maxf;

ML = 400; % number of slices
w_c = 0.85;
NS = 100;  % Number of Sinusoids

LF = capM*2*K1*(1:10)+1;  %359,159,239          % min length of LF should be capM*2*K1

std = 1e-1;
taus = [0 1.1+std*randn 2.2+std*randn]*TQ;
tausI = sort([taus(1) taus(2) taus(3) T(1)+taus(1) T(2)+taus(2) 2*T(1)+taus(1)]);
% display(tausI);

Frq = rand(1,NS)*w_c/2;
Amp = rand(1,NS)/(sqrt(NS)*2);
Phi = rand(1,NS)*2*pi;

input = zeros(1,ML*K);
inputN = zeros(1,ML*K1);
for k = 1:NS
  input = input + Amp(k)*sin(2*pi*Frq(k)*(0:ML*K-1)*TQ+Phi(k));
  inputN = inputN + Amp(k)*sin(2*pi*Frq(k)*(0:ML*K1-1)*TQ1+Phi(k));
end;

tauI = zeros(K,ML);
for p = 1:K
    tauI(p,:) = tausI(p)+(0:ML-1)*capT;
end;

x11 = zeros(K,ML);
for k = 1:NS
    x11 = x11 + Amp(k)*sin(2*pi*Frq(k)*tauI+Phi(k));
end;

tauPr = zeros(K,ML);
for p = 1:K
    tauPr(p,:) = -tausI(p)+(0:ML-1)*capT;
end;
x1p = zeros(K,ML);
for k = 1:NS
    x1p = x1p + Amp(k)*sin(2*pi*Frq(k)*tauPr+Phi(k));
end;

timeP = zeros(size(LF));
timeE = zeros(size(LF));
timeI = zeros(size(LF));
timeV = zeros(size(LF));
timePr = zeros(size(LF));
timeJ = zeros(size(LF));

MC_runs = 1;
for rrr = 1:length(LF)
display(rrr);
n = -(LF(rrr)-1)/2:1:(LF(rrr)-1)/2;

for tt = 1:MC_runs

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Filterbank Reconstruction of Bandlimited Signals from Nonuniform and
% Generalized Samples
% Authors: Y C Eldar and A V Oppenheim
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y = zeros(N,ML*K1);
for p = 1:N
    tau = taus(p)+(0:ML*M(p)-1)*T(p);
    x1 = zeros(1,ML*M(p));
    for k = 1:NS
        x1 = x1 + Amp(k)*sin(2*pi*Frq(k)*tau+Phi(k));
    end;
    LFE = M(p)*84*rrr+1;      % length of LF should be Multiple of LCM{M(p)}*2*K
    nE = -(LFE-1)/2:1:(LFE-1)/2;

    tic
    y1 = upsample(x1,K1);
    h = sinc((nE/K1)-(taus(p)/T(p))).*kaiser_mine1(LFE,18,-K1*(taus(p)/T(p)));

%%%%%%%%%%%%%%%%%%%%%%%%%%% Implementation 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%
%     G = M(setdiff((1:N),p));
%     F = taus(setdiff((1:N),p));
%     temp1 = G(1)-G(2);
%     temp2 = G(1)+G(2);
%     bb = zeros(2*(K-M(p))+1,M(p));
%     if temp1~=0
%         for l = 0:M(p)-1
%             c = 0.5*cos(pi*(l*T(p)*temp1+G(2)*F(2)-G(1)*F(1))/capT);
%             s = -0.5*sin(pi*(l*T(p)*temp1+G(2)*F(2)-G(1)*F(1))/capT);
%             bb(1+temp2-temp1,l+1) = 0.5*(c+1i*s);
%             bb(1+temp2+temp1,l+1) = conj(bb(1+temp2-temp1,l+1));
%         end;
%     else
%         bb(1+temp2,1:M(p)) = 0.5*cos(pi*G(1)*(F(2)-F(1))/capT);
%     end;
%     for l = 0:M(p)-1
%         c = -0.5*cos(pi*(l*T(p)*temp2-G(2)*F(2)-G(1)*F(1))/capT);
%         s = 0.5*sin(pi*(l*T(p)*temp2-G(2)*F(2)-G(1)*F(1))/capT);
%         bb(1+temp2-temp2,l+1) = 0.5*(c+1i*s);
%         bb(1+temp2+temp2,l+1) = conj(bb(1+temp2-temp2,l+1));
%     end;
%     aaa = ones(1,M(p));
%     bbvl = zeros(M(p),LFE);
%     for l = 1:M(p)
%         for q = 1:N
%             if q ~= p
%                 aaa(l) = aaa(l)*sin(pi*(taus(p)-taus(q)+(l-1)*T(p))/T(q));
%             end;
%         end;
%         bbv = zeros(2*(K-M(p))+1,LFE);
%         for v = -(K-M(p)):K-M(p)
%             bbv(K-M(p)+1+v,:) = bb(K-M(p)+1+v,l)*exp(1i*(pi/(K1*M(p)))*v*nE);
%         end;
%         bbvl(l,:) = sum(bbv,1)/aaa(l);
%     end;
%     bbb = zeros(M(p),LFE);
%     bbn = zeros(M(p),LFE);
%     for m = 1:M(p)
%         for l = 1:M(p)
%             bbb(l,:) = bbvl(l,:)*exp(1i*(2*pi/M(p))*(m-1)*(l-1));
%         end
%         bbb = sum(bbb,1);
%         bbn(m,:) = bbb.*exp(1i*(2*pi/M(p))*(m-1)*nE);
%     end;
%     bbn = sum(bbn,1);
%     bbn = bbn.*h/M(p);

%%%%%%%%%%%%%%%%%%%%%%%%%%% Implementation 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%
%     G = M(setdiff((1:N),p));
%     F = taus(setdiff((1:N),p));
%     temp1 = G(1)-G(2);
%     temp2 = G(1)+G(2);
%     bb = zeros(2*(K-M(p))+1,M(p));
%     if temp1~=0
%         for l = 0:M(p)-1
%             c = 0.5*cos(pi*(l*T(p)*temp1+G(2)*F(2)-G(1)*F(1))/capT);
%             s = -0.5*sin(pi*(l*T(p)*temp1+G(2)*F(2)-G(1)*F(1))/capT);
%             bb(1+temp2-temp1,l+1) = 0.5*(c+1i*s);
%             bb(1+temp2+temp1,l+1) = conj(bb(1+temp2-temp1,l+1));
%         end;
%     else
%         bb(1+temp2,1:M(p)) = 0.5*cos(pi*G(1)*(F(2)-F(1))/capT);
%     end;
%     for l = 0:M(p)-1
%         c = -0.5*cos(pi*(l*T(p)*temp2-G(2)*F(2)-G(1)*F(1))/capT);
%         s = 0.5*sin(pi*(l*T(p)*temp2-G(2)*F(2)-G(1)*F(1))/capT);
%         bb(1+temp2-temp2,l+1) = 0.5*(c+1i*s);
%         bb(1+temp2+temp2,l+1) = conj(bb(1+temp2-temp2,l+1));
%     end;
%     aaa = ones(1,M(p));
%     for l = 1:M(p)
%         for q = 1:N
%             if q ~= p
%                 aaa(l) = aaa(l)/sin(pi*(taus(p)-taus(q)+(l-1)*T(p))/T(q));
%             end;
%         end;
%     end;
%     AE = diag(aaa);
%     BE = bb.';
%
%     dim = K-M(p); w = -dim:1:dim;
%     FE = exp(1i*(pi/(K1*M(p))).*kron(nE,w'));
%
%     E1E = exp(1i*(2*pi/M(p)).*kron((0:M(p)-1),(0:M(p)-1)'))/M(p);
%
%     E2E = exp(1i*(2*pi/M(p)).*kron(nE,(0:M(p)-1)'));
%
%     temp = E1E*AE*BE*FE;
%     bbn = h.*sum(E2E.*temp,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%% Implementation 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%
    aaa = ones(1,M(p));
    bb = ones(M(p),LFE);
    for l = 1:M(p)
        for q = 1:N
            if q ~= p
                aaa(l) = aaa(l)*sin(pi*(taus(p)-taus(q)+(l-1)*T(p))/T(q));
                bb(l,:) = bb(l,:).*sin(pi*((nE*TQ1/M(p))-taus(q)+(l-1)*T(p))/T(q));
            end;
        end;
        bb(l,:) = bb(l,:)/aaa(l);
    end;

    bbb = zeros(M(p),LFE);
    bbn = zeros(M(p),LFE);
    for m = 1:M(p)
        for l = 1:M(p)
            bbb(l,:) = bb(l,:)*exp(1i*(2*pi/M(p))*(m-1)*(l-1));
        end
        bbb = sum(bbb,1);
        bbn(m,:) = bbb.*exp(1i*(2*pi/M(p))*(m-1)*nE);
    end;
    bbn = sum(bbn,1);
    bbn = bbn.*h/M(p);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    y1 = conv(y1,bbn);
    delay = (length(h)-1)/2;
    y(p,:) = y1(1+delay:M(p):end-delay);
    timeE(rrr) = timeE(rrr)+toc;
end;
y = real(sum(y,1));
x = inputN;
% y = y(160:end);
% x = x(160:end);
serE = 20*log10(norm(x,2)/norm(y-x,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A realization of Digital Filter Banks for Reconstruction of Uniformly
% sampled signals from nonuniform samples
% Authors: Itami, Watanabe, Nishihara
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
a = zeros(1,K);
for p = 1:K
    a(p) = 1;
    for q = 1:K
        if q ~= p
                a(p) = a(p)/sin(pi*(tausI(p)-tausI(q))/capT);
        end;
    end;
end;
% c = sin(pi*(tausI(0+1))/capT);
% s = cos(pi*(tausI(0+1))/capT);
% b(1,1) = 0.5*(c+1i*s);
% c = sin(pi*(tausI(1+1))/capT);
% s = cos(pi*(tausI(1+1))/capT);
% b(1,2) = 0.5*(c+1i*s);
% b(3,:) = conj(b(1,:));
% b(2,:) = 0;

b = zeros(2*K-1,K);
c = 0.25*(sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
s = 0.25*(cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
b(5,1) = 0.5*(c+1i*s);

c = 0.25*(sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(2+1)-tausI(0+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(0+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(0+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(0+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(0+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
s = 0.25*(cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(2+1)-tausI(0+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(0+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(0+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(0+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(0+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
b(5,2) = 0.5*(c+1i*s);

c = 0.25*(sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(0+1)-tausI(1+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(0+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(0+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(0+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(0+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(0+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(0+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
s = 0.25*(cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(0+1)-tausI(1+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(0+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(0+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(0+1)-tausI(1+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(0+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(0+1)-tausI(1+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(0+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
b(5,3) = 0.5*(c+1i*s);

c = 0.25*(sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)*cos(pi*(tausI(4+1)-tausI(0+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(0+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(0+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(0+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(0+1)+tausI(4+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(0+1)+tausI(4+1))/capT));
s = 0.25*(cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)*cos(pi*(tausI(4+1)-tausI(0+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(0+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(0+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(0+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(0+1)+tausI(4+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(0+1)+tausI(4+1))/capT));
b(5,4) = 0.5*(c+1i*s);

c = 0.25*(sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)*cos(pi*(tausI(0+1)-tausI(3+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(0+1)-tausI(3+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(0+1)-tausI(3+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(3+1)+tausI(0+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(0+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(0+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(0+1))/capT));
s = 0.25*(cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)*cos(pi*(tausI(0+1)-tausI(3+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(0+1)-tausI(3+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(0+1)-tausI(3+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(3+1)+tausI(0+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(0+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(0+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(0+1))/capT));
b(5,5) = 0.5*(c+1i*s);

c = 0.25*(sin(pi*(-tausI(0+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-0.5*sin(pi*(-tausI(0+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-0.5*cos(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-0.5*sin(pi*(-tausI(0+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(0+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-0.5*cos(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
s = 0.25*(cos(pi*(-tausI(0+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)+0.5*cos(pi*(-tausI(0+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-0.5*sin(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)+0.5*cos(pi*(-tausI(0+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*cos(pi*(-tausI(0+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-0.5*sin(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*cos(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
b(5,6) = 0.5*(c+1i*s);
b(7,:) = conj(b(5,:));

c = 0.125*(-sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)+cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
s = 0.125*(-cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
b(3,1) = 0.5*(c+1i*s);

c = 0.125*(-sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(0+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)+cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(0+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(0+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(0+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(0+1)+tausI(2+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
s = 0.125*(-cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(0+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(0+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(0+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(0+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(0+1)+tausI(2+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
b(3,2) = 0.5*(c+1i*s);

c = 0.125*(-sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(0+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)+cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(0+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(0+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(0+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(0+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(0+1)-tausI(1+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(1+1)+tausI(0+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(0+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
s = 0.125*(-cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(0+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(0+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(0+1)-tausI(1+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(0+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(0+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(0+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(1+1)+tausI(0+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(0+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
b(3,3) = 0.5*(c+1i*s);

c = 0.125*(-sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(0+1))/capT)+cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(0+1))/capT)-sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(0+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(0+1)+tausI(4+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(0+1)+tausI(4+1))/capT)+cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(4+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(0+1)+tausI(4+1))/capT));
s = 0.125*(-cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(0+1))/capT)-sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(0+1))/capT)-cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(0+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(0+1)+tausI(4+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(0+1)+tausI(4+1))/capT)-sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(4+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(0+1)+tausI(4+1))/capT));
b(3,4) = 0.5*(c+1i*s);

c = 0.125*(-sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(0+1)-tausI(3+1))/capT)+cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(0+1)-tausI(3+1))/capT)-sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(3+1)+tausI(0+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(0+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(0+1))/capT)+cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(0+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(0+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(0+1))/capT));
s = 0.125*(-cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(0+1)-tausI(3+1))/capT)-sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(0+1)-tausI(3+1))/capT)-cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(3+1)+tausI(0+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(0+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(0+1))/capT)-sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(0+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(0+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)+0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(0+1))/capT));
b(3,5) = 0.5*(c+1i*s);

c = 0.125*(-sin(pi*(-tausI(0+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)+cos(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-sin(pi*(-tausI(0+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(0+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+0.5*cos(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+cos(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*cos(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)-0.5*sin(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
s = 0.125*(-cos(pi*(-tausI(0+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-sin(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(4+1)-tausI(3+1))/capT)-cos(pi*(-tausI(0+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)-0.5*cos(pi*(-tausI(0+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+0.5*sin(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-sin(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(2+1)-tausI(1+1))/capT)+0.5*sin(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)+0.5*cos(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
b(3,6) = 0.5*(c+1i*s);
b(9,:) = conj(b(3,:));

c = 0.125*(0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
s = 0.125*(0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
b(1,1) = 0.5*(c+1i*s);

c = 0.125*(0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(0+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(0+1)+tausI(2+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
s = 0.125*(0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(0+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(0+1)+tausI(2+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
b(1,2) = 0.5*(c+1i*s);

c = 0.125*(0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(0+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(0+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(1+1)+tausI(0+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(0+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
s = 0.125*(0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(0+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(0+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(1+1)+tausI(0+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(0+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
b(1,3) = 0.5*(c+1i*s);

c = 0.125*(0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(0+1)+tausI(4+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(0+1)+tausI(4+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(4+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(0+1)+tausI(4+1))/capT));
s = 0.125*(0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(0+1)+tausI(4+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(0+1)+tausI(4+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(0+1)+tausI(4+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(0+1)+tausI(4+1))/capT));
b(1,4) = 0.5*(c+1i*s);

c = 0.125*(0.5*sin(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(0+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(0+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(0+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)-0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(0+1))/capT));
s = 0.125*(0.5*cos(pi*(-tausI(5+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(0+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(0+1))/capT)+0.5*sin(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(3+1)+tausI(0+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)-0.5*cos(pi*(-tausI(5+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(0+1))/capT));
b(1,5) = 0.5*(c+1i*s);

c = 0.125*(0.5*sin(pi*(-tausI(0+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-0.5*cos(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)-0.5*cos(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)-0.5*sin(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
s = 0.125*(0.5*cos(pi*(-tausI(0+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+0.5*sin(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*cos(pi*(tausI(3+1)+tausI(4+1))/capT)+0.5*sin(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT)*cos(pi*(tausI(1+1)+tausI(2+1))/capT)-0.5*cos(pi*(-tausI(0+1))/capT)*sin(pi*(tausI(1+1)+tausI(2+1))/capT)*sin(pi*(tausI(3+1)+tausI(4+1))/capT));
b(1,6) = 0.5*(c+1i*s);
b(11,:) = conj(b(1,:));
b = -b;

y1=upsample(x11.',K1).';

k = -(K-1):1:(K-1);
m = (0:1:(2*K1-1))';
F = exp(1i*(pi/K1).*kron(m,k));

y = zeros(K,size(y1,2));
for r = 1:K
    tempI = a(r)*y1(r,:);
%     tempI2 = F*b(:,r)*a(r);
    tempI2 = F*b(:,r);
%     y2 = tempI2*y1(r,:);
    y2 = tempI2*tempI;
%     y2 = F*b(:,r)*a(r)*y1(r,:);
    h = sinc((n/K1)-tausI(r)/capT).*kaiser_mine1(LF(rrr),18,-tausI(r)/TQ1);
    for i=1:2*K1
            h1 = upsample(downsample(h,2*K1,i-1),2*K1);
%             y2(i,:) = filter(h1,1,y2(i,:));
%             y2(i,:) = filter([zeros(1,i-1),1],1,y2(i,:));%,zeros(1,2*K-i)
            temp = filter([zeros(1,i-1),1],1,h1);
            y2(i,:) = filter(temp,1,y2(i,:));
    end;
    y(r,:) = sum(y2,1);
end;
y = -real(sum(y,1));
timeI(rrr) = timeI(rrr)+toc;
delayI = (LF(rrr)-1)/2;
x=inputN(1:end-delayI);
y=y(1+delayI:end);
% y = y(160:end);
% x = x(160:end);
serI = 20*log10(norm(x,2)/norm(y-x,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reconstruction of Nonuniformly Sampled Band-Limited Signals
% Using a Differentiator-Multiplier Cascade
% Authors: Stefan Tertinek and Christian Vogel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r = tausI-TQ*(0:K-1);
rr = r(mod((0:ML*K-1),K)+1);
x1 = reshape(x11,1,K*size(x11,2));
% x1 = x11;% LF = (0:11); Fs = 1;
% Differentiator Design
% figure();
% NFFT = 2^nextpow2(length(Hd)); % Next power of 2 from length of Hd
% HD = fftshift(fft(Hd,NFFT))/length(Hd);
% f = Fs*linspace(-1,1,NFFT);
% % Plot double-sided amplitude spectrum.
% plot(f,2*abs(HD(1:NFFT)))
% title('Double-Sided Amplitude Spectrum of Hd(n)')
% xlabel('Frequency (Hz)')
% ylabel('|HD(f)|')
tic
Hd = firpm(LF(rrr)-1,[0 w_c],[0 w_c*pi],'differentiator');
delayV = (LF(rrr)-1)/2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y1 = filter(Hd,1,x1);
x1 = filter([zeros(1,delayV),1],1,x1);
r2 = filter([zeros(1,delayV),1],1,rr);
e = y1.*r2;
y1 = x1-e;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y2 = filter(Hd,1,y1);
temp = filter([zeros(1,delayV),1],1,y2);
x1 = filter([zeros(1,2*delayV),1],1,x1);
r2 = filter([zeros(1,2*delayV),1],1,r2);
e1 = temp.*r2;
y2 = filter(Hd,1,y2);
e2 = 0.5*y2.*r2.^2;
y2 = x1-e1-e2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y3 = filter(Hd,1,y2);
temp = filter([zeros(1,2*delayV),1],1,y3);
x1 = filter([zeros(1,3*delayV),1],1,x1);
r2 = filter([zeros(1,3*delayV),1],1,r2);
e1 = temp.*r2;
y3 = filter(Hd,1,y3);
temp = filter([zeros(1,delayV),1],1,y3);
e2 = 0.5*temp.*r2.^2;
y3 = filter(Hd,1,y3);
e3 = (y3.*r2.^3)/6;
y3 = x1-e1-e2-e3;
timeV(rrr) = timeV(rrr)+toc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y = real(y3(1+6*delayV:end));
x = input(1:end-6*delayV);
% y = y(160:end);
% x = x(160:end);
serV = 20*log10(norm(x,2)/norm(y-x,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reconstruction of Band-Limited Periodic Nonuniformly Sampled Signals
% Through Multirate Filter Banks
% Ryan S Prendergast, Bernard C Levy, Paul J Hurst
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
H = zeros(K,LF(rrr));
for i=1:K
%     H(i,:) = sinc(n-tausI(i)).*conv(sinc(n-tausI(i)),kaiser(LF,10).','same');
    H(i,:) = sinc(n-tausI(i)).*kaiser_mine1(LF(rrr),18,-tausI(i));
end;
H = H(:,1:end-1);
EP = reshape(H.',K,length(H(1,:))/K,K);
capE = [];
for k=1:K
    temp = [];
    for i=1:K
        temp=[temp,toeplitz([EP(i,1,k),zeros(1,(length(H(1,:))/K)-1)],[EP(i,:,k),zeros(1,(length(H(1,:))/K)-1)])];
    end
    capE = [capE;temp];
end
d = ceil(size(capE,2)/(2*K));
P = kron(eye(K),[zeros(1,d),1,zeros(1,(size(capE,2)/K)-d-1)]);
R = P/capE;
% size(capE)
% size(zeros(LF-1,2*LF-2-K))
% size(R)
% size(zeros(K,LF-1))
% size(P)
% size(zeros(K,2*LF-2-K))
R = upsample(R.',K).';
rt = size(R,2)/K;
FR = zeros(K,rt);
for j=1:K
    for i=1:K
        temp = filter([zeros(1,K-i),1],1,R(i,(j-1)*rt+1:j*rt));
        FR(j,:) = FR(j,:)+temp;
    end
end
yb = upsample(x1p.',K).';
for i = 1:K
    yb(i,:) = filter(FR(i,:),1,yb(i,:));
end;
y = sum(yb,1);
timePr(rrr) = timePr(rrr)+toc;
delayPr = (size(FR,2))/2+K-1;
x=input(1:end-delayPr);
y=y(1+delayPr:end);
% y = y(160:end);
% x = x(160:end);
serPr = 20*log10(norm(x,2)/norm(y-x,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reconstruction of Periodically Nonuniformly Sampled Bandlimited Signals
% Using Time-Varying FIR Filters
% Authors: H. Johansson and Per Lowenborg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r = r.';
x1 = reshape(x11,1,size(x11,2)*K);
tic
w_o = w_c*pi*TQ;
hJ = zeros(K,LF(rrr));
C = zeros(1,LF(rrr));
Nt = (LF(rrr)-1)/2;
for i = 1:K
    C = -2*sin(w_o*(n-r(1+(mod(i-1-n,K)))'))./(pi*(n-r(1+(mod(i-1-n,K)))'));
    C(isnan(C)==1)=-2*w_o/pi;
    C = C.';
    S = zeros(LF(rrr),LF(rrr));
    for k = 1:LF(rrr)
        S(k,:) = sin(w_o*(-Nt+k-1-r(1+(mod(i-1-(-Nt+k-1),K)))-(n-r(1+(mod(i-1-n,K)))')))./(pi*(-Nt+k-1-r(1+(mod(i-1-(-Nt+k-1),K)))-(n-r(1+(mod(i-1-n,K)))')));
    end;
    S(isnan(S)==1)=w_o/pi;
    hJ(i,:) = -0.5*S\C;
end;

y1 = zeros(K,length(x1));
for j=1:K
    y1(j,:) = filter(hJ(j,:),1,x1);
    y1(j,:) = upsample(downsample(y1(j,:),K,j-1),K)/K;
    y1(j,:) = filter([zeros(1,j-1),1],1,y1(j,:));
end;
y = K*0.25*sum(y1,1);
timeJ(rrr) = timeJ(rrr)+toc;
delayJ = (size(hJ,2)-1)/2;
y = real(y(1+delayJ:end));
x = input(1:end-delayJ);
% y = y(160:end);
% x = x(160:end);
serJ = 20*log10(norm(x,2)/norm(y-x,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reconstruction of N-th order nonuniformly sampled bandlimited signals
% using digital filter banks
% Authors: S. K. Sindhi, K. M. M. Prabhu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xp = zeros(N,ML*K1);
for p = 1:N
    tau = taus(p)+(0:ML*M(p)-1)*T(p);
    x1 = zeros(1,ML*M(p));
    for k = 1:NS
        x1 = x1 + Amp(k)*sin(2*pi*Frq(k)*tau+Phi(k));
    end;
    tic
    m = (0:1:M(p)-1)'; lemda = 0:1:M(p)-1;
    W = exp(1i*(2*pi/M(p)).*kron(m,lemda)); % m*lemda

    aaa = ones(1,M(p));
    for l = 1:M(p)
        for q = 1:N
            if q ~= p
                    aaa(l) = aaa(l)/(sin(pi*M(q)*(taus(p)-taus(q)+(l-1)*T(p))/capT));
            end;
        end;
    end;
    A = diag(aaa); %display(A);

    G = M(setdiff((1:N),p));
    F = taus(setdiff((1:N),p));
    temp1 = G(1)-G(2);
    temp2 = G(1)+G(2);
    bb = zeros(2*(K-M(p))+1,M(p));
    if temp1~=0
        for l = 0:M(p)-1
            c = 0.5*cos(pi*(l*T(p)*temp1+G(2)*F(2)-G(1)*F(1))/capT);
            s = -0.5*sin(pi*(l*T(p)*temp1+G(2)*F(2)-G(1)*F(1))/capT);
            bb(1+temp2-temp1,l+1) = 0.5*(c+1i*s);
            bb(1+temp2+temp1,l+1) = conj(bb(1+temp2-temp1,l+1));
        end;
    else
        bb(1+temp2,1:M(p)) = 0.5*cos(pi*G(1)*(F(2)-F(1))/capT);
    end;
    for l = 0:M(p)-1
        c = -0.5*cos(pi*(l*T(p)*temp2-G(2)*F(2)-G(1)*F(1))/capT);
        s = 0.5*sin(pi*(l*T(p)*temp2-G(2)*F(2)-G(1)*F(1))/capT);
        bb(1+temp2-temp2,l+1) = 0.5*(c+1i*s);
        bb(1+temp2+temp2,l+1) = conj(bb(1+temp2-temp2,l+1));
    end;
    B = bb;% display(B);

    y1 = upsample(x1,K1);
    y1 = reshape(y1,M(p),length(y1)/M(p));

    if M(p)>1
        y1(2:end,:) = flipud(y1(2:end,:));
        for i = 1:M(p)-1
            y1(i+1,:) = filter([0,1],1,y1(i+1,:));
        end;
    end;

    dim = K-M(p); w = -dim:1:dim;

    xlemda = zeros(M(p),size(y1,2));
    for lemda = 0:M(p)-1

        rP = (lemda/M(p))+(0:1:(2*K1-1))';
        Fshift = exp(1i*(pi/K1).*kron(rP,w));   % r*w
        mtemp = A*W;
        mtemp = B*mtemp;
        mtemp = Fshift*mtemp;
        Htemp = mtemp*W(:,lemda+1);

        h = sinc((n*TQ1/T(p))+(lemda/K1)-(taus(p)/T(p))).*kaiser_mine1(LF(rrr),18,(lemda/M(p))-(taus(p)/TQ1));
        h1 = zeros(2*K1, length(h)+2*K1-1);
        for i = 2:2*K1
            h1(i,:) = [filter([zeros(1,i-1),1],1,upsample(downsample(h,2*K1,i-1),2*K1)) zeros(1,2*K1)]*Htemp(i);
        end;
        h1(1,:) = upsample(downsample(h,2*K1),2*K1)*Htemp(1);
        h1 = sum(h1,1);
        xlemda(lemda+1,:) = filter(h1,1,y1(lemda+1,:));
    end;
    xp(p,:) = sum(xlemda,1)/M(p);
    timeP(rrr) = timeP(rrr)+toc;
end;
y = sum(xp,1);
delayP = (length(h)-1)/2;
y = y(1+delayP:end);
x = inputN(1:end-delayP);
% y = y(160:end);
% x = x(160:end);
serP = 20*log10(norm(x,2)/norm(y-x,2));
end;
end;
timeP = timeP/MC_runs;
timeE = timeE/MC_runs;
timeI = timeI/MC_runs;
timeV = timeV/MC_runs;
timePr = timePr/MC_runs;
timeJ = timeJ/MC_runs;

figure();hold on;
plot(LF,timeJ,'kp-','LineWidth',2);
plot(LF,timePr,'ko-','LineWidth',2);
plot(LF,timeV,'ks-','LineWidth',2);
plot(LF,timeP,'kd-','LineWidth',2);
plot(LF,timeI,'k>-','LineWidth',2);
plot(LF,timeE,'k+-','LineWidth',2);
legend('Johansson','Prendergast','Tertinek','Proposed','Itami','Eldar');
xlabel('Filter length','fontsize',14,'fontweight','b');
ylabel('Time in seconds','fontsize',14,'fontweight','b');
grid on;box on;
set(gca,'fontsize',14,'fontweight','b')

% temp = (LF.^3);%+(K*LF*ML*K);
% temp = temp/max(temp);
% timePrI = timePr(rrr)*temp;
%
% temp = (LF.^3);%+(K*(LF+1)*ML*K);
% temp = temp/max(temp);
% timeJI = timeJ(rrr)*temp;
%
% temp = (LF.^2);%+(6*LF*ML*K);
% temp = temp/max(temp);
% timeVI = timeV(rrr)*temp;
% plot(LF,timePrI,'k--','LineWidth',2);
% plot(LF,timeJI,'k--','LineWidth',2);
% plot(LF,timeVI,'k--','LineWidth',2);

figure();hold on;
plot(LF,timeJ,'kp-','LineWidth',2);
temp = (LF.^3);%+(K*(LF+1)*ML*K);
temp = temp/max(temp);
timeJI = timeJ(rrr)*temp;
plot(LF,timeJI,'k--','LineWidth',2);
legend('Johansson','Ideal curve');
xlabel('Filter length','fontsize',14,'fontweight','b');
ylabel('Time in seconds','fontsize',14,'fontweight','b');
grid on;box on;
set(gca,'fontsize',14,'fontweight','b')

figure();hold on;
plot(LF,timePr,'ko-','LineWidth',2);
temp = (LF.^3);%+(K*LF*ML*K);
temp = temp/max(temp);
timePrI = timePr(rrr)*temp;
plot(LF,timePrI,'k--','LineWidth',2);
legend('Prendergast','Ideal curve');
xlabel('Filter length','fontsize',14,'fontweight','b');
ylabel('Time in seconds','fontsize',14,'fontweight','b');
grid on;box on;
set(gca,'fontsize',14,'fontweight','b')

figure();hold on;
plot(LF,timeV,'ks-','LineWidth',2);
temp = (LF.^2);%+(6*LF*ML*K);
temp = temp/max(temp);
timeVI = timeV(rrr)*temp;
plot(LF,timeVI,'k--','LineWidth',2);
legend('Tertinek','Ideal curve');
xlabel('Filter length','fontsize',14,'fontweight','b');
ylabel('Time in seconds','fontsize',14,'fontweight','b');
grid on;box on;
set(gca,'fontsize',14,'fontweight','b')

display(serI)
display(serV)
display(serPr)
display(serJ)
display(serP)
display(serE)

% LF = 84*(1:30)+1;      % length of LF should be Multiple of LCM{M(p)}*2*K
% timeJ = [timeJ,zeros(1,length(1:length(LF)-10))];
% for rrr = 11:length(LF)
% display(rrr);
% n = -(LF(rrr)-1)/2:1:(LF(rrr)-1)/2;
% for tt = 1:MC_runs
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Reconstruction of Periodically Nonuniformly Sampled Bandlimited Signals
% % Using Time-Varying FIR Filters
% % Authors: H. Johansson and Per Lowenborg
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% r = tausI-TQ*(0:K-1);
% r = r.';
% x1 = reshape(x11,1,size(x11,2)*K);
% tic
% w_o = w_c*pi*TQ;
% hJ = zeros(K,LF(rrr));
% C = zeros(1,LF(rrr));
% Nt = (LF(rrr)-1)/2;
% for i = 1:K
%     C = -2*sin(w_o*(n-r(1+(mod(i-1-n,K)))'))./(pi*(n-r(1+(mod(i-1-n,K)))'));
%     C(isnan(C)==1)=-2*w_o/pi;
%     C = C.';
%     S = zeros(LF(rrr),LF(rrr));
%     for k = 1:LF(rrr)
%         S(k,:) = sin(w_o*(-Nt+k-1-r(1+(mod(i-1-(-Nt+k-1),K)))-(n-r(1+(mod(i-1-n,K)))')))./(pi*(-Nt+k-1-r(1+(mod(i-1-(-Nt+k-1),K)))-(n-r(1+(mod(i-1-n,K)))')));
%     end;
%     S(isnan(S)==1)=w_o/pi;
%     hJ(i,:) = -0.5*S\C;
% end;
%
% y1 = zeros(K,length(x1));
% for j=1:K
%     y1(j,:) = filter(hJ(j,:),1,x1);
%     y1(j,:) = upsample(downsample(y1(j,:),K,j-1),K)/K;
%     y1(j,:) = filter([zeros(1,j-1),1],1,y1(j,:));
% end;
% y = K*0.25*sum(y1,1);
% timeJ(rrr) = timeJ(rrr)+toc;
% delayJ = (size(hJ,2)-1)/2;
% y = real(y(1+delayJ:end));
% x = input(1:end-delayJ);
% % y = y(160:end);
% % x = x(160:end);
% serJ = 20*log10(norm(x,2)/norm(y-x,2));
% end;
% end;
% timeJ(11:length(LF)) = timeJ(11:length(LF))/MC_runs;
% figure();hold on;
% plot(LF,timeJ,'kp-','LineWidth',2);
% temp = (LF.^3);%+(K*(LF+1)*ML*K);
% temp = temp/max(temp);
% timeJI = timeJ(rrr)*temp;
% plot(LF,timeJI,'k--','LineWidth',2);
% legend('Johansson','Prendergast','Tertinek','Proposed','Itami','Eldar');
% xlabel('Filter length','fontsize',14,'fontweight','b');
% ylabel('Time in seconds','fontsize',14,'fontweight','b');
% grid on;box on;
% set(gca,'fontsize',14,'fontweight','b')

% figure();
% subplot(2,1,1);
% plot(([x' y']));
% title('input / output signals');
% xlabel('sample');
% ylabel('signal value');
% grid on;
% subplot(2,1,2);
% plot((x'-y'));
% xlabel('time (sample)');
% ylabel('error value');
% grid on;

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
if strcmp(jobName, 'save_disparity')
  p = get_paths();
  C = cropCamera(getCameraParam('color'));
  outDir = p.ft_disparity_dir;
  imList = getImageSet('all');
  parfor i = 1:length(imList),
    args{i} = {imList{i}, C, outDir, [], []};
    saveDisparity(args{i}{:});
  end
end

if strcmp(jobName, 'save_color')
  p = get_paths();
  outDir = p.ft_image_dir;
  imList = getImageSet('all');
  parfor i = 1:length(imList),
    I = getImage(imList{i}, 'images');
    imwrite(I, fullfile_ext(outDir, imList{i}, 'png'), 'png');
  end
end

if strcmp(jobName, 'save_hha')
  p = get_paths();
  C = cropCamera(getCameraParam('color'));
  outDir = p.ft_hha_dir;
  imList = getImageSet('all');
  parfor i = 1:length(imList),
    args{i} = {imList{i}, C, outDir, [], []};
    saveHHA(args{i}{:});
  end
end

if strcmp(jobName, 'write_window_file')
  p = get_paths();
  c = benchmarkPaths();
  NYU_ROOT_DIR = c.dataDir;
  REGIONDIR = fullfile(p.output_dir, 'regions', 'release');
  SALT = 'release';
  MAX_BOXES = 2000;
  task = 'task-detection-with-cabinet';

  trainset = 'train';
  imdb = imdb_from_nyud2(c.dataDir, trainset, task, REGIONDIR, SALT, MAX_BOXES);
  imdb.roidb_func = @roidb_from_nyud2;
  window_file_train_color = rcnn_make_window_file(imdb, fullfile(p.wf_dir), 'color', p.ft_image_dir, 'png');
  window_file_train_hha = rcnn_make_window_file(imdb, fullfile(p.wf_dir), 'hha', p.ft_hha_dir, 'png');
  window_file_train_disparity = rcnn_make_window_file(imdb, fullfile(p.wf_dir), 'disparity', p.ft_disparity_dir, 'png');

  trainset = 'val';
  imdb = imdb_from_nyud2(NYU_ROOT_DIR, trainset, task, REGIONDIR, SALT, MAX_BOXES);
  imdb.roidb_func = @roidb_from_nyud2;
  window_file_val_color = rcnn_make_window_file(imdb, fullfile(p.wf_dir), 'color', p.ft_image_dir, 'png');
  window_file_val_hha = rcnn_make_window_file(imdb, fullfile(p.wf_dir), 'hha', p.ft_hha_dir, 'png');
  window_file_val_disparity = rcnn_make_window_file(imdb, fullfile(p.wf_dir), 'disparity', p.ft_disparity_dir, 'png');

  file_name = {'nyud2_finetune_solver.color', 'nyud2_finetune_solver.hha', 'nyud2_finetune_solver.disparity', ...
  'nyud2_finetune_train.color', 'nyud2_finetune_train.hha', 'nyud2_finetune_train.disparity', ...
  'nyud2_finetune_val.color', 'nyud2_finetune_val.hha', 'nyud2_finetune_val.disparity', ...
  'imagenet_deploy.prototxt'};

  args = {'PROTODIR', p.proto_dir, ...
    'WINDOW_FILE_TRAIN_COLOR', window_file_train_color, ...
    'WINDOW_FILE_TRAIN_HHA', window_file_train_hha, ...
    'WINDOW_FILE_TRAIN_DISPARITY', window_file_train_disparity, ...
    'WINDOW_FILE_VAL_COLOR', window_file_val_color, ...
    'WINDOW_FILE_VAL_HHA', window_file_val_hha, ...
    'WINDOW_FILE_VAL_DISPARITY', window_file_val_disparity, ...
    'MEAN_FILE_COLOR', fullfile_ext(p.mean_file_color, 'proto'), ...
    'MEAN_FILE_HHA', fullfile_ext(p.mean_file_hha, 'proto'), ...
    'MEAN_FILE_DISPARITY', fullfile_ext(p.mean_file_disparity, 'proto'), ...
    'N_CLASSES', sprintf('%d', imdb.num_classes+1), ...
    'SNAPSHOTDIR', p.snapshot_dir,...
  };
  make_protofiles(p.proto_dir, file_name, args);

  % Generate the caffe command
  fprintf('GLOG_logstdtoerr=1 caffe/build/tools/caffe.bin train -gpu 0 -iterations 30000 -solver %s -weights %s -model %s/imagenet_deploy.prototxt 2>&1 | tee %s\n\n\n', fullfile(p.proto_dir, 'nyud2_finetune_solver.color'), p.caffe_net, p.proto_dir, fullfile(p.proto_dir, 'color.log'));
  fprintf('GLOG_logstdtoerr=1 caffe/build/tools/caffe.bin train -gpu 0 -iterations 30000 -solver %s -weights %s -model %s/imagenet_deploy.prototxt 2>&1 | tee %s\n\n\n', fullfile(p.proto_dir, 'nyud2_finetune_solver.hha'), p.caffe_net, p.proto_dir, fullfile(p.proto_dir, 'hha.log'));

  % Control experiment for finetuning on the disparity only..
  % fprintf('GLOG_logstdtoerr=1 caffe/build/tools/caffe.bin train -gpu 1 -iterations 30000 -solver %s -weights %s -model %s/imagenet_deploy.prototxt 2>&1 | tee %s\n\n\n', fullfile(p.proto_dir, 'nyud2_finetune_solver.disparity'), p.caffe_net, p.proto_dir, fullfile(p.proto_dir, 'disparity.log'));
end

if strcmp(jobName, 'hha_cache_features')
  p = get_paths();
  c = benchmarkPaths();
  NYU_ROOT_DIR = c.dataDir;
  REGIONDIR = fullfile(p.output_dir, 'regions', 'release');
  SALT = 'release';
  MAX_BOXES = 2000;
  task = 'task-detection-with-cabinet';

  imset = 'all';
  imdb = imdb_from_nyud2(NYU_ROOT_DIR, imset, task, REGIONDIR, SALT, MAX_BOXES);
  imdb.roidb_func = @roidb_from_nyud2;

  p = get_paths();
  image_dir = fullfile(p.ft_hha_dir);
  image_ext = 'png';
  snapshot = 30000;
  net_file = fullfile_ext(p.snapshot_dir, sprintf('nyud2_finetune_hha_iter_%d', snapshot), 'caffemodel');
  feat_cache_dir = p.cnnF_cache_dir;
  net_def_file = fullfile('nyud2_finetuning', 'imagenet_hha_256_fc6.prototxt');
  mean_file = fullfile_ext(p.mean_file_hha, 'mat');
  cache_name = sprintf('hha_%d', snapshot);
  args = {};
  st = 1; sp = 1; e = 0; gpu_id = 1;
  args{1} = {'start', st, 'step', sp, 'end', e, ...
    'image_dir', image_dir, 'image_ext', image_ext, ...
    'feat_cache_dir', feat_cache_dir, ...
    'net_def_file', net_def_file, 'net_file', net_file, 'mean_file', mean_file, ...
    'cache_name', cache_name, 'gpu_id', gpu_id};
  rcnn_cache_features(imdb, args{1}{:});
end

if strcmp(jobName, 'disparity_cache_features')
  p = get_paths();
  c = benchmarkPaths();
  NYU_ROOT_DIR = c.dataDir;
  REGIONDIR = fullfile(p.output_dir, 'regions', 'release');
  SALT = 'release';
  MAX_BOXES = 2000;
  task = 'task-detection-with-cabinet';

  imset = 'all';
  imdb = imdb_from_nyud2(NYU_ROOT_DIR, imset, task, REGIONDIR, SALT, MAX_BOXES);
  imdb.roidb_func = @roidb_from_nyud2;

  p = get_paths();
  image_dir = fullfile(p.ft_disparity_dir);
  image_ext = 'png';
  snapshot = 30000;
  net_file = fullfile_ext(p.snapshot_dir, sprintf('nyud2_finetune_disparity_iter_%d', snapshot), 'caffemodel');
  feat_cache_dir = p.cnnF_cache_dir;
  net_def_file = fullfile('nyud2_finetuning', 'imagenet_disparity_256_fc6.prototxt');
  mean_file = fullfile_ext(p.mean_file_disparity, 'mat');
  cache_name = sprintf('disparity_%d', snapshot);
  args = {};
  st = 1; sp = 1; e = 0; gpu_id = 0;
  args{1} = {'start', st, 'step', sp, 'end', e, ...
    'image_dir', image_dir, 'image_ext', image_ext, ...
    'feat_cache_dir', feat_cache_dir, ...
    'net_def_file', net_def_file, 'net_file', net_file, 'mean_file', mean_file, ...
    'cache_name', cache_name, 'gpu_id', gpu_id};
  rcnn_cache_features(imdb, args{1}{:});
end

if strcmp(jobName, 'color_cache_features')
  p = get_paths();
  c = benchmarkPaths();
  NYU_ROOT_DIR = c.dataDir;
  REGIONDIR = fullfile(p.output_dir, 'regions', 'release');
  SALT = 'release';
  MAX_BOXES = 2000;
  task = 'task-detection-with-cabinet';

  imset = 'all';
  imdb = imdb_from_nyud2(NYU_ROOT_DIR, imset, task, REGIONDIR, SALT, MAX_BOXES);
  imdb.roidb_func = @roidb_from_nyud2;

  p = get_paths();
  image_dir = fullfile(p.ft_image_dir);
  image_ext = 'png';
  snapshot = 30000;
  net_file = fullfile_ext(p.snapshot_dir, sprintf('nyud2_finetune_color_iter_%d', snapshot), 'caffemodel');
  feat_cache_dir = p.cnnF_cache_dir;
  net_def_file = fullfile('nyud2_finetuning', 'imagenet_color_256_fc6.prototxt');
  mean_file = fullfile_ext(p.mean_file_color, 'mat');
  cache_name = sprintf('rgb_%d', snapshot);
  args = {};
  st = 1; sp = 1; e = 0; gpu_id = 0;
  args{1} = {'start', st, 'step', sp, 'end', e, ...
    'image_dir', image_dir, 'image_ext', image_ext, ...
    'feat_cache_dir', feat_cache_dir, ...
    'net_def_file', net_def_file, 'net_file', net_file, 'mean_file', mean_file, ...
    'cache_name', cache_name, 'gpu_id', gpu_id};
  rcnn_cache_features(imdb, args{1}{:});
end

if strcmp(jobName, 'rcnn_train'),
  % The final model
  res = rcnn_all('task-detection', 'rgb_hha', 'trainval', 'test');

  % Control experiments
  % res = rcnn_all('task-detection', 'rgb_hha', 'train', 'val');
  % res = rcnn_all('task-detection', 'hha', 'train', 'val');
  % res = rcnn_all('task-detection', 'rgb', 'train', 'val');
  % res = rcnn_all('task-detection', 'disparity', 'train', 'val');
end

if strcmp(jobName, 'proto_to_mat_mean')
  p = get_paths();
  image_mean = caffe('read_mean', fullfile_ext(p.mean_file_color, 'proto'));
  image_mean = permute(image_mean, [2 1 3]);
  save(fullfile_ext(p.mean_file_color, 'mat'), 'image_mean');

  image_mean = caffe('read_mean', fullfile_ext(p.mean_file_hha, 'proto'));
  image_mean = permute(image_mean, [2 1 3]);
  save(fullfile_ext(p.mean_file_hha, 'mat'), 'image_mean');

  image_mean = caffe('read_mean', fullfile_ext(p.mean_file_disparity, 'proto'));
  image_mean = permute(image_mean, [2 1 3]);
  save(fullfile_ext(p.mean_file_disparity, 'mat'), 'image_mean');
end

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function siLen = calcSupInfLength(structNum,planC)
% function siLen = calcSupInfLength(structNum,planC)
%
% This function computes the length in SI direction for the passed structNum.
%
% APA, 9/28/2017

if ~exist('planC','var')
    global planC
end

indexS = planC{end};

% Initialize min/max
maxZ = -inf;
minZ = inf;

% Loo over all transverse slices for z-values
for slc = 1:length(planC{indexS.structures}(structNum).contour)
    for seg = 1:length(planC{indexS.structures}(structNum).contour(slc).segments)
        if ~isempty(planC{indexS.structures}(structNum).contour(slc).segments(seg).points)
            zSegV = planC{indexS.structures}(structNum).contour(slc).segments(seg).points(:,3);
            minZ = min(minZ,zSegV(1));
            maxZ = max(maxZ,zSegV(1));
        end
    end
end

siLen = maxZ - minZ;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function varargout=choosenk(varargin)
%
% For calling details please see v_choosenk.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_choosenk.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: choosenk.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_choosenk(varargin{:});
else
    v_choosenk(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
%% Copyright (C) 2015 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @deftypefn  {Function File} {} print_usage ()
%% Wrapper for print_usage on Matlab.
%%
%% Matlab has no print_usage function.  Preparing
%% the Matlab package should rename this from
%% "my_print_usage.m" to "print_usage.m".
%%
%% @seealso{error}
%% @end deftypefn

function print_usage ()

  % if we are on Octave
  if (exist ('OCTAVE_VERSION', 'builtin'))
    %evalin ('caller', 'print_usage ();' );
    %return
  end

  H = dbstack;
  name = H(2).name;

  msgid = sprintf('%s:invalidCall', name);
  errmsg = sprintf('Invalid call to "%s": see documentation', name);

  %error('Invalid call to "%s": see documentation', name)
  ME = MException(msgid, errmsg);
  ME.throwAsCaller();

end

repo: octsympy


text:
RELATIONS.COOCCUR = 0;
RELATIONS.SUPPORT1 = -1; % row supports column
RELATIONS.SUPPORT2 = -2; % column supports row
RELATIONS.ATTACH = -3;
RELATIONS.SIDEBYSIDE = -4;
RELATIONS.SYMSUR_START = 10;

NODETYPES.COOCCUR = 2;
NODETYPES.SURROUND = 3;
NODETYPES.SUPPORT = 1;
NODETYPES.ROOM = 4;
NODETYPES.WALL = 5;

%% file locations
data_folder = ['..', filesep, '0-data'];
labelset_folder = [data_folder, filesep, 'labelsets'];

data_filename = [data_folder, filesep,'1-graphs',filesep,'bedroom_graph.mat']; % the dataset generated by step1.
labelsetname = [labelset_folder, filesep, 'bedroomlabelset.mat']; % coarse-grained labels
output_folder = [data_folder, filesep, '2-hierarchies'];
output_filename = [output_folder, filesep, 'bedroom_hierarchies.mat'];

repo: GRAINS


text:
function [ fileList, patientList ] = matRad_scanDicomImportFolder( patDir )
% matRad function to scan a folder for dicom data
%
% call
%   [ fileList, patientList ] = matRad_scanDicomImportFolder( patDir )
%
% input
%   patDir:         folder to be scanned
%
% output
%   fileList:       matlab struct with a list of dicom files including meta
%                   infomation (type, series number etc.)
%   patientList:    list of patients with dicom data in the folder
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

matRad_cfg = MatRad_Config.instance();

%% print current status of the import script
matRad_cfg.dispInfo('Dose series matched to the different plans are displayed and could be selected.\n');
matRad_cfg.dispInfo('Rechecking of correct matching procedure is recommended.\n');

global warnDlgDICOMtagShown;
warnDlgDICOMtagShown = false;

%% get all files in search directory

% dicom import needs image processing toolbox -> check if available
if ~license('checkout','image_toolbox')
    matRad_cfg.dispError('image processing toolbox and/or corresponding licence not available');
end

fileList = matRad_listAllFiles(patDir);

if ~isempty(fileList)
    %% check for dicom files and differentiate patients, types, and series
    numOfFiles = numel(fileList(:,1));
    h = waitbar(0,'Please wait...');
    %h.WindowStyle = 'Modal';
    steps = numOfFiles;
    for i = numOfFiles:-1:1
        waitbar((numOfFiles+1-i) / steps)
        try % try to get DicomInfo
            if verLessThan('matlab','9')
                info = dicominfo(fileList{i});
            else
                info = dicominfo(fileList{i},'UseDictionaryVR',true);
            end
        catch
            fileList(i,:) = [];
            matRad_progress(numOfFiles+1-i, numOfFiles);
            continue;
        end
        try
            fileList{i,2} = info.Modality;
        catch
            fileList{i,2} = NaN;
        end

        fileList = parseDicomTag(fileList,info,'PatientID',i,3);

        switch fileList{i,2}
            case 'CT'

               fileList = parseDicomTag(fileList,info,'SeriesInstanceUID',i,4);

            case {'RTPLAN','RTDOSE','RTSTRUCT'}

               fileList = parseDicomTag(fileList,info,'SOPInstanceUID',i,4);

           otherwise

              fileList = parseDicomTag(fileList,info,'SeriesInstanceUID',i,4);

        end

        fileList = parseDicomTag(fileList,info,'SeriesNumber',i,5,@seriesnum2str); %We want to make sure the series number is stored as string
        fileList = parseDicomTag(fileList,info,'FamilyName',i,6);
        fileList = parseDicomTag(fileList,info,'GivenName',i,7);
        fileList = parseDicomTag(fileList,info,'PatientBirthDate',i,8);

        try
            if strcmp(info.Modality,'CT')
                fileList{i,9} = num2str(info.PixelSpacing(1));
            else
                fileList{i,9} = NaN;
            end
        catch
            fileList{i,9} = NaN;
        end
        try
            if strcmp(info.Modality,'CT')
                fileList{i,10} = num2str(info.PixelSpacing(2));
            else
                fileList{i,10} = NaN;
            end
        catch
            fileList{i,10} = NaN;
        end
        try
            if strcmp(info.Modality,'CT')
                %usually the Attribute should be SliceThickness, but it
                %seems like some data uses "SpacingBetweenSlices" instead.
                if isfield(info,'SliceThickness') && ~isempty(info.SliceThickness)
                    fileList{i,11} = num2str(info.SliceThickness);
                elseif isfield(info,'SpacingBetweenSlices')
                    fileList{i,11} = num2str(info.SpacingBetweenSlices);
                else
                    matRad_cfg.dispError('Could not identify spacing between slices since neither ''SliceThickness'' nor ''SpacingBetweenSlices'' are specified');
                end
            else
                fileList{i,11} = NaN;
            end
        catch
            fileList{i,11} = NaN;
        end
        try
            if strcmp(info.Modality,'RTDOSE')
                dosetext_helper = strcat('Instance','_', num2str(info.InstanceNumber),'_', ...
                    info.DoseSummationType, '_', info.DoseType);
                fileList{i,12} = dosetext_helper;
            else
                fileList{i,12} = NaN;
            end
        catch
            fileList{i,12} = NaN;
        end
        % writing corresponding dose dist.
        try
            if strcmp(fileList{i,2},'RTPLAN')
                corrDose = [];
                numDose = length(fieldnames(info.ReferencedDoseSequence));
                for j = 1:numDose
                    fieldName = strcat('Item_',num2str(j));
                    corrDose{j} = info.ReferencedDoseSequence.(fieldName).ReferencedSOPInstanceUID;
                end
                fileList{i,13} = corrDose;
            else
                fileList{i,13} = {'NaN'};
            end

        catch
            fileList{i,13} = {'NaN'};
        end
        matRad_progress(numOfFiles+1-i, numOfFiles);

    end
    close(h)

    if ~isempty(fileList)
        patientList = unique(fileList(:,3));

        if isempty(patientList)
            msgbox('No patient found with DICOM CT _and_ RT structure set in patient directory!', 'Error','error');
        end
    else
        msgbox('No DICOM files found in patient directory!', 'Error','error');
        %h.WindowStyle = 'Modal';
        %error('No DICOM files found in patient directory');
    end
else
    msgbox('Search folder empty!', 'Error','error');

end

clear warnDlgDICOMtagShown;

end

function fileList = parseDicomTag(fileList,info,tag,row,column,parsefcn)

global warnDlgDICOMtagShown;

defaultPlaceHolder = '001';

if nargin < 6
    parsefcn = @(x) x;
end

try
   if isfield(info,tag)
      if ~isempty(info.(tag))
         fileList{row,column} = parsefcn(info.(tag));
      else
         fileList{row,column} = defaultPlaceHolder;
      end
   else
      fileList{row,column} = defaultPlaceHolder;
   end
catch
   fileList{row,column} = NaN;
end

if ~warnDlgDICOMtagShown && strcmp(fileList{row,column},defaultPlaceHolder) && (column == 3 || column == 4)

   dlgTitle    = 'Dicom Tag import';
   dlgQuestion = ['matRad_scanDicomImportFolder: Could not parse dicom tag: ' tag '. Using placeholder ' defaultPlaceHolder ' instead. Please check imported data carefully! Do you want to continue?'];
   answer      = questdlg(dlgQuestion,dlgTitle,'Yes','No', 'Yes');

   warnDlgDICOMtagShown = true;

   switch answer
      case 'No'
         matRad_cfg.dispError('Inconsistency in DICOM tags')
   end
end

end

function value = seriesnum2str(value)
    if isnumeric(value)
        value = num2str(value);
    end
end

repo: matRad


%% custom input
addpath('./CustomLayers/','./utils/')
% clear all % 如果更换模型，需重置静态函数(影响性能)，否则可以不用清理
cfg_file = 'cfg/yolov4-tiny.cfg';
weight_file = 'weights/yolov4-tiny.weights';
throushold = 0.1;
NMS = 0.4;

%% import all classes
fid = fopen('coco.names','r');
names = textscan(fid, '%s', 'Delimiter',{'   '});
fclose(fid);classesNames = categorical(names{1});
RGB = randi(255,length(classesNames),3);

%% 摄像头视频流识别
cap = webcam();
player = vision.DeployableVideoPlayer();
image = cap.snapshot();
step(player, image);

while player.isOpen()
    image = cap.snapshot();
    t1 = tic;
    outFeatures = yolov3v4Predict(cfg_file,weight_file,image);% M*(5+nc) ,为[x,y,w,h,Pobj,p1,p2,...,pn]
    fprintf('预测耗时：%.2f 秒\n',toc(t1));% yolov4大概0.4秒，yolov3大概0.2秒，yolov3-tiny大概0.06秒,yolov4-tiny大概0.07秒,yolov3-tiny-prn大概0.06秒

    %% 阈值过滤+NMS处理
    scores = outFeatures(:,5);
    outFeatures = outFeatures(scores>throushold,:);

    allBBoxes = outFeatures(:,1:4);
    allScores = outFeatures(:,5);
    [maxScores,indxs] = max(outFeatures(:,6:end),[],2);
    allScores = allScores.*maxScores;
    allLabels = classesNames(indxs);

    % NMS非极大值抑制
    if ~isempty(allBBoxes)
        [bboxes,scores,labels] = selectStrongestBboxMulticlass(allBBoxes,allScores,allLabels,...
            'RatioType','Min','OverlapThreshold',NMS);
        annotations = string(labels) + ": " + string(scores);
        [~,ids] = ismember(labels,classesNames);
        colors = RGB(ids,:);
        image = insertObjectAnnotation(image,...
            'rectangle',bboxes,cellstr(annotations),...
            'Color',colors,...
            'LineWidth',3);
    end
    step(player,image);
end
release(player);


text:
function printXferOrbitToTextbox(hXfrOrbitText, xferOrbitTextHeader, xfrOrbit, form, paddLen)
%printXferOrbitToTextbox Summary of this function goes here
%   Detailed explanation goes here

    hRule = getHRule();

    xferOrbitText = xferOrbitTextHeader;
    xferOrbitText{end+1} = hRule;
    xferOrbitText{end+1} = [paddStr('Semi-major Axis = ',paddLen), num2str(xfrOrbit(1), form), ' km'];
    xferOrbitText{end+1} = [paddStr('Eccentricity = ', paddLen), num2str(xfrOrbit(2))];
    xferOrbitText{end+1} = [paddStr('Inclination = ',paddLen), num2str(rad2deg(AngleZero2Pi(xfrOrbit(3))), form), ' deg'];
    xferOrbitText{end+1} = [paddStr('Right Ascension of AN = ',paddLen), num2str(rad2deg(AngleZero2Pi(xfrOrbit(4))), form), ' deg'];
    xferOrbitText{end+1} = [paddStr('Argument of Periapse = ',paddLen), num2str(rad2deg(AngleZero2Pi(xfrOrbit(5))), form), ' deg'];
    xferOrbitText{end+1} = '---------------------';
    xferOrbitText{end+1} = [paddStr('Departure True Anomaly = ',paddLen), num2str(rad2deg(AngleZero2Pi(xfrOrbit(6))), form), ' deg'];
    xferOrbitText{end+1} = [paddStr('Arrival True Anomaly = ',paddLen), num2str(rad2deg(AngleZero2Pi(xfrOrbit(7))), form), ' deg'];
    set(hXfrOrbitText,'String',xferOrbitText);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot