load speakerVerificationExampleData.mat ubm afe normFactors

enrollYourself = false;
if enrollYourself
    numToRecord = 5;
    ID = 'self';
    helperAddUser (afe.SampleRate, numToRecord, ID);
end

ads = audioDatastore(pwd);

[~,fileName] = cellfun(@(x)fileparts(x),ads.Files,'UniformOutput',false);
fileName = split(fileName,'_');
speaker = strcat(fileName(:,1));
ads.Labels = categorical(speaker);

if enrollYourself
    enrollLabel = ID;
else
    enrollLabel = 'BHm';
end

forEnrollment = ads.Labels==enrollLabel;
forEnrollment(find(forEnrollment==1,1)) = false;
adsEnroll = subset(ads,forEnrollment);
adsTest = subset(ads,~forEnrollment);

speakerGMM = helperEnroll(ubm,afe,normFactors,adsEnroll);

threshold = 0.7;
сброс (adsTest)
while hasdata (adsTest)
    fprintf ('Identity to confirm: %s\n', enrollLabel)
    [audioData, adsInfo] = read (adsTest);
    
    fprintf (' | Speaker identity: %s\n', string (adsInfo.Label))
    
    verificationStatus = helperVerify (audioData, afe, normFactors, speakerGMM, ubm, порог);

    if verificationStatus
        fprintf (' | Confirmed.\n');
    else
        fprintf (' | Imposter!\n');
    end
end

url = 'https://storage.googleapis.com/download.tensorflow.org/data/speech_commands_v0.01.tar.gz';

downloadFolder = tempdir;
datasetFolder = fullfile(downloadFolder,'google_speech');

if ~exist(datasetFolder,'dir')
    disp('Downloading Google speech commands data set (1.9 GB)...')
    untar(url,datasetFolder)
end

ads = audioDatastore(datasetFolder,"Includesubfolders",true,'LabelSource','folderNames');

ads = subset(ads,ads.Labels==categorical("stop"));

[~,fileName] = cellfun(@(x)fileparts(x),ads.Files,'UniformOutput',false);
fileName = split(fileName,'_');
speaker = strcat('a',fileName(:,1));
ads.Labels = categorical(speaker);

numSpeakersToEnroll = 10;
labelCount = countEachLabel (объявления);
forEnrollAndTestSet = labelCount {:, 1} (labelCount {:, 2} > = 3);
forEnroll = forEnrollAndTestSet (randi ([1, numel (forEnrollAndTestSet)], numSpeakersToEnroll, 1));
tf = ismember (ads.Labels, forEnroll);
adsEnrollAndValidate = подмножество (объявления, tf);
adsEnroll = splitEachLabel (adsEnrollAndValidate, 2);

adsTest = subset (ads, ismember (ads.Labels, forEnrollAndTestSet));
adsTest = subset (adsTest, ~ ismember (adsTest.Files, файлы));

forUBMTraining = ~ (ismember (ads.Files, файлы) | ismember (ads.Files, файлы));
adsTrainUBM = подмножество (объявления, forUBMTraining);

[audioData,audioInfo] = read(adsTrainUBM);
fs = audioInfo.SampleRate;

sound(audioData,fs)

reset(adsTrainUBM)

windowDuration = 0.04;
hopDuration = 0.01;
windowSamples = round(windowDuration*fs);
hopSamples = round(hopDuration*fs);
overlapSamples = windowSamples - hopSamples;

afe = audioFeatureExtractor( ...
    'SampleRate',fs, ...
    'Window',hann(windowSamples,'periodic'), ...
    'OverlapLength',overlapSamples, ...
    ...
    'mfcc',true);

audioData = audioData./max(abs(audioData));

detectSpeech(audioData,fs);

idx = detectSpeech(audioData,fs);
audioData = audioData(idx(1,1):idx(1,2));

features = extract(afe,audioData);
[numHops,numFeatures] = size(features)

