text:
function [x,fval,exitflag,info] = opti_csdp(f,A,b,lb,ub,sdcone,x0,opts)
%OPTI_CSDP Solve a SDP using CSDP [Primal Form]
%
%   min f'*x      subject to:     A*x <= b
%    x                            lb <= x <= ub
%                                 X = F1*x1 + F2*x2 + ... + Fn*xn - F0
%                                 X >= 0 [positive semidefinite]
%
%
%
%   x = opti_csdp(f,A,b,lb,ub) solves a LP where f is the objective
%   vector, A,b are the linear inequality constraints, and lb,ub are the
%   bounds.
%
%   x = opti_csdp(f,...,ub,sdcone) uses the cell array sdcone to pass
%   semidefinite constraints to the solver. Each cell is a sparse matrix of
%   the form [F0(:) F1(:) F2(:)...] (or in alternative notation, [C A0 A1 A2..])
%   where each F matrix has been converted to a column vector and concatenated.
%
%   x = opti_csdp(f,...,x0,opts) uses opts to pass optiset options to the
%   solver.
%
%   [x,fval,exitflag,info] = opti_csdp(...) returns the objective value at
%   the solution, together with the solver exitflag, and an information
%   structure.
%
%   THIS IS A WRAPPER FOR CSDP USING THE MEX INTERFACE
%   See supplied Eclipse Public License

%   Copyright (C) 2013 Jonathan Currie (I2C2)

t = tic;

% Handle missing arguments
if nargin < 8, opts = optiset; end
if nargin < 7, x0 = []; end
if nargin < 6, sdcone = []; end
if nargin < 5, ub = []; end
if nargin < 4, lb = []; end
if nargin < 3, error('You must supply at least 3 arguments to opti_csdp'); end

warn = strcmpi(opts.warnings,'all');

%Check sparsity
if(~isempty(A) && ~issparse(A))
    if(warn), optiwarn('opti:sparse','The A matrix should be sparse, correcting: [sparse(A)]'); end
    A = sparse(A);
end
if(~isempty(sdcone))
    ok = 1;
    if(iscell(sdcone))
        for i = 1:length(sdcone)
            if(~issparse(sdcone{i}))
                sdcone{i} = sparse(sdcone{i});
                ok = 0;
            end
        end
    elseif(~issparse(sdcone))
        sdcone = sparse(sdcone);
        ok = 0;
    end
    if(~ok && warn)
        optiwarn('opti:sparse','The Semidefinite Constraints should be sparse, correcting: [sparse(sdcone)]');
    end
end

%Addin csdp settings if specified
if(isfield(opts,'solverOpts') && ~isempty(opts.solverOpts))
    copts = csdpset(opts.solverOpts);
else
    copts = [];
end
%Add OPTI Options
copts.maxtime = opts.maxtime;
copts.maxiter = opts.maxiter;
copts.display = dispLevel(opts.display);

%MEX contains error checking
[x,fvals,exitflag,stats,X] = csdp(f, A, b, lb, ub, sdcone, x0, copts);

%Assign primal fval
fval = fvals.pval;

%Assign Outputs
info.Iterations = stats.iter;
info.Time = toc(t);
info.Algorithm = 'CSDP: Predictor-Corrector Primal-Dual SDP Solver';

switch(exitflag)
    case 0
        info.Status = 'CSDP Converged';
        exitflag = 1;
    case 1
        info.Status = 'Primal Infeasible';
        exitflag = -1;
    case 2
        info.Status = 'Dual Infeasible';
        exitflag = -1;
    case 3
        info.Status = 'Partial Success (Full Accuracy Not Achieved)';
        exitflag = 3;
    case 4
        info.Status = 'Exceeded Maximum Iterations';
        exitflag = 0;
    case 5
        info.Status = 'Stuck at Edge of Primal Feasibility';
        exitflag = -2;
    case 6
        info.Status = 'Stuck at Edge of Dual Infeasibility';
        exitflag = -2;
    case 7
        info.Status = 'Lack of Progress';
        exitflag = -3;
    case 8
        info.Status = 'X, Z or O was singular';
        exitflag = -3;
    case 9
        info.Status = 'Detected NaN or Inf';
        exitflag = -3;
    case 10
        info.Status = 'General easy_sdp() failure';
        exitflag = -3;
    case 11
        info.Status = 'Failed check on C (check symmetry)';
        exitflag = -3;
    case 12
        info.Status = 'Failed constraint check';
        exitflag = -3;
    case -27
        info.Status = 'Exceeded Maximum Time';
        exitflag = 0;
    case -50
        info.Status = 'User Exited';
        exitflag = -5;
    otherwise
        info.Status = [];
end

info.DualObjective = fvals.dval;
info.X = X;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function script_faster_rcnn_demo()
close all;
clc;
clear mex;
clear is_valid_handle; % to clear init_key
run(fullfile(fileparts(fileparts(mfilename('fullpath'))), 'startup'));
%% -------------------- CONFIG --------------------
opts.caffe_version          = 'caffe_faster_rcnn';
opts.gpu_id                 = auto_select_gpu;
active_caffe_mex(opts.gpu_id, opts.caffe_version);

opts.per_nms_topN           = 6000;
opts.nms_overlap_thres      = 0.7;
opts.after_nms_topN         = 300;
opts.use_gpu                = true;

opts.test_scales            = 600;

%% -------------------- INIT_MODEL --------------------
model_dir                   = fullfile(pwd, 'output', 'faster_rcnn_final', 'faster_rcnn_VOC0712_vgg_16layers'); %% VGG-16
%model_dir                   = fullfile(pwd, 'output', 'faster_rcnn_final', 'faster_rcnn_VOC0712_ZF'); %% ZF
proposal_detection_model    = load_proposal_detection_model(model_dir);

proposal_detection_model.conf_proposal.test_scales = opts.test_scales;
proposal_detection_model.conf_detection.test_scales = opts.test_scales;
if opts.use_gpu
    proposal_detection_model.conf_proposal.image_means = gpuArray(proposal_detection_model.conf_proposal.image_means);
    proposal_detection_model.conf_detection.image_means = gpuArray(proposal_detection_model.conf_detection.image_means);
end

% caffe.init_log(fullfile(pwd, 'caffe_log'));
% proposal net
rpn_net = caffe.Net(proposal_detection_model.proposal_net_def, 'test');
rpn_net.copy_from(proposal_detection_model.proposal_net);
% fast rcnn net
fast_rcnn_net = caffe.Net(proposal_detection_model.detection_net_def, 'test');
fast_rcnn_net.copy_from(proposal_detection_model.detection_net);

% set gpu/cpu
if opts.use_gpu
    caffe.set_mode_gpu();
else
    caffe.set_mode_cpu();
end

%% -------------------- WARM UP --------------------
% the first run will be slower; use an empty image to warm up

for j = 1:2 % we warm up 2 times
    im = uint8(ones(375, 500, 3)*128);
    if opts.use_gpu
        im = gpuArray(im);
    end
    [boxes, scores]             = proposal_im_detect(proposal_detection_model.conf_proposal, rpn_net, im);
    aboxes                      = boxes_filter([boxes, scores], opts.per_nms_topN, opts.nms_overlap_thres, opts.after_nms_topN, opts.use_gpu);
    if proposal_detection_model.is_share_feature
        [boxes, scores]             = fast_rcnn_conv_feat_detect(proposal_detection_model.conf_detection, fast_rcnn_net, im, ...
            rpn_net.blobs(proposal_detection_model.last_shared_output_blob_name), ...
            aboxes(:, 1:4), opts.after_nms_topN);
    else
        [boxes, scores]             = fast_rcnn_im_detect(proposal_detection_model.conf_detection, fast_rcnn_net, im, ...
            aboxes(:, 1:4), opts.after_nms_topN);
    end
end

%% -------------------- TESTING --------------------
im_names = {'001763.jpg', '004545.jpg', '000542.jpg', '000456.jpg', '001150.jpg'};
% these images can be downloaded with fetch_faster_rcnn_final_model.m

running_time = [];
for j = 1:length(im_names)

    im = imread(fullfile(pwd, im_names{j}));

    if opts.use_gpu
        im = gpuArray(im);
    end

    % test proposal
    th = tic();
    [boxes, scores]             = proposal_im_detect(proposal_detection_model.conf_proposal, rpn_net, im);
    t_proposal = toc(th);
    th = tic();
    aboxes                      = boxes_filter([boxes, scores], opts.per_nms_topN, opts.nms_overlap_thres, opts.after_nms_topN, opts.use_gpu);
    t_nms = toc(th);

    % test detection
    th = tic();
    if proposal_detection_model.is_share_feature
        [boxes, scores]             = fast_rcnn_conv_feat_detect(proposal_detection_model.conf_detection, fast_rcnn_net, im, ...
            rpn_net.blobs(proposal_detection_model.last_shared_output_blob_name), ...
            aboxes(:, 1:4), opts.after_nms_topN);
    else
        [boxes, scores]             = fast_rcnn_im_detect(proposal_detection_model.conf_detection, fast_rcnn_net, im, ...
            aboxes(:, 1:4), opts.after_nms_topN);
    end
    t_detection = toc(th);

    fprintf('%s (%dx%d): time %.3fs (resize+conv+proposal: %.3fs, nms+regionwise: %.3fs)\n', im_names{j}, ...
        size(im, 2), size(im, 1), t_proposal + t_nms + t_detection, t_proposal, t_nms+t_detection);
    running_time(end+1) = t_proposal + t_nms + t_detection;

    % visualize
    classes = proposal_detection_model.classes;
    boxes_cell = cell(length(classes), 1);
    thres = 0.6;
    for i = 1:length(boxes_cell)
        boxes_cell{i} = [boxes(:, (1+(i-1)*4):(i*4)), scores(:, i)];
        boxes_cell{i} = boxes_cell{i}(nms(boxes_cell{i}, 0.3), :);

        I = boxes_cell{i}(:, 5) >= thres;
        boxes_cell{i} = boxes_cell{i}(I, :);
    end
    figure(j);
    showboxes(im, boxes_cell, classes, 'voc');
    pause(0.1);
end
fprintf('mean time: %.3fs\n', mean(running_time));

caffe.reset_all();
clear mex;

end

function proposal_detection_model = load_proposal_detection_model(model_dir)
    ld                          = load(fullfile(model_dir, 'model'));
    proposal_detection_model    = ld.proposal_detection_model;
    clear ld;

    proposal_detection_model.proposal_net_def ...
                                = fullfile(model_dir, proposal_detection_model.proposal_net_def);
    proposal_detection_model.proposal_net ...
                                = fullfile(model_dir, proposal_detection_model.proposal_net);
    proposal_detection_model.detection_net_def ...
                                = fullfile(model_dir, proposal_detection_model.detection_net_def);
    proposal_detection_model.detection_net ...
                                = fullfile(model_dir, proposal_detection_model.detection_net);

end

function aboxes = boxes_filter(aboxes, per_nms_topN, nms_overlap_thres, after_nms_topN, use_gpu)
    % to speed up nms
    if per_nms_topN > 0
        aboxes = aboxes(1:min(length(aboxes), per_nms_topN), :);
    end
    % do nms
    if nms_overlap_thres > 0 && nms_overlap_thres < 1
        aboxes = aboxes(nms(aboxes, nms_overlap_thres, use_gpu), :);
    end
    if after_nms_topN > 0
        aboxes = aboxes(1:min(length(aboxes), after_nms_topN), :);
    end
end

author: ShaoqingRen


text:
function signal = toneBurst(sample_freq, signal_freq, num_cycles, varargin)
%TONEBURST Create an enveloped single frequency tone burst.
%
% DESCRIPTION:
%       toneBurst creates an enveloped single frequency tone burst for use
%       in ultrasound simulations. If an array is given for the optional
%       input 'SignalOffset', a matrix of tone bursts is created where each
%       row corresponds to a tone burst for particular value of the
%       'SignalOffset'. If a value for the optional input 'SignalLength' is
%       given, the tone burst/s are zero padded to this length (in
%       samples).
%
% USAGE:
%       signal = toneBurst(sample_freq, signal_freq, num_cycles)
%       signal = toneBurst(sample_freq, signal_freq, num_cycles, ...)
%
% INPUTS:
%       sample_freq     - sampling frequency [Hz]
%       signal_freq     - frequency of the tone burst signal [Hz]
%       num_cycles      - number of sinusoidal oscillations
%
% OPTIONAL INPUTS:
%       Optional 'string', value pairs that may be used to modify the
%       default computational settings.
%
%       'Envelope'      - envelope used to taper the tone burst, can be set to
%                         either 'Gaussian' (default) or 'Rectangular'
%       'Plot'          - Boolean controlling whether the created tone
%                         burst is plotted
%       'SignalLength'  - signal length in number of samples, if longer
%                         than the tone burst length, the signal is
%                         appended with zeros.
%       'SignalOffset'  - signal offset before the tone burst starts in
%                         number of samples
%
% OUTPUTS:
%       signal      - created tone burst
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 4th December 2009
%       last update - 21st December 2011
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also gaussian

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

% set usage defaults
num_req_input_variables = 3;
envelope = 'Gaussian';
signal_length = [];
signal_offset = 0;
plot_signal = false;

% replace with user defined values if provided
if nargin < num_req_input_variables
    error('Incorrect number of inputs');
elseif ~isempty(varargin)
    for input_index = 1:2:length(varargin)
        switch varargin{input_index}
            case 'Envelope'
                envelope = varargin{input_index + 1};
            case 'Plot'
                plot_signal = varargin{input_index + 1};
            case 'SignalOffset'
                signal_offset = varargin{input_index + 1};
                signal_offset = round(signal_offset);       % force integer
            case 'SignalLength'
                signal_length = varargin{input_index + 1};
                signal_length = round(signal_length);       % force integer
            otherwise
                error('Unknown optional input');
        end
    end
end

% calculate the temporal spacing
dt = 1/sample_freq;

% create the tone burst
tone_length = num_cycles/(signal_freq);
tone_t = 0:dt:tone_length;
tone_burst = sin(2*pi*signal_freq*tone_t);
tone_index = round(signal_offset);

% create the envelope
switch envelope
    case 'Gaussian'
        x_lim = 3;
        window_x = -x_lim:2*x_lim/(length(tone_burst)-1):x_lim;
        window = gaussian(window_x, 1, 0, 1);
    case 'Rectangular'
        window = ones(size(tone_burst));
    otherwise
        error(['Unknown envelope ' envelope]);
end

% apply the envelope
tone_burst = tone_burst.*window;

