ivs = speakerRecognition;

detectionErrorTradeoff(ivs)

targetSampleRate = ivs.SampleRate;
[adsEnroll,adsDev,adsCalibrate] = loadDataset(targetSampleRate);

enroll(ivs,adsEnroll,adsEnroll.Labels)

detectionErrorTradeoff(ivs,adsDev,adsDev.Labels)

reset(adsDev)
[audioIn,fileInfo] = read(adsDev);
targetLabel = fileInfo.Label;

[decision,rawScore] = verify(ivs,audioIn,targetLabel,'css')

enrolledLabels = categorical(ivs.EnrolledLabels.Properties.RowNames);
nontargetIdx = find(~ismember(enrolledLabels,targetLabel));
nontargetLabel = enrolledLabels(nontargetIdx(1));

[decision,rawScore] = verify(ivs,audioIn,nontargetLabel,'css')

calibrationIvecs = ivector(ivs,adsCalibrate);

targets = true(size(calibrationIvecs,2),size(calibrationIvecs,2));
calibrationLabels = adsCalibrate.Labels;
for ii = 1:size(calibrationIvecs,2)
    targets(:,ii) = ismember(calibrationLabels,calibrationLabels(ii));
end

targets = targets + diag(diag(targets)*nan);
[targetScores,nontargetScores] = scoreTargets(calibrationIvecs,calibrationIvecs,targets);

plotScoreDistributions(targetScores,nontargetScores,'Analyze','group')

tS = cat(1,targetScores{:});
ntS = cat(1,nontargetScores{:});
x = [tS;ntS];

yplus = (numel(tS) + 1)./(numel(tS) + 2);
yminus = 1./(numel(ntS) + 2);
y = [yplus*ones(numel(tS),1);yminus*ones(numel(ntS),1)];

init = [1,1];
AB = fminsearch(@(AB)logRegCost(y,x,AB),init);

[x,idx] = sort(x,'ascend');
trueLabel = [ones(numel(tS),1);zeros(numel(ntS),1)];
trueLabel = trueLabel(idx);

calibratedScores = calibrateScores(x,AB);

plot(x,trueLabel,'o')
hold on
plot(x,calibratedScores,'LineWidth',1.5)
grid on
xlabel('Raw Score')
ylabel('Calibrated Score')
hold off

reset(adsDev)
[audioIn,fileInfo] = read(adsDev);
targetLabel = fileInfo.Label;

[decision,rawScore] = verify(ivs,audioIn,targetLabel,'css')

calibratedScore = calibrateScores(rawScore,AB)

nontargetIdx = find(~ismember(enrolledLabels,targetLabel));
nontargetLabel = enrolledLabels(nontargetIdx(randperm(numel(nontargetIdx),1)));

[decision,rawScore] = verify(ivs,audioIn,nontargetLabel,'css')

calibratedScore = calibrateScores(rawScore,AB)

scoringMap = isotonicRegression(x,trueLabel);

plot(x,trueLabel,'o')
hold on
plot(scoringMap.Raw,scoringMap.Calibrated,'LineWidth',1.5)
grid on
xlabel('Raw Score')
ylabel('Calibrated Score')
hold off

reset(adsDev)
[audioIn,fileInfo] = read(adsDev);
targetLabel = fileInfo.Label;

[decision,rawScore] = verify(ivs,audioIn,targetLabel,'css')

calibratedScore = calibrateScores(rawScore,scoringMap)

nontargetIdx = find(~ismember(enrolledLabels,targetLabel));
nontargetLabel = enrolledLabels(nontargetIdx(1));

[decision,rawScore] = verify(ivs,audioIn,nontargetLabel,'css')

calibratedScore = calibrateScores(rawScore,scoringMap)

devIvecs = ivector(ivs,adsDev);

devLabels = adsDev.Labels;
targets = true(size(devIvecs,2),size(devIvecs,2));
for ii = 1:size(devIvecs,2)
    targets(:,ii) = ismember(devLabels,devLabels(ii));
end
targets = targets + diag(diag(targets)*nan);

[targetScores,nontargetScores] = scoreTargets(devIvecs,devIvecs,targets);