features = (features' - normFactors.Mean) ./ normFactors.Variance;

features = features - mean(features,'all');

featuresAll = {};
if ~isempty(ver('parallel'))
    numPar = 18;
else
    numPar = 1;
end

parfor ii = 1:numPar
    adsPart = partition(ads,numPar,ii);
    featuresPart = cell(0,numel(adsPart.Files));
    for iii = 1:numel(adsPart.Files)
        audioData = read(adsPart);
        featuresPart{iii} = helperFeatureExtraction(audioData,afe,[]);
    end
    featuresAll = [featuresAll,featuresPart];
end

allFeatures = cat(2,featuresAll{:});

normFactors.Mean = mean(allFeatures,2,'omitnan');
normFactors.STD = std(allFeatures,[],2,'omitnan');

numComponents =32;
alpha = таковые (1,numComponents )/numComponents;

mu = randn(numFeatures,numComponents);
sigma = rand(numFeatures,numComponents);
ubm = struct('ComponentProportion',alpha,'mu',mu,'sigma',sigma);

maxIter = 20;
targetLogLikelihood = 0;
tol = 0.5;
pastL = -inf; % initialization of previous log-likelihood

tic
for iter = 1:maxIter
    
    % EXPECTATION
    N = zeros(1,numComponents);
    F = zeros(numFeatures,numComponents);
    S = zeros(numFeatures,numComponents);
    L = 0;
    parfor ii = 1:numPar
        adsPart = partition(adsTrainUBM,numPar,ii);
        while hasdata(adsPart)
            audioData = read(adsPart);
            
            % Extract features
            features = helperFeatureExtraction(audioData,afe,normFactors);
 
            % Compute a posteriori log-likelihood
            logLikelihood = helperGMMLogLikelihood(features,ubm);

            % Compute a posteriori normalized probability
            logLikelihoodSum = helperLogSumExp(logLikelihood);
            gamma = exp(logLikelihood - logLikelihoodSum)';
            
            % Compute Baum-Welch statistics
            n = sum(gamma,1);
            f = features * gamma;
            s = (features.*features) * gamma;
            
            % Update the sufficient statistics over utterances
            N = N + n;
            F = F + f;
            S = S + s;
            
            % Update the log-likelihood
            L = L + sum(logLikelihoodSum);
        end
    end
    
    % Print current log-likelihood and stop if it meets criteria.
    L = L/numel(adsTrainUBM.Files);
    fprintf('\tIteration %d, Log-likelihood = %0.3f\n',iter,L)
    if L > targetLogLikelihood || abs(pastL - L) < tol
        break
    else
        pastL = L;
    end
    
    % MAXIMIZATION
    N = max(N,eps);
    ubm.ComponentProportion = max(N/sum(N),eps);
    ubm.ComponentProportion = ubm.ComponentProportion/sum(ubm.ComponentProportion);
    ubm.mu = bsxfun(@rdivide,F,N);
    ubm.sigma = max(bsxfun(@rdivide,S,N) - ubm.mu.^2,eps);
end

fprintf('UBM training completed in %0.2f seconds.\n',toc)

relevanceFactor = 16;

speakers = unique(adsEnroll.Labels);
numSpeakers = numel(speakers);

gmmCellArray = cell(numSpeakers,1);
tic
parfor ii = 1:numSpeakers
    % Subset the datastore to the speaker you are adapting.
    adsTrainSubset = subset(adsEnroll,adsEnroll.Labels==speakers(ii));
    
    N = zeros(1,numComponents);
    F = zeros(numFeatures,numComponents);
    S = zeros(numFeatures,numComponents);
    while hasdata(adsTrainSubset)
        audioData = read(adsTrainSubset);
        features = helperFeatureExtraction(audioData,afe,normFactors);
        [n,f,s,l] = helperExpectation(features,ubm);
        N = N + n;
        F = F + f;
        S = S + s;
    end
    
    % Determine the maximum likelihood
    gmm = helperMaximization(N,F,S);
    
    % Determine adaption coefficient
    alpha = N ./ (N + relevanceFactor);
    
    % Adapt the means
    gmm.mu = alpha.*gmm.mu + (1-alpha).*ubm.mu;
    
    % Adapt the variances
    gmm.sigma = alpha.*(S./N) + (1-alpha).*(ubm.sigma + ubm.mu.^2) - gmm.mu.^2;
    gmm.sigma = max(gmm.sigma,eps);
    
    % Adapt the weights
    gmm.ComponentProportion = alpha.*(N/sum(N)) + (1-alpha).*ubm.ComponentProportion;
    gmm.ComponentProportion = gmm.ComponentProportion./sum(gmm.ComponentProportion);

    gmmCellArray{ii} = gmm;
end
fprintf('Enrollment completed in %0.2f seconds.\n',toc)

for i = 1:numel(gmmCellArray)
    enrolledGMMs.(string(speakers(i))) = gmmCellArray{i};
end

speakers = unique(adsEnroll.Labels);
numSpeakers = numel(speakers);
llr = cell(numSpeakers,1);
tic
parfor speakerIdx = 1:numSpeakers
    localGMM = enrolledGMMs.(string(speakers(speakerIdx))); 
    adsTestSubset = subset(adsTest,adsTest.Labels==speakers(speakerIdx));
    llrPerSpeaker = zeros(numel(adsTestSubset.Files),1);
    for fileIdx = 1:numel(adsTestSubset.Files)
        audioData = read(adsTestSubset);
        [x,numFrames] = helperFeatureExtraction(audioData,afe,normFactors);
        
        logLikelihood = helperGMMLogLikelihood(x,localGMM);
        Lspeaker = helperLogSumExp(logLikelihood);
        
        logLikelihood = helperGMMLogLikelihood(x,ubm);
        Lubm = helperLogSumExp(logLikelihood);
        
        llrPerSpeaker(fileIdx) = mean(movmedian(Lspeaker - Lubm,3));
    end
    llr{speakerIdx} = llrPerSpeaker;
end
fprintf('False rejection rate computed in %0.2f seconds.\n',toc)

llr = cat(1,llr{:});

thresholds = -0.5:0.01:2.5;
FRR = mean(llr<thresholds);

plot(thresholds,FRR*100)
title('False Rejection Rate (FRR)')
xlabel('Threshold')
ylabel('Incorrectly Rejected (%)')
grid on

speakersTest = unique(adsTest.Labels);
llr = cell(numSpeakers,1);
tic
parfor speakerIdx = 1:numel(speakers)
    localGMM = enrolledGMMs.(string(speakers(speakerIdx)));
    adsTestSubset = subset(adsTest,adsTest.Labels~=speakers(speakerIdx));
    llrPerSpeaker = zeros(numel(adsTestSubset.Files),1);
    for fileIdx = 1:numel(adsTestSubset.Files)
        audioData = read(adsTestSubset);
        [x,numFrames] = helperFeatureExtraction(audioData,afe,normFactors);
        
        logLikelihood = helperGMMLogLikelihood(x,localGMM);
        Lspeaker = helperLogSumExp(logLikelihood);
        
        logLikelihood = helperGMMLogLikelihood(x,ubm);
        Lubm = helperLogSumExp(logLikelihood);
        
        llrPerSpeaker(fileIdx) = mean(movmedian(Lspeaker - Lubm,3));
    end
    llr{speakerIdx} = llrPerSpeaker;
end
fprintf('FAR computed in %0.2f seconds.\n',toc)

llr = cat(1,llr{:});

FAR = mean(llr>thresholds);

plot(thresholds,FAR*100)
title('False Acceptance Rate (FAR)')
xlabel('Threshold')
ylabel('Incorrectly Rejected (%)')
grid on

x1 = FAR*100;
y1 = FRR*100;
plot(x1,y1)
grid on
xlabel('False Acceptance Rate (%)')
ylabel('False Rejection Rate (%)')
title('Detection Error Tradeoff (DET) Curve')

[~,EERThresholdIdx] = min(abs(FAR - FRR));
EERThreshold = thresholds(EERThresholdIdx);
EER = mean([FAR(EERThresholdIdx),FRR(EERThresholdIdx)]);
plot(thresholds,FAR,'k', ...
     thresholds,FRR,'b', ...
     EERThreshold,EER,'ro','MarkerFaceColor','r')
title(sprintf('Equal Error Rate = %0.2f, Threshold = %0.2f',EER,EERThreshold))
xlabel('Threshold')
ylabel('Error Rate')
legend('False Acceptance Rate (FAR)','False Rejection Rate (FRR)','Equal Error Rate (EER)')
grid on

resave = false;
if повторно спасти
    сохранить'speakerVerificationExampleData.mat','ubm','afe','normFactors')
