function [MLE,HST] = HoopStats_EstimMixtureModel(Responses,Targets,IncUniform,nVonmComps,MuVonmComps)
% "HoopStats_EstimMixtureModel"
% R01.03
% By Sam Berens (sam.berens@york.ac.uk).
%
% Estimates a mixture model for circularly distributed data. Components
% within the model can include a uniform distribution and arbitrarily many
% target/fixed position von Mises distributions.
%
% MLE = HoopStats_EstimMixtureModel(Responses,Targets,IncUniform,nVonM,tVonM)
%
% INPUT:
% - "Responses" : an nx1 numeric denoting n angles (in radians) of
% circular responses;
% - "Targets" : an nxt numeric denoting n angles (in radians) of t
% targets distributions. If sum(isnan(Targets),2) ==
% size(Targets,1), then response to target mappings are
% assumed to be mutually exclusive;
% - "IncUniform" : a 1x1 logical specifying whether model should include
% a uniform distribution;
% - "nVonmComps" : a 1x1 numeric denoting the number of non-target von
% Mises distributions to be modelled;
% - "MuVonmComps" : [Optional] a cx1 numeric denoting the preferred
% direction (in radians) for each of c non-target von
% Mises distributions;
%
% OUTPUT:
% - "MLE" : The best fitting model;
% - "HST" : Iteration history;
%% Clear, define, and Shuffle RNG:
rng('shuffle');
global Responses_X Responses_Y ScatterPointRadius LineThickness TerminateOndNll Tolerance_MaxDeltaParams;
Responses_X = real(exp(Responses.*1i));
Responses_Y = imag(exp(Responses.*1i));
if ~exist('MuVonmComps','var')
MuVonmComps = nan(nVonmComps,1);
elseif numel(MuVonmComps) ~= nVonmComps
error('MuVonmComps should be an nx1 double where n = nVonmComps;');
elseif size(MuVonmComps,1) ~= nVonmComps
MuVonmComps = MuVonmComps(:);
end
%% Estimation options:
EstimatePriors = true;
PriorBreakIn = Inf;
TerminateOndNll = true;
Tolerance_dNll = 1e-3;
Tolerance_MaxDeltaParams = 1e-3;
MaxIterations = 10^3;
NStarts = 3;
WNLLWPOST = 1; % Weight NLL with posterior probability: 0 or 1;
%%% Plot appearance:
ShowFigure = true;
ScreenSize = get(0,'ScreenSize');
% FigureHight = round(ScreenSize(1,4).*0.95);
% FigureSize = [ScreenSize(3)-(FigureHight.*16/9),(ScreenSize(4)-FigureHight)+20,FigureHight.*(16/9).*0.85,FigureHight.*0.85];
FigureSize = [1,ScreenSize(1)+40,ScreenSize(3),ScreenSize(4)-120];
ScatterPointRadius = 30;
LineThickness = 3;
%% Start Loop:
try
UniqueStarts = struct;
for iStart = 1:1:NStarts
%%% Specify starting parametres in ParamStruct:
ParamStruct = struct;
ParamStruct(1,1).Iteration = 0;
for iTarget = 1:1:size(Targets,2)
TargetId = sprintf('T%i',iTarget-1);
ParamStruct(1,1).(TargetId) = struct;
ParamStruct(1,1).(TargetId).DistributionType = 'vonMises';
ParamStruct(1,1).(TargetId).Mu = 0;
ParamStruct(1,1).(TargetId).K = (rand(1,1)*1)+0.5;
ParamStruct(1,1).(TargetId).Prior = 1/(nVonmComps+double(IncUniform)+size(Targets,2));
end
if IncUniform
ParamStruct(1,1).C0 = struct;
ParamStruct(1,1).C0.DistributionType = 'Uniform';
ParamStruct(1,1).C0.Prior = 1/(nVonmComps+1+size(Targets,2));
%%% Here the p(Responses_i|C0) is always = 1/(2*pi);
%%% Only the prior is estimated;
%%% The posterior is proportunal to p(Responses_i|C0) * Prior;
end
StartMus = sort(rand(nVonmComps,1)*2*pi);
for iVonmComp = 1:1:nVonmComps
ComponentId = sprintf('C%i',iVonmComp);
ParamStruct(1,1).(ComponentId) = struct;
ParamStruct(1,1).(ComponentId).DistributionType = 'vonMises';
if isnan(MuVonmComps(iVonmComp,1))
ParamStruct(1,1).(ComponentId).Mu = StartMus(iVonmComp,1);
else
ParamStruct(1,1).(ComponentId).Mu = MuVonmComps(iVonmComp,1);
end
ParamStruct(1,1).(ComponentId).K = (rand(1,1)*1)+0.5;
ParamStruct(1,1).(ComponentId).Prior = 1/(nVonmComps+double(IncUniform)+size(Targets,2));
end
for iTarget = 1:1:size(Targets,2)
TargetId = sprintf('T%i',iTarget-1);
ParamStruct(1,1).(['W_',TargetId]) = nan(size(Responses,1),1);
end
if IncUniform
ParamStruct(1,1).W_C0 = nan(size(Responses,1),1);
end
for iVonmComp = 1:1:nVonmComps
ComponentId = sprintf('C%i',iVonmComp);
ParamStruct(1,1).(['W_',ComponentId]) = nan(size(Responses,1),1);
% W is for Posterior;
end
ParamStruct(1,1).MaxDeltaParams = NaN;
ParamStruct(1,1).NLL = NaN;
%% Get the numbers and types of components being fitted:
AFun_CellFun00 = @(s) strcmp(s(1,1),'C');
AFun_CellFun01 = @(s) s(1,2);
Fields = fieldnames(ParamStruct);
DistNames = Fields(cellfun(AFun_CellFun00,Fields));
DistNumsPre = cellfun(AFun_CellFun01,DistNames);
DistNums = nan(1,size(DistNumsPre,1));
for i_DistNumsPre = 1:1:size(DistNumsPre,1)
DistNums(1,i_DistNumsPre) = str2double(DistNumsPre(i_DistNumsPre,1));
end
%% Draw initial figure:
if ShowFigure
FigureH = figure;
set(FigureH,'Position',FigureSize);
set(FigureH,'color','white');
PlotLastEmStep(ParamStruct,FigureH);
pause(1);
end
%% Start EM:
Finished = false;
iIteration = 0;
CurrentState = struct;
while ~Finished
%%% Incrament iteration count and set variables:
iIteration = iIteration + 1;
%%% Targets:
for iTarget = 1:1:size(Targets,2)
TargetId = sprintf('T%i',iTarget-1);
CurrentState.(TargetId).K = ParamStruct(iIteration,1).(TargetId).K;
CurrentState.(TargetId).Prior = ParamStruct(iIteration,1).(TargetId).Prior;
end
if IncUniform
CurrentState.C0.Prior = ParamStruct(iIteration,1).C0.Prior;
CurrentState.C0.Pdents = ones(size(Responses,1),1) .* (1/(2*pi));
end
for iVonmComp = 1:1:nVonmComps
% Remember that the first entry of 'ParamStruct' (index (i) = 1)
% contains model starting parameters. As such, the data for
% iteration 'a' is contained in index i = a + 1. Below, we read
% in model values from the previous iteration (i = a);
ComponentId = sprintf('C%i',iVonmComp);
CurrentState.(ComponentId).Mu = ParamStruct(iIteration,1).(ComponentId).Mu;
CurrentState.(ComponentId).K = ParamStruct(iIteration,1).(ComponentId).K;
CurrentState.(ComponentId).Prior = ParamStruct(iIteration,1).(ComponentId).Prior;
end
%% E-Step:
%%% Fixed Targets:
for iTarget = 1:1:size(Targets,2)
TargetId = sprintf('T%i',iTarget-1);
Selector = ~isnan(Targets(:,iTarget));
Acc = wrapToPi(Responses(Selector) - Targets(Selector,iTarget));
[CurrentState.(TargetId).Pdents,~] = HoopStats_VonmFit(0,CurrentState.(TargetId).K,Acc);
end
%%% Componants:
for iVonmComp = 1:1:nVonmComps
ComponentId = sprintf('C%i',iVonmComp);
[CurrentState.(ComponentId).Pdents,~] = HoopStats_VonmFit(CurrentState.(ComponentId).Mu,CurrentState.(ComponentId).K,Responses);
end
%%% Make Denominator:
Denominator = zeros(size(Responses));
% Targets:
for iTarget = 1:1:size(Targets,2)
TargetId = sprintf('T%i',iTarget-1);
Selector = ~isnan(Targets(:,iTarget));
TermToAdd = zeros(size(Responses));
TermToAdd(Selector) = CurrentState.(TargetId).Prior .* CurrentState.(TargetId).Pdents;
Denominator = Denominator + TermToAdd;
end
% Components:
for nComp = DistNums % This includes the Uniform if present;
ComponentId = sprintf('C%i',nComp);
TermToAdd = CurrentState.(ComponentId).Prior .* CurrentState.(ComponentId).Pdents;
Denominator = Denominator + TermToAdd;
end
%%% Compute weights:
% Targets:
for iTarget = 1:1:size(Targets,2)
TargetId = sprintf('T%i',iTarget-1);
Selector = ~isnan(Targets(:,iTarget));
CurrentState.(TargetId).W = zeros(size(Responses));
CurrentState.(TargetId).W(Selector) = (CurrentState.(TargetId).Prior .* CurrentState.(TargetId).Pdents) ./ Denominator(Selector);
end
% Components:
for nComp = DistNums
ComponentId = sprintf('C%i',nComp);
CurrentState.(ComponentId).W = (CurrentState.(ComponentId).Prior .* CurrentState.(ComponentId).Pdents) ./ Denominator;
end
%% M-Step:
% Targets:
for iTarget = 1:1:size(Targets,2)
TargetId = sprintf('T%i',iTarget-1);
Selector = ~isnan(Targets(:,iTarget));
Acc = wrapToPi(Responses(Selector) - Targets(Selector,iTarget));
[~,CurrentState.(TargetId).K] = HoopStats_GetParams(Acc,CurrentState.(TargetId).W(Selector));
CurrentState.(TargetId).Prior = mean(CurrentState.(TargetId).W(Selector),1);
end
% Components:
for nComp = DistNums
ComponentId = sprintf('C%i',nComp);
if nComp ~= 0
if isnan(MuVonmComps(nComp,1))
[CurrentState.(ComponentId).Mu,CurrentState.(ComponentId).K] = HoopStats_GetParams(Responses,CurrentState.(ComponentId).W);
else
[~,CurrentState.(ComponentId).K] = HoopStats_GetParams(Responses,CurrentState.(ComponentId).W);
end
end
CurrentState.(ComponentId).Prior = mean(CurrentState.(ComponentId).W,1);
end
%% Update DeltaParams:
iEntry = 0;
DeltaParams = nan(size(Targets,2)+numel(DistNums),3);
for iTarget = 1:1:size(Targets,2)
iEntry = iEntry + 1;
TargetId = sprintf('T%i',iTarget-1);
DeltaParams(iEntry,2) = abs((1/CurrentState.(TargetId).K) - (1/ParamStruct(iIteration,1).(TargetId).K)); % Here, DeltaParams for K are difined as [1/K_i - 1/K_(i+1)];
DeltaParams(iEntry,3) = abs(CurrentState.(TargetId).Prior - ParamStruct(iIteration,1).(TargetId).Prior);
end
for nComp = DistNums
iEntry = iEntry + 1;
ComponentId = sprintf('C%i',nComp);
if nComp ~= 0
DeltaParams(iEntry,1) = abs(CurrentState.(ComponentId).Mu - ParamStruct(iIteration,1).(ComponentId).Mu);
DeltaParams(iEntry,2) = abs(CurrentState.(ComponentId).K - ParamStruct(iIteration,1).(ComponentId).K);
end
DeltaParams(iEntry,3) = abs(CurrentState.(ComponentId).Prior - ParamStruct(iIteration,1).(ComponentId).Prior);
end
MaxDeltaParamsExcludingPriors = DeltaParams(:,1:2);
MaxDeltaParamsExcludingPriors = nanmax(MaxDeltaParamsExcludingPriors(:));
MDPEP = MaxDeltaParamsExcludingPriors;
clear MaxDeltaParamsExcludingPriors;
MaxDeltaParams = nanmax(DeltaParams(:));
%% Calculate NLL:
% Here, likelihood is calculated to be the sum of the probability
% densities for each component (T&C) and is weighted by the
% posterior terms P(C*|Responses) [W*].
iDist = 0;
Ps = nan(size(Responses,1),size(Targets,2)+numel(DistNums));
for iTarget = 1:1:size(Targets,2)
iDist = iDist + 1;
TargetId = sprintf('T%i',iTarget-1);
Selector = ~isnan(Targets(:,iTarget));
Acc = wrapToPi(Responses(Selector) - Targets(Selector,iTarget));
[CurrentState.(TargetId).Pdents,~] = HoopStats_VonmFit(0,CurrentState.(TargetId).K,Acc);
Ps(Selector,iDist) = CurrentState.(TargetId).Pdents .* (CurrentState.(TargetId).W(Selector).^WNLLWPOST);
end
for iVonmComp = 1:1:nVonmComps
ComponentId = sprintf('C%i',iVonmComp);
[CurrentState.(ComponentId).Pdents,~] = HoopStats_VonmFit(CurrentState.(ComponentId).Mu,CurrentState.(ComponentId).K,Responses);
end
for nComp = DistNums
iDist = iDist + 1;
ComponentId = sprintf('C%i',nComp);
Ps(:,iDist) = CurrentState.(ComponentId).Pdents .* (CurrentState.(ComponentId).W.^WNLLWPOST);
end
NLL = -sum(log(nansum(Ps,2)),1);
%% Update ParamStruct:
ParamStruct(iIteration+1,1).Iteration = iIteration;
%%% Targets:
for iTarget = 1:1:size(Targets,2)
TargetId = sprintf('T%i',iTarget-1);
ParamStruct(iIteration+1,1).(TargetId).Mu = 0;
ParamStruct(iIteration+1,1).(TargetId).K = CurrentState.(TargetId).K;
ParamStruct(iIteration+1,1).(['W_',TargetId]) = CurrentState.(TargetId).W;
if EstimatePriors && (MDPEP < PriorBreakIn)
ParamStruct(iIteration+1,1).(TargetId).Prior = CurrentState.(TargetId).Prior;
else
ParamStruct(iIteration+1,1).(TargetId).Prior = ParamStruct(iIteration,1).(TargetId).Prior;
end
end
%%% Components:
for nComp = DistNums
ComponentId = sprintf('C%i',nComp);
if nComp ~= 0
ParamStruct(iIteration+1,1).(ComponentId).Mu = CurrentState.(ComponentId).Mu;
ParamStruct(iIteration+1,1).(ComponentId).K = CurrentState.(ComponentId).K;
end
ParamStruct(iIteration+1,1).(['W_',ComponentId]) = CurrentState.(ComponentId).W;
if EstimatePriors && (MDPEP < PriorBreakIn)
ParamStruct(iIteration+1,1).(ComponentId).Prior = CurrentState.(ComponentId).Prior;
else
ParamStruct(iIteration+1,1).(ComponentId).Prior = ParamStruct(iIteration,1).(ComponentId).Prior;
end
end
ParamStruct(iIteration+1,1).MaxDeltaParams = MaxDeltaParams;
ParamStruct(iIteration+1,1).NLL = NLL;
%% Plot:
if ShowFigure% && (nanmin([ParamStruct.NLL],[],2) == NLL)
PlotLastEmStep(ParamStruct,FigureH)
end
%% Terminate if MaxDetlaParams is less than tollerence:
if TerminateOndNll
dNll = ParamStruct(iIteration,1).NLL - ParamStruct(iIteration+1,1).NLL;
if abs(dNll) < Tolerance_dNll
Finished = true;
elseif (abs(dNll) < 0.01) && (MaxDeltaParams < 0.001)
Finished = true;
end
else
if EstimatePriors %#ok<UNRCH>
if MaxDeltaParams < Tolerance_MaxDeltaParams
Finished = true;
end
else
if MDPEP < Tolerance_MaxDeltaParams
Finished = true;
end
end
end
%% Terminate if reached Iteration limit.
if iIteration == MaxIterations
Finished = true;
end
end
%%% Record Result:
UniqueStarts(iStart,1).ParamStruct = ParamStruct;
UniqueStarts(iStart,1).NLL_Last = ParamStruct(end,1).NLL;
[UniqueStarts(iStart,1).NLL_Min,UniqueStarts(iStart,1).NLL_MinIndex] = nanmin([ParamStruct.NLL],[],2);
close all;
end
UsedBestNll = true;
if UsedBestNll
[~,BestStart_Index] = min([UniqueStarts.NLL_Min],[],2);
BestIte_Index = UniqueStarts(BestStart_Index,1).NLL_MinIndex;
else
[~,BestStart_Index] = min([UniqueStarts.NLL_Last],[],2); %#ok<UNRCH>
BestIte_Index = size(UniqueStarts(BestStart_Index,1).ParamStruct,1);
end
HST = UniqueStarts(BestStart_Index,1).ParamStruct(1:BestIte_Index,1);
MLE = HST(end,1);
close all;
FigureH = figure;
set(FigureH,'Position',FigureSize);
set(FigureH,'color','white');
PlotLastEmStep(HST,FigureH);
Notification_A01;
catch
Notification_A02;
warning('Mixture model failed to converge! %04d%02d%02d-%02d%02d%02d;',fix(clock));
end
return
function [] = PlotLastEmStep(ParamStruct,FigureH)
global Responses_X Responses_Y ScatterPointRadius LineThickness TerminateOndNll Tolerance_MaxDeltaParams;
%VectorTextOffset = 0.6 + (-0.1*1i);
%% Open figure:
set(0,'CurrentFigure',FigureH);
clf;
subplot(2,2,[1,3]);
DrawCircle;
hold on;
%% Get Colours and scatter:
Fields = fieldnames(ParamStruct);
M = [];
for iFields = 1:1:size(Fields,1)
if strcmp(Fields{iFields,1}(1,1),'W')
M = [M,ParamStruct(end,1).(Fields{iFields,1})]; %#ok<AGROW>
end
end
[PointColours_RGB,Stereotypes_RGB] = MakeColours(M);
scatter(Responses_X,Responses_Y,ScatterPointRadius,PointColours_RGB,'filled',...
'MarkerEdgeColor',[0,0,0],'LineWidth',1);
%%
axis equal;
axis off;
text(0.5,1,sprintf('Iteration: %i;',ParamStruct(end,1).Iteration),'FontSize',10);
iDist = 0;
Priors = nan(0,1);
DistNames = cell(0,1);
ParamsToDisplay = nan(0,2);
InfoContent = nan(0,1);
for iFields = 1:1:size(Fields,1)
if strcmp(Fields{iFields,1}(1,1),'T')
iDist = iDist + 1;
DistNames{iDist,1} = Fields{iFields,1};
Prior = ParamStruct(end,1).(DistNames{iDist,1}).Prior;
K = ParamStruct(end,1).(DistNames{iDist,1}).K;
Priors(iDist,1) = Prior;
ParamsToDisplay(iDist,2) = K;
InfoContent(iDist,1) = HoopStats_InfoContent(Prior,K);
elseif strcmp(Fields{iFields,1}(1,1),'C')
iDist = iDist + 1;
DistNum = str2double(Fields{iFields,1}(1,2));
DistNames{iDist,1} = Fields{iFields,1};
if DistNum == 0
%%% Uniform distrobution:
Prior = ParamStruct(end,1).C0.Prior;
ParamsToDisplay(iDist,1) = NaN;
ParamsToDisplay(iDist,2) = NaN;
Thetas = linspace(0,2*pi);
X = cos(Thetas).*Prior;
Y = sin(Thetas).*Prior;
plot(X,Y,'Color',Stereotypes_RGB(iDist,:),'LineWidth',LineThickness./2);
InfoContent(iDist,1) = Prior * log(2*pi);
else
%% vonM:
Mu = ParamStruct(end,1).(DistNames{iDist,1}).Mu;
K = ParamStruct(end,1).(DistNames{iDist,1}).K;
Mag = HoopStats_K2R(K);
Prior = ParamStruct(end,1).(DistNames{iDist,1}).Prior;
ParamsToDisplay(iDist,1) = Mu;
ParamsToDisplay(iDist,2) = K;
InfoContent(iDist,1) = HoopStats_InfoContent(Prior,K);
%%% Maths:
Vector = exp(Mu*1i) * Mag;
Vector_X = real(Vector);
Vector_Y = imag(Vector);
% TextPos = Vector * VectorTextOffset;
% Text_X = real(TextPos);
% Text_Y = imag(TextPos);
%%% Draw:
plot([0;Vector_X],[0;Vector_Y],...
'Color',Stereotypes_RGB(iDist,:),...
'LineStyle','-',...
'LineWidth',LineThickness);
% text(Text_X,Text_Y,...
% ['\mu = ',num2str(Mu),';',10,...
% '\it{k} = ',num2str(K),';']);
end
Priors(iDist,1) = Prior;
end
end
ParamsToDisplay = [ParamsToDisplay,Priors];
%%
subplot(2,2,2);
for iInfoContent = 1:1:size(InfoContent,1)
H = bar(iInfoContent,InfoContent(iInfoContent,1));
if iInfoContent == 1
hold on;
end
set(H,'FaceColor',Stereotypes_RGB(iInfoContent,:));
end
ParamsToDisplay(ParamsToDisplay==0) = NaN;
for iParamsToDisplay = 1:1:size(ParamsToDisplay,1)
Text_X = iParamsToDisplay;
Text_Y = -0.5;
Mu = ParamsToDisplay(iParamsToDisplay,1);
K = ParamsToDisplay(iParamsToDisplay,2);
Prior = ParamsToDisplay(iParamsToDisplay,3);
if isnan(Mu) && isnan(K)
text(Text_X,Text_Y,...
['\it{p} = ',sprintf('%0.2f',Prior),';'],...
'HorizontalAlignment','center');
elseif isnan(Mu)
text(Text_X,Text_Y,...
['\it{k} = ',sprintf('%0.2f',K),';',10,...
'\it{p} = ',sprintf('%0.2f',Prior),';'],...
'HorizontalAlignment','center');
else
text(Text_X,Text_Y,...
['\mu = ',sprintf('%0.2f',Mu),';',10,...
'\it{k} = ',sprintf('%0.2f',K),';',10,...
'\it{p} = ',sprintf('%0.2f',Prior),';'],...
'HorizontalAlignment','center');
end
end
title('Parameter estimates');
ylabel('Information content (nats)');
ylim([-1,2.2]);
set(gca, 'XTick',1:size(DistNames,1));
set(gca,'XTickLabel',DistNames');
plot(xlim,[log(2*pi),log(2*pi)],'LineStyle','--','LineWidth',LineThickness./4,'Color','k');
%%
subplot(2,2,4);
XaxisData = (1:size(ParamStruct,1))';
Y2Add = nan(0,1);
if size(ParamStruct,1) < 100
N2Add = 100 - size(ParamStruct,1);
XaxisData = (1:1:100)';
Y2Add = nan(N2Add,1);
end
[HFit_Ax,HFit_L1,HFit_L2] = plotyy(...
XaxisData,[[ParamStruct.NLL]';Y2Add],...
XaxisData,[log10([ParamStruct.MaxDeltaParams]');Y2Add]); %#ok<ASGLU>
title('Iteration history');
ylabel(HFit_Ax(1),'Negative loig-likelihood');
ylabel(HFit_Ax(2),'Log_1_0 max \Delta params');
xlabel('Iteration number');
if TerminateOndNll
% Nothing
else
HFit_Ax(1,2).YLim = [log10(Tolerance_MaxDeltaParams),2];
end
HFit_Ax(1,2).YTick = flip([0,-2,-4,-6,-8,-10]);
if size(ParamStruct,1) < 100
HFit_Ax(1,1).XLim = [0,100];
HFit_Ax(1,2).XLim = [0,100];
end
pause(0.05);
hold off;
return
function [] = DrawCircle()
global LineThickness;
Thetas = linspace(0,2*pi);
X = cos(Thetas);
Y = sin(Thetas);
plot(X,Y,'k-','LineWidth',LineThickness./3);
return
function [PointColours_RGB,Stereotypes_RGB,FSDisp] = MakeColours(M)
if sum(sum(isnan(M))) == numel(M)
M = ones(size(M)) ./ size(M,2);
else
M(isnan(M)) = 0;
end
nColours = size(M,2);
Range = 1;
Offset = 0;
%%
Stereotypes_Angles = mod((((Range/nColours).*((1:1:nColours)-1))+Offset),1)' .* (2*pi);
PointColours_Complex = M * exp(Stereotypes_Angles*1i);
PointColours_Angles = angle(PointColours_Complex);
PointColours_Angles(PointColours_Angles<0) = (2*pi) + PointColours_Angles(PointColours_Angles<0);
PointColours_Angles = PointColours_Angles ./ (2*pi);
PointColours_Mags = abs(PointColours_Complex);
PointColours_HSV = [PointColours_Angles,PointColours_Mags,ones(size(PointColours_Complex,1),1)];
PointColours_RGB = hsv2rgb(PointColours_HSV);
%%
Stereotypes_HSV = [Stereotypes_Angles./(2*pi),ones(size(Stereotypes_Angles)),ones(size(Stereotypes_Angles))];
Stereotypes_RGB = hsv2rgb(Stereotypes_HSV);
FSDisp = @() image(reshape(Stereotypes_RGB,size(Stereotypes_RGB,1),1,3));
return
function [] = Notification_A01()
Fs = 48000;
Vol = 0.4;
L1 = 0.1; % Seconds;
P1 = 1000;
L2 = 0.1; % Seconds;
P2 = 1200;
%
[T1] = MakeSound(Fs,L1,P1,Vol);
[T2] = MakeSound(Fs,L2,P2,Vol);
sound(T1,Fs);
pause(L1);
sound(T2,Fs);
pause(L2);
sound(T2,Fs);
pause(L1);
sound(T1,Fs);
return
function [] = Notification_A02()
Fs = 48000;
Vol = 0.4;
L1 = 0.1; % Seconds;
P1 = 1200;
L2 = 0.1; % Seconds;
P2 = 900;
%
[T1] = MakeSound(Fs,L1,P1,Vol);
[T2] = MakeSound(Fs,L2,P2,Vol);
sound(T1,Fs);
pause(L1);
sound(T2,Fs);
pause(L2);
sound(T2,Fs);
pause(L1);
sound(T2,Fs);
return
function [fT] = MakeSound(Fs,Length,Pitch,Vol)
FadeTime = 0.01; % Seconds;
fT = 0:(Pitch/Fs):(Pitch*Length);
fT = sin(fT.*2.*pi);
InEnv = 0:(1/(Fs*FadeTime)):1;
OutEnv = 1:-(1/(Fs*FadeTime)):0;
fT(1:numel(InEnv)) = fT(1:numel(InEnv)) .* InEnv;
fT(end-numel(InEnv)+1:end) = fT(end-numel(InEnv)+1:end) .* OutEnv;
fT = fT .* Vol;
return