ts = cat(1,targetScores{:});
nts = cat(1,nontargetScores{:});
scores = [ts;nts];
trueLabels = [true(numel(ts),1);false(numel(nts),1)];

calibratedScoresPlattScaling = calibrateScores(scores,AB);

calibratedScoresIsotonicRegression = calibrateScores(scores,scoringMap);

numBins = 10;

reliabilityDiagram(trueLabels,calibratedScoresPlattScaling,numBins)

reliabilityDiagram(trueLabels,calibratedScoresIsotonicRegression,numBins)

function [adsEnroll,adsDev,adsCalibrate] = loadDataset(targetSampleRate)
%LOADDATASET Load PTDB-TUG data set
% [adsEnroll,adsDev,adsCalibrate] = loadDataset(targetSampleteRate)
% downloads the PTDB-TUG data set, resamples it to the specified target
% sample rate and save the results in your current folder. The function
% then creates and returns three audioDatastore objects. The enrollment set
% includes two utterances per speaker. The calibrate set does not overlap
% with the other data sets.

% Copyright 2021 The MathWorks, Inc.

rng(0)
url = 'https://www2.spsc.tugraz.at/databases/PTDB-TUG/SPEECH_DATA_ZIPPED.zip';
downloadFolder = tempdir;
datasetFolder = fullfile(downloadFolder,'PTDB-TUG');

% Download and unzip the dataset if it doesn't already exist.
if ~isfolder(datasetFolder)
    disp('Downloading PTDB-TUG (3.9 G) ...')
    unzip(url,datasetFolder)
end

% Resample the dataset and save to current folder if it doesn't already
% exist.
if ~isfolder(fullfile(pwd,'MIC'))
    ads = audioDatastore([fullfile(datasetFolder,"SPEECH DATA","FEMALE","MIC"),fullfile(datasetFolder,"SPEECH DATA","MALE","MIC")], ...
        'IncludeSubfolders',true, ...
        'FileExtensions','.wav', ...
        'LabelSource','foldernames');
    reduceDataset = false;
    if reduceDataset
        ads = splitEachLabel(ads,10);
    end
    adsTransform = transform(ads,@(x,y)fileResampler(x,y,targetSampleRate),'IncludeInfo',true);
    writeall(adsTransform,pwd,'OutputFormat','flac','UseParallel',~isempty(ver('parallel')))
end

% Create a datastore that points to the resampled dataset. Use the folder
% names as the labels.
ads = audioDatastore(fullfile(pwd,'MIC'),'IncludeSubfolders',true,'LabelSource',"foldernames");

% Split the data set into enrollment, development, and calibration sets.
calibrationLabels = categorical(["M01","M03","M05","M7","M9","F01","F03","F05","F07","F09"]);

adsCalibrate = subset(ads,ismember(ads.Labels,calibrationLabels));

adsDev = subset(ads,~ismember(ads.Labels,calibrationLabels));

numToEnroll = 2;
[adsEnroll,adsDev] = splitEachLabel(adsDev,numToEnroll);

end

function [audioOut,adsInfo] = fileResampler(audioIn,adsInfo,targetSampleRate)
%FILERESAMPLER Resample audio files
% [audioOut,adsInfo] = fileResampler(audioIn,adsInfo,targetSampleRate)
% resamples the input audio to the target sample rate and updates the info
% passed through the datastore.

% Copyright 2021 The MathWorks, Inc.

arguments
    audioIn (:,1) {mustBeA(audioIn,{'single','double'})}
    adsInfo (1,1) {mustBeA(adsInfo,'struct')}
    targetSampleRate (1,1) {mustBeNumeric,mustBePositive}
end

% Isolate the original sample rate
originalSampleRate = adsInfo.SampleRate;

% Resample if necessary
if originalSampleRate ~= targetSampleRate
    audioOut = resample(audioIn,targetSampleRate,originalSampleRate);
    amax = max(abs(audioOut));
    if max(amax>1)
        audioOut = audioOut./amax;
    end
end

% Update the info passed through the datastore
adsInfo.SampleRate = targetSampleRate;

end