% force the ends to be zero by applying a second window
tone_burst = tone_burst.*(getWin(length(tone_burst), 'Tukey', 'Param', 0.05).');

% % calculate the expected FWHM in the frequency domain
% t_var = tone_length/(2*x_lim);
% w_var = 1/(4*pi^2*t_var);
% fw = 2 * sqrt(2 * log(2) * w_var)

% create the signal with the offset tone burst
if isempty(signal_length)
    signal = zeros(length(tone_index), max(signal_offset(:)) + length(tone_burst));
else
    signal = zeros(length(tone_index), signal_length);
end
for offset = 1:length(tone_index)
    signal(offset, tone_index(offset) + 1:tone_index(offset) + length(tone_burst)) = tone_burst;
end

% plot the signal if required
if plot_signal

    % compute suitable axis scaling
    time_axis = (0:length(signal)-1)*dt;
    [t_sc, scale, prefix] = scaleSI(max(time_axis(:)));

    % create figure
    figure;
    if numel(signal_offset) == 1
        plot(time_axis*scale, signal, 'k-');
    else
        plot(time_axis*scale, signal + repmat(2*max(signal(:))*(1:length(signal_offset)).', 1, length(time_axis)), 'k-');
    end
    xlabel(['Time [' prefix 's]']);
    ylabel('Signal Amplitude');
    axis tight;

end

repo: Aero-matlab


text:
% DEMCLASSIFICATIONSVARGPLVM2
% COPYRIGHT Andreas C. Damianou, 2012
% SEEALSO: demClassificationSvargplvm, demClassificationSvargplvm3,  demClassificationSvargplvm4
% VARGPLVM

% Like demClassificationSvargplvm but we train
% with RANDOM subsets of the training data and in the end we take the mean.

% Fix seeds

if ~exist('onlyTest'), onlyTest = false; end
if ~exist('doPredictions'),    doPredictions = false; end
if ~exist('testDataPerClass'), testDataPerClass = -1; end %%% All data
if ~exist('printPlots'), printPlots = true; end

dataPerClass = -1; %%%% ALL data

addpath(genpath('../../bc-vargplvm/matlab'));
addpath(genpath('../../vargplvm/matlab'));

% This script initialises the options structure 'globalOpt'.
svargplvm_init;

if onlyTest
    load(['demOilSvargplvm' num2str(experimentNo)]);
    globalOpt = model.globalOpt;
end

globalOpt.dataSetNames = {'oilData', 'oilLabels'};
globalOpt.dataSetName = 'oil';
globalOpt.dataToKeep = -1; %%

if exist('timeStamps')
    [globalOpt, YtrAll] = bc_loadData(globalOpt, timeStamps);
else
    [globalOpt, YtrAll] = bc_loadData(globalOpt);
end
Yall{1} = YtrAll.Y;
lblsTrain = YtrAll.lbls;
labelsTrain = transformLabels(lblsTrain);
% ---- We want labels with -1 and 1, not 0 and 1 (anti-correlation)
YtrAll.lbls(YtrAll.lbls == 0) = -1;
Yall{2} = YtrAll.lbls;

globalOptTest = globalOpt;
globalOptTest.dataPerClass = testDataPerClass;
globalOptTest.dataSetName = 'oilTest';
if exist('timeStamps')
    [globalOptTest, YtsAll] = bc_loadData(globalOptTest, timeStamps);
else
    [globalOptTest, YtsAll] = bc_loadData(globalOptTest);
end
lblsTest = YtsAll.lbls;
labelsTest = transformLabels(lblsTest);
Yts{1} = YtsAll.Y;
YtsAll.lbls(YtsAll.lbls == 0) = -1;
Yts{2} = YtsAll.lbls;

if isfield(globalOpt, 'normaliseData') && globalOpt.normaliseData
    Yall{1} = utils_normaliseData(Yall{1});
    Yts{1} = utils_normaliseData(Yts{1});
end

%%%%%%%%%%
indTemp = randperm(size(Yall{1},1));
indTrial = indTemp(1:dataToKeep);
Yall{1} = Yall{1}(indTrial,:);
Yall{2} = Yall{2}(indTrial,:);
lblsTrain = lblsTrain(indTrial,:);
labelsTrain = labelsTrain(indTrial);
%%%%%%%%%%%

if onlyTest
    model = svargplvmRestorePrunedModel(model, Yall);
end

numberOfDatasets = length(Yall);

% globalOpt.baseKern = {'rbfardjit', 'linard2'};

%globalOpt.baseKern = {'rbfardjit', 'rbfardjit'};

globalOpt.indPoints = min(globalOpt.indPoints, size(Yall{1},1));

%-- Load datasets
for i=1:numberOfDatasets
    Y = Yall{i};
    dims{i} = size(Y,2);
    N{i} = size(Y,1);
    indTr = globalOpt.indTr;
    if indTr == -1
        indTr = 1:N{i};
    end
    if ~exist('Yts')
        indTs = setdiff(1:size(Y,1), indTr);
        Yts{i} = Y(indTs,:);
    end
    Ytr{i} = Y(indTr,:);

    t{i} = linspace(0, 2*pi, size(Y, 1)+1)'; t{i} = t{i}(1:end-1, 1);
    timeStampsTraining{i} = t{i}(indTr,1); %timeStampsTest = t(indTs,1);
end

for i=2:numberOfDatasets
    if N{i} ~= N{i-1}
        error('The number of observations in each dataset must be the same!');
    end
end

%%
if ~onlyTest
    % Free up some memory
    clear('Y')

    options = svargplvmOptions(Ytr, globalOpt, labelsTrain);

    if ~isempty(globalOpt.dynamicsConstrainType)
        for i=1:numberOfDatasets
            % Set up dynamcis (i.e. back-constraints) model
            optionsDyn{i}.type = 'vargpTime';
            optionsDyn{i}.inverseWidth=30;
            %   optionsDyn.vardistCovars = vardistCovarsMult;
            optionsDyn{i}.initX = globalOpt.initX;
            optionsDyn{i}.constrainType = globalOpt.dynamicsConstrainType;

            if exist('timeStampsTraining')
                optionsDyn{i}.t = timeStampsTraining;
            end
            if exist('labelsTrain') && ~isempty(labelsTrain)
                optionsDyn{i}.labels = labelsTrain;
            end
        end
    else
        optionsDyn= [];
    end

    model = svargplvmModelCreate(Ytr, globalOpt, options, optionsDyn);
    if exist('diaryFile')
        model.diaryFile = diaryFile;
    end

    model.globalOpt = globalOpt;
    model.options = options;

    %-- Define what level of parallelism to use (w.r.t submodels or/and w.r.t
    % datapoints).
    %{
fprintf('# Parallel computations w.r.t the submodels!\n');
model.parallel = 1;
model = svargplvmPropagateField(model,'parallel', 1);
%
fprintf('# Parallel computations w.r.t the datapoints!\n');
model.vardist.parallel = 1;
for i=1:model.numModels
    model.comp{i}.vardist.parallel = 1;
end
    %}

            %%%% TEMP
    if exist('whiteVar')
        fprintf('aaa\n\n\n')
        model.dynamics.kern.comp{2}.variance = whiteVar;
    end
    %%%%

    % Force kernel computations
    params = svargplvmExtractParam(model);
    model = svargplvmExpandParam(model, params);

    %%
    fprintf('# Median of vardist. covars: %d \n',median(median(model.vardist.covars)));
    fprintf('# Min of vardist. covars: %d \n',min(min(model.vardist.covars)));
    fprintf('# Max of vardist. covars: %d \n',max(max(model.vardist.covars)));

    model = svargplvmOptimiseModel(model);

    svargplvmShowScales(model)
    figure
end

%%
    capName = model.globalOpt.dataSetName;
    capName(1) = upper(capName(1));
    errors = fgplvmNearestNeighbour(model.comp{1}, lblsTrain);
    model2 = vargplvmReduceModel(model.comp{1}, 2);
    errors2 = fgplvmNearestNeighbour(model2, lblsTrain);
    fprintf('# Visualisation errors in all dims/in 2D:  %d / %d\n', errors, errors2)
    if printPlots
        vargplvmPrintPlot(model2, lblsTrain, [capName 'Vargplvm'], model.globalOpt.experimentNo);

        %%
        % For 3D plots:
        % labels = model.options{1}.labels; dims = [1 3 5];
        % plot3k({model.X(:,dims(1)) model.X(:,dims(2)) model.X(:,dims(3))}, 'ColorData', labels, 'Marker', {'x',6});
    end

%%
obsMod = 1; % one of the involved sub-models (the one for which we have the data)
infMod = setdiff(1:2, obsMod);

sharedDims = svargplvmFindSharedDims(model);

%---------------------------- PREDICTIONS ---------------
if ~doPredictions
    return
end

%%

if ~exist('testOnTraining')
    testOnTraining=0;
end

%--------------------%
svargplvmPredictions %---- Script returning: ZpredMuAll and mini(the indices for NN)
%--------------------%

%--- For classification we are interested in labels
ZpredMuAll(ZpredMuAll > 0) = 1;
ZpredMuAll(ZpredMuAll < 0) = -1;
NNpred = Ytr{infMod}(mini,:);
realLabels = lblsTest(testInd,:);

% Find the error in the labels ------------
[labelErrors, labelSharingErrors, wrongIndices] = findLabelErrors(realLabels, ZpredMuAll);
[labelErrorsNN, void, wrongIndicesNN] = findLabelErrors(realLabels,NNpred);

fprintf('# Gplvm label errors: %d\n', labelErrors);
fprintf('# Gplvm label sharing errors: %d\n', labelSharingErrors);
fprintf('# NN label errors:    %d\n', labelErrorsNN);

% Confusion matrix
cmat_a = confMatrix( XpredAll, model.X, labelsTrain, labelsTest, 1);
cmat_rel = cmat_a ./ repmat(sum(cmat_a, 2), 1, length(unique(labelsTrain)));
cmat_rel

author: SheffieldML

repo: vargplvm


text:
% PLOT THE COLLISION EVENT DESCRIPTION FIGURES
function [currentFigure,figureSet] = GetFigure_objectCollision(SIM,objectIndex,DATA,currentFigure)
% This function is designed to move through the collsion event history and
% generate figure for each of the collisions.

figureSet = [];
% CHECK COLLISIONS OCCURED
if ~isfield(DATA,'uniqueCollisions') || isempty(DATA.uniqueCollisions)
    fprintf('[%s]\t...No collision data available.\n',SIM.phase);
    return
end

for collisionNumber = 1:numel(DATA.uniqueCollisions)
    [currentFigure,figureSet(collisionNumber)] = BuildObjectCollisionFigure(...
        SIM,...
        objectIndex,...
        DATA,...
        currentFigure,...
        DATA.uniqueCollisions(collisionNumber));
end

% ASSEMBLE TABBED FIGURE
windowHandle = GetTabbedFigure(figureSet,'OpenMAS Collision Overview');
set(windowHandle,'Position',DATA.figureProperties.windowSettings);        % Maximise the figure in the tab
savefig(windowHandle,[SIM.outputPath,'collision-overview']);              % Save the output figure
end

% Build the collision figure for the event
function [currentFigure,figureHandle] = BuildObjectCollisionFigure(SIM,objectIndex,DATA,currentFigure,collisionEvent)

% Generate the figure
figureHandle = figure('Name',sprintf('[t=%.0fs] %s & %s',collisionEvent.time,collisionEvent.name_A,collisionEvent.name_B));
set(figureHandle,'Position', DATA.figureProperties.windowSettings);        % [x y width height]
set(figureHandle,'Color',DATA.figureProperties.figureColor);               % Background colour
ax = axes(figureHandle);
hold on;
view([70 25]);
titlestr = sprintf('Collision Event between %s and %s at t=%s',...         % Declare title string for figure
           collisionEvent.name_A,collisionEvent.name_B,num2str(collisionEvent.time));

% /////////////////////////// ENTITY DATA /////////////////////////////////
% EVENT, OBJECT AND META DATA
META_A = SIM.OBJECTS(SIM.globalIDvector == collisionEvent.objectID_A);
META_B = SIM.OBJECTS(SIM.globalIDvector == collisionEvent.objectID_B);     % The associated META datas
object_A = objectIndex{SIM.globalIDvector == collisionEvent.objectID_A};
object_B = objectIndex{SIM.globalIDvector == collisionEvent.objectID_B};   % The associated object structures
geometry_A = object_A.GEOMETRY;
geometry_B = object_B.GEOMETRY;
[objectAStates] = OMAS_getTrajectoryData_mex(DATA.globalTrajectories,SIM.globalIDvector,META_A.objectID,inf);
[objectBStates] = OMAS_getTrajectoryData_mex(DATA.globalTrajectories,SIM.globalIDvector,META_B.objectID,inf);

% We need to plot:
% - The object geometries
% - The object hitboxes

% ////////////////////// PLOT THE OBJECT GEOMETRY /////////////////////////
% SUBSTITUTE A'S GEOMETRY IF NECESSARY
if size(object_A.GEOMETRY.vertices,1) < 1
    [geometry_A] = OMAS_graphics.defineSphere(zeros(3,1),META_A.radius,10);
end
% GET A'S ROTATIONS
R_A = OMAS_geometry.quaternionToRotationMatrix(collisionEvent.state_A(7:10));
% OBJECT A'S GEOMETRY
patch(ax,...
    'Vertices',(geometry_A.vertices*R_A + collisionEvent.state_A(1:3)'),...
    'Faces',geometry_A.faces,...
    'FaceColor',META_A.colour,...
    'EdgeColor','k',...
    'EdgeAlpha',DATA.figureProperties.edgeAlpha,...
    'FaceAlpha',DATA.figureProperties.faceAlpha,...
    'FaceLighting','gouraud',...
    'LineWidth',DATA.figureProperties.lineWidth);

% SUBSTITUTE B'S GEOMETRY IF NECESSARY
if size(object_B.GEOMETRY.vertices,1) < 1
    [geometry_B] = OMAS_graphics.defineSphere(zeros(3,1),META_B.radius,10);
end
% GET B'S ROTATIONS
R_B = OMAS_geometry.quaternionToRotationMatrix(collisionEvent.state_B(7:10));
% OBJECT B'S GEOMETRY
patch(ax,...
    'Vertices',(geometry_B.vertices*R_B + collisionEvent.state_B(1:3)'),...
    'Faces',geometry_B.faces,...
    'FaceColor',META_B.colour,...
    'EdgeColor','k',...
    'EdgeAlpha',DATA.figureProperties.edgeAlpha,...
    'FaceAlpha',DATA.figureProperties.faceAlpha,...
    'FaceLighting','gouraud',...
    'LineWidth',DATA.figureProperties.lineWidth);

% ////////////////////// PLOT THE OBJECT'S HITBOX /////////////////////////
% Get the hit-box geometry
[hitBoxGeometry_A] = OMAS_graphics.getHitBoxGeometry(object_A.GetGLOBAL(),object_A.GEOMETRY);
% DEFINE A's HITBOX
patch(ax,...
    'Vertices',(hitBoxGeometry_A.vertices + collisionEvent.state_A(1:3)'),...
    'Faces',hitBoxGeometry_A.faces,...
    'FaceColor','r',...
    'EdgeColor','k',...
    'EdgeAlpha',DATA.figureProperties.edgeAlpha*0.8,...
    'FaceAlpha',DATA.figureProperties.faceAlpha*0.8,...
    'FaceLighting','gouraud',...
    'LineWidth',DATA.figureProperties.lineWidth);

% Get the hit-box geometry
[hitBoxGeometry_B] = OMAS_graphics.getHitBoxGeometry(object_B.GetGLOBAL(),object_B.GEOMETRY);
% DEFINE B's HITBOX
patch(ax,...
    'Vertices',(hitBoxGeometry_B.vertices + collisionEvent.state_B(1:3)'),...
    'Faces',hitBoxGeometry_B.faces,...
    'FaceColor','r',...
    'EdgeColor','k',...
    'EdgeAlpha',DATA.figureProperties.edgeAlpha*0.8,...
    'FaceAlpha',DATA.figureProperties.faceAlpha*0.8,...
    'FaceLighting','gouraud',...
    'LineWidth',DATA.figureProperties.lineWidth);

% ////////////////////// PLOT THE TRAJECTORY TRAILS ///////////////////////
plot3(ax,objectAStates(1,:),objectAStates(2,:),objectAStates(3,:),...
    'LineStyle',DATA.figureProperties.lineStyle,...
    'LineWidth',DATA.figureProperties.lineWidth,...
    'Color',META_A.colour);
plot3(ax,objectBStates(1,:),objectBStates(2,:),objectBStates(3,:),...
    'LineStyle',DATA.figureProperties.lineStyle,...
    'LineWidth',DATA.figureProperties.lineWidth,...
    'Color',META_B.colour);

% ////////////////////// PLOT THE VELOCITY VECTORS ////////////////////////
% PLOT THE OBJECT VELOCITY VECTORS
q = quiver3(collisionEvent.state_A(1),...
    collisionEvent.state_A(2),...
    collisionEvent.state_A(3),...
    collisionEvent.state_A(4),...
    collisionEvent.state_A(5),...
    collisionEvent.state_A(6),'r');
q.AutoScaleFactor = 1;
q.LineWidth = DATA.figureProperties.lineWidth;
q = quiver3(collisionEvent.state_B(1),...
    collisionEvent.state_B(2),...
    collisionEvent.state_B(3),...
    collisionEvent.state_B(4),...
    collisionEvent.state_B(5),...
    collisionEvent.state_B(6),'g');
q.AutoScaleFactor = 1;
q.LineWidth = DATA.figureProperties.lineWidth;

% ///////////////////////// ADD LABEL ANNOTATIONS /////////////////////////
annotationText = sprintf('    %s [ID-%d]',META_A.name,META_A.objectID);
text(collisionEvent.state_A(1),collisionEvent.state_A(2),collisionEvent.state_A(3),annotationText);
annotationText = sprintf('    %s [ID-%d]',META_B.name,META_B.objectID);
text(collisionEvent.state_B(1),collisionEvent.state_B(2),collisionEvent.state_B(3),annotationText);

% DEFINE THE PLOT LAYOUT
% Title
title(ax,titlestr,...
    'Interpreter',DATA.figureProperties.interpreter,...
    'fontname',DATA.figureProperties.fontName,...
    'Fontweight',DATA.figureProperties.fontWeight,...
    'FontSize',DATA.figureProperties.titleFontSize,...
    'FontSmoothing','on');
% X-label
xlabel(ax,'x(m)',...
    'Interpreter',DATA.figureProperties.interpreter,...
    'fontname',DATA.figureProperties.fontName,...
    'Fontweight',DATA.figureProperties.fontWeight,...
    'FontSize',DATA.figureProperties.axisFontSize,...
    'FontSmoothing','on');
% Y-Label
ylabel(ax,'y(m)',...
    'Interpreter',DATA.figureProperties.interpreter,...
    'fontname',DATA.figureProperties.fontName,...
    'Fontweight',DATA.figureProperties.fontWeight,...
    'FontSize',DATA.figureProperties.axisFontSize,...
    'FontSmoothing','on');
% Z-Label
zlabel(ax,'z(m)',...
    'Interpreter',DATA.figureProperties.interpreter,...
    'fontname',DATA.figureProperties.fontName,...
    'Fontweight',DATA.figureProperties.fontWeight,...
    'FontSize',DATA.figureProperties.axisFontSize,...
    'FontSmoothing','on');
% Axes
set(ax,...
    'TickLabelInterpreter',DATA.figureProperties.interpreter,...
    'fontname',DATA.figureProperties.fontName,...
    'Fontweight',DATA.figureProperties.fontWeight,...
    'FontSize',DATA.figureProperties.axisFontSize,...
    'FontSmoothing','on',...
    'Color',DATA.figureProperties.axesColor,...
    'GridLineStyle','--',...
    'GridAlpha',0.25,...
    'GridColor','k');
axisLimit = (DATA.figureProperties.objectMaximalRadii + DATA.figureProperties.maxAbsPosition);
xlim(ax,[-axisLimit,axisLimit]);
ylim(ax,[-axisLimit,axisLimit]);
zlim(ax,[-axisLimit,axisLimit]);
axis vis3d equal;
grid on;  box on;
hold off;

% SAVE THE OUTPUT FIGURE
figurePath = [SIM.outputPath,...
    sprintf('collision-event-id-%d-%s-id-%d-%s.fig',...
    collisionEvent.objectID_A,collisionEvent.name_A,...
    collisionEvent.objectID_B,collisionEvent.name_B)];

savefig(figureHandle,figurePath);
% ITERATE PLOT
currentFigure = currentFigure + 1;
end

repo: OpenMAS


text:
function exportfig(varargin)
%EXPORTFIG  Export a figure to Encapsulated Postscript.
%   EXPORTFIG(H, FILENAME) writes the figure H to FILENAME.  H is
%   a figure handle and FILENAME is a string that specifies the
%   name of the output file.
%
%   EXPORTFIG(...,PARAM1,VAL1,PARAM2,VAL2,...) specifies
%   parameters that control various characteristics of the output
%   file.
%
%   Format Paramter:
%     'Format'  one of the strings 'eps','eps2','jpeg','png','preview'
%          specifies the output format. Defaults to 'eps'.
%          The output format 'preview' does not generate an output
%          file but instead creates a new figure window with a
%          preview of the exported figure. In this case the
%          FILENAME parameter is ignored.
%
%     'Preview' one of the strings 'none', 'tiff'
%          specifies a preview for EPS files. Defaults to 'none'.
%
%   Size Parameters:
%     'Width'   a positive scalar
%          specifies the width in the figure's PaperUnits
%     'Height'  a positive scalar
%          specifies the height in the figure's PaperUnits
%
%     Specifying only one dimension sets the other dimension
%     so that the exported aspect ratio is the same as the
%     figure's current aspect ratio.
%     If neither dimension is specified the size defaults to
%     the width and height from the figure's PaperPosition.
%
%   Rendering Parameters:
%     'Color'     one of the strings 'bw', 'gray', 'cmyk'
%         'bw'    specifies that lines and text are exported in
%                 black and all other objects in grayscale
%         'gray'  specifies that all objects are exported in grayscale
%         'cmyk'  specifies that all objects are exported in color
%                 using the CMYK color space
%     'Renderer'  one of the strings 'painters', 'zbuffer', 'opengl'
%         specifies the renderer to use
%     'Resolution'   a positive scalar
%         specifies the resolution in dots-per-inch.
%
%     The default color setting is 'bw'.
%
%   Font Parameters:
%     'FontMode'     one of the strings 'scaled', 'fixed'
%     'FontSize'     a positive scalar
%          in 'scaled' mode multiplies with the font size of each
%          text object to obtain the exported font size
%          in 'fixed' mode specifies the font size of all text
%          objects in points
%     'FontEncoding' one of the strings 'latin1', 'adobe'
%          specifies the character encoding of the font
%
%     If FontMode is 'scaled' but FontSize is not specified then a
%     scaling factor is computed from the ratio of the size of the
%     exported figure to the size of the actual figure. The minimum
%     font size allowed after scaling is 5 points.
%     If FontMode is 'fixed' but FontSize is not specified then the
%     exported font sizes of all text objects is 7 points.
%
%     The default 'FontMode' setting is 'scaled'.
%
%   Line Width Parameters:
%     'LineMode'     one of the strings 'scaled', 'fixed'
%     'LineWidth'    a positive scalar
%          the semantics of LineMode and LineWidth are exactly the
%          same as FontMode and FontSize, except that they apply
%          to line widths instead of font sizes. The minumum line
%          width allowed after scaling is 0.5 points.
%          If LineMode is 'fixed' but LineWidth is not specified
%          then the exported line width of all line objects is 1
%          point.
%
%   Examples:
%     exportfig(gcf,'fig1.eps','height',3);
%       Exports the current figure to the file named 'fig1.eps' with
%       a height of 3 inches (assuming the figure's PaperUnits is
%       inches) and an aspect ratio the same as the figure's aspect
%       ratio on screen.
%
%     exportfig(gcf, 'fig2.eps', 'FontMode', 'fixed',...
%                'FontSize', 10, 'color', 'cmyk' );
%       Exports the current figure to 'fig2.eps' in color with all
%       text in 10 point fonts. The size of the exported figure is
%       the figure's PaperPostion width and height.

if (nargin < 2)
  error('Too few input arguments');
end

% exportfig(H, filename, ...)
H = varargin{1};
if ~ishandle(H) | ~strcmp(get(H,'type'), 'figure')
  error('First argument must be a handle to a figure.');
end
filename = varargin{2};
if ~ischar(filename)
  error('Second argument must be a string.');
end
paramPairs = varargin(3:end);

% Do some validity checking on param-value pairs
if (rem(length(paramPairs),2) ~= 0)
  error(['Invalid input syntax. Optional parameters and values' ...
	 ' must be in pairs.']);
end

format = 'eps';
preview = 'none';
width = -1;
height = -1;
color = 'bw';
fontsize = -1;
fontmode='scaled';
linewidth = -1;
linemode=[];
fontencoding = 'latin1';
renderer = [];
resolution = [];

% Process param-value pairs
args = {};
for k = 1:2:length(paramPairs)
  param = lower(paramPairs{k});
  if (~ischar(param))
    error('Optional parameter names must be strings');
  end
  value = paramPairs{k+1};

  switch (param)
   case 'format'
    format = value;
    if (~strcmp(format,{'eps','eps2','jpeg','png','preview'}))
      error(['Format must be ''eps'', ''eps2'', ''jpeg'', ''png'' or' ...
	     ' ''preview''.']);
    end
   case 'preview'
    preview = value;
    if (~strcmp(preview,{'none','tiff'}))
      error('Preview must be ''none'' or ''tiff''.');
    end
   case 'width'
    width = LocalToNum(value);
    if(~LocalIsPositiveScalar(width))
      error('Width must be a numeric scalar > 0');
    end
   case 'height'
    height = LocalToNum(value);
    if(~LocalIsPositiveScalar(height))
      error('Height must be a numeric scalar > 0');
    end
   case 'color'
    color = lower(value);
    if (~strcmp(color,{'bw','gray','cmyk'}))
      error('Color must be ''bw'', ''gray'' or ''cmyk''.');
    end
   case 'fontmode'
    fontmode = lower(value);
    if (~strcmp(fontmode,{'scaled','fixed'}))
      error('FontMode must be ''scaled'' or ''fixed''.');
    end
   case 'fontsize'
    fontsize = LocalToNum(value);
    if(~LocalIsPositiveScalar(fontsize))
      error('FontSize must be a numeric scalar > 0');
    end
   case 'fontencoding'
    fontencoding = lower(value);
    if (~strcmp(fontencoding,{'latin1','adobe'}))
      error('FontEncoding must be ''latin1'' or ''adobe''.');
    end
   case 'linemode'
    linemode = lower(value);
    if (~strcmp(linemode,{'scaled','fixed'}))
      error('LineMode must be ''scaled'' or ''fixed''.');
    end
   case 'linewidth'
    linewidth = LocalToNum(value);
    if(~LocalIsPositiveScalar(linewidth))
      error('LineWidth must be a numeric scalar > 0');
    end
   case 'renderer'
    renderer = lower(value);
    if (~strcmp(renderer,{'painters','zbuffer','opengl'}))
      error('Renderer must be ''painters'', ''zbuffer'' or ''opengl''.');
    end
   case 'resolution'
    resolution = LocalToNum(value);
    if ~(isnumeric(value) & (prod(size(value)) == 1) & (value >= 0));
      error('Resolution must be a numeric scalar >= 0');
    end
   otherwise
    error(['Unrecognized option ' param '.']);
  end
end

allLines  = findall(H, 'type', 'line');
allText   = findall(H, 'type', 'text');
allAxes   = findall(H, 'type', 'axes');
allImages = findall(H, 'type', 'image');
allLights = findall(H, 'type', 'light');
allPatch  = findall(H, 'type', 'patch');
allSurf   = findall(H, 'type', 'surface');
allRect   = findall(H, 'type', 'rectangle');
allFont   = [allText; allAxes];
allColor  = [allLines; allText; allAxes; allLights];
allMarker = [allLines; allPatch; allSurf];
allEdge   = [allPatch; allSurf];
allCData  = [allImages; allPatch; allSurf];

old.objs = {};
old.prop = {};
old.values = {};

% Process format and preview parameter
showPreview = strcmp(format,'preview');
if showPreview
  format = 'png';
  filename = [tempName '.png'];
end
if strncmp(format,'eps',3) & ~strcmp(preview,'none')
  args = {args{:}, ['-' preview]};
end

hadError = 0;
try
  % Process size parameters
  paperPos = get(H, 'PaperPosition');
  old = LocalPushOldData(old, H, 'PaperPosition', paperPos);
  figureUnits = get(H, 'Units');
  set(H, 'Units', get(H,'PaperUnits'));
  figurePos = get(H, 'Position');
  aspectRatio = figurePos(3)/figurePos(4);
  set(H, 'Units', figureUnits);
  if (width == -1) & (height == -1)
    width = paperPos(3);
    height = paperPos(4);
  elseif (width == -1)
    width = height * aspectRatio;
  elseif (height == -1)
    height = width / aspectRatio;
  end
  set(H, 'PaperPosition', [0 0 width height]);
  paperPosMode = get(H, 'PaperPositionMode');
  old = LocalPushOldData(old, H, 'PaperPositionMode', paperPosMode);
  set(H, 'PaperPositionMode', 'manual');

  % Process rendering parameters
  switch (color)
   case {'bw', 'gray'}
    if ~strcmp(color,'bw') & strncmp(format,'eps',3)
      format = [format 'c'];
    end
    args = {args{:}, ['-d' format]};

    %compute and set gray colormap
    oldcmap = get(H,'Colormap');
    newgrays = 0.30*oldcmap(:,1) + 0.59*oldcmap(:,2) + 0.11*oldcmap(:,3);
    newcmap = [newgrays newgrays newgrays];
    old = LocalPushOldData(old, H, 'Colormap', oldcmap);
    set(H, 'Colormap', newcmap);

    %compute and set ColorSpec and CData properties
    old = LocalUpdateColors(allColor, 'color', old);
    old = LocalUpdateColors(allAxes, 'xcolor', old);
    old = LocalUpdateColors(allAxes, 'ycolor', old);
    old = LocalUpdateColors(allAxes, 'zcolor', old);
    old = LocalUpdateColors(allMarker, 'MarkerEdgeColor', old);
    old = LocalUpdateColors(allMarker, 'MarkerFaceColor', old);
    old = LocalUpdateColors(allEdge, 'EdgeColor', old);
    old = LocalUpdateColors(allEdge, 'FaceColor', old);
    old = LocalUpdateColors(allCData, 'CData', old);

   case 'cmyk'
    if strncmp(format,'eps',3)
      format = [format 'c'];
      args = {args{:}, ['-d' format], '-cmyk'};
    else
      args = {args{:}, ['-d' format]};
    end
   otherwise
    error('Invalid Color parameter');
  end
  if (~isempty(renderer))
    args = {args{:}, ['-' renderer]};
  end
  if (~isempty(resolution)) | ~strncmp(format,'eps',3)
    if isempty(resolution)
      resolution = 0;
    end
    args = {args{:}, ['-r' int2str(resolution)]};
  end

  % Process font parameters
  if (~isempty(fontmode))
    oldfonts = LocalGetAsCell(allFont,'FontSize');
    switch (fontmode)
     case 'fixed'
      oldfontunits = LocalGetAsCell(allFont,'FontUnits');
      old = LocalPushOldData(old, allFont, {'FontUnits'}, oldfontunits);
      set(allFont,'FontUnits','points');
      if (fontsize == -1)
	set(allFont,'FontSize',7);
      else
	set(allFont,'FontSize',fontsize);
      end
     case 'scaled'
      if (fontsize == -1)
	wscale = width/figurePos(3);
	hscale = height/figurePos(4);
	scale = min(wscale, hscale);
      else
	scale = fontsize;
      end
      newfonts = LocalScale(oldfonts,scale,5);
      set(allFont,{'FontSize'},newfonts);
     otherwise
      error('Invalid FontMode parameter');
    end
    % make sure we push the size after the units
    old = LocalPushOldData(old, allFont, {'FontSize'}, oldfonts);
  end
  if strcmp(fontencoding,'adobe') & strncmp(format,'eps',3)
    args = {args{:}, '-adobecset'};
  end

  % Process linewidth parameters
  if (~isempty(linemode))
    oldlines = LocalGetAsCell(allMarker,'LineWidth');
    old = LocalPushOldData(old, allMarker, {'LineWidth'}, oldlines);
    switch (linemode)
     case 'fixed'
      if (linewidth == -1)
	set(allMarker,'LineWidth',1);
      else
	set(allMarker,'LineWidth',linewidth);
      end
     case 'scaled'
      if (linewidth == -1)
	wscale = width/figurePos(3);
	hscale = height/figurePos(4);
	scale = min(wscale, hscale);
      else
	scale = linewidth;
      end
      newlines = LocalScale(oldlines, scale, 0.5);
      set(allMarker,{'LineWidth'},newlines);
     otherwise
      error('Invalid LineMode parameter');
    end
  end

  % Export
  print(H, filename, args{:});

catch
  hadError = 1;
end

% Restore figure settings
for n=1:length(old.objs)
  set(old.objs{n}, old.prop{n}, old.values{n});
end

if hadError
  error(deblank(lasterr));
end

% Show preview if requested
if showPreview
  X = imread(filename,'png');
  delete(filename);
  f = figure( 'Name', 'Preview', ...
	      'Menubar', 'none', ...
	      'NumberTitle', 'off', ...
	      'Visible', 'off');
  image(X);
  axis image;
  ax = findobj(f, 'type', 'axes');
  set(ax, 'Units', get(H,'PaperUnits'), ...
	  'Position', [0 0 width height], ...
	  'Visible', 'off');
  set(ax, 'Units', 'pixels');
  axesPos = get(ax,'Position');
  figPos = get(f,'Position');
  rootSize = get(0,'ScreenSize');
  figPos(3:4) = axesPos(3:4);
  if figPos(1) + figPos(3) > rootSize(3)
    figPos(1) = rootSize(3) - figPos(3) - 50;
  end
  if figPos(2) + figPos(4) > rootSize(4)
    figPos(2) = rootSize(4) - figPos(4) - 50;
  end
  set(f, 'Position',figPos, ...
	 'Visible', 'on');
end

%
%  Local Functions
%

function outData = LocalPushOldData(inData, objs, prop, values)
outData.objs = {inData.objs{:}, objs};
outData.prop = {inData.prop{:}, prop};
outData.values = {inData.values{:}, values};

function cellArray = LocalGetAsCell(fig,prop);
cellArray = get(fig,prop);
if (~isempty(cellArray)) & (~iscell(cellArray))
  cellArray = {cellArray};
end

function newArray = LocalScale(inArray, scale, minValue)
n = length(inArray);
newArray = cell(n,1);
for k=1:n
  newArray{k} = max(minValue,scale*inArray{k}(1));
end

function newArray = LocalMapToGray(inArray);
n = length(inArray);
newArray = cell(n,1);
for k=1:n
  color = inArray{k};
  if (~isempty(color))
    if ischar(color)
      switch color(1)
       case 'y'
	color = [1 1 0];
       case 'm'
	color = [1 0 1];
       case 'c'
	color = [0 1 1];
       case 'r'
	color = [1 0 0];
       case 'g'
	color = [0 1 0];
       case 'b'
	color = [0 0 1];
       case 'w'
	color = [1 1 1];
       case 'k'
	color = [0 0 0];
       otherwise
	newArray{k} = color;
      end
    end
    if ~ischar(color)
      color = 0.30*color(1) + 0.59*color(2) + 0.11*color(3);
    end
  end
  if isempty(color) | ischar(color)
    newArray{k} = color;
  else
    newArray{k} = [color color color];
  end
end

function newArray = LocalMapCData(inArray);
n = length(inArray);
newArray = cell(n,1);
for k=1:n
  color = inArray{k};
  if (ndims(color) == 3) & isa(color,'double')
    gray = 0.30*color(:,:,1) + 0.59*color(:,:,2) + 0.11*color(:,:,3);
    color(:,:,1) = gray;
    color(:,:,2) = gray;
    color(:,:,3) = gray;
  end
  newArray{k} = color;
end

function outData = LocalUpdateColors(inArray, prop, inData)
value = LocalGetAsCell(inArray,prop);
outData.objs = {inData.objs{:}, inArray};
outData.prop = {inData.prop{:}, {prop}};
outData.values = {inData.values{:}, value};
if (~isempty(value))
  if strcmp(prop,'CData')
    value = LocalMapCData(value);
  else
    value = LocalMapToGray(value);
  end
  set(inArray,{prop},value);
end

function bool = LocalIsPositiveScalar(value)
bool = isnumeric(value) & ...
       prod(size(value)) == 1 & ...
       value > 0;

function value = LocalToNum(value)
if ischar(value)
  value = str2num(value);
end

author: johannesgerer

repo: jburkardt-m


text:
% =========================================================================
% (c) 2018 Ronald Nissel, ronald.nissel@gmail.com
% =========================================================================
% Allows to reproduce Figure 2, 3, 4 and 5 of "Doubly-Selective Channel
% Estimation in FBMC-OQAM and OFDM Systems", R. Nissel, et al, IEEE VTC
% Fall, 2018. In particular, this script simulates an FBMC and OFDM
% transmission over a doubly-selective channel, including doubly-selective
% MMSE channel estimation and interference cancellation. Note that, to
% truly reproduce the figures, the lines 41-46 must be uncommented.

clear; close all;
addpath('./Theory');

%% Parameters
% Simulation
M_SNR_dB                  = [10:5:40];              % Signal-to-Noise Ratio in dB
NrRepetitions             = 25;                     % Number of Monte Carlo repetition (different channel realizations)
ZeroThresholdSparse       = 8;                      % Set some matrix values, smaller than "10^(-ZeroThresholdSparse)", to zero.
PlotIterationStepsSNRdB   = 35;                     % Plot the BER over iteration step for an SNR of 35dB.

% FBMC and OFDM parameters
L                         = 12*2;                   % Number of Subcarriers, one resource block consists of 12 subcarriers (and 0.5ms in time)
F                         = 15e3;                   % Subcarrier Spacing in Hz
SamplingRate              = F*12*2;                 % Sampling rate in (Samples/s)
NrSubframes               = 1;                      % Number of subframes. One subframe requires 1ms for F=15kHz.
QAM_ModulationOrder       = 256;                    % QAM signal constellation order, 4, 16, 64, 256, 1024,...

% Channel estimation parameters
PilotToDataPowerOffset    = 2;                      % Pilot to data power offset for OFDM. In FBMC data spreading, the power offset is twice this number.
PilotToDataPowerOffsetAux = 4.685;                  % Pilot to data power offset for FBMC, auxiliary method.
NrIterations              = 4;                      % Number of iterations for interference cancellation scheme.

% Channel
Velocity_kmh              = 500;                    % Velocity in km/h. Note that [mph]*1.6=[kmh] and [m/s]*3.6=[kmh]
PowerDelayProfile         = 'VehicularA';           % Channel model, either string or vector: 'Flat', 'AWGN', 'PedestrianA', 'PedestrianB', 'VehicularA', 'VehicularB', 'ExtendedPedestrianA', 'ExtendedPedestrianB', or 'TDL-A_xxns','TDL-B_xxns','TDL-C_xxns' (with xx the RMS delay spread in ns, e.g. 'TDL-A_30ns'), or [1 0 0.2] (Self-defined power delay profile which depends on the sampling rate)

% ###########################################################################
% % In the paper:
% M_SNR_dB                = [10:2:40];
% PlotIterationStepsSNRdB = 32;
% NrRepetitions           = 1000;
% SamplingRate            = F*14*14;
% NrSubframes             = 2;
% ###########################################################################

%% FBMC Object
FBMC = Modulation.FBMC(...
    L,...                               % Number subcarriers
    30*NrSubframes,...                  % Number FBMC symbols
    F,...                               % Subcarrier spacing (Hz)
    SamplingRate,...                    % Sampling rate (Samples/s)
    0,...                               % Intermediate frequency first subcarrier (Hz)
    false,...                           % Transmit real valued signal
    'Hermite-OQAM',...                  % Prototype filter (Hermite, PHYDYAS, RRC) and OQAM or QAM,
    8, ...                              % Overlapping factor (also determines oversampling in the frequency domain)
    0, ...                              % Initial phase shift
    true ...                            % Polyphase implementation
    );

%% OFDM Object (Add zeroes to the OFDM signal so that it fits the FBMC signal)
ZeroGuardTimeLength = ((FBMC.Nr.SamplesTotal-(round((1/15e3/14)*SamplingRate)+round(SamplingRate/15e3))*14*NrSubframes)/2)/SamplingRate;
OFDM = Modulation.OFDM(...
    L,...                           % Number Subcarriers
    14*NrSubframes,...              % Number OFDM Symbols
    F,...                           % Subcarrier spacing (Hz)
    SamplingRate,...                % Sampling rate (Samples/s)
    0,...                           % Intermediate frequency first subcarrier (Hz)
    false,...                       % Transmitreal valued signal
    1/15e3/14, ...                  % Cyclic prefix length (s)
    ZeroGuardTimeLength ...         % Zero guard length (s)
    );

%% Check Number of Samples
if  OFDM.Nr.SamplesTotal~=FBMC.Nr.SamplesTotal
   error('Total number of samples must be the same for OFDM and FBMC.');
end
N = OFDM.Nr.SamplesTotal;

%% PAM and QAM Object
PAM = Modulation.SignalConstellation(sqrt(QAM_ModulationOrder),'PAM');
QAM = Modulation.SignalConstellation(QAM_ModulationOrder,'QAM');

%% Pilot Matrices, 0="Data", 1=Pilot
PilotMatrix_OFDM = zeros(OFDM.Nr.Subcarriers,14);
PilotMatrix_OFDM(2:2*6:end,2:2*3.5:end)=1;
PilotMatrix_OFDM(5:2*6:end,6:2*3.5:end)=1;
PilotMatrix_OFDM(8:2*6:end,2:2*3.5:end)=1;
PilotMatrix_OFDM(11:2*6:end,6:2*3.5:end)=1;
PilotMatrix_OFDM = repmat(PilotMatrix_OFDM,[1 NrSubframes]);

PilotMatrix_FBMC =  zeros(FBMC.Nr.Subcarriers,30);
PilotMatrix_FBMC(2:12:end,3:16:end) = 1;
PilotMatrix_FBMC(5:12:end,11:16:end) = 1;
PilotMatrix_FBMC(8:12:end,3+1:16:end) = 1;
PilotMatrix_FBMC(11:12:end,11+1:16:end) = 1;
PilotMatrix_FBMC = repmat(PilotMatrix_FBMC,[1 NrSubframes]);

AuxilaryPilotMatrix_FBMC = PilotMatrix_FBMC;
[a,b] = find(PilotMatrix_FBMC);
for i_pilot = 1:length(a)
    AuxilaryPilotMatrix_FBMC(a(i_pilot)+1,b(i_pilot))=-1;
    AuxilaryPilotMatrix_FBMC(a(i_pilot)-1,b(i_pilot))=-1;
    AuxilaryPilotMatrix_FBMC(a(i_pilot),b(i_pilot)+1)=-1;
    AuxilaryPilotMatrix_FBMC(a(i_pilot),b(i_pilot)-1)=-1;
end

%% Cancel Imaginary Interference At Pilot Position Object (Precoding Matrix)
AuxiliaryMethod = ChannelEstimation.ImaginaryInterferenceCancellationAtPilotPosition(...
    'Auxiliary', ...                                    % Cancellation method
    AuxilaryPilotMatrix_FBMC, ...                       % PilotAndAuxiliaryMatrix
    FBMC.GetFBMCMatrix, ...                             % Imaginary interference matrix
    28, ...                                             % Cancel 28 closest interferers
    PilotToDataPowerOffsetAux ...                       % Pilot to data power offset
    );
CodingMethod = ChannelEstimation.ImaginaryInterferenceCancellationAtPilotPosition(...
    'Coding', ...                                       % Cancellation method
    PilotMatrix_FBMC, ...                               % PilotMatrix
    FBMC.GetFBMCMatrix, ...                             % Imaginary interference matrix
    20, ...                                             % Cancel 28 closest interferers
    2*PilotToDataPowerOffset ...                        % Pilot to data power offset
    );

NrPilotSymbols = sum(PilotMatrix_OFDM(:)==1);
NrDataSymbols_OFDM = sum(PilotMatrix_OFDM(:)==0);

PilotMapping_OFDM = zeros(numel(PilotMatrix_OFDM));
PilotMapping_OFDM(PilotMatrix_OFDM(:)==1,1:NrPilotSymbols) = sqrt(PilotToDataPowerOffset)*eye(NrPilotSymbols);
PilotMapping_OFDM(PilotMatrix_OFDM(:)==0,NrPilotSymbols+1:end) = eye(NrDataSymbols_OFDM);
PilotMapping_OFDM = PilotMapping_OFDM/sqrt(mean(diag(PilotMapping_OFDM*PilotMapping_OFDM')));
DataPowerReduction_OFDM = numel(PilotMatrix_OFDM)/(NrPilotSymbols*PilotToDataPowerOffset+NrDataSymbols_OFDM);

Kappa_Aux  = AuxiliaryMethod.PilotToDataPowerOffset*AuxiliaryMethod.DataPowerReduction;
Kappa_Cod  = CodingMethod.PilotToDataPowerOffset*CodingMethod.DataPowerReduction;
Kappa_OFDM = PilotToDataPowerOffset*DataPowerReduction_OFDM;

%% Evaluate only the BER near the center. i.e., ignore edges
ConsideredTimeFrequencyPositions_FBMC = zeros(size(PilotMatrix_FBMC));
ConsideredTimeFrequencyPositions_FBMC(5:end-4,11:end-10)= 1;

ConsideredTimeFrequencyPositions_OFDM = zeros(size(PilotMatrix_OFDM));
ConsideredTimeFrequencyPositions_OFDM(5:end-4,6:end-5)= 1;

for i_lk = 1: AuxiliaryMethod.NrDataSymbols
    xD_Temp = zeros(AuxiliaryMethod.NrDataSymbols,1);
    xD_Temp(i_lk) = 1;
    x_FBMC_Aux = AuxiliaryMethod.PrecodingMatrix*[zeros(NrPilotSymbols,1);xD_Temp];
    ConsideredDataPositions_FBMC_Aux(i_lk) = sum(abs(x_FBMC_Aux(1==ConsideredTimeFrequencyPositions_FBMC.*not(AuxilaryPilotMatrix_FBMC))))>AuxiliaryMethod.DataPowerReduction*0.9;
end
for i_lk = 1: CodingMethod.NrDataSymbols
    xD_Temp = zeros(CodingMethod.NrDataSymbols,1);
    xD_Temp(i_lk) = 1;
    x_FBMC_Cod = CodingMethod.PrecodingMatrix*[zeros(NrPilotSymbols,1);xD_Temp];
    ConsideredDataPositions_FBMC_Cod(i_lk) =not(any(x_FBMC_Cod(not(ConsideredTimeFrequencyPositions_FBMC))));
end
for i_lk = 1: NrDataSymbols_OFDM
    xD_Temp = zeros(NrDataSymbols_OFDM,1);
    xD_Temp(i_lk) = 1;
    x_OFDM = PilotMapping_OFDM*[zeros(NrPilotSymbols,1);xD_Temp];
    ConsideredDataPositions_OFDM(i_lk) = sum(abs(x_OFDM(1==ConsideredTimeFrequencyPositions_OFDM.*not(PilotMatrix_OFDM))))>DataPowerReduction_OFDM*0.9;
end

ConsideredBits_FBMC_Aux = reshape(repmat(ConsideredDataPositions_FBMC_Aux,log2(QAM_ModulationOrder)/2,1),[],1);
ConsideredBits_FBMC_Cod = reshape(repmat(ConsideredDataPositions_FBMC_Cod,log2(QAM_ModulationOrder)/2,1),[],1);
ConsideredBits_OFDM     = reshape(repmat(ConsideredDataPositions_OFDM,log2(QAM_ModulationOrder),1),[],1);

%% Channel Model Object
ChannelModel = Channel.FastFading(...
    SamplingRate,...                                   % Sampling rate (Samples/s)
    PowerDelayProfile,...                              % Power delay profile, either string or vector: 'Flat', 'AWGN', 'PedestrianA', 'PedestrianB', 'VehicularA', 'VehicularB', 'ExtendedPedestrianA', 'ExtendedPedestrianB', or 'TDL-A_xxns','TDL-B_xxns','TDL-C_xxns' (with xx the RMS delay spread in ns, e.g. 'TDL-A_30ns'), or [1 0 0.2] (Self-defined power delay profile which depends on the sampling rate)
    N,...                                              % Number of total samples
    Velocity_kmh/3.6*2.5e9/2.998e8,...                 % Maximum Doppler shift: Velocity_kmh/3.6*CarrierFrequency/2.998e8
    'Jakes',...                                        % Which Doppler model: 'Jakes', 'Uniform', 'Discrete-Jakes', 'Discrete-Uniform'. For "Discrete-", we assume a discrete Doppler spectrum to improve the simulation time. This only works accuratly if the number of samples and the velocity is sufficiently large
    200,...                                            % Number of paths for the WSSUS process. Only relevant for a 'Jakes' and 'Uniform' Doppler spectrum
    1,...                                              % Number of transmit antennas
    1,...                                              % Number of receive antennas
    1 ...                                              % Gives a warning if the predefined delay taps of the channel do not fit the sampling rate. This is usually not much of a problem if they are approximatly the same.
    );
R_vecH = ChannelModel.GetCorrelationMatrix;

%% Precalculate Transmit and Receive Matrices
G_FBMC = FBMC.GetTXMatrix;
Q_FBMC = (FBMC.GetRXMatrix)';

G_OFDM = OFDM.GetTXMatrix;
Q_OFDM = (OFDM.GetRXMatrix)';

GP_FBMC = G_FBMC(:,PilotMatrix_FBMC(:)==1);
GP_OFDM = G_OFDM(:,PilotMatrix_OFDM(:)==1);

QP_FBMC = Q_FBMC(:,PilotMatrix_FBMC(:)==1);
QP_OFDM = Q_OFDM(:,PilotMatrix_OFDM(:)==1);

G_Aux = G_FBMC*AuxiliaryMethod.PrecodingMatrix;
G_Cod = G_FBMC*CodingMethod.PrecodingMatrix;
G_OFDM_PilotMapping = G_OFDM*PilotMapping_OFDM;

%% Calculate Correlation Matrices
disp('Calculate Correlation Matrix of Pilot Estimates (no Noise, no Interference)...');
R_hP_FBMC  = nan(NrPilotSymbols,NrPilotSymbols);
R_hP_OFDM = nan(NrPilotSymbols,NrPilotSymbols);
for j_pilot = 1:NrPilotSymbols
    R_hP_FBMC(:,j_pilot) = sum((QP_FBMC'*reshape(R_vecH*kron(GP_FBMC(:,j_pilot).',QP_FBMC(:,j_pilot)')',N,N)).*(GP_FBMC.'),2);
    R_hP_OFDM(:,j_pilot) = sum((QP_OFDM'*reshape(R_vecH*kron(GP_OFDM(:,j_pilot).',QP_OFDM(:,j_pilot)')',N,N)).*(GP_OFDM.'),2);
end

disp('Calculate Correlation Matrix of Pilot Estimates (no Interference)...');
R_hP_est_noNoise_FBMC_Aux = R_hP_FBMC;
R_hP_est_noNoise_FBMC_Cod = R_hP_FBMC;
R_hP_est_noNoise_OFDM = R_hP_OFDM;
for i_pilots = 1: NrPilotSymbols
    % FBMC Auxiliary Method, Similar as Equation (13) but computationally more efficient
    Temp = kron(sparse(eye(N)),QP_FBMC(:,i_pilots)')/sqrt(Kappa_Aux);
    R_hP_est_noNoise_FBMC_Aux(i_pilots,i_pilots)=abs(sum(sum((G_Aux.'*(Temp*R_vecH*Temp')).*G_Aux',2)));

    % FBMC Coding
    Temp = kron(sparse(eye(N)),QP_FBMC(:,i_pilots)')/sqrt(Kappa_Cod);
    R_hP_est_noNoise_FBMC_Cod(i_pilots,i_pilots)=abs(sum(sum((G_Cod.'*(Temp*R_vecH*Temp')).*G_Cod',2)));

    % OFDM
    Temp = kron(sparse(eye(N)),QP_OFDM(:,i_pilots)')/sqrt(Kappa_OFDM);
    R_hP_est_noNoise_OFDM(i_pilots,i_pilots)=abs(sum(sum((G_OFDM_PilotMapping.'*(Temp*R_vecH*Temp')).*G_OFDM_PilotMapping',2)));
end

disp('Calculate Correlation Matrix of Pilot Estimates...');
R_hP_est_FBMC_Aux = repmat(R_hP_est_noNoise_FBMC_Aux,[1 1 length(M_SNR_dB)]);
R_hP_est_FBMC_Cod = repmat(R_hP_est_noNoise_FBMC_Cod,[1 1 length(M_SNR_dB)]);
R_hP_est_OFDM = repmat(R_hP_est_noNoise_OFDM,[1 1 length(M_SNR_dB)]);
for i_SNR = 1:length(M_SNR_dB)
    SNR_dB = M_SNR_dB(i_SNR);
    Pn_time = SamplingRate/(F*L)*10^(-SNR_dB/10);
    for i_pilots = 1: NrPilotSymbols
        R_hP_est_FBMC_Aux(i_pilots,i_pilots,i_SNR)=R_hP_est_noNoise_FBMC_Aux(i_pilots,i_pilots)+Pn_time*QP_FBMC(:,i_pilots)'*QP_FBMC(:,i_pilots)/(Kappa_Aux);
        R_hP_est_FBMC_Cod(i_pilots,i_pilots,i_SNR)=R_hP_est_noNoise_FBMC_Cod(i_pilots,i_pilots)+Pn_time*QP_FBMC(:,i_pilots)'*QP_FBMC(:,i_pilots)/(Kappa_Cod);
        R_hP_est_OFDM(i_pilots,i_pilots,i_SNR)=R_hP_est_noNoise_OFDM(i_pilots,i_pilots)+Pn_time*QP_OFDM(:,i_pilots)'*QP_OFDM(:,i_pilots)/Kappa_OFDM;
    end
end

R_hP_est_noInterference_FBMC_Aux = R_hP_est_FBMC_Aux-repmat((R_hP_est_noNoise_FBMC_Aux-R_hP_FBMC),[1 1 length(M_SNR_dB)]);
R_hP_est_noInterference_FBMC_Cod = R_hP_est_FBMC_Cod-repmat((R_hP_est_noNoise_FBMC_Cod-R_hP_FBMC),[1 1 length(M_SNR_dB)]);
R_hP_est_noInterference_OFDM = R_hP_est_OFDM-repmat((R_hP_est_noNoise_OFDM-R_hP_OFDM),[1 1 length(M_SNR_dB)]);

disp('Calculate Correlation Matrix between Transmission-Matrix D and Pilot Estimates...');
R_Dij_hP_FBMC = sparse(size(G_FBMC,2)^2,NrPilotSymbols);
R_Dij_hP_OFDM = sparse(size(G_OFDM,2)^2,NrPilotSymbols);
for i_pilot = 1: NrPilotSymbols
    R_Dij_hP_FBMC_Temp = reshape(Q_FBMC'*reshape(R_vecH*kron(GP_FBMC(:,i_pilot).',QP_FBMC(:,i_pilot)')',N,N)*G_FBMC,[],1);
    R_Dij_hP_OFDM_Temp = reshape(Q_OFDM'*reshape(R_vecH*kron(GP_OFDM(:,i_pilot).',QP_OFDM(:,i_pilot)')',N,N)*G_OFDM,[],1);

    R_Dij_hP_FBMC_Temp(abs(R_Dij_hP_FBMC_Temp)<10^(-ZeroThresholdSparse))=0;
    R_Dij_hP_OFDM_Temp(abs(R_Dij_hP_OFDM_Temp)<10^(-ZeroThresholdSparse))=0;

    R_Dij_hP_FBMC(:,i_pilot) = R_Dij_hP_FBMC_Temp;
    R_Dij_hP_OFDM(:,i_pilot) = R_Dij_hP_OFDM_Temp;
end
clear R_vecH;

%% Calculate SIR at Pilot Positions
SIR_P_FBMC_Aux_dB = 10*log10(trace(abs(R_hP_FBMC))./trace(abs(R_hP_est_noNoise_FBMC_Aux-R_hP_FBMC)));
SIR_P_FBMC_Cod_dB = 10*log10(trace(abs(R_hP_FBMC))./trace(abs(R_hP_est_noNoise_FBMC_Cod-R_hP_FBMC)));
SIR_P_OFDM_dB     = 10*log10(trace(abs(R_hP_OFDM))./trace(abs(R_hP_est_noNoise_OFDM-R_hP_OFDM)));

%% Calculate MMSE Estimation Matrix Using Correlation
disp('Calculate MMSE Solution ...');
W_MMSE_FBMC_Aux = sparse(size(G_FBMC,2)*size(G_FBMC,2)*NrPilotSymbols,length(M_SNR_dB));
W_MMSE_FBMC_Cod = sparse(size(G_FBMC,2)*size(G_FBMC,2)*NrPilotSymbols,length(M_SNR_dB));
W_MMSE_OFDM     = sparse(size(G_OFDM,2)*size(G_OFDM,2)*NrPilotSymbols,length(M_SNR_dB));
for i_SNR = 1:length(M_SNR_dB)
    W_MMSE_FBMC_Aux_Temp = reshape(R_Dij_hP_FBMC*pinv(R_hP_est_FBMC_Aux(:,:,i_SNR)),size(G_FBMC,2)*size(G_FBMC,2)*NrPilotSymbols,1);
    W_MMSE_FBMC_Cod_Temp = reshape(R_Dij_hP_FBMC*pinv(R_hP_est_FBMC_Cod(:,:,i_SNR)),size(G_FBMC,2)*size(G_FBMC,2)*NrPilotSymbols,1);
    W_MMSE_OFDM_Temp     = reshape(R_Dij_hP_OFDM*pinv(R_hP_est_OFDM(:,:,i_SNR)),size(G_OFDM,2)*size(G_OFDM,2)*NrPilotSymbols,1);

    W_MMSE_FBMC_Aux_Temp(abs(W_MMSE_FBMC_Aux_Temp)<10^(-ZeroThresholdSparse))=0;
    W_MMSE_FBMC_Cod_Temp(abs(W_MMSE_FBMC_Cod_Temp)<10^(-ZeroThresholdSparse))=0;
    W_MMSE_OFDM_Temp(abs(W_MMSE_OFDM_Temp)<10^(-ZeroThresholdSparse))=0;

    W_MMSE_FBMC_Aux(:,i_SNR) = W_MMSE_FBMC_Aux_Temp;
    W_MMSE_FBMC_Cod(:,i_SNR) = W_MMSE_FBMC_Cod_Temp;
    W_MMSE_OFDM(:,i_SNR)     = W_MMSE_OFDM_Temp;
end

%% Calculate MMSE Estimation if no Interference is present at the pilot positions
disp('Calculate MMSE Solution (no Interference)...');
W_MMSE_noInterference_FBMC_Aux = sparse(size(G_FBMC,2)*size(G_FBMC,2)*NrPilotSymbols,length(M_SNR_dB));
W_MMSE_noInterference_FBMC_Cod = sparse(size(G_FBMC,2)*size(G_FBMC,2)*NrPilotSymbols,length(M_SNR_dB));
W_MMSE_noInterference_OFDM     = sparse(size(G_OFDM,2)*size(G_OFDM,2)*NrPilotSymbols,length(M_SNR_dB));
for i_SNR = 1:length(M_SNR_dB)
    W_MMSE_noInterference_FBMC_Aux_Temp = reshape(R_Dij_hP_FBMC*pinv(R_hP_est_noInterference_FBMC_Aux(:,:,i_SNR)),size(G_FBMC,2)*size(G_FBMC,2)*NrPilotSymbols,1);
    W_MMSE_noInterference_FBMC_Cod_Temp = reshape(R_Dij_hP_FBMC*pinv(R_hP_est_noInterference_FBMC_Cod(:,:,i_SNR)),size(G_FBMC,2)*size(G_FBMC,2)*NrPilotSymbols,1);
    W_MMSE_noInterference_OFDM_Temp = reshape(R_Dij_hP_OFDM*pinv(R_hP_est_noInterference_OFDM(:,:,i_SNR)),size(G_OFDM,2)*size(G_OFDM,2)*NrPilotSymbols,1);

    W_MMSE_noInterference_FBMC_Aux_Temp(abs(W_MMSE_noInterference_FBMC_Aux_Temp)<10^(-ZeroThresholdSparse))=0;
    W_MMSE_noInterference_FBMC_Cod_Temp(abs(W_MMSE_noInterference_FBMC_Cod_Temp)<10^(-ZeroThresholdSparse))=0;
    W_MMSE_noInterference_OFDM_Temp(abs(W_MMSE_noInterference_OFDM_Temp)<10^(-ZeroThresholdSparse))=0;

    W_MMSE_noInterference_FBMC_Aux(:,i_SNR) = W_MMSE_noInterference_FBMC_Aux_Temp;
    W_MMSE_noInterference_FBMC_Cod(:,i_SNR) = W_MMSE_noInterference_FBMC_Cod_Temp;
    W_MMSE_noInterference_OFDM(:,i_SNR)     = W_MMSE_noInterference_OFDM_Temp;
end

%% Theoretical BEP for a Doubly Flat Rayleigh Channel
M_SNR_dB_morePoints = min(M_SNR_dB):1:max(M_SNR_dB);
BitErrorProbability = BitErrorProbabilityDoublyFlatRayleigh(M_SNR_dB_morePoints,QAM.SymbolMapping,QAM.BitMapping);

%% Preallocate for Parfor
BER_FBMC_Aux_PerfectCSI_InterferenceCancellation        = nan(length(M_SNR_dB),NrRepetitions,NrIterations);
BER_FBMC_Aux_PerfectCSI_InterferenceCancellation_NoEdge = nan(length(M_SNR_dB),NrRepetitions,NrIterations);
BER_FBMC_Cod_PerfectCSI_InterferenceCancellation        = nan(length(M_SNR_dB),NrRepetitions,NrIterations);
BER_FBMC_Cod_PerfectCSI_InterferenceCancellation_NoEdge = nan(length(M_SNR_dB),NrRepetitions,NrIterations);
BER_OFDM_PerfectCSI_InterferenceCancellation            = nan(length(M_SNR_dB),NrRepetitions,NrIterations);
BER_OFDM_PerfectCSI_InterferenceCancellation_NoEdge     = nan(length(M_SNR_dB),NrRepetitions,NrIterations);
BER_FBMC_Aux_OneTapEqualizer_PerfectCSI                 = nan(length(M_SNR_dB),NrRepetitions);
BER_FBMC_Aux_OneTapEqualizer_PerfectCSI_NoEdge          = nan(length(M_SNR_dB),NrRepetitions);
BER_FBMC_Cod_OneTapEqualizer_PerfectCSI                 = nan(length(M_SNR_dB),NrRepetitions);
BER_FBMC_Cod_OneTapEqualizer_PerfectCSI_NoEdge          = nan(length(M_SNR_dB),NrRepetitions);
BER_OFDM_OneTapEqualizer_PerfectCSI                     = nan(length(M_SNR_dB),NrRepetitions);
BER_OFDM_OneTapEqualizer_PerfectCSI_NoEdge              = nan(length(M_SNR_dB),NrRepetitions);
BER_FBMC_Aux_InterferenceCancellation                   = nan(length(M_SNR_dB),NrRepetitions,NrIterations);
BER_FBMC_Cod_InterferenceCancellation                   = nan(length(M_SNR_dB),NrRepetitions,NrIterations);
BER_OFDM_InterferenceCancellation                       = nan(length(M_SNR_dB),NrRepetitions,NrIterations);
BER_FBMC_Aux_InterferenceCancellation_NoEdge            = nan(length(M_SNR_dB),NrRepetitions,NrIterations);
BER_FBMC_Cod_InterferenceCancellation_NoEdge            = nan(length(M_SNR_dB),NrRepetitions,NrIterations);
BER_OFDM_InterferenceCancellation_NoEdge                = nan(length(M_SNR_dB),NrRepetitions,NrIterations);
BER_FBMC_Aux_OneTapEqualizer                            = nan(length(M_SNR_dB),NrRepetitions);
BER_FBMC_Aux_OneTapEqualizer_NoEdge                     = nan(length(M_SNR_dB),NrRepetitions);
BER_FBMC_Cod_OneTapEqualizer                            = nan(length(M_SNR_dB),NrRepetitions);
BER_FBMC_Cod_OneTapEqualizer_NoEdge                     = nan(length(M_SNR_dB),NrRepetitions);
BER_OFDM_OneTapEqualizer                                = nan(length(M_SNR_dB),NrRepetitions);
BER_OFDM_OneTapEqualizer_NoEdge                         = nan(length(M_SNR_dB),NrRepetitions);

%% Start Simulation
tic
disp('Monte Carlo Simulation ...');
for i_rep = 1:NrRepetitions
    %% Update Channel
    ChannelModel.NewRealization;

    %% Binary Data
    BinaryDataStream_FBMC_Aux = randi([0 1],AuxiliaryMethod.NrDataSymbols*log2(PAM.ModulationOrder),1);
    BinaryDataStream_FBMC_Cod = randi([0 1],CodingMethod.NrDataSymbols*log2(PAM.ModulationOrder),1);
    BinaryDataStream_OFDM     = randi([0 1],NrDataSymbols_OFDM*log2(QAM.ModulationOrder),1);

    %% Data Symbols
    xD_FBMC_Aux = PAM.Bit2Symbol(BinaryDataStream_FBMC_Aux);
    xD_FBMC_Cod = PAM.Bit2Symbol(BinaryDataStream_FBMC_Cod);
    xD_OFDM     = QAM.Bit2Symbol(BinaryDataStream_OFDM);

    %% Pilot Symbols
    xP_FBMC = PAM.SymbolMapping(randi(PAM.ModulationOrder,AuxiliaryMethod.NrPilotSymbols,1));
    xP_FBMC = xP_FBMC./abs(xP_FBMC);
    xP_OFDM = QAM.SymbolMapping(randi(QAM.ModulationOrder,NrPilotSymbols,1));
    xP_OFDM = xP_OFDM./abs(xP_OFDM);

    %% Transmitted Data Symbols (Map bin to symbol)
    x_FBMC_Aux = AuxiliaryMethod.PrecodingMatrix*[xP_FBMC;xD_FBMC_Aux];
    x_FBMC_Cod = CodingMethod.PrecodingMatrix*[xP_FBMC;xD_FBMC_Cod];
    x_OFDM     = PilotMapping_OFDM*[xP_OFDM;xD_OFDM];

    %% Transmitted Signal (time domain)
    s_FBMC_Aux = G_FBMC*x_FBMC_Aux(:); % Same as "FBMC.Modulation(x_FBMC_Aux)" which is computationally more efficient. But G_FBMC is consistent with the paper.
    s_FBMC_Cod = G_FBMC*x_FBMC_Cod(:);
    s_OFDM     = G_OFDM*x_OFDM(:);

    %% Channel
    ConvolutionMatrix = ChannelModel.GetConvolutionMatrix{1};

    r_FBMC_Aux_noNoise = ConvolutionMatrix*s_FBMC_Aux;
    r_FBMC_Cod_noNoise = ConvolutionMatrix*s_FBMC_Cod;
    r_OFDM_noNoise     = ConvolutionMatrix*s_OFDM;

    %% Transmission Matrix
    D_FBMC = Q_FBMC'*ConvolutionMatrix*G_FBMC;
    D_OFDM = Q_OFDM'*ConvolutionMatrix*G_OFDM;

    %% One-Tap Channel (for perfect channel knowledge)
    h_FBMC = diag(D_FBMC);
    h_OFDM = diag(D_OFDM);

    for i_SNR = 1:length(M_SNR_dB)
        %% Add Noise
        SNR_dB  = M_SNR_dB(i_SNR);
        Pn_time = SamplingRate/(F*L)*10^(-SNR_dB/10);
        noise   = sqrt(Pn_time/2)*(randn(size(s_OFDM))+1j*randn(size(s_OFDM)));

        r_FBMC_Aux = r_FBMC_Aux_noNoise+noise;
        r_FBMC_Cod = r_FBMC_Cod_noNoise+noise;
        r_OFDM     = r_OFDM_noNoise+noise;

        %% Demodulate FBMC signal
        y_FBMC_Aux         = Q_FBMC'*r_FBMC_Aux; % Same as "FBMC.Demodulation(r_FBMC_Aux)"
        y_FBMC_Cod         = Q_FBMC'*r_FBMC_Cod; % Same as "FBMC.Demodulation(r_FBMC_Cod)"
        y_FBMC_Cod_PostCod = CodingMethod.PrecodingMatrix'*y_FBMC_Cod;
        y_OFDM             = Q_OFDM'*r_OFDM; % Same as "OFDM.Demodulation(r_OFDM)"

        %% Channel Estimation at Pilot Position
        hP_est_FBMC_Aux = y_FBMC_Aux(PilotMatrix_FBMC==1)./xP_FBMC/sqrt(Kappa_Aux);
        hP_est_FBMC_Cod = y_FBMC_Cod(PilotMatrix_FBMC==1)./xP_FBMC/sqrt(Kappa_Cod);
        hP_est_OFDM     = y_OFDM(PilotMatrix_OFDM==1)./xP_OFDM/sqrt(Kappa_OFDM);

        %% Estimate Transmit Matix
        D_FBMC_est_Aux = sum(bsxfun(@times,...
            reshape(full(W_MMSE_FBMC_Aux(:,i_SNR)),size(G_FBMC,2),size(G_FBMC,2),NrPilotSymbols),...
            reshape(hP_est_FBMC_Aux,1,1,[])),3);
        D_FBMC_est_Cod = sum(bsxfun(@times,...
            reshape(full(W_MMSE_FBMC_Cod(:,i_SNR)),size(G_FBMC,2),size(G_FBMC,2),NrPilotSymbols),...
            reshape(hP_est_FBMC_Cod,1,1,[])),3);
        D_OFDM_est = sum(bsxfun(@times,...
            reshape(full(W_MMSE_OFDM(:,i_SNR)),size(G_OFDM,2),size(G_OFDM,2),NrPilotSymbols),...
            reshape(hP_est_OFDM,1,1,[])),3);

        %% One-Tap Equalizer
        h_est_FBMC_Aux = diag(D_FBMC_est_Aux);
        x_est_OneTapEqualizer_FBMC_Aux = y_FBMC_Aux./h_est_FBMC_Aux;
        xD_est_OneTapEqualizer_FBMC_Aux = real(x_est_OneTapEqualizer_FBMC_Aux(AuxilaryPilotMatrix_FBMC(:)==0)./sqrt(AuxiliaryMethod.DataPowerReduction));
        DetectedBitStream_OneTapEqualizer_FBMC_Aux = PAM.Symbol2Bit(xD_est_OneTapEqualizer_FBMC_Aux);
        BER_FBMC_Aux_OneTapEqualizer(i_SNR,i_rep) = mean(BinaryDataStream_FBMC_Aux~=DetectedBitStream_OneTapEqualizer_FBMC_Aux);
        BER_FBMC_Aux_OneTapEqualizer_NoEdge(i_SNR,i_rep) = mean(BinaryDataStream_FBMC_Aux(ConsideredBits_FBMC_Aux)~=DetectedBitStream_OneTapEqualizer_FBMC_Aux(ConsideredBits_FBMC_Aux));

        h_est_FBMC_Cod = diag(D_FBMC_est_Cod);
        x_est_OneTapEqualizer_FBMC_Cod = CodingMethod.PrecodingMatrix'*(y_FBMC_Cod./h_est_FBMC_Cod);
        xD_est_OneTapEqualizer_FBMC_Cod = real(x_est_OneTapEqualizer_FBMC_Cod(NrPilotSymbols+1:end))/CodingMethod.DataPowerReduction;
        DetectedBitStream_OneTapEqualizer_FBMC_Cod = PAM.Symbol2Bit(xD_est_OneTapEqualizer_FBMC_Cod);
        BER_FBMC_Cod_OneTapEqualizer(i_SNR,i_rep) = mean(BinaryDataStream_FBMC_Cod~=DetectedBitStream_OneTapEqualizer_FBMC_Cod);
        BER_FBMC_Cod_OneTapEqualizer_NoEdge(i_SNR,i_rep) = mean(BinaryDataStream_FBMC_Cod(ConsideredBits_FBMC_Cod)~=DetectedBitStream_OneTapEqualizer_FBMC_Cod(ConsideredBits_FBMC_Cod));

        h_est_OFDM = diag(D_OFDM_est);
        x_est_OneTapEqualizer_OFDM = y_OFDM./h_est_OFDM;
        xD_est_OneTapEqualizer_OFDM = x_est_OneTapEqualizer_OFDM(PilotMatrix_OFDM(:)==0)./sqrt(DataPowerReduction_OFDM);
        DetectedBitStream_OneTapEqualizer_OFDM = QAM.Symbol2Bit(xD_est_OneTapEqualizer_OFDM);
        BER_OFDM_OneTapEqualizer(i_SNR,i_rep) = mean(BinaryDataStream_OFDM~=DetectedBitStream_OneTapEqualizer_OFDM);
        BER_OFDM_OneTapEqualizer_NoEdge(i_SNR,i_rep) = mean(BinaryDataStream_OFDM(ConsideredBits_OFDM)~=DetectedBitStream_OneTapEqualizer_OFDM(ConsideredBits_OFDM));

        %% One-Tap Equalizer, Perfect Channel Knowledge
        x_est_OneTapEqualizer_FBMC_Aux_PerfectCSI = y_FBMC_Aux./h_FBMC;
        xD_est_OneTapEqualizer_FBMC_Aux_PerfectCSI = real(x_est_OneTapEqualizer_FBMC_Aux_PerfectCSI(AuxilaryPilotMatrix_FBMC(:)==0)./sqrt(AuxiliaryMethod.DataPowerReduction));
        DetectedBitStream_OneTapEqualizer_FBMC_Aux_PerfectCSI = PAM.Symbol2Bit(xD_est_OneTapEqualizer_FBMC_Aux_PerfectCSI);
        BER_FBMC_Aux_OneTapEqualizer_PerfectCSI(i_SNR,i_rep) = mean(BinaryDataStream_FBMC_Aux~=DetectedBitStream_OneTapEqualizer_FBMC_Aux_PerfectCSI);
        BER_FBMC_Aux_OneTapEqualizer_PerfectCSI_NoEdge(i_SNR,i_rep) = mean(BinaryDataStream_FBMC_Aux(ConsideredBits_FBMC_Aux)~=DetectedBitStream_OneTapEqualizer_FBMC_Aux_PerfectCSI(ConsideredBits_FBMC_Aux));

        x_est_OneTapEqualizer_FBMC_Cod_PerfectCSI = CodingMethod.PrecodingMatrix'*(y_FBMC_Cod./h_FBMC);
        xD_est_OneTapEqualizer_FBMC_Cod_PerfectCSI = real(x_est_OneTapEqualizer_FBMC_Cod_PerfectCSI(NrPilotSymbols+1:end))/CodingMethod.DataPowerReduction;
        DetectedBitStream_OneTapEqualizer_FBMC_Cod_PerfectCSI = PAM.Symbol2Bit(xD_est_OneTapEqualizer_FBMC_Cod_PerfectCSI);
        BER_FBMC_Cod_OneTapEqualizer_PerfectCSI(i_SNR,i_rep) = mean(BinaryDataStream_FBMC_Cod~=DetectedBitStream_OneTapEqualizer_FBMC_Cod_PerfectCSI);
        BER_FBMC_Cod_OneTapEqualizer_PerfectCSI_NoEdge(i_SNR,i_rep) = mean(BinaryDataStream_FBMC_Cod(ConsideredBits_FBMC_Cod)~=DetectedBitStream_OneTapEqualizer_FBMC_Cod_PerfectCSI(ConsideredBits_FBMC_Cod));

        x_est_OneTapEqualizer_OFDM_PerfectCSI = y_OFDM./h_OFDM;
        xD_est_OneTapEqualizer_OFDM_PerfectCSI = x_est_OneTapEqualizer_OFDM_PerfectCSI(PilotMatrix_OFDM(:)==0)./sqrt(DataPowerReduction_OFDM);
        DetectedBitStream_OneTapEqualizer_OFDM_PerfectCSI = QAM.Symbol2Bit(xD_est_OneTapEqualizer_OFDM_PerfectCSI);
        BER_OFDM_OneTapEqualizer_PerfectCSI(i_SNR,i_rep) = mean(BinaryDataStream_OFDM~=DetectedBitStream_OneTapEqualizer_OFDM_PerfectCSI);
        BER_OFDM_OneTapEqualizer_PerfectCSI_NoEdge(i_SNR,i_rep) = mean(BinaryDataStream_OFDM(ConsideredBits_OFDM)~=DetectedBitStream_OneTapEqualizer_OFDM_PerfectCSI(ConsideredBits_OFDM));

        %% Improved Channel Estimation and Data Detection
        xD_est_FBMC_Aux_Temp = xD_est_OneTapEqualizer_FBMC_Aux; % initialize with one tap estimates
        xD_est_FBMC_Cod_Temp = xD_est_OneTapEqualizer_FBMC_Cod; % initialize with one tap estimates
        xD_est_OFDM_Temp     = xD_est_OneTapEqualizer_OFDM; % initialize with one tap estimates
        xD_est_FBMC_Aux_PerfectCSI_Temp = xD_est_OneTapEqualizer_FBMC_Aux_PerfectCSI; % initialize with one tap estimates
        xD_est_FBMC_Cod_PerfectCSI_Temp = xD_est_OneTapEqualizer_FBMC_Cod_PerfectCSI; % initialize with one tap estimates
        xD_est_OFDM_PerfectCSI_Temp     = xD_est_OneTapEqualizer_OFDM_PerfectCSI; % initialize with one tap estimates
        D_FBMC_est_Aux_Temp  = D_FBMC_est_Aux;
        D_FBMC_est_Cod_Temp  = D_FBMC_est_Cod;
        D_OFDM_est_Temp      = D_OFDM_est;
        h_est_FBMC_Aux_Temp  = h_est_FBMC_Aux;
        h_est_FBMC_Cod_Temp  = h_est_FBMC_Cod;
        h_est_OFDM_Temp      = h_est_OFDM;
        for i_iteration = 1:NrIterations
            y_FBMC_Aux_InterferenceCancellation = (y_FBMC_Aux(:) - (D_FBMC_est_Aux_Temp-diag(h_est_FBMC_Aux_Temp))*AuxiliaryMethod.PrecodingMatrix*[xP_FBMC;PAM.SymbolQuantization(xD_est_FBMC_Aux_Temp)]);
            y_FBMC_Cod_InterferenceCancellation = (y_FBMC_Cod - (D_FBMC_est_Cod_Temp-diag(h_est_FBMC_Cod_Temp))*CodingMethod.PrecodingMatrix*[xP_FBMC;PAM.SymbolQuantization(xD_est_FBMC_Cod_Temp)]);
            y_OFDM_InterferenceCancellation     = (y_OFDM - (D_OFDM_est_Temp-diag(h_est_OFDM_Temp))*PilotMapping_OFDM*[xP_OFDM;QAM.SymbolQuantization(xD_est_OFDM_Temp)]);

            % New Channel Estimates at Pilot Positions
            hP_est_FBMC_Aux_Temp = y_FBMC_Aux_InterferenceCancellation(PilotMatrix_FBMC==1)./xP_FBMC/sqrt(Kappa_Aux);
            hP_est_FBMC_Cod_Temp = y_FBMC_Cod_InterferenceCancellation(PilotMatrix_FBMC==1)./xP_FBMC/sqrt(Kappa_Cod);
            hP_est_OFDM_Temp     = y_OFDM_InterferenceCancellation(PilotMatrix_OFDM==1)./xP_OFDM/sqrt(Kappa_OFDM);

            % Improved Channel Estimation
            if i_iteration<=NrIterations/2
                D_FBMC_est_Aux_Temp = sum(bsxfun(@times,...
                    reshape(full(W_MMSE_FBMC_Aux(:,i_SNR)),size(G_FBMC,2),size(G_FBMC,2),NrPilotSymbols),...
                    reshape(hP_est_FBMC_Aux_Temp,1,1,[])),3);
                D_FBMC_est_Cod_Temp = sum(bsxfun(@times,...
                    reshape(full(W_MMSE_FBMC_Cod(:,i_SNR)),size(G_FBMC,2),size(G_FBMC,2),NrPilotSymbols),...
                    reshape(hP_est_FBMC_Cod_Temp,1,1,[])),3);
                D_OFDM_est_Temp = sum(bsxfun(@times,...
                    reshape(full(W_MMSE_OFDM(:,i_SNR)),size(G_OFDM,2),size(G_OFDM,2),NrPilotSymbols),...
                    reshape(hP_est_OFDM_Temp,1,1,[])),3);
            else
                D_FBMC_est_Aux_Temp = sum(bsxfun(@times,...
                    reshape(full(W_MMSE_noInterference_FBMC_Aux(:,i_SNR)),size(G_FBMC,2),size(G_FBMC,2),NrPilotSymbols),...
                    reshape(hP_est_FBMC_Aux_Temp,1,1,[])),3);
                D_FBMC_est_Cod_Temp = sum(bsxfun(@times,...
                    reshape(full(W_MMSE_noInterference_FBMC_Cod(:,i_SNR)),size(G_FBMC,2),size(G_FBMC,2),NrPilotSymbols),...
                    reshape(hP_est_FBMC_Cod_Temp,1,1,[])),3);
                D_OFDM_est_Temp = sum(bsxfun(@times,...
                    reshape(full(W_MMSE_noInterference_OFDM(:,i_SNR)),size(G_OFDM,2),size(G_OFDM,2),NrPilotSymbols),...
                    reshape(hP_est_OFDM_Temp,1,1,[])),3);
            end

            % One-Tap Channel
            h_est_FBMC_Aux_Temp = diag(D_FBMC_est_Aux_Temp);
            h_est_FBMC_Cod_Temp = diag(D_FBMC_est_Cod_Temp);
            h_est_OFDM_Temp     = diag(D_OFDM_est_Temp);

            x_est_FBMC_Aux_Temp = y_FBMC_Aux_InterferenceCancellation(:)./h_est_FBMC_Aux_Temp;
            x_est_FBMC_Cod_Temp = CodingMethod.PrecodingMatrix'*(y_FBMC_Cod_InterferenceCancellation./h_est_FBMC_Cod_Temp);
            x_est_OFDM_Temp     =  y_OFDM_InterferenceCancellation./h_est_OFDM_Temp;

            xD_est_FBMC_Aux_Temp = real(x_est_FBMC_Aux_Temp(AuxilaryPilotMatrix_FBMC(:)==0)./sqrt(AuxiliaryMethod.DataPowerReduction));
            xD_est_FBMC_Cod_Temp = real(x_est_FBMC_Cod_Temp(NrPilotSymbols+1:end))/CodingMethod.DataPowerReduction;
            xD_est_OFDM_Temp     = x_est_OFDM_Temp(PilotMatrix_OFDM(:)==0)./sqrt(DataPowerReduction_OFDM);

            DetectedBitStream_FBMC_Aux_Temp = PAM.Symbol2Bit(xD_est_FBMC_Aux_Temp);
            DetectedBitStream_FBMC_Cod_Temp = PAM.Symbol2Bit(xD_est_FBMC_Cod_Temp);
            DetectedBitStream_OFDM_Temp     = QAM.Symbol2Bit(xD_est_OFDM_Temp);

            BER_FBMC_Aux_InterferenceCancellation(i_SNR,i_rep,i_iteration) = mean(BinaryDataStream_FBMC_Aux~=DetectedBitStream_FBMC_Aux_Temp);
            BER_FBMC_Cod_InterferenceCancellation(i_SNR,i_rep,i_iteration) = mean(BinaryDataStream_FBMC_Cod~=DetectedBitStream_FBMC_Cod_Temp);
            BER_OFDM_InterferenceCancellation(i_SNR,i_rep,i_iteration)     = mean(BinaryDataStream_OFDM~=DetectedBitStream_OFDM_Temp);

            BER_FBMC_Aux_InterferenceCancellation_NoEdge(i_SNR,i_rep,i_iteration) = mean(BinaryDataStream_FBMC_Aux(ConsideredBits_FBMC_Aux)~=DetectedBitStream_FBMC_Aux_Temp(ConsideredBits_FBMC_Aux));
            BER_FBMC_Cod_InterferenceCancellation_NoEdge(i_SNR,i_rep,i_iteration) = mean(BinaryDataStream_FBMC_Cod(ConsideredBits_FBMC_Cod)~=DetectedBitStream_FBMC_Cod_Temp(ConsideredBits_FBMC_Cod));
            BER_OFDM_InterferenceCancellation_NoEdge(i_SNR,i_rep,i_iteration)     = mean(BinaryDataStream_OFDM(ConsideredBits_OFDM)~=DetectedBitStream_OFDM_Temp(ConsideredBits_OFDM));

            % Perfect Channel Knoweledge
            y_FBMC_Aux_InterferenceCancellation_PerfectCSI = (y_FBMC_Aux(:) - (D_FBMC-diag(h_FBMC))*AuxiliaryMethod.PrecodingMatrix*[xP_FBMC;PAM.SymbolQuantization(xD_est_FBMC_Aux_PerfectCSI_Temp)]);
            y_FBMC_Cod_InterferenceCancellation_PerfectCSI = (y_FBMC_Cod - (D_FBMC-diag(h_FBMC))*CodingMethod.PrecodingMatrix*[xP_FBMC;PAM.SymbolQuantization(xD_est_FBMC_Cod_PerfectCSI_Temp)]);
            y_OFDM_InterferenceCancellation_PerfectCSI     = (y_OFDM - (D_OFDM-diag(h_OFDM))*PilotMapping_OFDM*[xP_OFDM;QAM.SymbolQuantization(xD_est_OFDM_PerfectCSI_Temp)]);

            x_est_FBMC_Aux_PerfectCSI_Temp = y_FBMC_Aux_InterferenceCancellation_PerfectCSI./h_FBMC;
            xD_est_FBMC_Aux_PerfectCSI_Temp = real(x_est_FBMC_Aux_PerfectCSI_Temp(AuxilaryPilotMatrix_FBMC(:)==0)./sqrt(AuxiliaryMethod.DataPowerReduction));
            DetectedBitStream_FBMC_Aux_PerfectCSI_Temp = PAM.Symbol2Bit(xD_est_FBMC_Aux_PerfectCSI_Temp);
            BER_FBMC_Aux_PerfectCSI_InterferenceCancellation(i_SNR,i_rep,i_iteration) = mean(BinaryDataStream_FBMC_Aux~=DetectedBitStream_FBMC_Aux_PerfectCSI_Temp);
            BER_FBMC_Aux_PerfectCSI_InterferenceCancellation_NoEdge(i_SNR,i_rep,i_iteration) = mean(BinaryDataStream_FBMC_Aux(ConsideredBits_FBMC_Aux)~=DetectedBitStream_FBMC_Aux_PerfectCSI_Temp(ConsideredBits_FBMC_Aux));

            x_est_FBMC_Cod_PerfectCSI_Temp = CodingMethod.PrecodingMatrix'*(y_FBMC_Cod_InterferenceCancellation_PerfectCSI./h_FBMC);
            xD_est_FBMC_Cod_PerfectCSI_Temp = real(x_est_FBMC_Cod_PerfectCSI_Temp(NrPilotSymbols+1:end))/CodingMethod.DataPowerReduction;
            DetectedBitStream_FBMC_Cod_PerfectCSI_Temp = PAM.Symbol2Bit(xD_est_FBMC_Cod_PerfectCSI_Temp);
            BER_FBMC_Cod_PerfectCSI_InterferenceCancellation(i_SNR,i_rep,i_iteration) = mean(BinaryDataStream_FBMC_Cod~=DetectedBitStream_FBMC_Cod_PerfectCSI_Temp);
            BER_FBMC_Cod_PerfectCSI_InterferenceCancellation_NoEdge(i_SNR,i_rep,i_iteration) = mean(BinaryDataStream_FBMC_Cod(ConsideredBits_FBMC_Cod)~=DetectedBitStream_FBMC_Cod_PerfectCSI_Temp(ConsideredBits_FBMC_Cod));

            x_est_OFDM_PerfectCSI_Temp = y_OFDM_InterferenceCancellation_PerfectCSI./h_OFDM;
            xD_est_OFDM_PerfectCSI_Temp = x_est_OFDM_PerfectCSI_Temp(PilotMatrix_OFDM(:)==0)./sqrt(DataPowerReduction_OFDM);
            DetectedBitStream_OFDM_PerfectCSI_Temp = QAM.Symbol2Bit(xD_est_OFDM_PerfectCSI_Temp);
            BER_OFDM_PerfectCSI_InterferenceCancellation(i_SNR,i_rep,i_iteration) = mean(BinaryDataStream_OFDM~=DetectedBitStream_OFDM_PerfectCSI_Temp);
            BER_OFDM_PerfectCSI_InterferenceCancellation_NoEdge(i_SNR,i_rep,i_iteration) = mean(BinaryDataStream_OFDM(ConsideredBits_OFDM)~=DetectedBitStream_OFDM_PerfectCSI_Temp(ConsideredBits_OFDM));

        end

    end
    TimeNeededSoFar = toc;
    disp([int2str(i_rep/NrRepetitions*100) '% Completed! Time Left: ' int2str(TimeNeededSoFar/i_rep*(NrRepetitions-i_rep)/60) 'min, corresponding to approx. '  int2str(TimeNeededSoFar/i_rep*(NrRepetitions-i_rep)/3600) 'hour']);

    %% Plot results
    % OFDM
    figure(2);
    Markersize = 4;
    hold off;
    semilogy(M_SNR_dB_morePoints, BitErrorProbability,'Color',[1 1 1]*0.75);
    hold on;
    semilogy(M_SNR_dB, nanmean(BER_OFDM_PerfectCSI_InterferenceCancellation(:,:,end),2),'-x black','Markersize',Markersize);
    semilogy(M_SNR_dB, nanmean(BER_OFDM_InterferenceCancellation(:,:,end),2),'-s magenta','Markersize',Markersize);
    semilogy(M_SNR_dB, nanmean(BER_OFDM_InterferenceCancellation_NoEdge(:,:,end),2),'-o blue','Markersize',Markersize);
    semilogy(M_SNR_dB, nanmean(BER_OFDM_OneTapEqualizer_PerfectCSI,2),'-x','Color',[1 1 0]*0.7,'Markersize',Markersize);
    semilogy(M_SNR_dB, nanmean(BER_OFDM_OneTapEqualizer,2),'-s red','Markersize',Markersize);
    ylim([10^-2 0.5]);
    title(['OFDM, Realization ' int2str(i_rep) '/'  int2str(NrRepetitions)])
    legend({'Doubly-Flat Theory','Cancellation (Perfect CSI)','Cancellation','Cancellation (no Edges)','One-Tap (Perfect CSI)','One-Tap'});
    ylabel('Bit Error Ratio');
    xlabel('Signal-to-Noise Ratio [dB]');

    % FBMC, Auxiliary Method
    figure(3);
    hold off;
    semilogy(M_SNR_dB_morePoints, BitErrorProbability,'Color',[1 1 1]*0.75);
    hold on;
    semilogy(M_SNR_dB, nanmean(BER_FBMC_Aux_PerfectCSI_InterferenceCancellation(:,:,end),2),'-x black','Markersize',Markersize);
    semilogy(M_SNR_dB, nanmean(BER_FBMC_Aux_InterferenceCancellation(:,:,end),2),'-s magenta','Markersize',Markersize);
    semilogy(M_SNR_dB, nanmean(BER_FBMC_Aux_InterferenceCancellation_NoEdge(:,:,end),2),'-o blue','Markersize',Markersize);
    semilogy(M_SNR_dB, nanmean(BER_FBMC_Aux_OneTapEqualizer_PerfectCSI,2),'-x','Color',[1 1 0]*0.7,'Markersize',Markersize);
    semilogy(M_SNR_dB, nanmean(BER_FBMC_Aux_OneTapEqualizer,2),'-s red','Markersize',Markersize);
    ylim([10^-2 0.5]);
    semilogy([PlotIterationStepsSNRdB PlotIterationStepsSNRdB],[10^-2 10^-1],'Color',[1 1 1]*0.5,'Linewidth',1);
    title(['FBMC Auxiliary Symbols, Realization ' int2str(i_rep) '/'  int2str(NrRepetitions)])
    legend({'Doubly-Flat Theory','Cancellation (Perfect CSI)','Cancellation','Cancellation (no Edges)','One-Tap (Perfect CSI)','One-Tap'});
    ylabel('Bit Error Ratio');
    xlabel('Signal-to-Noise Ratio [dB]');

    % FBMC, Data Spreading Method
    figure(4);
    hold off;
    semilogy(M_SNR_dB_morePoints, BitErrorProbability,'Color',[1 1 1]*0.75);
    hold on;
    semilogy(M_SNR_dB, nanmean(BER_FBMC_Cod_PerfectCSI_InterferenceCancellation(:,:,end),2),'-x black','Markersize',Markersize);
    semilogy(M_SNR_dB, nanmean(BER_FBMC_Cod_InterferenceCancellation(:,:,end),2),'-s magenta','Markersize',Markersize);
    semilogy(M_SNR_dB, nanmean(BER_FBMC_Cod_InterferenceCancellation_NoEdge(:,:,end),2),'-o blue','Markersize',Markersize);
    semilogy(M_SNR_dB, nanmean(BER_FBMC_Cod_OneTapEqualizer_PerfectCSI,2),'-x','Color',[1 1 0]*0.7,'Markersize',Markersize);
    semilogy(M_SNR_dB, nanmean(BER_FBMC_Cod_OneTapEqualizer,2),'-s red','Markersize',Markersize);
    ylim([10^-2 0.5]);
    title(['FBMC Data Spreading, Realization ' int2str(i_rep) '/'  int2str(NrRepetitions)])
    legend({'Doubly-Flat Theory','Cancellation (Perfect CSI)','Cancellation','Cancellation (no Edges)','One-Tap (Perfect CSI)','One-Tap'});
    ylabel('Bit Error Ratio');
    xlabel('Signal-to-Noise Ratio [dB]');

    % FBMC, Auxiliary Method, BER over Interation
    figure(5);
    hold off;
    semilogy(0:NrIterations, repmat(BitErrorProbability(find(PlotIterationStepsSNRdB==M_SNR_dB_morePoints)),NrIterations+1,1),'Color',[1 1 1]*0.75);
    hold on;
    Index = find(PlotIterationStepsSNRdB==M_SNR_dB);
    semilogy(0:NrIterations, [nanmean(BER_FBMC_Aux_OneTapEqualizer_PerfectCSI(Index,:),2);squeeze(nanmean(BER_FBMC_Aux_PerfectCSI_InterferenceCancellation(Index,:,:),2))],'-x black','Markersize',Markersize);
    semilogy(0:NrIterations, [nanmean(BER_FBMC_Aux_OneTapEqualizer(Index,:),2);squeeze(nanmean(BER_FBMC_Aux_InterferenceCancellation(Index,:,:),2))],'-s magenta','Markersize',Markersize);
    semilogy(0:NrIterations, [nanmean(BER_FBMC_Aux_OneTapEqualizer_NoEdge(Index,:),2);squeeze(nanmean(BER_FBMC_Aux_InterferenceCancellation_NoEdge(Index,:,:),2))],'-o blue','Markersize',Markersize);
    semilogy(0:NrIterations, repmat(nanmean(BER_FBMC_Aux_OneTapEqualizer_PerfectCSI(Index,:),2),NrIterations+1,1),'-x','Color',[1 1 0]*0.7,'Markersize',Markersize);
    semilogy(0:NrIterations, repmat(nanmean(BER_FBMC_Aux_OneTapEqualizer(Index,:),2),NrIterations+1,1),'-s red','Markersize',Markersize);
    title(['FBMC Auxiliary Symbols, Realization ' int2str(i_rep) '/'  int2str(NrRepetitions)])
    legend({'Doubly-Flat Theory','Cancellation (Perfect CSI)','Cancellation','Cancellation (no Edges)','One-Tap (Perfect CSI)','One-Tap'});
    set(gca, 'XTick',0:NrIterations);
    ylabel('Bit Error Ratio');
    xlabel('Iteration Step i');

    pause(0.01);

end

%% Plot Additional Information
fprintf('=============================\n');
fprintf('========= Data Rate =========\n');
fprintf('OFDM       |%7.2f Mbit/s  | \n', length(BinaryDataStream_OFDM)     / (OFDM.PHY.TimeSpacing*OFDM.Nr.MCSymbols)/1e6   );
fprintf('FBMC, Aux. |%7.2f Mbit/s  | \n', length(BinaryDataStream_FBMC_Aux) / (OFDM.PHY.TimeSpacing*OFDM.Nr.MCSymbols)/1e6   );
fprintf('FBMC, Cod. |%7.2f Mbit/s  | \n', length(BinaryDataStream_FBMC_Cod) / (OFDM.PHY.TimeSpacing*OFDM.Nr.MCSymbols)/1e6   );
fprintf('=============================\n');

% The power is normalized so that the average transmit power is one
fprintf('================================================\n');
fprintf('============== Relative SNR Shift ==============\n');
fprintf('================================================\n');
fprintf('           |    SNR    |  Data SNR  | Pilot SNR |\n');
fprintf('OFDM       |   %2.1fdB   |   %2.1fdB   |   %2.1fdB   |\n', 0, 10*log10(DataPowerReduction_OFDM), 10*log10(Kappa_OFDM))
fprintf('FBMC, Aux. |   %2.1fdB   |   %2.1fdB   |   %2.1fdB   |\n', 0, 10*log10(AuxiliaryMethod.DataPowerReduction), 10*log10(Kappa_Aux/2))
fprintf('FBMC, Cod. |   %2.1fdB   |   %2.1fdB   |   %2.1fdB   |\n', 0, 10*log10(CodingMethod.DataPowerReduction), 10*log10(Kappa_Cod/2))
fprintf('================================================\n');

author: rnissel

repo: Channel-Estimation


text:
function xs=delif(x,t);
% DELIF
% DELIF(x,t) deletes the rows of x for which t=1
% x: NxK matrix, t: Nx1 matrix of 0's and 1's
nt=repmat(1-t,1,rows(x'));
nx=x.*nt;
xs=nx(nx~=0);
col=rows(x');
xs=reshape(xs,rows(xs)/col,col);

author: Sable

repo: mcbench-benchmarks


text:
function result = ball_unit_f1_nd ( func, n )

%*****************************************************************************80
%
%% BALL_UNIT_F1_ND approximates an integral inside a unit ball in ND.
%
%  Integration region:
%
%    Points X(1:N) such that:
%
%      Sum ( X(1:N)**2 ) <= 1.
%
%  Discussion:
%
%    An (N+1)*2**N point 5-th degree formula is used, Stroud number SN:5-6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 November 2004
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Arthur H Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971.
%
%  Parameters:
%
%    Input, external FUNC, the name of the user supplied
%    function which evaluates F at the N-vector X, of the form
%      function value = func ( n, x )
%
%    Input, integer N, the dimension of the space.
%
%    Output, real RESULT, the approximate integral of the function.
%
  u2 = ( 1.0E+00 - 2.0E+00 * sqrt ( 1.0E+00 / ( n + 4 ) ) ) / ( n + 2 );
  u = sqrt ( u2 );
  x(1:n) = - u;

  w = 1.0E+00 / ( ( n + 1 ) * 2^n );

  quad = 0.0E+00;
  ihi = 2^n;

  for i = 1 : ihi

    itemp = i - 1;

    for j = 1 : n

      if ( mod ( itemp, 2 ) == 1 )
        x(j) = - abs ( x(j) );
      else
        x(j) = abs ( x(j) );
      end

      itemp = floor ( itemp / 2 );

    end

    quad = quad + w * feval ( func, n, x );

  end

  temp = sqrt ( n + 4 );

  t = sqrt ( 2.0E+00 * ( n + 1 ) / ( n + 2 ) ) / ( n * temp );

  y = ( 1.0E+00 + 2.0E+00 / ( n * temp ) ) / ( n + 2 );
  v = sqrt ( y - t );
  u = sqrt ( y + ( n - 1 ) * t );

  khi = 2^n;

  for i = 1 : n

    x(1:n) = - v;

    x(i) = - u;

    for k = 1 : khi

      ktemp = k - 1;

      for j = 1 : n

        if ( mod ( ktemp, 2 ) == 1 )
          x(j) = - abs ( x(j) );
        else
          x(j) = abs ( x(j) );
        end

        ktemp = floor ( ktemp / 2 );

      end

      quad = quad + w * feval ( func, n, x );

    end

    x(i) = - v;

  end

  volume = ball_unit_volume_nd ( n );
  result = quad * volume;

  return
end


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  June 2013
% ------------------------------------------------------------------------
%  Code obtained from:
%  https://imatge.upc.edu/web/resources/supervised-evaluation-image-segmentation
% ------------------------------------------------------------------------
% This file is part of the SEISM package presented in:
%    Jordi Pont-Tuset, Ferran Marques,
%    "Measures and Meta-Measures for the Supervised Evaluation of Image Segmentation,"
%    Computer Vision and Pattern Recognition (CVPR), 2013.
% If you use this code, please consider citing the paper.
% ------------------------------------------------------------------------
%
% This function returns the root directory where the package is installed.
% Before running the script 'install.m', this directory should be updated
%
% ------------------------------------------------------------------------

function mcg_root_dir = mcg_root_dir()
   mcg_root_dir = fullfile(fileparts(mfilename('fullpath'))); %pwd(), 'regions');
end

author: s-gupta

repo: rcnn-depth

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


text:
classdef AdaptiveChol < handle
   properties (Constant)
      mex = AdaptiveChol.mexSelector();
   end

   methods (Static)
      function func = mexSelector()
         if isarm()
            func = @AdaptiveCholArmMex;
         else
            func = @AdaptiveCholMex;
         end
      end
   end

   properties
      % the constraint matrix
      A
      w = NaN
      cholTol = 1e-4

      % private
      uid
      lastChol = 0; % 1 = double, 2 = ddouble, 3 = qdouble
   end

   methods (Static)
      function o = loadobj(s)
         s.uid = AdaptiveChol.mex('new', uint64(randi(2^32-1,'uint32')), s.A);
         if ~any(isnan(s.w))
            w = s.w; s.w = NaN;
            s.factorize(w);
         end
         o = s;
      end
   end

   methods
      function o = AdaptiveChol(A, cholTol)
         if nargin <= 1, cholTol = 1e-4; end;

         o.A = A;
         o.cholTol = cholTol;
         if isobject(A), A = A.x; end
         o.uid = AdaptiveChol.mex('new', uint64(randi(2^32-1,'uint32')), A);
      end

      function b = saveobj(a)
         b = a;
         b.uid = [];
      end

      function delete(o)
         if ~isempty(o.uid)
            AdaptiveChol.mex('delete', o.uid);
         end
      end

      function err = factorize(o, w, offset)
         if nargin <= 2, offset = 0.0; end
         if isvector(w), w = diag(sparse(w)); end

         o.w = w;
         err = +Inf;

         okay = AdaptiveChol.mex('factorize', o.uid, double(w), offset);
         o.lastChol = 1;
         if okay, err = o.cholAccuracy(); end
         if err < o.cholTol, return; end

         okay = AdaptiveChol.mex('factorize', o.uid, ddouble.toMex(w), offset);
         o.lastChol = 2;
         if okay, err = o.cholAccuracy(); end
         if err < o.cholTol, return; end

         okay = AdaptiveChol.mex('factorize', o.uid, qdouble.toMex(w), offset);
         o.lastChol = 3;
         if okay
            err = o.cholAccuracy();
         else
            err = +Inf;
         end
      end

      function ls = leverageScore(o, JLDim)
         % Warning: This compute (W A' (AWA')^-1 A)_ii
         % This is not exactly leverageScore unless W is diagonal.

         % tau = A' L^{-1} zeta
         tau = AdaptiveChol.mex('halfProj', o.uid, JLDim);
         if (o.lastChol == 2)
            tau = ddouble(tau);
         elseif (o.lastChol == 3)
            tau = qdouble(tau);
         end

         % ls = avg_j (W tau_j)_i (tau_j)_i
         ls = sum((o.w * tau) .* tau,2) / JLDim;
      end

      function err = cholAccuracy(o)
         err = abs(sum(o.leverageScore(1)) - size(o.A, 1));
      end

      function r = diagonal(o)
         assert(o.lastChol, 'factorize must be called before diagonal.');

         r = AdaptiveChol.mex('diagonal', o.uid);

         if (o.lastChol == 2)
            r = ddouble(r);
         elseif (o.lastChol == 3)
            r = qdouble(r);
         end
      end

      function y = solve(o, b, w, iter)
         assert(logical(o.lastChol), 'factorize must be called before solve.');
         assert(nargin >= 2);
         if nargin <= 2 || isempty(w), w = o.w; end
         if nargin <= 3, iter = 3; end

         className = class(b);
         if issparse(b), b = full(b); end
         if isobject(b), b = b.x; end
         if isobject(w), w = w.x; end

         y = AdaptiveChol.mex('solve', o.uid, b, w, iter);

         if (strcmp(className, 'ddouble')) %#ok<*STISA>
            y = ddouble(y);
         elseif (strcmp(className, 'qdouble'))
            y = qdouble(y);
         end
      end

   end
end

author: opencobra

repo: cobratoolbox