end

function helperAddUser(fs,numToRecord,ID)
% Create an audio device reader to read from your audio device
deviceReader = audioDeviceReader('SampleRate',fs);

% Initialize variables
numRecordings = 1;
audioIn = [];

% Record the requested number
while numRecordings <= numToRecord
    fprintf('Say "stop" once (recording %i of %i) ...',numRecordings,numToRecord)
    tic
    while toc<2
        audioIn = [audioIn;deviceReader()];
    end
    fprintf('complete.\n')
    idx = detectSpeech(audioIn,fs);
    if isempty(idx)
        fprintf('Speech not detected. Try again.\n')
    else
        audiowrite(sprintf('%s_%i.flac',ID,numRecordings),audioIn,fs)
        numRecordings = numRecordings+1;
    end
    pause(0.2)
    audioIn = [];
end

% Release the device
release(deviceReader)
end

function speakerGMM = helperEnroll(ubm,afe,normFactors,adsEnroll)
% Initialization
numComponents = numel(ubm.ComponentProportion);
numFeatures = size(ubm.mu,1);
N = zeros(1,numComponents);
F = zeros(numFeatures,numComponents);
S = zeros(numFeatures,numComponents);
NumFrames = 0;

while hasdata(adsEnroll)
    % Read from the enrollment datastore
    audioData = read(adsEnroll);

    % 1. Extract the features and apply feature normalization
    [features,numFrames] = helperFeatureExtraction(audioData,afe,normFactors);
    
    % 2. Calculate the a posteriori probability. Use it to determine the
    % sufficient statistics (the count, and the first and second moments)
    [n,f,s] = helperExpectation(features,ubm);
    
    % 3. Update the sufficient statistics
    N = N + n;
    F = F + f;
    S = S + s;
    NumFrames = NumFrames + numFrames;