function [targetScores,nontargetScores] = scoreTargets(e,t,targetMap,nvargs)
%SCORETARGETS Score i-vector pairs
% [targetScores,nontargetScores] = scoreTargets(e,t,targetMap) exhaustively
% scores i-vectors in e against i-vectors in t. Specify e as an M-by-N
% matrix, where M corresponds to the i-vector dimension, and N corresponds
% to the number of i-vectors in e. Specify t as an M-by-P matrix, where P
% corresponds to the number of i-vectors in t. Specify targetMap as a
% P-by-N numeric matrix that maps which i-vectors in e and t are target
% pairs (derived from the same speaker) and which i-vectors in e and t are
% non-target pairs (derived from different speakers). Values in targetMap
% set to NaN are discarded. The outputs, targetScores and nontargetScores,
% are N-element cell arrays. Each cell contains a vector of scores between
% the i-vector in e and either all the targets or nontargets in t.
%
% [targetScores,nontargetScores] =
% scoreTargets(e,t,targetMap,'NormFactorsSe',NFSe,'NormFactorsSt',NFSt)
% normalizes the scores by the specified normalization statistics contained
% in structs NFSe and NFSt. If unspecified, no normalization is applied.

% Copyright 2021 The MathWorks, Inc.

arguments
    e (:,:) {mustBeA(e,{'single','double'})}
    t (:,:) {mustBeA(t,{'single','double'})}
    targetMap (:,:)
    nvargs.NormFactorsSe = [];
    nvargs.NormFactorsSt = [];
end

% Score the i-vector pairs
scores = cosineSimilarityScore(e,t);

