text:
% [INPUT]
% ds = A structure representing the dataset.
% sn = A string representing the serial number of the result file.
% temp = A string representing the full path to the Excel spreadsheet used as template for the result file.
% out = A string representing the full path to the Excel spreadsheet to which the results are written, eventually replacing the previous ones.
% bw = An integer [21,252] representing the dimension of each rolling window (optional, default=252).
% sst = A float (0.0,0.1] representing the statistical significance threshold for the linear Granger-causality test (optional, default=0.05).
% rp = A boolean indicating whether to use robust p-values for the linear Granger-causality test (optional, default=false).
% k = A float (0.00,0.20] representing the Granger-causality threshold for non-causal relationships (optional, default=0.06).
% analyze = A boolean that indicates whether to analyse the results and display plots (optional, default=false).
%
% [OUTPUT]
% result = A structure representing the original dataset inclusive of intermediate and final calculations.
% stopped = A boolean that indicates whether the process has been stopped through user input.

function [result,stopped] = run_connectedness(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('ds',@(x)validateattributes(x,{'struct'},{'nonempty'}));
        ip.addRequired('sn',@(x)validateattributes(x,{'char'},{'nonempty' 'size' [1 NaN]}));
        ip.addRequired('temp',@(x)validateattributes(x,{'char'},{'nonempty' 'size' [1 NaN]}));
        ip.addRequired('out',@(x)validateattributes(x,{'char'},{'nonempty' 'size' [1 NaN]}));
        ip.addOptional('bw',252,@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 21 '<=' 252 'scalar'}));
        ip.addOptional('sst',0.05,@(x)validateattributes(x,{'double'},{'real' 'finite' '>' 0 '<=' 0.1 'scalar'}));
        ip.addOptional('rp',false,@(x)validateattributes(x,{'logical'},{'scalar'}));
        ip.addOptional('k',0.06,@(x)validateattributes(x,{'double'},{'real' 'finite' '>' 0 '<=' 0.20 'scalar'}));
        ip.addOptional('analyze',false,@(x)validateattributes(x,{'logical'},{'scalar'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    ds = validate_dataset(ipr.ds,'Connectedness');
    sn = ipr.sn;
    temp = validate_template(ipr.temp);
    out = validate_output(ipr.out);
    bw = ipr.bw;
    sst = ipr.sst;
    rp = ipr.rp;
    k = ipr.k;
    analyze = ipr.analyze;

    nargoutchk(1,2);

    [result,stopped] = run_connectedness_internal(ds,sn,temp,out,bw,sst,rp,k,analyze);

end

function [result,stopped] = run_connectedness_internal(ds,sn,temp,out,bw,sst,rp,k,analyze)

    result = [];
    stopped = false;
    e = [];

    ds = initialize(ds,sn,bw,sst,rp,k);
    t = ds.T;

    bar = waitbar(0,'Initializing connectedness measures...','CreateCancelBtn',@(src,event)setappdata(gcbf(),'Stop',true));
    setappdata(bar,'Stop',false);
    cleanup = onCleanup(@()delete(bar));

    pause(1);
    waitbar(0,bar,'Calculating connectedness measures...');
    pause(1);

    try

        windows = extract_rolling_windows(ds.Returns,ds.BW);

        futures(1:t) = parallel.FevalFuture;
        futures_max = 0;
        futures_results = cell(t,1);

        for i = 1:t
            futures(i) = parfeval(@main_loop,1,windows{i},ds.SST,ds.RP,ds.GroupDelimiters);
        end

        for i = 1:t
            if (getappdata(bar,'Stop'))
                stopped = true;
                break;
            end

            [future_index,value] = fetchNext(futures);
            futures_results{future_index} = value;

            futures_max = max([future_index futures_max]);
            waitbar((futures_max - 1) / t,bar);

            if (getappdata(bar,'Stop'))
                stopped = true;
                break;
            end
        end

    catch e
    end

    try
        cancel(futures);
    catch
    end

    if (~isempty(e))
        delete(bar);
        rethrow(e);
    end

    if (stopped)
        delete(bar);
        return;
    end

    pause(1);
    waitbar(1,bar,'Finalizing connectedness measures...');
    pause(1);

    try
        ds = finalize(ds,futures_results);
    catch e
        delete(bar);
        rethrow(e);
    end

    pause(1);
    waitbar(1,bar,'Writing connectedness measures...');
    pause(1);

    try
        write_results(ds,temp,out);
        delete(bar);
    catch e
        delete(bar);
        rethrow(e);
    end

    if (analyze)
        analyse_result(ds);
    end

    result = ds;

end

%% PROCESS

function ds = initialize(ds,sn,bw,sst,rp,k)

    n = ds.N;
    t = ds.T;

    ds.Result = 'Connectedness';
    ds.ResultDate = now(); %#ok<TNOW1>
    ds.ResultAnalysis = @(ds)analyse_result(ds);
    ds.ResultSerial = sn;

    ds.BW = bw;
    ds.K = k;
    ds.RP = rp;
    ds.SST = sst;

    if (ds.RP)
        all_label = [' (SST=' num2str(ds.SST) ', K=' num2str(ds.K) ', R)'];
    else
        all_label = [' (SST=' num2str(ds.SST) ', K=' num2str(ds.K) ')'];
    end

    ds.LabelsCentralities = {'Betweenness Centrality' 'Closeness Centrality' 'Degree Centrality' 'Eigenvector Centrality' 'Katz Centrality' 'Clustering Coefficient'};

    ds.LabelsIndicatorsSimple = {'DCI' 'CIO' 'CIOO'};
    ds.LabelsIndicators = {['DCI' all_label] ['CIO' all_label] ['CIOO' all_label]};

    ds.LabelsSheetsSimple = {'Indicators' 'Average Adjacency Matrix' 'Average Centrality Measures'};
    ds.LabelsSheets = {['Indicators' all_label] 'Average Adjacency Matrix' 'Average Centrality Measures'};

    ds.AdjacencyMatrices = cell(t,1);
    ds.BetweennessCentralities = NaN(t,n);
    ds.ClosenessCentralities = NaN(t,n);
    ds.DegreeCentralities = NaN(t,n);
    ds.EigenvectorCentralities = NaN(t,n);
    ds.KatzCentralities = NaN(t,n);
    ds.ClusteringCoefficients = NaN(t,n);
    ds.Degrees = NaN(t,n);
    ds.DegreesIn = NaN(t,n);
    ds.DegreesOut = NaN(t,n);

    ds.Indicators = NaN(t,numel(ds.LabelsIndicators));

    ds.AverageAdjacencyMatrix = NaN(n);
    ds.AverageBetweennessCentralities = NaN(1,n);
    ds.AverageClosenessCentralities = NaN(1,n);
    ds.AverageDegreeCentralities = NaN(1,n);
    ds.AverageEigenvectorCentralities = NaN(1,n);
    ds.AverageKatzCentralities = NaN(1,n);
    ds.AverageClusteringCoefficients = NaN(1,n);
    ds.AverageDegreesIn = NaN(1,n);
    ds.AverageDegreesOut = NaN(1,n);
    ds.AverageDegrees = NaN(1,n);

    if (ds.Groups == 0)
        ds.ComparisonReferences = {'Indicators' 1:2 strcat({'CO-'},ds.LabelsIndicatorsSimple)};
    else
        ds.ComparisonReferences = {'Indicators' 1:3 strcat({'CO-'},ds.LabelsIndicatorsSimple)};
    end

end

function window_results = main_loop(r,sst,rp,gd)

    window_results = struct();

    am = causal_adjacency(r,sst,rp);
    window_results.AdjacencyMatrix = am;

    [dci,cio,cioo] = connectedness_metrics(am,gd);
    window_results.DCI = dci;
    window_results.ConnectionsInOut = cio;
    window_results.ConnectionsInOutOther = cioo;

    [bc,cc,dc,ec,kc,clc,deg,deg_in,deg_out] = network_centralities(am);
    window_results.BetweennessCentralities = bc;
    window_results.ClosenessCentralities = cc;
    window_results.DegreeCentralities = dc;
    window_results.EigenvectorCentralities = ec;
    window_results.KatzCentralities = kc;
    window_results.ClusteringCoefficients = clc;
    window_results.Degrees = deg;
    window_results.DegreesIn = deg_in;
    window_results.DegreesOut = deg_out;

end

function ds = finalize(ds,results)

    t = ds.T;

    for i = 1:t
        result = results{i};

        ds.AdjacencyMatrices{i} = result.AdjacencyMatrix;
        ds.BetweennessCentralities(i,:) = result.BetweennessCentralities;
        ds.ClosenessCentralities(i,:) = result.ClosenessCentralities;
        ds.DegreeCentralities(i,:) = result.DegreeCentralities;
        ds.EigenvectorCentralities(i,:) = result.EigenvectorCentralities;
        ds.KatzCentralities(i,:) = result.KatzCentralities;
        ds.ClusteringCoefficients(i,:) = result.ClusteringCoefficients;
        ds.Degrees(i,:) = result.Degrees;
        ds.DegreesIn(i,:) = result.DegreesIn;
        ds.DegreesOut(i,:) = result.DegreesOut;

        ds.Indicators(i,:) = [result.DCI result.ConnectionsInOut result.ConnectionsInOutOther];
    end

    am = sum(cat(3,ds.AdjacencyMatrices{:}),3) ./ numel(ds.AdjacencyMatrices);
    am_threshold = mean(mean(am));
    am(am < am_threshold) = 0;
    am(am >= am_threshold) = 1;
    ds.AverageAdjacencyMatrix = am;

    [bc,cc,dc,ec,kc,clc,deg,deg_in,deg_out] = network_centralities(am);
    ds.AverageBetweennessCentralities = bc;
    ds.AverageClosenessCentralities = cc;
    ds.AverageDegreeCentralities = dc;
    ds.AverageEigenvectorCentralities = ec;
    ds.AverageKatzCentralities = kc;
    ds.AverageClusteringCoefficients = clc;
    ds.AverageDegrees = deg;
    ds.AverageDegreesIn = deg_in;
    ds.AverageDegreesOut = deg_out;

end

function write_results(ds,temp,out)

    [out_path,~,~] = fileparts(out);

    try
        if (exist(out_path,'dir') ~= 7)
            mkdir(out_path);
        end

        if (exist(out,'file') == 2)
            delete(out);
        end
    catch
        error('A system I/O error occurred while writing the results.');
    end

    copy_result = copyfile(temp,out,'f');

    if (copy_result == 0)
        error('The output file could not be created from the template file.');
    end

    firm_names = ds.FirmNames';

    vars = [ds.DatesStr num2cell(ds.Indicators)];
    labels = [{'Date'} ds.LabelsIndicatorsSimple];
    tab = cell2table(vars,'VariableNames',labels);
    writetable(tab,out,'FileType','spreadsheet','Sheet',ds.LabelsSheetsSimple{1},'WriteRowNames',true);

    vars = [firm_names num2cell(ds.AverageAdjacencyMatrix)];
    labels = {'Firms' ds.FirmNames{:,:}};
    tab = cell2table(vars,'VariableNames',labels);
    writetable(tab,out,'FileType','spreadsheet','Sheet',ds.LabelsSheetsSimple{2},'WriteRowNames',true);

    vars = [firm_names num2cell(ds.AverageBetweennessCentralities') num2cell(ds.AverageClosenessCentralities') num2cell(ds.AverageDegreeCentralities') num2cell(ds.AverageEigenvectorCentralities') num2cell(ds.AverageKatzCentralities') num2cell(ds.AverageClusteringCoefficients')];
    labels = [{'Firms'} strrep(ds.LabelsCentralities,' ','')];
    tab = cell2table(vars,'VariableNames',labels);
    writetable(tab,out,'FileType','spreadsheet','Sheet',ds.LabelsSheetsSimple{3},'WriteRowNames',true);

    worksheets_batch(out,ds.LabelsSheetsSimple,ds.LabelsSheets);

end

%% PLOTTING

function analyse_result(ds)

    safe_plot(@(id)plot_indicators(ds,id));
    safe_plot(@(id)plot_network(ds,id));
    safe_plot(@(id)plot_adjacency_matrix(ds,id));
    safe_plot(@(id)plot_centralities(ds,id));

end

function plot_indicators(ds,id)

    dci = smooth_data(ds.Indicators(:,1));
    cio = smooth_data(ds.Indicators(:,2));
    cioo = smooth_data(ds.Indicators(:,3));

    connections_max = max(max([cio cioo])) * 1.1;

    threshold_indices = dci >= ds.K;
    threshold = NaN(ds.T,1);
    threshold(threshold_indices) = connections_max;

    if (ds.RP)
        label = [' (SST=' num2str(ds.SST) ', K=' num2str(ds.K) ', R)'];
    else
        label = [' (SST=' num2str(ds.SST) ', K=' num2str(ds.K) ')'];
    end

    f = figure('Name','Connectedness Measures > Indicators','Units','normalized','Position',[100 100 0.85 0.85],'Tag',id);

    sub_1 = subplot(2,1,1);
    p1 = plot(sub_1,ds.DatesNum,dci);
    hold on;
        p2 = plot(sub_1,ds.DatesNum,repmat(ds.K,[ds.T 1]),'Color',[1.000 0.400 0.400]);
    hold off;
    set(sub_1,'XLim',[ds.DatesNum(1) ds.DatesNum(end)],'XTickLabelRotation',45);
    set(sub_1,'XGrid','on','YGrid','on');
    legend(sub_1,[p1 p2],'Indicator','Threshold','Location','eastoutside');
    title(sub_1,['DCI' label]);

    sub_2 = subplot(2,1,2);
    a1 = area(sub_2,ds.DatesNum,threshold,'EdgeColor','none','FaceColor',[1.000 0.400 0.400]);
    hold on;
        a2 = area(sub_2,ds.DatesNum,cio,'EdgeColor','none','FaceColor','b');
        if (ds.Groups == 0)
            a3 = area(sub_2,ds.DatesNum,NaN(ds.T,1),'EdgeColor','none','FaceColor',[0.678 0.922 1]);
        else
            a3 = area(sub_2,ds.DatesNum,cioo,'EdgeColor','none','FaceColor',[0.678 0.922 1]);
        end
    hold off;
    set(sub_2,'XLim',[ds.DatesNum(1) ds.DatesNum(end)],'XTickLabelRotation',45,'YLim',[0 connections_max]);
    legend(sub_2,[a2 a3 a1],'CIO','CIOO','Threshold Exceeded','Location','eastoutside');
    title(sub_2,['Connections' label]);

    if (ds.MonthlyTicks)
        date_ticks([sub_1 sub_2],'x','mm/yyyy','KeepLimits','KeepTicks');
    else
        date_ticks([sub_1 sub_2],'x','yyyy','KeepLimits');
    end

    sub_1_position = get(sub_1,'Position');
    sub_2_position = get(sub_2,'Position');
    set(sub_1,'Position',[sub_2_position(1) sub_1_position(2) sub_2_position(3) sub_2_position(4)]);

    figure_title(f,'Indicators');

    maximize_figure(f);

end

function plot_network(ds,id)

    if (ds.Groups == 0)
        group_colors = repmat(lines(1),ds.N,1);
    else
        group_colors = zeros(ds.N,3);
        group_delimiters_len = length(ds.GroupDelimiters);
        group_lines = lines(ds.Groups);

        for i = 1:group_delimiters_len
            group_delimiter = ds.GroupDelimiters(i);

            if (i == 1)
                group_colors(1:group_delimiter,:) = repmat(group_lines(i,:),group_delimiter,1);
            else
                group_delimiter_prev = ds.GroupDelimiters(i-1) + 1;
                group_colors(group_delimiter_prev:group_delimiter,:) = repmat(group_lines(i,:),group_delimiter - group_delimiter_prev + 1,1);
            end

            if (i == group_delimiters_len)
                group_colors(group_delimiter+1:end,:) = repmat(group_lines(i+1,:),ds.N - group_delimiter,1);
            end
        end
    end

    weights = mean(ds.Degrees,1,'omitnan');
    weights = weights ./ mean(weights);
    weights = (weights - min(weights)) ./ (max(weights) - min(weights));
    weights = (weights .* 3.75) + 0.25;

    theta = linspace(0,(2 * pi),(ds.N + 1)).';
    theta(end) = [];
    xy = [cos(theta) sin(theta)];
    [i,j] = find(ds.AverageAdjacencyMatrix);
    [~,order] = sort(max(i,j));
    i = i(order);
    j = j(order);
    x = [xy(i,1) xy(j,1)].';
    y = [xy(i,2) xy(j,2)].';

    f = figure('Name','Connectedness Measures > Network Graph','Units','normalized','Position',[100 100 0.85 0.85],'Tag',id);

    sub = subplot(100,1,10:100);

    hold on;
        for i = 1:size(x,2)
            index = ismember(xy,[x(1,i) y(1,i)],'rows');
            plot(sub,x(:,i),y(:,i),'Color',group_colors(index,:));
        end
    hold off;

    if (ds.Groups == 0)
        hold on;
            for i = 1:size(xy,1)
                line(xy(i,1),xy(i,2),'Color',group_colors(i,:),'LineStyle','none','Marker','.','MarkerSize',(35 + (15 * weights(i))));
            end
        hold off;
    else
        d_inc = ds.GroupDelimiters + 1;

        lines_ref = NaN(ds.Groups,1);
        lines_off = 1;

        hold on;
            for i = 1:size(xy,1)
                group_color = group_colors(i,:);
                line(xy(i,1),xy(i,2),'Color',group_color,'LineStyle','none','Marker','.','MarkerSize',(35 + (15 * weights(i))));

                if ((i == 1) || any(d_inc == i))
                    lines_ref(lines_off) = line(xy(i,1),xy(i,2),'Color',group_color,'LineStyle','none','Marker','.','MarkerSize',35);
                    lines_off = lines_off + 1;
                end
            end
        hold off;

        legend(sub,lines_ref,ds.GroupShortNames,'Units','normalized','Position',[0.725 0.131 0.040 0.076]);
    end

    axis(sub,[-1 1 -1 1]);
    axis('equal','off');

    labels = text((xy(:,1) .* 1.075), (xy(:,2) .* 1.075),ds.FirmNames,'FontSize',10);
    set(labels,{'Rotation'},num2cell(theta * (180 / pi())));

    figure_title(f,'Network Graph');

    maximize_figure(f);

end

function plot_adjacency_matrix(ds,id)

    am = ds.AverageAdjacencyMatrix;
    am(logical(eye(ds.N))) = 0.5;
    am = padarray(am,[1 1],'post');

    off = ds.N + 0.5;

    f = figure('Name','Connectedness Measures > Average Adjacency Matrix','Units','normalized','Position',[100 100 0.85 0.85],'Tag',id);

    pcolor(am);
    colormap([1 1 1; 0.65 0.65 0.65; 0.749 0.862 0.933]);
    axis image;

    ax = gca();
    set(ax,'TickLength',[0 0]);
    set(ax,'XAxisLocation','top','XTick',1.5:off,'XTickLabels',ds.FirmNames,'XTickLabelRotation',45);
    set(ax,'YDir','reverse','YTick',1.5:off,'YTickLabels',ds.FirmNames,'YTickLabelRotation',45);

    figure_title(f,'Average Adjacency Matrix');

    maximize_figure(f);

end

function plot_centralities(ds,id)

    seq = 1:ds.N;

    [bc,order] = sort(ds.AverageBetweennessCentralities);
    bc_names = ds.FirmNames(order);
    [cc,order] = sort(ds.AverageClosenessCentralities);
    cc_names = ds.FirmNames(order);
    [dc,order] = sort(ds.AverageDegreeCentralities);
    dc_names = ds.FirmNames(order);
    [ec,order] = sort(ds.AverageEigenvectorCentralities);
    ec_names = ds.FirmNames(order);
    [kc,order] = sort(ds.AverageKatzCentralities);
    kc_names = ds.FirmNames(order);
    [clc,order] = sort(ds.AverageClusteringCoefficients);
    clc_names = ds.FirmNames(order);

    f = figure('Name','Connectedness Measures > Average Centrality Measures','Units','normalized','Position',[100 100 0.85 0.85],'Tag',id);

    sub_1 = subplot(2,3,1);
    bar(sub_1,seq,bc,'FaceColor',[0.749 0.862 0.933]);
    set(sub_1,'XTickLabel',bc_names);
    title(ds.LabelsCentralities{1});

    sub_2 = subplot(2,3,2);
    bar(sub_2,seq,cc,'FaceColor',[0.749 0.862 0.933]);
    set(sub_2,'XTickLabel',cc_names);
    title(ds.LabelsCentralities{2});

    sub_3 = subplot(2,3,3);
    bar(sub_3,seq,dc,'FaceColor',[0.749 0.862 0.933]);
    set(sub_3,'XTickLabel',dc_names);
    title(ds.LabelsCentralities{3});

    sub_4 = subplot(2,3,4);
    bar(sub_4,seq,ec,'FaceColor',[0.749 0.862 0.933]);
    set(sub_4,'XTickLabel',ec_names);
    title(ds.LabelsCentralities{4});

    sub_5 = subplot(2,3,5);
    bar(sub_5,seq,kc,'FaceColor',[0.749 0.862 0.933]);
    set(sub_5,'XTickLabel',kc_names);
    title(ds.LabelsCentralities{5});

    sub_6 = subplot(2,3,6);
    bar(sub_6,seq,clc,'FaceColor',[0.749 0.862 0.933]);
    set(sub_6,'XTickLabel',clc_names);
    title(ds.LabelsCentralities{6});

    set([sub_1 sub_2 sub_3 sub_4 sub_5 sub_6],'XLim',[0 (ds.N + 1)],'XTick',seq,'XTickLabelRotation',90);
    set([sub_1 sub_2 sub_3 sub_4 sub_5 sub_6],'YGrid','on');

    figure_title(f,'Average Centrality Measures');

    maximize_figure(f);

end

%% VALIDATION

function out = validate_output(out)

    [path,name,extension] = fileparts(out);

    if (~strcmpi(extension,'.xlsx'))
        out = fullfile(path,[name extension '.xlsx']);
    end

end

function temp = validate_template(temp)

    sheets = {'Indicators' 'Average Adjacency Matrix' 'Average Centrality Measures'};
    file_sheets = validate_xls(temp,'T');

    if (~all(ismember(sheets,file_sheets)))
        error(['The template must contain the following sheets: ' sheets{1} sprintf(', %s',sheets{2:end}) '.']);
    end

    worksheets_batch(temp,sheets);

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
function conf = SFS_config()
%SFS_CONFIG returns a struct with the default configuration settings
%
%   Usage: conf = SFS_config
%
%   Output parameters:
%       conf    - struct containing all configuration variables
%
%   SFS_CONFIG() creates the struct conf containing the default
%   configuration values. If you want to create other entries, please set
%   them in your script (e.g. conf.fs = 48000) and pass the conf struct to
%   the desired function as last input (e.g. tapering_window(x0,conf)).
%
%   Please don't edit this function to change your configuration settings!
%
%   See also: SFS_start

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ===================================
nargmin = 0;
nargmax = 0;
narginchk(nargmin,nargmax);

%% ===== Table of Contents ===============================================
%
% - Misc
% - Audio
% - Delayline
% - Sound Field Synthesis (SFS)
%   * Dimensionality
%   * Driving functions
%   * Impulse responses
%   * 2.5D
%   * Tapering
% - Sound Field Simulations
% - Secondary Sources
% - Wave Field Synthesis (WFS)
%   * Pre-equalization
% - Spectral Division Method (SDM)
% - Near-Field Compensated Higher Order Ambisonics (NFC-HOA)
% - Local Sound Field Synthesis (LSFS)
% - Binaural Reproduction
%   * Headphone compensation
%   * SoundScape Renderer
% - Plotting
% - References

%% ===== Misc ============================================================
conf.tmpdir = '/tmp/sfs'; % string
% Debugging level. We are supporting two levels:
%   0 - normal mode
%   1 - debug modus, showing interim results and plots
conf.debug = 0; % 0 or 1
% Show a progress bar in selected loops (for example sound_field_mono). This can
% be useful, if you are using a high number of secondary sources.
conf.showprogress = false; % boolean

%% ===== Audio ===========================================================
% Audio and signal processing settings
%
% Samplingrate
conf.fs = 44100; % / Hz
% Speed of sound
conf.c = 343; % / m/s

%% ===== Delayline =======================================================
% Delaying of time signals. This can be critical as very often the wanted delays
% are given as fractions of samples. This configuration section handles how
% those delays should be handled. Beside choosing the actual delayline filter,
% the signal can also be resampled before delaying. As the default setting,
% integer only delays are used by rounding to the next larger integer delay.
% If you want to use a fractional delayline, a setting with a high accuracy is:
%   conf.delayline.resampling = 'pm';
%   conf.delayline.resamplingfactor = 8;
%   conf.delayline.resamplingorder = 64;
%   conf.delayline.filter = 'lagrange';
%   conf.delayline.filterorder = 9;
% Note, that the necessary interpolation accuracy highly depends on the
% actual use case and parametrisation, compare Winter and Spors (2016).
%
% Delayline filter
%   'integer'       - round to nearest integer delay (default)
%   'zoh'           - round to next larger integer delay
%   'lagrange'      - lagrange interpolator (FIR Filter)
%   'least_squares' - least squares FIR interpolation filter
%   'thiran'        - Thiran's allpass IIR filter
%   'farrow'        - use the Farrow structure (to be implemented)
conf.delayline.filter = 'integer';  % string
% Resample signal
%   'none'   - no resampling (default)
%   'matlab' - use matlab's resample() function
%   'pm'     - use Parks-McClellan-Method to compute resample filter (firpm)
conf.delayline.resampling = 'none'; % / string
% Oversamplingfactor factor >= 1
% This should be in the order of (1/stepsize of fractional delays)
conf.delayline.resamplingfactor = 100; % / 1
% Order of Parks-McClellan resample filter (only for 'pm')
% This results in a filter length of resamplingfactor*resamplingorder
conf.delayline.resamplingorder = 64;
% Order of delayline filter (only for Lagrange, Least-Squares & Thiran)
conf.delayline.filterorder = 0;  % / 1
% Number of parallel filters in Farrow structure
% (only for 'farrow');
conf.delayline.filternumber = 1; % / 1

%% ===== Sound Field Synthesis (SFS) =====================================
% Common sound field synthesis settings
%
% === Dimensionality ===
% Dimensionality of the secondary sources and the sound field synthesis driving
% functions:
% '2D'    - line sources as secondary sources, arranged in a circle, line, ...
% '2.5D'  - point sources as secondary sources, arranged in a circle, line, ...
% '3D'    - point sources as secondary sources, arranged in a sphere, plane, ...
conf.dimension = '2.5D'; % string
%
% === Driving functions ===
% Implementation of driving functions. For the default ones use 'default'. These
% functions are described at https://sfs.readthedocs.io/. For possible other
% flags have a look into the driving functions as they can be quite specific.
% Most users can safely use the 'default' flag here.
conf.driving_functions = 'default'; % string
%
% === Impulse responses ===
% Length of impulse responses used in the time domain driving functions
% and for the creation of the binaural simulations.
% Don't worry, SFS checks for you if conf.N is large enough.
conf.N = 2048; % samples
%
% === 2.5D ===
% The amplitude will be correct at the point xref for 2.5D
% synthesis.
conf.xref = [0 0 0]; % / m
%
% === Tapering ===
% The truncation of the loudspeaker array leads to diffraction of the
% synthesized sound field. It has been shown that the truncation can be
% described by cylindrical waves originating from the edges of the array, see
% Sect. 8.3.2 in Born, Wolf (1999) for the general principle and Sect. 3.2 in
% Wierstorf (2014) for how it relates to WFS. Therefore a good method to reduce
% artifacts due to the diffraction edge waves is to fade out the amplitude of
% the driving function at the edges of the array. This method is called tapering
% and implemented using a Hann window in the SFS Toolbox.
conf.usetapwin = true; % boolean
% Size of the tapering window
conf.tapwinlen = 0.3; % / percent of array length, 0..1
%
% === Time Domain Implementation ===
% Adjust the starting time in time domain driving functions.
% This can be set to
%   'system'   - the first secondary source will be active at t=0
%   'source'   - the virtual source will be active at t=0
% Setting it to 'system' is most convenient when simulating single sources as
% you will always see activity in the sound field for t>0. Setting it to
% 'source' helps you to simulate different sources as you can time align them
% easily. Note, that for virtual sources outside of the array this can mean you
% will see no activity inside the listening area until the time has passed, that
% the virtual source needs from its position until the nearest secondary source.
% Also note, using 'source' for systems with unbounded listening areas -- e.g.
% linear arrays -- focussed virtual sources may not be placed arbitrarily
% far from the secondary sources.
conf.t0 = 'system'; % string
% Bandpass filter applied in sound_field_imp()
conf.usebandpass = false; % boolean
conf.bandpassflow = 10; % / Hz
conf.bandpassfhigh = 20000; % / Hz
%
% === Modal Weighting ===
% Additional weighting of the modal coefficients by window function used
% for Near-Field Compensated Higher Order Ambisonics (NFC-HOA) and Local Sound
% Field Synthesis using Spatial Bandwidth Limitation (LSFS-SBL).
%
% Window type. Available windows are:
%   'rect'                     - all coefficients are weighted by 1.0
%   'max-rE'                   - 2D max-rE weighting
%   'kaiser', 'kaiser-bessel'  - Kaiser aka. Kaiser-Bessel window
%   'tukey'                    - modified Tukey (tapered cosine) window
conf.modal_window = 'rect';  % string
% Scalar parameter for window, if applicable. Effect for distinct window:
%   'rect'    - no effect
%   'max-rE'  - no effect
%   'kaiser'  - [0,inf]. trade-off between main-lobe width and side-lobe levels.
%               0.0 results in the rectangular window and the smallest main-lobe
%               width. infinity results in a dirac impulse.
%   'tukey'   - [0,1]. width of cosine tapering relative to the modal order.
%               0.0 results in the rectangular window, 1.0 results in a modified
%               Hann-window
conf.modal_window_parameter = 0.0;  % float

%% ===== Sound Field Simulations =========================================
% Simulations of monochromatic or time domain sound field
%
% xyz-resolution for sound field simulations, this value is applied along every
% desired dimension, except if only one point is desired
conf.resolution = 300; % / samples
% Phase of omega of sound field (change this value to create monochromatic sound
% fields with different phases, for example this can be useful to create a movie)
conf.phase = 0; % / rad

% ===== Secondary Sources ================================================
% Settings of the used loudspeaker array
%
% Number of secondary sources
conf.secondary_sources.number = 64; % integer
% Diameter/Length of secondary source array
conf.secondary_sources.size = 3; % / m
% Center of array, X0
conf.secondary_sources.center = [0 0 0]; % / m
% Array geometry
% Possible values are: 'line', 'box', 'rounded-box', 'circle', 'sphere', 'custom'
conf.secondary_sources.geometry = 'circle'; % string
% Exclusive for 'rounded-box' array geometry. Defines the bending radius for
% the corners of the smoothed box
conf.secondary_sources.corner_radius = 0.0; % / m
% Vector containing custom secondary source positions and directions.
% This is used if geometry = 'custom' is specified.
% conf.secondary_sources.x0 = [x0; y0; z0; nx0; ny0; nz0; weight];
% Or it could also be a SOFA struct or file name, in this case the positions are
% extracted from the provided SOFA file.
conf.secondary_sources.x0 = []; % / m
% For spherical arrays, available grids are:
%   'equally_spaced_points' - Sphere with equal distance between grid points
%   'gauss'                 - Gauss grid
%   'fabian'                - grid of 3D HRTF measurement, available at
%                             https://doi.org/10.14279/depositonce-5718
%
% Note, that 'equally_spaced_points' and 'fabian' are precomputed grids that
% will be automatically downloaded and cached on your disk. All available number
% of secondary sources for those grids can be seen at:
% https://github.com/sfstoolbox/data/tree/master/spherical_grids
% An exception are Gauss grids, which are available via 'gauss' and will be
% calculated on the fly allowing very high number of secondary sources.
conf.secondary_sources.grid = 'equally_spaced_points'; % string
% Grid for linear arrays: Allows for logarithmically spaced arrays. It defines
% the ratio of the spacing of the loudspeakers at the extremities of the array
% and the spacing of the loudspeakers at the center of the array. 1.0 for
% equi-distant sampling.
conf.secondary_sources.logspread = 1.0; % float

%% ===== Wave Field Synthesis (WFS) ======================================
% Settings for WFS, see Spors et al. (2008) for an introduction
%
% === Pre-Equalization ===
% WFS can be implemented very efficiently using a delay-line with different
% amplitudes and convolving the whole signal once with the so called
% pre-equalization filter, see Spors and Ahrens (2010a).
% Use WFS preequalization-filter
conf.wfs.usehpre = true; % boolean
% FIR or IIR pre-equalization filter
% NOTE: only FIR is working under octave at the moment
conf.wfs.hpretype = 'FIR'; % 'FIR' or 'IIR'
% Lower frequency limit of preequalization filter
% If we have a finite length (<10m) of the secondary source distribution we
% will have a 3dB increase at the very low frequencies and the pre-equalization
% filter should only start above those frequencies, see Sect. 7.2 in Spors and
% Ahrens (2010b).
conf.wfs.hpreflow = 50; % / Hz
% Upper frequency limit of preequalization filter
% If we have aliasing in our system we only want to use
% the pre-equalization filter until the aliasing frequency, because of the
% energy the aliasing is adding to the spectrum above this frequency, see
% Sect. 6.3 in Spors and Ahrens (2010a).
conf.wfs.hprefhigh = 1200; % / Hz
% IIR bandwidth for the Lagrange interpolation region
conf.wfs.hpreBandwidth_in_Oct = 2; % / octaves
% desired IIR filter order
conf.wfs.hpreIIRorder = 4; % integer
% desired FIR filter order, results in N+1 taps
conf.wfs.hpreFIRorder = 128; % even integer

%% ===== Spectral Division Method (SDM) ==================================
% Settings for SDM, see Ahrens and Spors (2010) for an introduction
%
% Use the evanescent part of the driving function for SDM
conf.sdm.withev = true; % boolean

%% ===== Near-Field Compensated Higher Order Ambisonics (NFC-HOA) ========
% Settings for NFC-HOA, see Ahrens (2012) for an introduction
%
% Highest order used with NFC-HOA. If this is set to [], band-limited NFC-HOA is
% used and the order is set by nfchoa_order() which returns the highest order
% for which no aliasing occurs.
conf.nfchoa.order = []; % integer

%% ===== Local Wave Field Synthesis (LWFS) ===============================
% Settings for Local WFS
%
% === Local Wave Field Synthesis using Virtual Secondary Sources (LWFS-VSS)
% See Spors and Ahrens (2010b) for an introduction.
%
% Method the virtual secondary sources should be driven
conf.localwfs_vss.method = 'wfs'; % 'wfs' or 'nfchoa'
% WFS settings for virtual secondary sources
conf.localwfs_vss.wfs = conf.wfs;
% Tapering of virtual secondary sources (only applied for WFS)
conf.localwfs_vss.usetapwin = false; % boolean
conf.localwfs_vss.tapwinlen = 0.5; % 0..1
% NFC-HOA settings for virtual secondary sources
conf.localwfs_vss.nfchoa = conf.nfchoa;
% Virtual secondary sources (see also: conf.secondary_sources)
conf.localwfs_vss.size = 0.4; % / m
conf.localwfs_vss.center = [0, 0, 0]; % / m
conf.localwfs_vss.geometry = 'circular'; % string
conf.localwfs_vss.number = 56; % integer
conf.localwfs_vss.grid = 'equally_spaced_points'; % string
% Driving function for virtual secondary sources
conf.localwfs_vss.driving_functions = 'default'; % string
% Linear VSS distribution: rotate the distribution orthogonal to the progation
% direction of the desired sound source.
% Circular VSS distribution: truncate the distribution to a circular arc
% which satisfies the secondary source selection criterions (source normal
% aligns with propagation directions of desired sound source).
conf.localwfs_vss.consider_target_field = true; % boolean
% VSS distribution is further truncated if parts of it cannot be correctly
% reproduced, because they lie outside the area which is surrounded by the real
% loudspeakers (secondary sources)
conf.localwfs_vss.consider_secondary_sources = true; % boolean
%
% === Local Wave Field Synthesis using Spatial Bandwidth Limitation (LWFS-SBL)
% See Hahn et al. (2016) for an introduction.
% The centre of the local synthesis region is set by conf.xref
%
% Maximum modal order aka. spatial bandwidth of desired sound field. If left
% empty, the value is set by nfchoa_order(), which may suboptimal depending on
% the geometry, e.g. number of secondary sources and shape of the secondary
% source distribution.
conf.localwfs_sbl.order = []; % integer
% Due to stability issues for the time-domain implementation of synthesis
% of a point source, conventional WFS has to be used for the low frequencies.
% fc defines the crossover frequency between the WFS and LWFS-SBL. If left
% empty, this frequency is estimated by aliasing_frequency().
conf.localwfs_sbl.fc = []; % float
% The spatially bandwidth-limited sound field is converted into plane wave
% decomposition which is then synthesised using conventional WFS for each
% individual plane wave. Npw defines the number of plane waves with their
% directions distributed equi-angularly on the unit circle. If left empty,
% it is estimated based on the sampling frequency and size of the secondary
% source distribution.
conf.localwfs_sbl.Npw = []; % integer

%% ===== Binaural reproduction ===========================================
% Settings regarding all the stuff with impulse responses from the SFS_ir and
% SFS_binaural_synthesis folders
%
% Use interpolation to get the desired HRTF or BRIR for binaural simulation.
% If this is disabled, the HRTF/BRIR returned by a nearest neighbour search is
% used instead.
conf.ir.useinterpolation = true; % boolean
% You can choose the way the points for interpolation are selected. Depending on
% the geometry of the measured HRTF/BRIR data set and the interpolation method,
% a different number of HRTFs are selected for interpolation. See validation
% script test_interpolation_point_selection.m for examples. Available methods:
%   'nearestneighbour'  - Interpolation between nearest neighbours. This only
%                         works for interpolation points on a circle in the
%                         horizontal plane.
%   'delaunay'          - Interpolation between surrounding points according to
%                         Delaunay triangulation. This only works for
%                         interpolation points on a sphere.
%   'voronoi'           - For spherical Voronoi interpolation, all coordinates
%                         are projected on the unit sphere. Voronoi regions
%                         (including their area) on the sphere are calculated
%                         for all measured positions included in the HRTF/BRIR
%                         data set. The query position is added to the point
%                         cloud and the regions are calculated, again.
%                         The weights result from the area "stolen" by the query
%                         position from each of the other coordinates
conf.ir.interpolationpointselection = 'nearestneighbour';
% You can choose between the following interpolation methods:
%   'simple'      - Interpolation in the time domain performed samplewise. This
%                   does not heed the times of arrival of the impulse responses.
%   'freqdomain'  - Interpolation in the frequency domain performed separately
%                   for magnitude and phase.
%                   This method cannot work properly if there is too much noise in
%                   the phase information at low frequencies which is often the
%                   case for measured HRTFs. Low frequencies can be corrected
%                   according to theory, see e.g. the corrected KEMAR HRTFs published
%                   at https://github.com/spatialaudio/lf-corrected-kemar-hrtfs.
%                   The implementation of this method suffers from circular shifting,
%                   see test_interpolation_methods.m in the validation folder. For
%                   typical HRIRs with leading and trailing zeros, the error is
%                   negligible.
%   'timedomain'  - Interpolation in the time domain with cross-correlation for
%                   estimation of time of arrival (TOA) differences.
conf.ir.interpolationmethod = 'simple';
%
% If you have HRIRs in the form of the SimpleFreeFieldHRIR SOFA convention, zeros
% are padded at the beginning of every impulse response corresponding to their
% measurement distance. If you know that your measured HRIRs already have a
% given pre-delay, add the pre-delay here and accordingly less zero padding will
% be applied. In this case you can lose samples from the beginning of the
% impulse response. If you are not sure, choose a value of 0.
conf.ir.hrirpredelay = 0; % / samples
%
% === Headphone compensation ===
% Headphone compensation
conf.ir.usehcomp = false; % boolean
% Headphone compensation file for left and right ear.
conf.ir.hcompfile = 'data/headphone_compensation/QU_KEMAR_AKGK601_hcomp.wav'; % string
%
% === SoundScape Renderer ===
% To use a dynamic binaural simulation together with the SoundScape Renderer
% (SSR) and a headtracker, brs sets can be created. If these sets should be
% used in BRS mode of the SSR, the angles have to be:
% conf.ir.brsangles = 0:1:359;
% If the brs set should be used as IRs for the SSR, the angles have to be:
% conf.ir.brsangles = 360:-1:1;
conf.ir.brsangles = 0:1:359; % / degree

%% ===== Plotting ========================================================
% Plot the results (sound fields etc.) directly
conf.plot.useplot = false; % boolean
% Normalize the sound field for plotting
conf.plot.usenormalisation = true; % boolean
% Normalisation method. Available methods are:
%   'auto'      - 'center' if center of sound field > 0.3, otherwise 'max'
%   'center'    - center of sound field == 1
%   'max'       - max of sound field == 1
conf.plot.normalisation = 'auto'; % string
% Plot mode (uses the GraphDefaults function). Available modes are:
%   'monitor'   - displays the plot on the monitor
%   'paper'     - eps output in conf.plot.outfile
%   'png'       - png output in conf.plot.outfile
conf.plot.mode = 'monitor'; % string
% Plot amplitudes in dB (e.g. sound field plots)
conf.plot.usedb = false; % boolean
% caxis settings (leave blank, if you would use the default values of the given
% plot function)
conf.plot.caxis = []; % [min max]
% Default colormap to use
% The Toolbox comes with two own color maps, if you set 'default' or 'moreland'
% you will get a blue/red-colormap after
% http://www.kennethmoreland.com/color-advice/
% If you set 'gray' or 'grey' you will get a colormap ranging from white to
% black. In addition you can add every other map you can specify in
% Matlab/Octave. For example to get the Matlab default colormap set 'parula'.
conf.plot.colormap = 'default'; % string
% Plot loudspeakers in the sound field plots
conf.plot.loudspeakers = true; % boolean
% Use real loudspeakers symbols (otherwise crosses are used)
conf.plot.realloudspeakers = false; % boolean
% Size of the loudspeaker
conf.plot.lssize = 0.16; % m
% Size of the plot
conf.plot.size_unit = 'px'; % 'px','cm','inches'
conf.plot.size = [540 404]; % [xsize ysize]
% Resolution of plot in dpi
conf.plot.resolution = 150; % integer
% Additional plot command
conf.plot.cmd = ''; % string
% output of plot (file or screen)
conf.plot.usefile = false; % boolean
% File name, if this is provided with as *.png or *.eps file, the figure is
% plotted to the regarding file
conf.plot.file = ''; % string

%% ===== References ======================================================
%
% Ahrens (2012) - "Analytic Methods of Sound Field Synthesis", Springer,
% https://doi.org/10.1007/978-3-642-25743-8
%
% Ahrens and Spors (2010) - "Sound Field Reproduction Using Planar and Linear
% Arrays of Loudspeakers", Transactions on Audio, Speech, and Language
% Processing, vol. 18, no. 8, pp. 2038-2050,
% https://doi.org/10.1109/TASL.2010.2041106
%
% Born and Wolf (1999) - "Principles of Optics", Cambridge University Press,
% 7th edition, ISBN 9780521642224
%
% Hahn, Winter, Spors (2016) - "Local Wave Field Synthesis by Spatial
% Band-limitation in the Circular/Spherical Harmonics Domain", 140th
% Convention of the Audio Engineering Society, Paper 9596,
% http://www.aes.org/e-lib/browse.cfm?elib=18294
%
% Spors and Ahrens (2010a) - "Analysis and Improvement of Pre-Equalization in
% 2.5-Dimensional Wave Field Synthesis", 128th Convention of the Audio
% Engineering Society, Paper 8121,
% http://www.aes.org/e-lib/browse.cfm?elib=15418
%
% Spors and Ahrens (2010b) - "Local Sound Field Synthesis by Virtual Secondary
% Sources", 40th Conference of the Audio Engineering Society, Paper 6-3,
% http://www.aes.org/e-lib/browse.cfm?elib=15561
%
% Spors, Rabenstein, Ahrens (2008) - "The Theory of Wave Field Synthesis
% Revisited," 124th Convention of the Audio Engineering Society, Paper 7358,
% http://www.aes.org/e-lib/browse.cfm?elib=14488
%
% Wierstorf (2014) - "Perceptual Assessment of Sound Field Synthesis",
% TU Berlin, https://doi.org/10.14279/depositonce-4310
%
% Winter and Spors (2016) - "On fractional delay interpolation for local wave
% field synthesis", 24th European Signal Processing Conference (EUSIPCO),
% pp. 2415-2419, https://doi.org/10.1109/EUSIPCO.2016.7760682

author: sfstoolbox

repo: sfs-matlab


text:
function Rx=x_rot(phi)
%X_ROT Matrix creating a rotation around the x axis by an angle phi.
%   phi: angle in radians.

Rx = [1 0 0; 0 cos(phi) -sin(phi);0 sin(phi) cos(phi)];

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
%sim_dMul.m
%Jamie Near, 2020.
%
% USAGE:
% d_out = sim_dMul(d_in,factor)
%
% DESCRIPTION:
% Multiply a density matrix by a scalar.  This function is necessary becuase
% the density matrix is a cell array, so cannot be scaled using simple
% multiplication.
%
% INPUTS:
% d_in      = input density matrix to be multiplied.
% factor    = scalar multiplcation factor.
%
% OUTPUTS:
% d_out     = output, result of d_in * factor.

function d_out = sim_dMul(d_in,factor)

d_out=cell(size(d_in));
for m=1:length(d_in)
    d_out{m}=d_in{m}.*factor;
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function g = gscale(f, varargin)
%GSCALE Scales the intensity of the input image.
%   G = GSCALE(F, 'full8') scales the intensities of F to the full
%   8-bit intensity range [0, 255].  This is the default if there is
%   only one input argument.
%
%   G = GSCALE(F, 'full16') scales the intensities of F to the full
%   16-bit intensity range [0, 65535].
%
%   G = GSCALE(F, 'minmax', LOW, HIGH) scales the intensities of F to
%   the range [LOW, HIGH]. These values must be provided, and they
%   must be in the range [0, 1], independently of the class of the
%   input. GSCALE performs any necessary scaling. If the input is of
%   class double, and its values are not in the range [0, 1], then
%   GSCALE scales it to this range before processing.
%
%   The class of the output is the same as the class of the input.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/11/21 14:36:09 $

if length(varargin) == 0 % If only one argument it must be f.
   method = 'full8';
else
   method = varargin{1};
end

if strcmp(class(f), 'double') & (max(f(:)) > 1 | min(f(:)) < 0)
   f = mat2gray(f);
end

% Perform the specified scaling.
switch method
case 'full8'
   g = im2uint8(mat2gray(double(f)));
case 'full16'
   g = im2uint16(mat2gray(double(f)));
case 'minmax'
   low = varargin{2}; high = varargin{3};
   if low > 1 | low < 0 | high > 1 | high < 0
      error('Parameters low and high must be in the range [0, 1].')
   end
   if strcmp(class(f), 'double')
      low_in = min(f(:));
      high_in = max(f(:));
   elseif strcmp(class(f), 'uint8')
      low_in = double(min(f(:)))./255;
      high_in = double(max(f(:)))./255;
   elseif strcmp(class(f), 'uint16')
      low_in = double(min(f(:)))./65535;
      high_in = double(max(f(:)))./65535;
   end
   % imadjust automatically matches the class of the input.
   g = imadjust(f, [low_in high_in], [low high]);
otherwise
   error('Unknown method.')
end

repo: weiyanmin


text:
function classifier = ocsvmp(options,data)
% {svmp} trains a One-Class SVM classifier on the primal.
%
%      classifier = ocsvmp(options,data)
%
%      Note: This function is only a wrapper to 'lapsvmp', where the
%            contribute of unlabeled data is discarded. See 'lapsvmp'
%            for the description of all required parameters. All parameters
%            related to unlabeled data are obviously not checked.
%
% Author: Stefano Melacci (2012), Salvatore Frandina (2012)
%         mela@dii.unisi.it, salvatore.frandina@gmail.com

options.Hinge=1;
options.gamma_I=0;
options.UseBias=1;
data.Y(data.Y==-1)=0;
classifier=lapsvmp(options,data);
classifier.name='ocsvmp';

author: viggin

repo: domain-adaptation-toolbox


text:
function nFrames = numFrames(view, scan)
% Number of time samples (frames) for the view/scan combination
%
%   nFrames = numFrames(view, scan)
%
% Accesses dataTYPES to get this info.
%
% global dataTYPES;
%
% Example:
%   numFrames(INPLANE{1},1)
if notDefined('scan'), scan = view.curScan; end
nFrames = viewGet(view,'nFrames',scan);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% most of the time. they both this and pcasvd takes about the same time
% when smallest = 0, this should be faster but it's not always
% when it's 1, due to null(U1'), it prob slows it down but not always
function [U,d] = pcaeig(X,n,smallest)
% X : nt x n, where nt examples of feature vecs of size n
% n returns n dimensions
% smallest = 1 when we want the n evecs assoc with the smallest
% evals
% U is evecs, d is diagonals
% Copyright (c) 2013, Vipin Vijayan.
if size(X,1) < size(X,2)
    [U,D] = eig(X*X');
    d = diag(D);
    [~,ix] = sort(d,'descend');
    U = U(:,ix);
    d = d(ix);
else
    [V,D1] = eig(X'*X);
    d1 = diag(D1);
    [~,ix] = sort(d1,'descend');
    V = V(:,ix);
    d1 = d1(ix);
    U1 = X*V; % go from evecs of X'*X to evecs of X*X'
    for ni = 1:size(U1,2), U1(:,ni) = U1(:,ni)./norm(U1(:,ni)); end
    if ~smallest
        U = U1;
        d = d1;
    else
        U2 = null(U1'); % this runs svd on U1'
        U = [U1 U2];
        d = [d1 ; zeros(size(U2,2),1)];
    end
end
if ~smallest,
    tol = eps(class(X));
    n = min(n, sum(d > tol));
    U = U(:,1:n);
    d = d(1:n);
else
    n = min(n,size(U,2));
    U = U(:,end:-1:end-n+1);
    d = max(0,d);
    d = d(end:-1:end-n+1);
end
% D = diag(d);
end

% very slow with large dimensions
function [U,d] = pcaeigslow(X,n,smallest)
[U,D] = eig(X*X');
d = diag(D);
[~,ix] = sort(d,'descend');
U = U(:,ix);
d = d(ix);
if ~smallest,
    tol = eps(class(X));
    n = min(n, sum(diag(D) > tol));
    U = U(:,1:n);
    d = d(1:n);
else
    n = min(n,size(U,2));
    U = U(:,end:-1:end-n+1);
    d = abs(diag(D));
    d = d(end:-1:end-n+1);
end
% D = diag(d);
end

author: Sable

repo: mcbench-benchmarks


text:
function polys = meshFacePolygons(varargin)
%MESHFACEPOLYGONS Returns the set of polygons that constitutes a mesh.
%
%   POLYGONS = meshFacePolygons(V, F)
%   POLYGONS = meshFacePolygons(MESH)
%
%   Example
%     [v f] = createCubeOctahedron;
%     polygons = meshFacePolygons(v, f);
%     areas = polygonArea3d(polygons);
%     sum(areas)
%     ans =
%         18.9282
%
%   See also
%     meshes3d, meshFace, polygonArea3d

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2013-08-20,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2013 INRA - Cepia Software Platform.

% extract vertices and faces
[v, f] = parseMeshData(varargin{:});

% number of faces
if iscell(f)
    nFaces = length(f);
else
    nFaces = size(f, 1);
end

% allocate cell array for result
polys = cell(nFaces, 1);

% compute polygon corresponding to each face
if iscell(f)
    for i = 1:nFaces
        polys{i} = v(f{i}, :);
    end
else
    for i = 1:nFaces
        polys{i} = v(f(i,:), :);
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function liga(pattern)
%*************************************************
%DATE: 13/8/2007 (created); 	13/8/2007 (modified)
%FUNCTION:
%INPUTS:
%OUTPUTS:
%muwn.gu@gmail.com
%*************************************************

global layerMatrix;

segLengthArray = extractPattern (pattern);

for segIndex = 1: size(segLengthArray, 2)

	segFrom = round ( 96 * segLengthArray(segIndex).from);
	if (segFrom == 0)
		segFrom = 1;
	end
	segEnd	= round ( 96 * segLengthArray(segIndex).end);

	for pointIndex = segFrom:segEnd 			% horzontal sweep

		prIndex = find(layerMatrix(:,pointIndex) == 9); % find the photo resister sublines in current pointIndex
		layerMatrix (prIndex, pointIndex) = 0;

	end
end

function [segLengthArray] = extractPattern(pattern)
%*************************************************
%DATE: 11/8/2007 (created); 	11/8/2007 (modified)
%FUNCTION:	extract segment number and length array from pattern
%INPUTS:	pattern string
%OUTPUTS:	segLengthArray
%*************************************************
%'  -- -- -'
%'1234567890
segIndex = 1;
patternLength = length (pattern);

if (pattern(1, 1) == ' ')
segLengthArray(1).from = 0;
end

for charIndex =  1:patternLength - 1

	if (pattern(1, charIndex) == '-') & (pattern(1, charIndex+1) == ' ')
		segLengthArray(segIndex).from = (charIndex)/(patternLength);
	end

	if (pattern(1, charIndex) == ' ') & (pattern(1, charIndex+1) == '-')
		segLengthArray(segIndex).end = (charIndex)/(patternLength);
		segIndex =  segIndex + 1;
	end
end

	if (pattern(1, patternLength) == ' ')
		segLengthArray(segIndex).end = 1;
	end


text:
function [b bbar sigeps]=panel3prior(Xibar,Xbar,yi,y,N,q)

% first obtain b, the mean group estimator
% estimate the first term
term1=sparse(q,q);
for ii=1:N
term1=term1+Xibar{ii,1}'*Xibar{ii,1};
end

% estimate the second term
term2=sparse(q,1);
for ii=1:N
term2=term2+Xibar{ii,1}'*yi(:,:,ii);
end

% obtain b
b=term1\term2;

% obtain bbar
bbar=repmat(b,N,1);

% obtain sigma_epsilon, the common residual variance term
eps=y-Xbar*bbar;
sigeps=var(eps);

author: european-central-bank

repo: BEAR-toolbox


text:
close all; clear all; clc;

audio = audioload();

file = 'text.txt';
fid  = fopen(file, 'r');
text = fread(fid,'*char')';
fclose(fid);

msg = dsss_dec(audio.data, 8*length(text));
err = BER(text,msg);
nc  = NC(text, msg);

fprintf('Text: %s\n', msg);
fprintf('BER : %d\n', err);
fprintf('NC  : %d\n', nc);

author: ktekeli

repo: audio-steganography-algorithms


text:
function ret = nonlinear(chrom,sizepop)

for i=1:sizepop
    x=fmincon(inline('-5*sin(x(1))*sin(x(2))*sin(x(3))*sin(x(4))*sin(x(5))-sin(5*x(1))*sin(5*x(2))*sin(5*x(3))*sin(5*x(4))*sin(5*x(5))'),chrom(i,:)',[],[],[],[],[0 0 0 0 0],[2.8274 2.8274 2.8274 2.8274 2.8274]);
    ret(i,:)=x';
end

author: HuangCongQing


text:
function range = p06_range ( )

%*****************************************************************************80
%
%% P06_RANGE returns an interval bounding the root for problem 6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real RANGE(2), the minimum and maximum values of
%    an interval containing the root.
%
  range(1) =  0.00001;
  range(2) = 20.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Lmk = lmkJacobians(Lmk)

% LMKJACOBIANS Compute Jacobians for projection onto the manifold.
%   LMKJACOBIANS computes all the Jacobians of the landmarks in structure
%   array Lmk for the projection of the landmark errors into the landmark
%   manifolds. The computed Jacobians are stored in Lmk.state.M. Its
%   expression depends on the landmark type.
%
%   Landmark types are specifyed by Lmk.type. In the case of Euclidean
%   points 'eucPnt', details follow.
%
%   The landmark's position manifold is defined trivially, so that
%       dp = [dpx dpy dpz]'
%   and the composition is just a sum,
%       p+ = p + dp.
%   In such case, the Jacobian
%       M = d p+ / d dp
%   is the identity matrix. To save some computations, this function
%   returns just M = 1 (scalar 'one').
%
%   For details on the Jacobians of other landmark parametrizations, please
%   refer to the code by editing the file lmkJacobians.m.
%
%   See also FRMJACOBIANS.

% Copyright 2015-     Joan Sola @ IRI-UPC-CSIC.

for lmk = [Lmk([Lmk.used]).lmk]
    switch Lmk(lmk).type
        case 'eucPnt'
            Lmk(lmk).state.M = 1; % trivial Jac
        case 'hmgPnt'
            [~,~,H_dh] = composeHmgPnt(Lmk(lmk).state.x, zeros(3,1));
            Lmk(lmk).state.M = H_dh;
        otherwise
            error('??? Unknown landmark type ''%s'' or Jacobian not implemented.',Lmk.type)
    end
end

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function sig=sigmerge(x1,x2,ratio);
%SIGMERGE Add two signals with given energy ratio in dB.
%	SIG=SIGMERGE(X1,X2,RATIO) adds two signals so that a given
%	energy ratio expressed in deciBels is satisfied.
%
%	X1, X2 : input signals.
%	RATIO  : Energy ratio in deciBels	(default : 0 dB).
%	X      : output signal.
%	X= X1+H*X2, such that 10*log(Energy(X1)/Energy(H*X2))=RATIO
%
%	Example :
%	 sig=fmlin(64,0.01,0.05,1); noise=hilbert(randn(64,1));
%	 SNR=15; x=sigmerge(sig,noise,SNR);
%	 Esig=mean(abs(sig).^2); Enoise=mean(abs(x-sig).^2);
%	 10*log10(Esig/Enoise)

%	F. Auger, July 1995.
%	Copyright (c) 1996 by CNRS (France).
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin<2)
 error('At least two parameters are required');
elseif nargin==2,
 ratio=0;
end;

[x1row,x1col] = size(x1);
[x2row,x2col] = size(x2);

if (x1col~=1)|(x2col~=1),
 error('X1 and X2 must have only one column');
elseif (x1row~=x2row),
 error('X1 and X2 must have the same number of rows');
elseif (length(ratio)~=1),
 error('RATIO must be a scalar');
elseif (ratio==inf),
 sig = x1;
else
 Ex1=mean(abs(x1).^2);
 Ex2=mean(abs(x2).^2);
 h=sqrt(Ex1/(Ex2*10^(ratio/10)));
 sig=x1+h*x2;
end;

author: HeLiangHIT


text:
function [Y,Z,dy,dz] = define_grid_coordinates(Ly,Lz,Ny,Nz,HubHt)
dy = Ly/(Ny-1); % Spacing along Y axis
dz = Lz/(Nz-1); % Spacing along Z axis
if isequal(mod(Ny,2),0)
    iky = flip([(-Ny/2:-1) (1:Ny/2)]);
else
    iky = flip(-floor(Ny/2):ceil(Ny/2-1));
end

if isequal(mod(Nz,2),0)
    ikz = [(-Nz/2:-1) (1:Nz/2)];
else
    ikz = -floor(Nz/2):ceil(Nz/2-1);
end

% define Y and Z coordinates for each grid-point
[Y,Z] = ndgrid(iky*dy,(ikz*dz + HubHt));
Y = Y';
Z = Z';
end

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


text:
function [xin,iin,cin] = faces_in_polygon(V,F,X,Y,Vin)
  % FACES_IN_POLYGON test whether faces of mesh are inside a given polygon
  %
  % [xin,iin,cin] = faces_in_polygon(V,F,X,Y)
  % [xin,iin,cin] = faces_in_polygon(V,F,X,Y)
  % [xin,iin,cin] = faces_in_polygon(V,F,X,Y,Vin)
  %
  % Inputs:
  %   V  #V by dim list of vertex positions
  %   F  #F by 3 list of face indices
  %   X  #P by 1 list of polygon x-coordinates
  %   Y  #P by 1 list of polygon y-coordinates
  %   Optional
  %     Vin  #V list of flags overriding whether each position in V is defined
  %       as inside polygon formed by (X,Y), default is to use inpolygon
  % Outputs:
  %   xin #F list of flags revealing whether query face indices are all in (X,Y)
  %   iin #F list of flags revealing whether query faces have at least one
  %     index in (X,Y)
  %   cin #F list of flags revealing whether query face barycenters are in (X,Y)
  %
  % See in_mesh, inpolygon
  %

  % NaNs are not recognized as they are in inpolygon
  assert(~any(isnan(X)));
  assert(~any(isnan(Y)));

  dim = size(V,2);
  % only works in 2D
  assert(dim == 2);

  if ~exist('exclusive','var')
    exclusive = true;
  end

  if ~exist('Vin','var')
    % first determine mesh points strictly in or on polygon
    Vin = inpolygon(V(:,1),V(:,2),X,Y);
  end

  % find faces that are inside because all vertices are inside
  [~,xin] = limit_faces(F,Vin,true);
  [~,iin] = limit_faces(F,Vin,false);

  % Note: if you can assume that if a triangle's corners are all inside (X,Y)
  % then also it's centroid is inside, then you could speed this up by only
  % testing triangles not in xin

  % barycenters for all triangles
  B = barycenter(V,F);

  % only do this if asked to
  if nargout >= 3
    % determine if each triangle barycenter is in polygon (X,Y)
    cin = inpolygon(B(:,1),B(:,2),X,Y);
  end

end

author: alecjacobson

repo: gptoolbox


text:
function summarize(X,Time,Spec,vintage)
%summarize Display the detail table for data entering the DFM
%
% Description:
%   Display the detail table for the nowcast, decomposing nowcast changes
%   into news and impacts for released data series.

% Print results to command window
fprintf('\n\n\n');
fprintf('Table 2: Data Summary \n');

[T,N] = size(X);
%fprintf('Vintage: %s \n',datestr(datenum(vintage,'yyyy-mm-dd'),'mmmm dd, yyyy'));
fprintf('N = %4d data series \n',N);
fprintf('T = %4d observations from %10s to %10s \n',T,datestr(Time(1),'yyyy-mm-dd'),datestr(Time(end),'yyyy-mm-dd'));

fprintf('%30s | %17s    %12s    %10s    %8s    %8s    %8s    %8s \n',...
        'Data Series','Observations','Units','Frequency','Mean','Std. Dev.','Min','Max');
fprintf([repmat('-',1,130) '\n']);

for i = 1:N

    % time indexes for which there are observed values for series i
    t_obs = ~isnan(X(:,i));

    data_series = Spec.SeriesName{i};
    if length(data_series) > 30
        data_series = [data_series(1:27) '...'];
    end
    series_id   = Spec.SeriesID{i};
    if length(series_id) > 28
        series_id = [series_id(1:25) '...'];
    end
    series_id   = ['[' series_id ']'];
    num_obs     = length(X(t_obs,i));
    freq        = Spec.Frequency{i};
    if strcmp('m',freq)
        format_date = 'mmm yyyy';
        frequency   = 'Monthly';
    elseif strcmp('q',freq)
        format_date = 'QQ yyyy';
        frequency   = 'Quarterly';
    end
    units       = Spec.Units{i};
    transform   = Spec.Transformation{i};
    % display transformed units
    if strfind(units,'Index')
        units_transformed = 'Index';
    elseif strcmp('chg',transform)
        if strfind(units,'%')
            units_transformed = 'Ppt. change';
        else
            units_transformed = 'Level change';
        end
    elseif strcmp('pch',transform) && strcmp('m',freq)
        units_transformed = 'MoM%';
    elseif strcmp('pca',transform) && strcmp('q',freq)
        units_transformed = 'QoQ% AR';
    else
        if length([units ' (' transform ')']) > 12
            units_transformed = [units(1:6) ' (' transform ')'];
        else
            units_transformed = [units ' (' transform ')'];
        end
    end
    t_obs_start = find(t_obs,1,'first');
    t_obs_end   = find(t_obs,1,'last');
    obs_date_start  = datestr(Time(t_obs_start),format_date);
    obs_date_end    = datestr(Time(t_obs_end  ),format_date);
    date_range  = [obs_date_start '-' obs_date_end];
    y           = X(t_obs,i);
    d           = Time(t_obs);
    mean_series = mean(y,'omitnan');
    stdv_series = std(y,'omitnan');
    [min_series,t_min] = min(y);
    [max_series,t_max] = max(y);
    min_date = datestr(d(t_min),format_date);
    max_date = datestr(d(t_max),format_date);
    fprintf('%30s | %17d    %12s    %10s    %8.1f    %8.1f    %8.1f    %8.1f \n',...
            data_series,num_obs,units_transformed,frequency,mean_series,stdv_series,min_series,max_series);
    fprintf('%30s | %17s    %12s    %10s    %8s    %8s    %8s    %8s \n',...
            series_id,date_range,'','','','',min_date,max_date);

end

fprintf('\n\n\n');

end

author: FRBNY-TimeSeriesAnalysis

repo: Nowcasting


text:
function orientation = histogramOrientation(hist_data)

%initialize output
orientation = [];

try
    %check to see if patch is in the shape of "vertical" rectangles :)
    if  size(hist_data.XData,1)==4  && size(hist_data.XData, 2) > 1 && ...
            all(hist_data.XData(1,:)==hist_data.XData(2,:)) && ...
            all(hist_data.XData(3,:)==hist_data.XData(4,:)) && ...
            all(hist_data.YData(1,:)==hist_data.YData(4,:)) && ...
            all(hist_data.YData(2,:)==hist_data.YData(3,:));
        orientation = 'v';
        %check to see if patch is in the shape of "horizontal" rectangles :)
    elseif size(hist_data.YData,1)==4 && size(hist_data.YData, 2) > 1 && ...
            all(hist_data.YData(1,:)==hist_data.YData(2,:)) && ...
            all(hist_data.YData(3,:)==hist_data.YData(4,:)) && ...
            all(hist_data.XData(1,:)==hist_data.XData(4,:)) && ...
            all(hist_data.XData(2,:)==hist_data.XData(3,:));
        orientation = 'h';
    end
end

author: plotly


text:
function all_boxes = selective_search(image_filenames, output_filename)

addpath('Dependencies');

if(~exist('anigauss'))
    mex Dependencies/anigaussm/anigauss_mex.c Dependencies/anigaussm/anigauss.c -output anigauss
end

if(~exist('mexCountWordsIndex'))
    mex Dependencies/mexCountWordsIndex.cpp
end

if(~exist('mexFelzenSegmentIndex'))
    mex Dependencies/FelzenSegment/mexFelzenSegmentIndex.cpp -output mexFelzenSegmentIndex;
end

colorTypes = {'Hsv', 'Lab', 'RGI', 'H', 'Intensity'};
colorType = colorTypes{1}; % Single color space for demo

% Here you specify which similarity functions to use in merging
simFunctionHandles = {@SSSimColourTextureSizeFillOrig, @SSSimTextureSizeFill, @SSSimBoxFillOrig, @SSSimSize};
simFunctionHandles = simFunctionHandles(1:2); % Two different merging strategies

% Thresholds for the Felzenszwalb and Huttenlocher segmentation algorithm.
% Note that by default, we set minSize = k, and sigma = 0.8.
k = 200; % controls size of segments of initial segmentation.
minSize = k;
sigma = 0.8;

% Process all images.
all_boxes = {};
for i=1:length(image_filenames)
    im = imread(image_filenames{i});
    [boxes blobIndIm blobBoxes hierarchy] = Image2HierarchicalGrouping(im, sigma, k, minSize, colorType, simFunctionHandles);
    all_boxes{i} = BoxRemoveDuplicates(boxes);
end

if nargin > 1
    all_boxes
    save(output_filename, 'all_boxes', '-v7');
end

author: sergeyk


text:
function D = driving_function_mono_wfs_vss(x0,xv,srcv,Dv,f,conf)
%DRIVING_FUNCTION_MONO_WFS_VSS driving signal for a given set of virtual
%secondary sources and their driving signals
%
%   Usage: D = driving_function_mono_wfs_vss(x0,xv,srcv,Dv,f,conf)
%
%   Input parameters:
%       x0          - position, direction, and weights of the real secondary
%                     sources / m [nx7]
%       xv          - position, direction, and weights of the virtual secondary
%                     sources / m [mx7]
%       srcv        - type of virtual secondary sources [mx7]
%                         'pw' - plane wave (xv(:,1:3) defines the direction of
%                                the plane waves in this case)
%                         'fs' - focused source (xv(:,1:6) defines the position
%                                and orientation of the focused sources in this
%                                case)
%       Dv          - driving functions of virtual secondary sources [mx1]
%       f           - frequency / Hz
%       conf        - optional configuration struct (see SFS_config)
%
%   Output parameters:
%       D           - driving function [nx1]
%
%   See also: driving_function_mono_wfs, driving_function_mono_wfs_fs
%
%   References:
%       Spors and Ahrens (2010) - "Local Sound Field Synthesis by Virtual
%       Secondary Sources", in 40th Conference of the Audio Engineering Society,
%       Paper 6-3, http://www.aes.org/e-lib/browse.cfm?elib=15561

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 6;
nargmax = 6;
narginchk(nargmin,nargmax);
isargvector(Dv);
isargpositivescalar(f);
isargchar(srcv);
isargsecondarysource(x0,xv);
isargstruct(conf);

%% ===== Configuration ==================================================
dimension = conf.dimension;

%% ===== Computation ====================================================
% Secondary source selection and driving function to synthesise a single virtual
% secondary source
switch srcv
case 'fs'
    ssd_select = @(X0,XS) secondary_source_selection(X0,XS(1:6),'fs');
    driv = @(X0,XS) driving_function_mono_wfs_fs(X0(:,1:3),X0(:,4:6),XS,f,conf);

    if strcmp('2.5D',dimension) || strcmp('3D',dimension)
        % === Focussed Point Sink ===
        conf.driving_functions = 'default';
    elseif strcmp('2D',dimension)
        % === Focussed Line Sink ===
        % We have to use the driving function setting directly, because in
        % opposite to the case of a non-focused source where 'ps' and 'ls' are
        % available as source types, for a focused source only 'fs' is
        % available. Have a look at driving_function_mono_wfs_fs() for details
        % on the implemented focused source types.
        conf.driving_functions = 'line_sink';
    else
        error('%s: %s is not a known source type.',upper(mfilename),dimension);
    end
case 'pw'
    ssd_select = @(X0,XS) secondary_source_selection(X0,XS(1:3),'pw');
    driv = @(X0,XS) driving_function_mono_wfs_pw(X0(:,1:3),X0(:,4:6),XS,f,conf);
end

% Get driving signals
Nv = size(xv,1);
N0 = size(x0,1);
Dmatrix = zeros(N0,Nv);

% it's ok to have zero secondary sources selected for pwd
warning('off','SFS:x0');
for idx=1:Nv
    [x0s, xdx] = ssd_select(x0,xv(idx,:));
    if (~isempty(x0s))
        % Virtual secondary source position
        xs = repmat(xv(idx,1:3),[size(x0s,1) 1]);
        % Optional tapering
        wtap = tapering_window(x0s,conf);
        Dmatrix(xdx,idx) = driv(x0s,xs) .* wtap;
    end
end
warning('on','SFS:x0');

D = Dmatrix*(Dv(:).*xv(:,7));

author: sfstoolbox


text:
% Cho, K., Raiko, T., Ilin, A., and Karhunen, J.
% A Two-stage Pretraining Algorithm for Deep Boltzmann Machines
% http://users.ics.aalto.fi/kcho/papers/nips12workshop.pdf

% add the path of RBM code
addpath('..');

% load MNIST
load 'mnist_14x14.mat';

% shuffle the training data
perm_idx = randperm (size(X,1));
X = X(perm_idx, :);
X_labels = X_labels(perm_idx);

X_labels = X_labels + 1;
X_test_labels = X_test_labels + 1;

% structure of DBM: 784 - 500 - 500
layers = [size(X,2), 500, 500];
n_layers = length(layers);

pretrain = 0;
centering = 1;

if pretrain
    % pretraining (stage 1)
    Xp = X;
    Rs = cell(n_layers, 1);
    Qpre = cell(n_layers, 1);
    Qpre_mask = zeros(n_layers, 1);

    for l = 1:n_layers
        if mod(l, 2) == 0
            continue;
        end

        if l+1 > n_layers
            break;
        end

        R = default_rbm (size(Xp, 2), layers(l+2));

        R.learning.lrate = 1e-2;
        R.learning.weight_decay = 1e-4;

        R.learning.persistent_cd = 0;

        R.fast.use = 0;
        R.fast.lrate = 1e-2;

        R.parallel_tempering.use = 0;

        R.adaptive_lrate.use = 1;
        R.adaptive_lrate.lrate_ub = R.learning.lrate;

        R.enhanced_grad.use = 1;
        R.learning.minibatch_sz = 128;

        % max. 100 epochs
        R.iteration.n_epochs = 100;

        if R.data.binary
            mH = mean(Xp, 1)';
            R.vbias = min(max(log(mH./(1 - mH)), -4), 4);
            R.fast.vbias = min(max(log(mH./(1 - mH)), -4), 4);
        else
            R.vbias = mean(Xp, 1)';
            R.fast.vbias = mean(Xp, 1)';
        end

        % set the stopping criterion
        R.stop.criterion = 0;
        R.stop.recon_error.tolerate_count = 1000;

        % save the intermediate data after every epoch
        R.hook.per_epoch = {@save_intermediate, {sprintf('rbm_%d.mat', l)}};
        R.hook.per_update = {};

        % print learining process
        R.verbose = 0;
        R.debug.do_display = 0;
        R.debug.display_interval = 10;
        R.debug.display_fid = 1;
        R.debug.display_function = @visualize_rbm;

        % train RBM
        fprintf(1, 'Training RBM\n');
        tic;
        R = train_rbm (R, Xp);
        fprintf(1, 'Training is done after %f seconds\n', toc);

        Rs{l} = R;

        Xp = rbm_get_hidden(Xp, R);

        Qpre{l+2} = Xp;
        Qpre_mask(l+2) = 1;
    end
end

[D] = default_dbm (layers);

if pretrain
    % pretraining (stage 2)
    D.learning.persistent_cd = 0;
    D.learning.lrate = 1e-2;
    D.learning.lrate0 = 1000;

    D.enhanced_grad.use = 0;
    D.adaptive_lrate.use = 1;
    D.adaptive_lrate.lrate_ub = D.learning.lrate;

    D.iteration.n_epochs = 200;
    D.hook.per_epoch = {@save_intermediate, {'dbm_pre.mat'}};

    fprintf(1, 'Training DBM\n');
    tic;
    [D] = dbm (D, X, 1, Qpre, Qpre_mask);
    fprintf(1, 'Training is done after %f seconds\n', toc);
end

% finetuning (stage 3)
D.learning.persistent_cd = 1;
if pretrain
    D.learning.lrate = 1e-4;
    D.learning.lrate0 = 5000;
else
    D.learning.lrate = 1e-3;
    D.learning.lrate0 = 5000;
end

D.enhanced_grad.use = 0;
D.adaptive_lrate.use = 1;
D.adaptive_lrate.lrate_ub = D.learning.lrate;

D.iteration.n_epochs = 100;
D.hook.per_epoch = {@save_intermediate, {'dbm_mnist.mat'}};

mH = mean(X, 1)';
D.biases{1} = min(max(log(mH./(1 - mH)), -4), 4);

if centering
    D = set_dbm_centers(D);
end

fprintf(1, 'Finetuning DBM\n');
tic;
[D] = dbm (D, X);
fprintf(1, 'Finetuning is done after %f seconds\n', toc);

% classification
[Q_mf] = dbm_get_hidden(X, D, 30, 1e-6, D.mf.reg);

perm_idx = randperm (size(X,1));

n_all = size(X, 1);
n_train = ceil(n_all * 3 / 4);
n_valid = floor(n_all /4);

X_valid = X(perm_idx(n_train+1:end), :);
X_valid_labels = X_labels(perm_idx(n_train+1:end));
X = X(perm_idx(1:n_train), :);
X_labels = X_labels(perm_idx(1:n_train));

for l = 1:n_layers
    Q_train{l} = Q_mf{l}(perm_idx(1:n_train), :);
    Q_valid{l} = Q_mf{l}(perm_idx(n_train+1:end), :);
end

[Q_test] = dbm_get_hidden(X_test, D, 30, 1e-6, D.mf.reg);

M = default_mlp ([layers, 10]);
M = set_mlp_dbm (M);

M.output.binary = 1;
M.hidden.use_tanh = 0;

M.valid_min_epochs = 10;
M.dropout.use = 0;

M.hook.per_epoch = {@save_intermediate, {'mlp_dbm_mnist.mat'}};

M.learning.lrate = 1e-2;
M.learning.lrate0 = 5000;
%M.learning.momentum = 0.9;
%M.learning.weight_decay = 0.0001;
M.learning.minibatch_sz = 128; M.adagrad.use = 1;
M.adagrad.epsilon = 1e-8;

M.noise.drop = 0;
M.noise.level = 0;

M.iteration.n_epochs = 100;

for l = 1:n_layers
    M.biases{l} = D.biases{l};
    if centering
        if l > 1
            M.biases{l} = M.biases{l} - M.W{l-1}' * D.centering.centers{l-1};
        end
        if l < n_layers
            M.biases{l} = M.biases{l} - M.W{l} * D.centering.centers{l+1};
        end
    end
    if l < n_layers
        M.W{l} = D.W{l};
        M.dbm.W{l} = D.W{l};
    end
end

fprintf(1, 'Training MLP\n');
tic;
M = mlp_dbm (M, X, Q_train, X_labels, X_valid, Q_valid, X_valid_labels, 0.1);
fprintf(1, 'Training is done after %f seconds\n', toc);

[pred] = mlp_classify (M, X_test, Q_test);
n_correct = sum(X_test_labels == pred);

fprintf(2, 'Correctly classified test samples: %d/%d\n', n_correct, size(X_test, 1));

author: kyunghyuncho

repo: deepmat


text:
function gray = gray_salt_and_pepper ( gray, level )

%*****************************************************************************80
%
%% GRAY_SALT_AND_PEPPER adds salt-and-pepper noise to a grayscale image.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, uint8 GRAY(:,:), the image.
%
%    Input, real LEVEL, the level of noise to add, between 0.0 (none)
%    and 1.0 (all).
%
%    Output, uint8 GRAY(:,:), the image with added noise.  A fraction of
%    about LEVEL of the pixels have been reset to 0 or 255.
%
  [ m, n ] = size ( gray );

  r = rand ( m, n );

  i0 = find ( r <= level / 2 );
  gray ( i0 ) = 0;

  i255 = find ( 1 - level / 2 <= r );
  gray ( i255 ) = 255;

  return
end


text:
function s = i4_to_s_roman ( i )

%*****************************************************************************80
%
%% I4_TO_S_ROMAN converts an integer to a string of Roman numerals.
%
%  Example:
%
%         I  S
%
%        -2  -II <-- Not a Roman numeral
%        -1  -I  <-- Not a Roman numeral
%         0   0  <-- Not a Roman numeral
%         1   I
%         2   II
%         3   III
%         4   IV
%         5   V
%        10   X
%        20   XX
%        30   XXX
%        40   XL
%        50   L
%        60   LX
%        70   LXX
%        80   LXXX
%        90   XC
%       100   C
%       500   D
%      1000   M
%      4999   MMMMCMLXLIX
%
%  Discussion:
%
%    To generate numbers greater than 4999, the numeral 'V' had a bar
%    above it, representing a value of 5000, a barred 'X' represented
%    10,000 and so on.
%
%    In the subtractive representation of 4 by 'IV', 9 by 'IX' and so on,
%    'I' can only subtract from 'V' or 'X',
%    'X' can only subtract from 'L' or 'C',
%    'C' can only subtract from 'D' or 'M'.
%    Under these rules, 1999 cannot be written IMM!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, an integer to be converted.  If the integer
%    has absolute value greater than 4999, the string '?' will be returned.
%    If the integer is 0, then the string '0' will be returned.  If
%    the integer is negative, then a minus sign will precede it, even
%    though this has nothing to do with Roman numerals.
%
%    Output, string S, the representation of the integer
%    as a Roman numeral.
%
  s = '';

  if ( 4999 < abs ( i ) )
    s = '?';
    s = char ( s );
    return
  end

  if ( i == 0 )
    s = '0';
    s = char ( s );
    return
  end

  if ( i <= 0 )
    s = '-';
    i = -i;
  end

  while ( 0 < i )

    if ( 1000 <= i )
      s = [ s 'M' ];
      i = i - 1000;
    elseif ( 900 <= i )
      s = [ s 'CM' ];
      i = i - 900;
    elseif ( 500 <= i )
      s = [ s 'D' ];
      i = i - 500;
    elseif ( 400 <= i )
      s = [ s 'CD' ];
      i = i - 400;
    elseif ( 100 <= i )
      s = [ s 'C' ];
      i = i - 100;
    elseif ( 90 <= i )
      s = [ s 'XC' ];
      i = i - 90;
    elseif ( 50 <= i )
      s = [ s 'L' ];
      i = i - 50;
    elseif ( 40 <= i )
      s = [ s 'XL' ];
      i = i - 40;
    elseif ( 10 <= i )
      s = [ s 'X' ];
      i = i - 10;
    elseif ( 9 <= i )
      s = [ s 'IX' ];
      i = i - 9;
    elseif ( 5 <= i )
      s = [ s 'V' ];
      i = i - 5;
    elseif ( 4 <= i )
      s = [ s 'IV' ];
      i = i - 4;
    else
      s = [ s 'I' ];
      i = i - 1;
    end

  end

  s = char ( s );

  return
end


text:
function obj = setMainSettings(obj)
	% Alter settings of object, in general don't need to use at the moment.
	% Biafra Ahanonu
	% branched from controllerAnalysis: 2014.08.01 [16:09:16]
	% inputs
		%
	% outputs
		%

	% changelog
		% 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		%

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	% h= uicontrol('Style','Text','String','hello: ','Position',[20 40 100 50],'BackgroundColor','white') ;
	% k = uicontrol('Style', 'popup','String', {'jet','hsv'},'Position', [100 40 200 50]);

	% the class should contain the available options as a property
	propertySettings = obj.settingOptions;

	propertyList = fieldnames(propertySettings);
	nPropertiesToChange = size(propertyList);
	% pulldownStr = '[choice] = pulldown(''settings''';

	% add current property to the top of the list
	for propertyNo = 1:nPropertiesToChange
		property = char(propertyList(propertyNo));
		propertyOptions = propertySettings.(property);
		currentProperty = obj.(property);
		propertySettings.(property) = cat(2,currentProperty,propertyOptions);
		propertySettings.(property)
	end

	uiListHandles = {};
	uiTextHandles = {};
	uiXIncrement = 0.05;
	figure(1337)
	for propertyNo = 1:nPropertiesToChange
		property = char(propertyList(propertyNo));
		uiTextHandles{propertyNo} = uicontrol('Style','Text','String',[property ': '],'Units','normalized','Position',[0.0 1-uiXIncrement*propertyNo 0.1 0.05],'BackgroundColor','white') ;
		uiListHandles{propertyNo} = uicontrol('Style', 'popup','String', propertySettings.(property),'Units','normalized','Position', [0.1 1-uiXIncrement*propertyNo 0.3 0.05]);
	end
	display('make choices then press enter...')
	pause

	for propertyNo = 1:nPropertiesToChange
		property = char(propertyList(propertyNo));
		uiListHandleData = get(uiListHandles{propertyNo});
		% uiListHandleData
		propertyOptions = propertySettings.(property);
		% obj.(property) = propertyOptions{choice(propertyNo)};
		% uiListHandleData.Value
		obj.(property) = propertyOptions{uiListHandleData.Value};
	end

	close(1337)
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
%**************************************************************
%* mex interface to Andy Liaw et al.'s C code (used in R package randomForest)
%* Added by Abhishek Jaiantilal ( abhishek.jaiantilal@colorado.edu )
%* License: GPLv2
%* Version: 0.02
%
% Calls Regression Random Forest
% A wrapper matlab file that calls the mex file
% This does prediction given the data and the model file
%**************************************************************

function [Y_hat,prediction_per_tree, nodes] = regRF_predict(X,model,extra_options)
    %function Y_hat = regRF_predict(X,model)
    %requires 2 arguments
    %X: data matrix
    %model: generated via regRF_train function
	if nargin<2
		error('need atleast 2 parameters,X matrix and model');
    end
    if exist('extra_options','var')
        if isfield(extra_options,'predict_all')
            predict_all = extra_options.predict_all;
        end
        if isfield(extra_options,'nodes')
            nodes = extra_options.nodes;
        end
    end

    if ~exist('predict_all','var'); predict_all=0;end

    if ~exist('nodes','var'); nodes=0;end

    if isfield(model, 'categorical_feature')
        % have to map prediction array to the correct categories
        for i=1:size(X,2)
            if model.categorical_feature(i)
                tmp_uniques_in_feature = model.orig_uniques_in_feature{i};
                tmp_mapped_uniques_in_feature = model.mapped_uniques_in_feature{i};
                X_loc = X(:,i); %cannot change the original array which may cause chained change of categories to something totally wrong
                for j=1:length(tmp_uniques_in_feature)
                    indices_to_change = find( X(:,i) == tmp_uniques_in_feature(j) );
                    X_loc(indices_to_change) = tmp_mapped_uniques_in_feature(j);
                end
                X(:,i) = X_loc;
            end
        end
        ncat = model.ncat;
    else
        ncat = ones(1,size(X,2));
    end

    maxcat = max(ncat);
	[Y_hat,prediction_per_tree, nodes] = mexRF_predict(X',model.lDau,model.rDau,model.nodestatus,model.nrnodes,model.upper,model.avnode,model.mbest,model.ndtree,model.ntree,predict_all, nodes, int32(ncat), maxcat );

    if ~isempty(find(model.coef)) %for bias corr
        Y_hat = model.coef(1) + model.coef(2)*Y_hat;
    end
	clear mexRF_predict

author: s-gupta

repo: rcnn-depth

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


text:
function h=spreadop(f,coef)
%SPREADOP  Spreading operator
%   Usage: h=spreadop(f,c);
%
%   `spreadop(f,c)` applies the operator with spreading function *c* to the
%   input *f*. *c* must be square.
%
%   `spreadop(f,c)` computes the following for *c* of size $L \times L$:
%
%   ..         L-1 L-1
%     h(l+1) = sum sum c(m+1,n+1)*exp(2*pi*i*l*m/L)*f(l-n+1)
%              n=0 m=0
%
%   .. math:: h\left(l+1\right)=\sum_{n=0}^{L-1}\sum_{m=0}^{L-1}c\left(m+1,n+1\right)e^{2{\pi}ilm/L}f\left(l-n+1\right)
%
%   where $l=0,\ldots,L-1$ and $l-n$ is computed modulo *L*.
%
%   The combined symbol of two spreading operators can be found by using
%   `tconv`. Consider two symbols *c1* and *c2* and define *f1* and *f2* by::
%
%     h  = tconv(c1,c2)
%     f1 = spreadop(spreadop(f,c2),c1);
%     f2 = spreadop(f,h);
%
%   then *f1* and *f2* are equal.
%
%   See also:  tconv, spreadfun, spreadinv, spreadadj
%
%   References: feko98

%   AUTHOR: Peter L. Søndergaard
%   TESTING: TEST_SPREAD
%   REFERENCE: REF_SPREADOP

complainif_argnonotinrange(nargin,2,2,mfilename);

if ndims(coef)>2 || size(coef,1)~=size(coef,2)
    error('Input symbol coef must be a square matrix.');
end;

L=size(coef,1);

% Change f to correct shape.
[f,Ls,W,wasrow,remembershape]=comp_sigreshape_pre(f,'DGT',0);

f=postpad(f,L);

h=zeros(L,W);

if issparse(coef) && nnz(coef)<L

    % The symbol is so sparse that the straighforward definition is
    % the fastest way to apply it.

    [mr,nr,cv]=find(coef);
    h=zeros(L,W);

    % We need mr and nr to be zero-indexed
    mr=mr-1;
    nr=nr-1;

    % This is the basic idea of the routine below
    %for ii=1:length(mr)
    %  for l=0:L-1
    %	h(l+1,:)=h(l+1,:)+cv(ii)*exp(-2*pi*i*l*mr(ii)/L)*f(mod(l-nr(ii),L)+1,:);
    %  end;
    %end;

    l=(-2*pi*i*(0:L-1)/L).';
    for ii=1:length(mr)
        bigmod=repmat(exp(l*mr(ii)),1,W);
        h=h+cv(ii)*(bigmod.*circshift(f,nr(ii)));
    end;

else

  % This version only touches coef one column at a time, and it suited
  % if coef is sparse.

  if issparse(coef)
    for n=0:L-1
      % The 'full' below is required for Matlab compatibility, as
      % Matlab refuses to do an ifft of a sparse matrix.
      cf=ifft(full(coef(:,n+1)))*L;
      modind=mod((0:L-1).'-n,L)+1;
      h=h+repmat(cf,1,W).*f(modind,:);
    end;
  else
    for n=0:L-1
      cf=ifft(coef(:,n+1))*L;
      modind=mod((0:L-1).'-n,L)+1;
      h=h+repmat(cf,1,W).*f(modind,:);
    end;
  end;

end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function p22_default ( dim_num )

%*****************************************************************************80
%
%% P22_DEFAULT sets default values for problem 22.
%
%  Discussion:
%
%    If a problem uses vector parameters, then the spatial dimension
%    DIM_NUM is needed as input, so that the vector parameters can be
%    properly dimensioned.
%
%    Every problem keeps a count of the number of function calls.  Calling
%    this routine causes that count to be reset to 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the problem.
%
  i4 = 0;
  p22_i4 ( 'D', '*', i4 );

  r8 = 0.0;
  p22_r8 ( 'D', '*', r8 );

  i4vec = [];
  p22_i4vec ( 'D', '*', dim_num, i4vec );

  return
end


text:
function res = GALA_clustering(lJcov,J1, S, distance, A)

    Nl = S.Nl;
    Nd = length(lJcov)/Nl;

    CJcov = GALA_find_localmin(lJcov,Nl,Nd,A,S.threshold);

    for p=1:Nl
        tCJcov=[];
        tCJcov.clvi=CJcov.clvi(p,:);
        tCJcov.maxs=CJcov.maxs{p};
        mxymcor = 0;
        while mxymcor<S.mincorr
            oldsize = length(tCJcov.maxs);
            cursize = length(tCJcov.maxs);
            for i=1:oldsize
                sxycor = corrcoef(J1(tCJcov.clvi{i}+(p-1)*Nd,:)');
                sxycor(sxycor==1)=NaN;
                sxymcor = nanmean(nanmean(sxycor));
                if length(tCJcov.clvi{i})>S.maxsize || sxymcor<S.mincorr
                    cltd = tCJcov.clvi{i};
    %                 ntd = ceil(length(cltd)/20);
                    ntd = 2;
                    Jcovpi = J1(cltd+(p-1)*Nd,:);
                    d = distance(cltd,cltd);
                    Jcd = ones(length(cltd),length(cltd))-corrcoef(Jcovpi');
                    dd = d+S.distratio1*Jcd;
                    Y = squareform(dd);
                    Z = linkage(Y,'complete');
                    CL = cluster(Z,'maxclust',ntd);
                    clm = CL(cltd==tCJcov.maxs(i));
                    tCJcov.clvi{i} = cltd(CL==clm);
                    cld = setdiff(1:max(CL),clm);
                    for j=cld
                        ind = find(CL==j);
                        tCJcov.clvi{cursize+1} = cltd(ind);
                        dd = distance(cltd(ind),cltd(ind));
                        [nu md] = min(sum(dd));
                        tCJcov.maxs(cursize+1) = cltd(ind(md));
                        cursize = cursize+1;
                    end
                end
            end
            i=1;
            while i<=cursize
                if length(tCJcov.clvi{i})<5
                    tCJcov.clvi(i) = [];
                    tCJcov.maxs(i) = [];
                    cursize = length(tCJcov.maxs);
                else
                    i=i+1;
                end
            end
            xymcor=[];
            for i=1:cursize
                xycor = corrcoef(J1(tCJcov.clvi{i}+(p-1)*Nd,:)');
                xycor(xycor==1)=NaN;
                xymcor(i) = nanmean(nanmean(xycor));
            end
            mxymcor = min(xymcor);
        end
        for i=1:max(length(CJcov.maxs{p}),cursize)
            if i<=cursize
                CJcov.clvi{p,i}=tCJcov.clvi{i};
                CJcov.maxs{p}(i)=tCJcov.maxs(i);
            else
                CJcov.clvi{p,i}=[];
                CJcov.maxs{p}(i)=[];
            end
        end
    end

    maxsa = spm_cat(CJcov.maxs);

    ldist = distance(maxsa,maxsa);
    lmaxs = [];
    smaxs = [];
    clind = [];
    for p=1:Nl
        lmaxs = horzcat(lmaxs,CJcov.maxs{p}+Nd*(p-1));
        smaxs(end+1:end+length(CJcov.maxs{p}),:) =...
            horzcat(repmat(p,[length(CJcov.maxs{p}) 1]),CJcov.maxs{p}');
        clind(end+1:end+length(CJcov.maxs{p}),:) =...
            horzcat(repmat(p,[length(CJcov.maxs{p}) 1]),(1:length(CJcov.maxs{p}))');

    end
    lmaxs = lmaxs';

    Jcor = ones(length(lmaxs),length(lmaxs))-corrcoef(J1(lmaxs,:)');
    dd = ldist+S.distratio2*Jcor;

    Y = squareform(dd);
    Z = linkage(Y,S.linkmeth);
    if isfield(S.cluster,'maxclust')
        CL = cluster(Z,'maxclust',S.cluster.maxclust.maxclustsize);
    elseif isfield(S.cluster,'cutoff')
        CL = cluster(Z,'cutoff',S.cluster.cutoff.cutoffthresh,'criterion','distance');
    end

    % remove multiple maxs for one subject in one cluster
    % and merge appropriate basins
    rs=[];
    for cl=1:max(CL)
       list = smaxs(CL==cl,:);
       ind = find(CL==cl);
       lm = lJcov(lmaxs(ind(1)));
       ind_lm = 1;
       for i=2:size(list,1)
           if list(i,1)==list(i-1,1)
               if lJcov(lmaxs(ind(i)))>lm
                   rs = [rs;ind(ind_lm)];
                   CJcov.clvi{clind(ind(i),1),clind(ind(i),2)}=...
                       horzcat(CJcov.clvi{clind(ind(i),1),clind(ind(i),2)},...
                       CJcov.clvi{clind(ind(ind_lm),1),clind(ind(ind_lm),2)});
                   CJcov.clvi{clind(ind(ind_lm),1),clind(ind(ind_lm),2)}=[];
                   lm = lJcov(lmaxs(ind(i)));
                   ind_lm = i;
               else
                   rs = [rs;ind(i)];
                   CJcov.clvi{clind(ind(ind_lm),1),clind(ind(ind_lm),2)}=...
                       horzcat(CJcov.clvi{clind(ind(i),1),clind(ind(i),2)},...
                       CJcov.clvi{clind(ind(ind_lm),1),clind(ind(ind_lm),2)});
                   CJcov.clvi{clind(ind(i),1),clind(ind(i),2)}=[];
               end
           else
               lm = lJcov(lmaxs(ind(i)));
               ind_lm = i;
           end
       end
    end

    CL(rs)=[];
    lmaxs(rs)=[];
    smaxs(rs,:)=[];
    clind(rs,:)=[];
    ldist(rs,:) = [];
    ldist(:,rs) = [];

    CLi = unique(CL);

thr = S.Nl-S.similarity;

dd=[];list=[];ncl=0;pclvi=[];pmaxs=[];tc=0;
for cl=1:length(CLi)
    if length(lmaxs(CL==CLi(cl)))>=thr
        ncl=ncl+1;
        ind = find(CL==CLi(cl));
        dd = ldist(ind,ind);
        [nu, md] = min(sum(dd));
        list{CLi(cl)} = smaxs(ind,:);
        mvd = list{CLi(cl)}(md,2);
        no = find(~ismember(1:Nl,list{CLi(cl)}(:,1)));
        iclvi = CJcov.clvi{clind(ind(1),1),clind(ind(1),2)};
        uclvi = CJcov.clvi{clind(ind(1),1),clind(ind(1),2)};
        wi = 1;
        while isempty(iclvi)
            wi=wi+1;
            iclvi = CJcov.clvi{clind(ind(wi),1),clind(ind(wi),2)};
            no = sort([no clind(ind(wi-1))]);
            ind(wi-1)=[];
        end
        pclvi{clind(ind(1),1),ncl} = iclvi;
        pmaxs{clind(ind(1),1)}(ncl) = smaxs(ind(1),2);
        for i=2:length(ind)
            pclvi{clind(ind(i),1),ncl} = CJcov.clvi{clind(ind(i),1),clind(ind(i),2)};
            pmaxs{clind(ind(i),1)}(ncl) = smaxs(ind(i),2);
            if isempty(pclvi{clind(ind(i),1),ncl})
                pclvi{clind(ind(i),1),ncl} = iclvi;
                pmaxs{clind(ind(i),1)}(ncl) = smaxs(ind(1),2);
            end
            iclvin = intersect(iclvi,pclvi{clind(ind(i),1),ncl});
            uclvin = union(uclvi,pclvi{clind(ind(i),1),ncl});
        end
        oth = setdiff(1:Nl,list{CLi(cl)}(:,1));
        for i=1:length(oth)
            if length(iclvi)>15
                pclvi{oth(i),ncl} = iclvi;
                pmaxs{oth(i)}(ncl) = smaxs(ind(1),2);
            else
                dd = distance(uclvin,uclvin);
                [nu, md] = min(sum(dd));
                pmaxs{oth(i)}(ncl) = uclvin(md);
                pclvi{oth(i),ncl} = find(distance(uclvin(md),:)<=2);
            end
        end

    end

end

res.pclvi = pclvi;
res.pmaxs = pmaxs;

%     figure;
%     ahp = [];
%     for p=1:Nl
%
%         srcs_disp1 = zeros(Nd,1);
%         for cl=1:size(pclvi,2)
%             srcs_disp1(pclvi{p,cl})=1;
%         end
%
%         ahp(p)=subplot(1,2,p);
%
%         cla; axis off
%
%         fig1 = patch('vertices',vert,'faces',face,'FaceVertexCData',srcs_disp1);
%
%         set(fig1,'FaceColor',[.5 .5 .5],'EdgeColor','none');
%         shading interp
%         lighting gouraud
%         %camlight
%         zoom off
%         lightangle(0,270);lightangle(270,0),lightangle(90,0),lightangle(0,45),lightangle(0,135);
%         material([.1 .1 .4 .5 .4]);
%         %view(140,15);
%         caxis([0 1])
%         colormap(jet);
%         hold on;
%         scatter3(vert(pmaxs{p},1),vert(pmaxs{p},2),vert(pmaxs{p},3),15,'r','fill');hold on;
%
%     end
%     hlink = linkprop(ahp, {'CameraPosition','CameraUpVector'});
%     key = 'graphics_linkprop';
%     setappdata(ahp(1),key,hlink);

end

author: spm