end
% Create the Gaussian mixture model that maximizes the expectation
speakerGMM = helperMaximization(N,F,S);

% Adapt the UBM to create the speaker model. Use a relevance factor of 16,
% as proposed in [2]
relevanceFactor = 16;

% Determine adaption coefficient
alpha = N ./ (N + relevanceFactor);

% Adapt the means
speakerGMM.mu = alpha.*speakerGMM.mu + (1-alpha).*ubm.mu;

% Adapt the variances
speakerGMM.sigma = alpha.*(S./N) + (1-alpha).*(ubm.sigma + ubm.mu.^2) - speakerGMM.mu.^2;
speakerGMM.sigma = max(speakerGMM.sigma,eps);

% Adapt the weights
speakerGMM.ComponentProportion = alpha.*(N/sum(N)) + (1-alpha).*ubm.ComponentProportion;
speakerGMM.ComponentProportion = speakerGMM.ComponentProportion./sum(speakerGMM.ComponentProportion);
end

function verificationStatus = helperVerify(audioData,afe,normFactors,speakerGMM,ubm,threshold)
    % Extract features
    x = helperFeatureExtraction(audioData,afe,normFactors);
    
    % Determine the log-likelihood the audio came from the GMM adapted to
    % the speaker
    post = helperGMMLogLikelihood(x,speakerGMM);
    Lspeaker = helperLogSumExp(post);
    
    % Determine the log-likelihood the audio came form the GMM fit to all
    % speakers
    post = helperGMMLogLikelihood(x,ubm);
    Lubm = helperLogSumExp(post);
    
    % Calculate the ratio for all frames. Apply a moving median filter
    % to remove outliers, and then take the mean across the frames
    llr = mean(movmedian(Lspeaker - Lubm,3));

    if llr > threshold
        verificationStatus = true;
    else
        verificationStatus = false;
    end
end

function [features,numFrames] = helperFeatureExtraction(audioData,afe,normFactors)
    % Normalize
    audioData = audioData/max(abs(audioData(:)));
    
    % Protect against NaNs
    audioData(isnan(audioData)) = 0;
    
    % Isolate speech segment
    % The dataset used in this example has one word per audioData, if more
    % than one is speech section is detected, just use the longest
    % detected.
    idx = detectSpeech(audioData,afe.SampleRate);
    if size(idx,1)>1
        [~,seg] = max(idx(:,2) - idx(:,1));
    else
        seg = 1;
    end
    audioData = audioData(idx(seg,1):idx(seg,2));
    
    % Feature extraction
    features = extract(afe,audioData);

    % Feature normalization
    if ~isempty(normFactors)
        features = (features-normFactors.Mean')./normFactors.STD';
    end
    features = features';
    
    % Cepstral mean subtraction (for channel noise)
    if ~isempty(normFactors)
        features = features - mean(features,'all');
    end
    
    numFrames = size(features,2);
end

function y = helperLogSumExp(x)
% Calculate the log-sum-exponent while avoiding overflow
a = max(x,[],1);
y = a + sum(exp(bsxfun(@minus,x,a)),1);
end

function [N,F,S,L] = helperExpectation(features,gmm)

post = helperGMMLogLikelihood(features,gmm);

% Sum the likelihood over the frames
L = helperLogSumExp(post);

% Compute the sufficient statistics
gamma = exp(post-L)';

N = sum(gamma,1);
F = features * gamma;
S = (features.*features) * gamma;
L = sum(L);
end

function gmm = helperMaximization(N,F,S)
    N = max(N,eps);
    gmm.ComponentProportion = max(N/sum(N),eps);
    gmm.mu = bsxfun(@rdivide,F,N);
    gmm.sigma = max(bsxfun(@rdivide,S,N) - gmm.mu.^2,eps);
end

function L = helperGMMLogLikelihood(x,gmm)
    xMinusMu = repmat(x,1,1,numel(gmm.ComponentProportion)) - permute(gmm.mu,[1,3,2]);
    permuteSigma = permute(gmm.sigma,[1,3,2]);
    
    Lunweighted = -0.5*(sum(log(permuteSigma),1) + sum(bsxfun(@times,xMinusMu,(bsxfun(@rdivide,xMinusMu,permuteSigma))),1) + size(gmm.mu,1)*log(2*pi));

    temp = squeeze(permute(Lunweighted,[1,3,2]));
    if size(temp,1)==1
        % If there is only one frame, the trailing singleton dimension was
        % removed in the permute. This accounts for that edge case
        temp = temp';
    end
    L = bsxfun(@plus,temp,log(gmm.ComponentProportion)');
end