% Apply as-norm1 if normalization factors supplied.
if ~isempty(nvargs.NormFactorsSe) && ~isempty(nvargs.NormFactorsSt)
    scores = 0.5*( (scores - nvargs.NormFactorsSe.mu)./nvargs.NormFactorsSe.std + (scores - nvargs.NormFactorsSt.mu')./nvargs.NormFactorsSt.std' );
end

% Separate the scores into targets and non-targets
targetScores = cell(size(targetMap,2),1);
nontargetScores = cell(size(targetMap,2),1);
removeIndex = isnan(targetMap);
for ii = 1:size(targetMap,2)
    localScores = scores(:,ii);
    localMap = targetMap(:,ii);
    localScores(removeIndex(:,ii)) = [];
    localMap(removeIndex(:,ii)) = [];

    targetScores{ii} = localScores(logical(localMap));
    nontargetScores{ii} = localScores(~localMap);
end
end

function scores = cosineSimilarityScore(a,b)
%COSINESIMILARITYSCORE Cosine similarity score
% scores = cosineSimilarityScore(a,b) scores matrix of i-vectors, a,
% against matrix of i-vectors b. Specify a as an M-by-N matrix of
% i-vectors. Specify b as an M-by-P matrix of i-vectors. scores is returned
% as a P-by-N matrix, where columns corresponds the i-vectors in a
% and rows corresponds to the i-vectors in b and the elements of the array
% are the cosine similarity scores between them.

% Copyright 2021 The MathWorks, Inc.

arguments
    a (:,:) {mustBeA(a,{'single','double'})}
    b (:,:) {mustBeA(b,{'single','double'})}
end

scores = squeeze(sum(a.*reshape(b,size(b,1),1,[]),1)./(vecnorm(a).*reshape(vecnorm(b),1,1,[])));
scores = scores';
end

function plotScoreDistributions(targetScores,nontargetScores,nvargs)
%PLOTSCOREDISTRIBUTIONS Plot target and non-target score distributions
% plotScoreDistribution(targetScores,nontargetScores) plots empirical
% estimations of the distribution for target scores and nontarget scores.
% Specify targetScores and nontargetScores as cell arrays where each
% element contains a vector of speaker-specific scores.
%
% plotScoreDistrubtions(targetScores,nontargetScores,'Analyze',ANALYZE)
% specifies the scope for analysis as either 'label' or 'group'. If ANALYZE
% is set to 'label', then a score distribution plot is created for each
% label. If ANALYZE is set to 'group', then a score distribution plot is
% created for the entire group by combining scores across speakers. If
% unspecified, ANALYZE defaults to 'group'.

% Copyright 2021 The MathWorks, Inc.

arguments
    targetScores (1,:) cell
    nontargetScores (1,:) cell
    nvargs.Analyze (1,:) char {mustBeMember(nvargs.Analyze,{'label','group'})} = 'group'
end

% Combine all scores to determine good bins for analyzing both the target
% and non-target scores together.
allScores = cat(1,targetScores{:},nontargetScores{:});
[~,edges] = histcounts(allScores);

% Determine the center of each bin for plotting purposes.
centers = movmedian(edges(:),2,'Endpoints','discard');

if strcmpi(nvargs.Analyze,'group')
     % Plot the score distributions for the group.

    targetScoresBinCounts = histcounts(cat(1,targetScores{:}),edges);
    targetScoresBinProb = targetScoresBinCounts(:)./sum(targetScoresBinCounts);
    nontargetScoresBinCounts = histcounts(cat(1,nontargetScores{:}),edges);
    nontargetScoresBinProb = nontargetScoresBinCounts(:)./sum(nontargetScoresBinCounts);

    figure
    plot(centers,[targetScoresBinProb,nontargetScoresBinProb])
    title("Score Distributions")
    xlabel('Score')
    ylabel('Probability')
    legend(["target","non-target"],'Location','northwest')
    axis tight

else

    % Create a tiled layout and plot the score distributions for each speaker.

    N = numel(targetScores);
    tiledlayout(N,1)
    for ii = 1:N
        targetScoresBinCounts = histcounts(targetScores{ii},edges);
        targetScoresBinProb = targetScoresBinCounts(:)./sum(targetScoresBinCounts);
        nontargetScoresBinCounts = histcounts(nontargetScores{ii},edges);
        nontargetScoresBinProb = nontargetScoresBinCounts(:)./sum(nontargetScoresBinCounts);
        nexttile
        hold on
        plot(centers,[targetScoresBinProb,nontargetScoresBinProb])
        title("Score Distribution for Speaker " + string(ii))
        xlabel('Score')
        ylabel('Probability')
        legend(["target","non-target"],'Location','northwest')
        axis tight
    end
end
end

function y = calibrateScores(score,scoreMapping)
%CALIBRATESCORES Calibrate scores
% y = calibrateScores(score,scoreMapping) maps the raw scores to calibrated
% scores, y, using the score mappinging information in scoreMapping.
% Specify score as a vector or matrix of raw scores. Specify score mapping
% as either struct or a two-element vector. If scoreMapping is specified as
% a struct, then it should have two fields: Raw and Calibrated, that
% together form a score mapping. If scoreMapping is specified as a vector,
% then the elements are used as the coefficients in the logistic function.
% y is returned as vector or matrix the same size as the raw scores.

% Copyright 2021 The MathWorks, Inc.

arguments
    score (:,:) {mustBeA(score,{'single','double'})}
    scoreMapping
end

if isstruct(scoreMapping)
    % Calibration using isotonic regression

    rawScore = scoreMapping.Raw;
    interpretedScore = scoreMapping.Calibrated;

    n = numel(score);

    % Find the index of the raw score in the mapping closest to the score provided.
    idx = zeros(n,1);
    for ii = 1:n
        [~,idx(ii)] = min(abs(score(ii)-rawScore));
    end

    % Get the calibrated score.
    y = interpretedScore(idx);

else

    % Calibration using logistic regression
    y = logistic(score,scoreMapping);

end
end

function reliabilityDiagram(targets,predictions,numBins)
%RELIABILITYDIAGRAM Plot reliability diagram
% reliabilityDiagram(targets,predictions) plots a reliability diagram for
% targets and predictions. Specify targets an M-by-1 logical vector.
% Specify predictions as an M-by-1 numeric vector.
%
% reliabilityDiagram(targets,predictions,numBins) specifies the number of
% bins for the reliability diagram. If unspecified, numBins defaults to 10.

% Copyright 2021 The MathWorks, Inc.

arguments
    targets (:,1) {mustBeA(targets,{'logical'})}
    predictions (:,1) {mustBeA(predictions,{'single','double'})}
    numBins (1,1) {mustBePositive,mustBeInteger} = 10;
end

% Bin the predictions into the requested number of bins. Count the number of
% predictions per bin.
[predictionsPerBin,~,predictionsInBin] = histcounts(predictions,numBins);

% Determine the mean of the predictions in the bin.
meanPredictions = accumarray(predictionsInBin,predictions)./predictionsPerBin(:);

% Determine the mean of the targets per bin. This is the fraction of
% positives--the number of targets in the bin over the total number of
% predictions in the bin.
meanTargets = accumarray(predictionsInBin,targets)./predictionsPerBin(:);

plot([0,1],[0,1])
hold on
plot(meanPredictions,meanTargets,'o')
legend('Ideal Calibration','Location','best')
xlabel('Mean Predicted Value')
ylabel('Fraction of Positives')
title('Reliability Diagram')
grid on
hold off

end

function cost = logRegCost(y,f,iparams)
%LOGREGCOST Logistic regression cost
% cost = logRegCost(y,f,iparams) calculates the cost of the logistic
% function given truth y, prediction f, and logistic params iparams.
% Specify y and f as column vectors. Specify iparams as a two-element row
% vector in the form [A,B], where A and B are the model parameters:
%
%                1
% p(x) = ------------------
%         1 + e^(-A*f - B)
%

% Copyright 2021 The MathWorks, Inc.

arguments
    y (:,1) {mustBeA(y,{'single','double'})}
    f (:,1) {mustBeA(f,{'single','double'})}
    iparams (1,2) {mustBeA(iparams,{'single','double'})}
end
p = logistic(f,iparams);
cost = -sum(y.*log(p) + (1-y).*log(1-p));
end

function p = logistic(f,iparams)
%LOGISTIC Logistic function
% p = logistic(f,iparams) applies the general logistic function to input f
% with parameters iparams. Specify f as a numeric array. Specify iparams as
% a two-element vector. p is returned as the same size as f.

% Copyright 2021 The MathWorks, Inc.

arguments
    f
    iparams = [1 0];
end
p = 1./(1+exp(-iparams(1).*f - iparams(2)));
end

function scoreMapping = isotonicRegression(x,y)
%ISOTONICREGRESSION Isotonic regression
% scoreMapping = isotonicRegression(x,y) fits a line yhat to data y under
% the monotonicity constraint that x(i)>x(j) -> yhat(i)>=yhat(j). That is,
% the values in yhat are monotontically non-decreasing with respect to x.
% The output, scoreMapping, is a struct containing the changepoints of yhat
% and the corresponding raw score in x.

% Copyright 2021, The MathWorks, Inc.

N = numel(x);

% Sort points in ascending order of x.
[x,idx] = sort(x(:),'ascend'); 
y = y(idx);

% Initialize fitted values to the given values.
m = y;

% Initialize blocks, one per point. These will merge and the number of
% blocks will reduce as the algorithm proceeds.
blockMap = 1:N;
w = ones(size(m));

while true

    diffs = diff(m);
    
    if all(diffs >= 0)

        % If all blocks are monotonic, end the loop.
        break;

    else

        % Find all positive changepoints. These are the beginnings of each
        % block.
        blockStartIndex = diffs>0;

        % Create group indices for each unique block.
        blockIndices = cumsum([1;blockStartIndex]);

        % Calculate the mean of each block and update the weights for the
        % blocks. We're merging all the points in the blocks here.
        m = accumarray(blockIndices,w.*m);
        w = accumarray(blockIndices,w);
        m = m ./ w;

        % Map which block corresponds to which index.
        blockMap = blockIndices(blockMap);

    end
end

% Broadcast merged blocks out to original points.
m = m(blockMap);

% Find the changepoints
changepoints = find(diff(m)>0);
changepoints = [changepoints;changepoints+1];
changepoints = sort(changepoints);

% Remove all points that aren't changepoints.
a = m(changepoints);
b = x(changepoints);

scoreMapping = struct('Raw',b,'Calibrated',a);
end