Верификация динамика с использованием i-векторов

Верификация типа «динамик» или аутентификация - это задача подтверждения того, что тождества динамика являются тем, кем они считаются. Верификация спикера была активной исследовательской областью на протяжении многих лет. Ранним прорывом в эффективности было использование смешанной гауссовской модели и универсальной фоновой модели (GMM-UBM) [1] на акустических функциях (обычно mfcc). Для получения примера смотрите Динамик Верификации Использование Смешанной гауссовской модели. Одна из основных трудностей систем GMM-UBM заключается в межсессионной изменчивости. Совместный факторный анализ (JFA) был предложен, чтобы компенсировать эту изменчивость, отдельно моделируя переменность между динамиками и изменчивость канала или сеанса [2] [3]. Однако [4] обнаружил, что канальные факторы в JFA также содержат информацию о динамиках, и предложил объединить пространства канала и динамика в общее пространство переменности. Затем интерсессионную изменчивость компенсировали с помощью процедур бэкэнда, таких как линейный дискриминантный анализ (LDA) и ковариационная нормализация внутри класса (WCCN), с последующей оценкой, такой как счет подобия косинуса. [5] предложил заменить оценку подобия косинуса вероятностной моделью LDA (PLDA). [11] и [12] предложили способ Гауссианизации i-векторов и, следовательно, сделать Гауссовы предположения в PLDA, называемые G-PLDA или упрощенной PLDA. В то время как i-векторы были первоначально предложены для верификации диктора, они были применены ко многим проблемам, таким как распознавание языка, диаризация диктора, распознавание эмоций, оценка возраста и анти-спуфинг [10]. Недавно были предложены методы глубокого обучения, чтобы заменить i-векторы на d-векторы или x-векторы [8] [6].

Использование i-векторной системы

Audio Toolbox обеспечивает ivectorSystem который инкапсулирует возможность обучать i-векторную систему, регистрировать динамики или другие звуковые метки, оценивать систему на предмет порога принятия решений и идентифицировать или проверить динамики или другие аудио метки. См. ivectorSystem для примеров использования этой функции и ее применения в нескольких приложениях.

Чтобы узнать больше о том, как работает i-векторная система, продолжите с примером.

Разрабатывайте i-векторную систему

В этом примере вы разрабатываете стандартную систему i-векторов для верификации динамика, которая использует бэкэнд LDA-WCCN с оценкой косинусоидальности или оценкой G-PLDA.

На протяжении всего примера вы найдете живые элементы управления настраиваемыми параметрами. Изменение элементов управления не приводит к повторному запуску примера. Если вы меняете элемент управления, необходимо перезапустить пример.

Управление набором данных

Этот пример использует базу данных отслеживания тангажа из Технологического университета Граца (PTDB-TUG) [7]. Набор данных состоит из 20 носителей английского языка, читающих 2342 фонетически богатых предложения из корпуса TIMIT. Загрузите и извлеките набор данных. В зависимости от вашей системы, загрузка и извлечение набора данных может занять приблизительно 1,5 часов.

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

if ~exist(datasetFolder,'dir')
    disp('Downloading PTDB-TUG (3.9 G) ...')
    unzip(url,datasetFolder)
end
Downloading PTDB-TUG (3.9 G) ...

Создайте audioDatastore объект, который указывает на набор данных. Набор данных первоначально предназначался для использования в обучении и оценке отслеживания основного тона и включает показания ларингографика и решения о базовом тангаже. Используйте только оригинальные аудиозаписи.

ads = audioDatastore([fullfile(datasetFolder,"SPEECH DATA","FEMALE","MIC"),fullfile(datasetFolder,"SPEECH DATA","MALE","MIC")], ...
                     'IncludeSubfolders',true, ...
                     'FileExtensions','.wav');
fileNames = ads.Files;

Имена файлов содержат идентификаторы динамиков. Декодируйте имена файлов, чтобы задать метки на audioDatastore объект.

speakerIDs = extractBetween(fileNames,'mic_','_');
ads.Labels = categorical(speakerIDs);
countEachLabel(ads)
ans=20×2 table
    Label    Count
    _____    _____

     F01      236 
     F02      236 
     F03      236 
     F04      236 
     F05      236 
     F06      236 
     F07      236 
     F08      234 
     F09      236 
     F10      236 
     M01      236 
     M02      236 
     M03      236 
     M04      236 
     M05      236 
     M06      236 
      ⋮

Разделите audioDatastore объект в наборы для обучения, оценки и тестирования. В набор обучающих данных 16 динамиков. Оценочный набор содержит четыре динамика и далее разделен на набор регистрации и набор для оценки компромисса ошибок обнаружения обученной i-векторной системы и тестовый набор.

developmentLabels = categorical(["M01","M02","M03","M04","M06","M07","M08","M09","F01","F02","F03","F04","F06","F07","F08","F09"]);
evaluationLabels = categorical(["M05","M010","F05","F010"]);

adsTrain = subset(ads,ismember(ads.Labels,developmentLabels));

adsEvaluate = subset(ads,ismember(ads.Labels,evaluationLabels));
numFilesPerSpeakerForEnrollment = 3;
[adsEnroll, adsTest, adsDET] = splitEachLabel (adsEvaluate, numFilesPerSpeakerForEnrollment, 2);

Отображение распределений меток полученной audioDatastore объекты.

countEachLabel(adsTrain)
ans=16×2 table
    Label    Count
    _____    _____

     F01      236 
     F02      236 
     F03      236 
     F04      236 
     F06      236 
     F07      236 
     F08      234 
     F09      236 
     M01      236 
     M02      236 
     M03      236 
     M04      236 
     M06      236 
     M07      236 
     M08      236 
     M09      236 

countEachLabel(adsEnroll)
ans=2×2 table
    Label    Count
    _____    _____

     F05       3  
     M05       3  

countEachLabel(adsDET)
ans=2×2 table
    Label    Count
    _____    _____

     F05      231 
     M05      231 

countEachLabel(adsTest)
ans=2×2 table
    Label    Count
    _____    _____

     F05       2  
     M05       2  

Прочтите аудио файла с обучающими данными набора, послушайте его и постройте график. Сбросьте datastore.

[audio,audioInfo] = read(adsTrain);
fs = audioInfo.SampleRate;

t = (0:size(audio,1)-1)/fs;
sound(audio,fs)
plot(t,audio)
xlabel('Time (s)')
ylabel('Amplitude')
axis([0 t(end) -1 1])
title('Sample Utterance from Training Set')

reset(adsTrain)

Можно уменьшить набор данных и количество параметров, используемых в этом примере для ускорения выполнения за счет эффективности. В целом сокращение набора данных является хорошей практикой для разработки и отладки.

speedUpExample = false;
if speedUpExample
    adsTrain = splitEachLabel (adsTrain, 30);
    adsDET = splitEachLabel (adsDET, 21);
end

Редукция данных

Создайте audioFeatureExtractor Целью является извлечение 20 MFCC, 20 дельта-MFCC и 20 дельта-дельта MFCC. Используйте дельта-окно длиной 9. Извлечь функции из 25 мс окон Ханна с 10 мс хмеля.

numCoeffs = 20;
deltaWindowLength = 9;

windowDuration = 0.025;
hopDuration = 0.01;

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

afe = audioFeatureExtractor  (...
    'SampleRate', фс ,...
    'Window', hann (windowSamples,'periodic'), ...
    'OverlapLength', overlapSamples ,...
    ...
    'mfcc'Правда, ...
    'mfccDelta'Правда, ...
    'mfccDeltaDelta'Правда;
setExtractorParams (afe,'mfcc','DeltaWindowLength', deltaWindowLength,'NumCoeffs', numCoeffs)

Извлеките функции из аудио, считанного из обучающего datastore. Функции возвращаются как numHops-by- numFeatures матрица.

features = extract(afe,audio);
[numHops,numFeatures] = size(features)
numHops = 797
numFeatures = 60

Обучение

Обучение i-векторной системы является в вычислительном отношении дорогим и длительным. Если у вас есть Parallel Computing Toolbox™, можно распределить работу по нескольким ядрам, чтобы ускорить пример. Определите оптимальное количество разделов для вашей системы. Если у вас нет Parallel Computing Toolbox™, используйте один раздел.

if ~isempty(ver('parallel')) && ~speedUpExample
    pool = gcp;
    numPar = numpartitions(adsTrain,pool);
else
    numPar = 1;
end
Starting parallel pool (parpool) using the 'local' profile ...
Connected to the parallel pool (number of workers: 6).

Коэффициенты нормализации функций

Используйте функцию helper, helperFeatureExtraction, чтобы извлечь все функции из набора данных. Функция helperFeatureExtraction извлекает MFCC из областей речи в аудио. Обнаружение речи выполняется detectSpeech функция.

featuresAll = {};
tic
parfor ii = 1:numPar
    adsPart = partition(adsTrain,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{:});
fprintf('Feature extraction from training set complete (%0.0f seconds).',toc)
Feature extraction from training set complete (58 seconds).

Вычислите глобальное среднее и стандартное отклонение каждой функции. Вы будете использовать их в будущих вызовах helperFeatureExtraction функция для нормализации функций.

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

Универсальная фоновая модель (UBM)

Инициализируйте Смешанную гауссовскую модель (GMM), которая будет универсальной фоновой моделью (UBM) в i-векторной системе. Веса компонентов инициализируются как равномерно распределенные. Системы, обученные на наборе данных TIMIT, обычно содержат около 2048 компонентов.

numComponents = 64;
if speedUpExample
    numComponents = 32;
end
alpha = таковые (1,numComponents )/numComponents;
mu = randn (numFeatures, numComponents);
vari = rand (numFeatures, numComponents) + eps;
ubm = struct ('ComponentProportion', альфа,'mu', му,'sigma', vari);

Обучите UBM с помощью алгоритма максимизации ожиданий (EM).

maxIter = 10;
if speedUpExample
    maxIter = 2;
end
tic
for iter = 1:maxIter
    tic
    % EXPECTATION
    N = нули (1, numComponents);
    F = нули (numFeatures, numComponents);
    S = нули (numFeatures, numComponents);
    L = 0;
    parfor ii = 1:numPar
        adsPart = раздел (adsTrain, numPar, ii);
        while hasdata (adsPart)
            audioData = read (adsPart);
            
            % Extract features
            Y = helperFeatureExtraction (audioData, afe, normFactors);
 
            % Compute a posteriori log-liklihood
            логарифмическая правдоподобность = helperGMMLogLikelihood (Y, ubm);

            % Compute a posteriori normalized probability
            amax = max (логарифмическая правдоподобность, [], 1);
            logLikelihoodSum = amax + журнал (сумма (exp (logLikelihood-amax), 1));
            gamma = exp (логарифмическая правдоподобность - logLikelihoodSum) ';
            
            % Compute Baum-Welch statistics
            n = сумма (гамма, 1);
            f = Y * гамма;
            s = (Y. * Y) * гамма;
            
            % 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
    fprintf ('Training UBM: %d/%d complete (%0.0f seconds), Log-likelihood = %0.0f\n', iter, maxIter, toc, L)
    
    % MAXIMIZATION
    N = max (N, eps);
    ubm.ComponentProportion = max (N/sum (N), eps);
    ubm.ComponentПропорция = ubm. ComponentProportion/sum (ubm.ComponentProportion);
    ubm.mu = F./N;
    ubm.sigma = max (S./N - ubm.mu ^ 2, eps);
end
Training UBM: 1/10 complete (59 seconds), Log-likelihood = -162907120
Training UBM: 2/10 complete (53 seconds), Log-likelihood = -82282814
Training UBM: 3/10 complete (54 seconds), Log-likelihood = -78667384
Training UBM: 4/10 complete (55 seconds), Log-likelihood = -77041863
Training UBM: 5/10 complete (54 seconds), Log-likelihood = -76338342
Training UBM: 6/10 complete (52 seconds), Log-likelihood = -75958218
Training UBM: 7/10 complete (52 seconds), Log-likelihood = -75724712
Training UBM: 8/10 complete (53 seconds), Log-likelihood = -75561701
Training UBM: 9/10 complete (54 seconds), Log-likelihood = -75417170
Training UBM: 10/10 complete (55 seconds), Log-likelihood = -75275185

Вычисление статистики Баума-Уэлча

Статистика Баума-Уэлча N (нулевой порядок) и F (первый порядок) статистики, используемой в алгоритме EM, вычисленной с использованием конечного UBM.

Nc(s)=tγt(c)

Fc(s)=tγt(c)Yt

  • Yt - вектор функции в момент времени t.

  • s{s1,s2,...,sN}, где N - количество дикторов. Для целей обучения общего пространства переменности каждый звуковой файл рассматривается как отдельный динамик (принадлежит ли он физическому одному динамику).

  • γt(c) - апостериорная вероятность того, что компонент UBM c учитывает вектор функции Yt.

Вычислите нули и статистику Баума-Уэлча первого порядка по набору обучающих данных.

numSpeakers = numel(adsTrain.Files);
Nc = {};
Fc = {};

tic
parfor ii = 1:numPar
    adsPart = partition(adsTrain,numPar,ii);
    numFiles = numel(adsPart.Files);
    
    Npart = cell(1,numFiles);
    Fpart = cell(1,numFiles);
    for jj = 1:numFiles
        audioData = read(adsPart);
        
        % Extract features
        Y = helperFeatureExtraction(audioData,afe,normFactors);
        
        % Compute a posteriori log-likelihood
        logLikelihood = helperGMMLogLikelihood(Y,ubm);
        
        % Compute a posteriori normalized probability
        amax = max(logLikelihood,[],1);
        logLikelihoodSum = amax + log(sum(exp(logLikelihood-amax),1));
        gamma = exp(logLikelihood - logLikelihoodSum)';
        
        % Compute Baum-Welch statistics
        n = sum(gamma,1);
        f = Y * gamma;
        
        Npart{jj} = reshape(n,1,1,numComponents);
        Fpart{jj} = reshape(f,numFeatures,1,numComponents);
    end
    Nc = [Nc,Npart];
    Fc = [Fc,Fpart];
end
fprintf('Baum-Welch statistics completed (%0.0f seconds).\n',toc)
Baum-Welch statistics completed (54 seconds).

Разверните статистику в матрицы и центр F(s), как описано в [3], таким образом, что

  • N(s) является CF×CF диагональная матрица, чьи блоки Nc(s)I(c=1,...C).

  • F(s) является CF×1 супервектор, полученный путем конкатенации Fc(s)  (c=1,...C).

  • C количество компонентов в UBM.

  • F - количество функций в векторе функций.

N = Nc;
F = Fc;
muc = reshape(ubm.mu,numFeatures,1,[]);
for s = 1:numSpeakers
    N{s} = repelem(reshape(Nc{s},1,[]),numFeatures);
    F{s} = reshape(Fc{s} - Nc{s}.*muc,[],1);
end

Поскольку этот пример принимает диагональную ковариационную матрицу для UBM, N являются также диагональными матрицами и сохраняются как векторы для эффективных расчетов.

Общее пространство переменностей

В модели i-вектора идеальный супервектор динамика состоит из независимого от динамика компонента и зависящего от динамика компонента. Зависящий от динамика компонент состоит из модели общего пространства переменностей и i-вектора динамика.

M=m+Tw

  • M является ли оратор высказыванием супервектора

  • m является независимым от громкой связи и канала супервектором, который может быть взят за супервектор UBM.

  • T является низкоуровневой прямоугольной матрицей и представляет общий подпространство переменности.

  • w - i-вектор для динамика

Размерность i-вектора, w, обычно намного ниже, чем C F-мерный супервектор громкоговорителя, что делает i-вектор или i-векторы намного более компактным и отслеживаемым представлением.

Чтобы обучить общее пространство переменностей, T, сначала случайным образом инициализируйте Tзатем выполните эти шаги итерационно [3]:

  1. Вычислите апостериорное распределение скрытой переменной.

lT(s)=I+T×Σ-1×N(s)×T

2. Накапливайте статистику по динамикам.

Κ=sF(s)×(lT-1(s)×T×Σ-1×F(s))

Ac=sNc(s)lT-1(s)

3. Обновите общее пространство переменностей.

Tc=Ac-1×Κ

T=[T1T2TC]

[3] предлагает инициализировать Σ по отклонению UBM, а затем по обновлению Σ согласно уравнению:

Σ=(sN(s))-1((sS(s))-diag(Κ×T))

где S (s) - центрированная статистика Баума-Уэлча второго порядка. Однако обновлениеΣ часто падает на практике, поскольку это имеет незначительный эффект. Этот пример не обновляется Σ.

Создайте переменную sigma.

Sigma = ubm.sigma(:);

Задайте размерность общего пространства переменностей. Типичное значение, используемое для набора данных TIMIT, составляет 1000.

numTdim = 32;
if speedUpExample
    numTdim = 16;
end

Инициализация T и матрицу тождеств и предварительное выделение массивов ячеек.

T = randn(numel(ubm.sigma),numTdim);
T = T/norm(T);

I = eye(numTdim);

Ey = cell(numSpeakers,1);
Eyy = cell(numSpeakers,1);
Linv = cell(numSpeakers,1);

Установите количество итераций для обучения. Типичное значение, сообщаемое как 20.

numIterations = 5;

Запустите цикл обучения.

for iterIdx = 1:numIterations
    tic
    
    % 1. Calculate the posterior distribution of the hidden variable
    TtimesInverseSSdiag = (T./Sigma)';
    parfor s = 1:numSpeakers
        L = (I + TtimesInverseSSdiag.*N{s}*T);
        Linv{s} = pinv(L);
        Ey{s} = Linv{s}*TtimesInverseSSdiag*F{s};
        Eyy{s} = Linv{s} + Ey{s}*Ey{s}';
    end
    
    % 2. Accumlate statistics across the speakers
    Eymat = cat(2,Ey{:});
    FFmat = cat(2,F{:});
    Kt = FFmat*Eymat';
    K = mat2cell(Kt',numTdim,repelem(numFeatures,numComponents));
    
    newT = cell(numComponents,1);
    for c = 1:numComponents
        AcLocal = zeros(numTdim);
        for s = 1:numSpeakers
            AcLocal = AcLocal + Nc{s}(:,:,c)*Eyy{s};
        end
        
    % 3. Update the Total Variability Space
        newT{c} = (pinv(AcLocal)*K{c})';
    end
    T = cat(1,newT{:});

    fprintf('Training Total Variability Space: %d/%d complete (%0.0f seconds).\n',iterIdx,numIterations,toc)
end
Training Total Variability Space: 1/5 complete (2 seconds).
Training Total Variability Space: 2/5 complete (2 seconds).
Training Total Variability Space: 3/5 complete (2 seconds).
Training Total Variability Space: 4/5 complete (1 seconds).
Training Total Variability Space: 5/5 complete (1 seconds).

i-векторное выделение

После вычисления общего пространства переменностей можно вычислить i-векторы как [4]:

w=(I+TΣ-1NT)TΣ-1F

На данной точке вы все еще рассматриваете каждый обучающий файл как отдельный динамик. Однако на следующем шаге, когда вы обучаете проекционную матрицу, чтобы уменьшить размерность и увеличить межзвуковые различия, i-векторы должны быть помечены соответствующими отдельными идентификаторами динамика.

Создайте массив ячеек, где каждый элемент массива ячеек содержит матрицу i-векторов через файлы для конкретного динамика.

speakers = unique(adsTrain.Labels);
numSpeakers = numel(speakers);
ivectorPerSpeaker = cell(numSpeakers,1);
TS = T./Sigma;
TSi = TS';
ubmMu = ubm.mu;
tic
parfor speakerIdx = 1:numSpeakers
    
    % Subset the datastore to the speaker you are adapting.
    adsPart = subset(adsTrain,adsTrain.Labels==speakers(speakerIdx));
    numFiles = numel(adsPart.Files);
    
    ivectorPerFile = zeros(numTdim,numFiles);
    for fileIdx = 1:numFiles
        audioData = read(adsPart);
        
        % Extract features
        Y = helperFeatureExtraction(audioData,afe,normFactors);
        
        % Compute a posteriori log-likelihood
        logLikelihood = helperGMMLogLikelihood(Y,ubm);
        
        % Compute a posteriori normalized probability
        amax = max(logLikelihood,[],1);
        logLikelihoodSum = amax + log(sum(exp(logLikelihood-amax),1));
        gamma = exp(logLikelihood - logLikelihoodSum)';
        
        % Compute Baum-Welch statistics
        n = sum(gamma,1);
        f = Y * gamma - n.*(ubmMu);

        ivectorPerFile(:,fileIdx) = pinv(I + (TS.*repelem(n(:),numFeatures))' * T) * TSi * f(:);
    end
    ivectorPerSpeaker{speakerIdx} = ivectorPerFile;
end
fprintf('I-vectors extracted from training set (%0.0f seconds).\n',toc)
I-vectors extracted from training set (60 seconds).

Матрица проекций

Для i-векторов было предложено много различных бэкендов. Наиболее простым и все еще хорошо работающим является комбинация линейного дискриминантного анализа (LDA) и ковариационной нормализации внутри класса (WCCN).

Создайте матрицу обучающих векторов и карту, указывающую, какой i-вектор соответствует какому динамику. Инициализируйте проекционную матрицу как тождества матрицу.

w = ivectorPerSpeaker;
utterancePerSpeaker = cellfun(@(x)size(x,2),w);

ivectorsTrain = cat(2,w{:});
projectionMatrix = eye(size(w{1},1));

LDA пытается минимизировать отклонение внутри внутрикласса и максимизировать отклонение между динамиками. Он может быть вычислен как изложено в [4]:

Учитывая:

Sb=s=1S(ws-w)(ws-w)

Sw=s=1S1nsi=1ns(wis-ws)(wis-ws)

где

  • ws=(1ns)i=1nswis - среднее значение i-векторов для каждого динамика.

  • w=1Ns=1Si=1nswis - средний вектор i для всех динамиков.

  • ns количество высказываний для каждого оратора.

Решить собственное уравнение для лучших собственных векторов:

Sbv=λSwv

Лучшими собственными векторами являются те, у которых самые высокие собственные значения.

performLDA = true;
if performLDA
    tic
    numEigenvectors = 16;

    Sw = нули (размер (projectionMatrix, 1));
    Sb = нули (размер (projectionMatrix, 1));
    wbar = среднее (кошка (2, w {:}), 2);
    for ii = 1: numel (w)
        ws = w {ii};
        wsbar = среднее (ws, 2);
        Sb = Sb + (wsbar - wbar) * (wsbar - wbar) ';
        Sw = Sw + co (ws ', 1);
    end
    
    [A, ~] = eigs (Sb, Sw, numEigenvectors);
    A = (A./vecnorm (A)) ';

    ivectorsTrain = A * ivectorsTrain;
    
    w = mat2cell (ivectorsTrain, size (ivectorsTrain, 1), adterancePerSpeaker);
    
    projectionMatrix = A * projectionMatrix;
    
    fprintf ('LDA projection matrix calculated (%0.2f seconds).', toc)
end
LDA projection matrix calculated (0.22 seconds).

WCCN пытается масштабировать i-векторное пространство обратно к ковариации в классе, так что направления высокой изменчивости внутри динамика отменяются в i-векторных сравнениях [9].

Учитывая ковариационную матрицу внутри класса:

W=1Ss=1S1nsi=1ns(wis-ws)(wis-ws)

где

  • ws=(1ns)i=1nswis - среднее значение i-векторов для каждого динамика.

  • ns количество высказываний для каждого оратора.

Решение для B с использованием разложения Холецкого:

W-1=BB

performWCCN = true;
if performWCCN
    tic
    alpha = 0.9;
    
    W = нули (размер (projectionMatrix, 1));
    for ii = 1: numel (w)
        W = W + cov (w {ii} ', 1);
    end
    W = W/numel (w);
    
    W = (1 - альфа) * W + альфа * глаз (размер (W, 1));
    
    B = хол (pinv (W),'lower');
    
    projectionMatrix = B * projectionMatrix;
    
    fprintf ('WCCN projection matrix calculated (%0.4f seconds).', toc)
end
WCCN projection matrix calculated (0.0063 seconds).

Этап обучения завершен. Теперь можно использовать универсальную фоновую модель (UBM), общее пространство переменности (T) и проекционную матрицу, чтобы зарегистрировать и проверить динамики.

Обучите модель G-PLDA

Примените проекционную матрицу к набору train.

ivectors = cellfun(@(x)projectionMatrix*x,ivectorPerSpeaker,'UniformOutput',false);

Этот алгоритм, реализованный в этом примере, является PLDA Гауссова, как изложено в [13]. В PLDA по Гауссу i-вектор представлен следующим уравнением:

ϕij=μ+Vyi+εij

yiΝ(0,Ι)

εijΝ(0,Λ-1)

где μ- глобальное среднее значение i-векторов, Λ - полная матрица точности шумового термина εij, и V - матрица загрузки фактора, также известная как собственные представления.

Задайте количество собственных значений для использования. Обычно числа находятся между 10 и 400.

numEigenVoices = 16;

Определите количество несвязанных людей, количество размерностей в векторах функций и количество высказываний на говорящего.

K = numel(ivectors);
D = size(ivectors{1},1);
utterancePerSpeaker = cellfun(@(x)size(x,2),ivectors);

Найдите общее количество выборок и центрируйте i-векторы.

N=i=1Kni

μ=1Ni,jϕi,j

φij=ϕij-μ

ivectorsMatrix = cat(2,ivectors{:});
N = size(ivectorsMatrix,2);
mu = mean(ivectorsMatrix,2);

ivectorsMatrix = ivectorsMatrix - mu;

Определите отбеливающую матрицу из обучающих i-векторов и затем отбельте i-векторы. Задайте либо отбеливание ZCA, либо отбеливание PCA, либо отсутствие отбеливания.

whiteningType = 'ZCA';

if strcmpi (beteningType,'ZCA')
    S = cov (ivectorsMatrix ');
    [~, sD, sV] = svd (S);
    W = diag (1 ./( sqrt (diag (sD)) + eps)) * sV ';
    ivectorsMatrix = W * ivectorsMatrix;
elseif strcmpi (beteningType,'PCA')
    S = cov (ivectorsMatrix ');
    [sV, sD] = eig (S);
    W = diag (1 ./( sqrt (diag (sD)) + eps)) * sV ';
    ivectorsMatrix = W * ivectorsMatrix;
else
    W = глаз (размер (ivectorsMatrix, 1));
end

Примените нормализацию длины, а затем преобразуйте обучающую матрицу i-вектора назад в массив ячеек.

ivectorsMatrix = ivectorsMatrix./vecnorm(ivectorsMatrix);

Вычислите глобальный момент второго порядка как

S=ijφijφijT

S = ivectorsMatrix*ivectorsMatrix';

Преобразуйте обучающую матрицу i-вектора назад в массив ячеек.

ivectors = mat2cell(ivectorsMatrix,D,utterancePerSpeaker);

Отсортируйте людей по количеству выборок и затем сгруппируйте i-векторы по количеству высказываний на динамик. Предварительный расчет момента первого порядка i-th person as

fi=j=1niφij

uniqueLengths = unique(utterancePerSpeaker);
numUniqueLengths = numel(uniqueLengths);

speakerIdx = 1;
f = zeros(D,K);
for uniqueLengthIdx = 1:numUniqueLengths
    idx = find(utterancePerSpeaker==uniqueLengths(uniqueLengthIdx));
    temp = {};
    for speakerIdxWithinUniqueLength = 1:numel(idx)
        rho = ivectors(idx(speakerIdxWithinUniqueLength));
        temp = [temp;rho]; %#ok<AGROW>

        f(:,speakerIdx) = sum(rho{:},2);
        speakerIdx = speakerIdx+1;
    end
    ivectorsSorted{uniqueLengthIdx} = temp; %#ok<SAGROW> 
end

Инициализируйте матрицу собственных значений, V и член обратного отклонения шума, Λ.

V = randn(D,numEigenVoices);
Lambda = pinv(S/N);

Задайте количество итераций для алгоритма EM и необходимость применения минимального расхождения.

numIter = 5;
minimumDivergence = true;

Обучите модель G-PLDA с помощью алгоритма EM, описанного в [13].

for iter = 1:numIter
    % EXPECTATION
    gamma = zeros(numEigenVoices,numEigenVoices);
    EyTotal = zeros(numEigenVoices,K);
    R = zeros(numEigenVoices,numEigenVoices);
    
    idx = 1;
    for lengthIndex = 1:numUniqueLengths
        ivectorLength = uniqueLengths(lengthIndex);
        
        % Isolate i-vectors of the same given length
        iv = ivectorsSorted{lengthIndex};
        
        % Calculate M
        M = pinv(ivectorLength*(V'*(Lambda*V)) + eye(numEigenVoices)); % Equation (A.7) in [13]
        
        % Loop over each speaker for the current i-vector length
        for speakerIndex = 1:numel(iv)
            
            % First moment of latent variable for V
            Ey = M*V'*Lambda*f(:,idx); % Equation (A.8) in [13]
            
            % Calculate second moment.
            Eyy = Ey * Ey';
            
            % Update Ryy 
            R = R + ivectorLength*(M + Eyy); % Equation (A.13) in [13]
            
            % Append EyTotal
            EyTotal(:,idx) = Ey;
            idx = idx + 1;
            
            % If using minimum divergence, update gamma.
            if minimumDivergence
                gamma = gamma + (M + Eyy); % Equation (A.18) in [13]
            end
        end
    end
    
    % Calculate T
    TT = EyTotal*f'; % Equation (A.12) in [13]
    
    % MAXIMIZATION
    V = TT'*pinv(R); % Equation (A.16) in [13]
    Lambda = pinv((S - V*TT)/N); % Equation (A.17) in [13]

    % MINIMUM DIVERGENCE
    if minimumDivergence
        gamma = gamma/K; % Equation (A.18) in [13]
        V = V*chol(gamma,'lower');% Equation (A.22) in [13]
    end
end

После обучения модели G-PLDA можно использовать ее для вычисления счета на основе логарифмического коэффициента логарифмической правдоподобности, как описано в [14]. Учитывая два i-вектора, которые были центрированы, отбелены и нормированы к длине, счет вычисляется как:

score(w1,wt)=[w1TwtT][Σ+VVTVVTVVTΣ+VVT][w1wt]-w1T[Σ+VVT]-1w1-wtT[Σ+VVT]-1wt+C

где w1 и wt являются регистрационными и тестовыми i-векторами, Σ - матрица отклонений члена шума, V - матрица собственных значений. C термин является факторизованным и может быть исключен на практике.

speakerIdx = 2;
adteranceIdx = 1;
w1 = ivectors {speakerIdx} (:, reteranceIdx);

speakerIdx = 1;
adteranceIdx = 10;
wt = ivectors {speakerIdx} (:, retteranceIdx);

VVt = V * V ';
Sigma Plus VV t = pinv (Lambda) + VV t;

term1 = pinv ([SigmaPlusVVt VVt; VVt SigmaPlusVVt]);
term2 = pinv (SigmaPlusVVt);

w1wt = [w1; wt];
счет = w1wt '* term1 * w1wt - w1' * term2 * w1 - wt '* term2 * wt
score = 52.4507

В pratice тестовые i-векторы, и в зависимости от вашей системы, ивекторы регистрации, не используются в обучении модели G-PLDA. В следующем разделе оценки для регистрации и верификации используются ранее не просмотренные данные. Поддерживающая функция gpldaScore инкапсулирует шаги оценки выше и дополнительно выполняет центрирование, отбеливание и нормализацию. Сохраните обученную модель G-PLDA как struct для использования с поддерживающей функцией gpldaScore.

gpldaModel = struct('mu',mu, ...
                    'WhiteningMatrix',W, ...
                    'EigenVoices',V, ...
                    'Sigma',pinv(Lambda));

Зарегистрироваться

Зарегистрируйте новых докладчиков, которые не были в наборе обучающих данных.

Создайте i-векторы для каждого файла для каждого динамика в наборе регистрации с помощью этой последовательности шагов:

  1. Редукция данных

  2. Статистика Baum-Welch: Определите нули и статистику первого порядка

  3. i-векторная экстракция

  4. Межсессионная компенсация

Затем усредните i-векторы между файлами, чтобы создать модель i-вектора для динамика. Повторите для каждого динамика.

speakers = unique(adsEnroll.Labels);
numSpeakers = numel(speakers);
enrolledSpeakersByIdx = cell(numSpeakers,1);
tic
parfor speakerIdx = 1:numSpeakers
    % Subset the datastore to the speaker you are adapting.
    adsPart = subset(adsEnroll,adsEnroll.Labels==speakers(speakerIdx));
    numFiles = numel(adsPart.Files);
    
    ivectorMat = zeros(size(projectionMatrix,1),numFiles);
    for fileIdx = 1:numFiles
        audioData = read(adsPart);
        
        % Extract features
        Y = helperFeatureExtraction(audioData,afe,normFactors);
        
        % Compute a posteriori log-likelihood
        logLikelihood = helperGMMLogLikelihood(Y,ubm);
        
        % Compute a posteriori normalized probability
        amax = max(logLikelihood,[],1);
        logLikelihoodSum = amax + log(sum(exp(logLikelihood-amax),1));
        gamma = exp(logLikelihood - logLikelihoodSum)';
        
        % Compute Baum-Welch statistics
        n = sum(gamma,1);
        f = Y * gamma - n.*(ubmMu);
        
        %i-vector Extraction
        w = pinv(I + (TS.*repelem(n(:),numFeatures))' * T) * TSi * f(:);

        % Intersession Compensation
        w = projectionMatrix*w;

        ivectorMat(:,fileIdx) = w;
    end
    % i-vector model
    enrolledSpeakersByIdx{speakerIdx} = mean(ivectorMat,2);
end
fprintf('Speakers enrolled (%0.0f seconds).\n',toc)
Speakers enrolled (0 seconds).

В целях бухгалтерии преобразуйте массив ячеек i-векторов в структуру с идентификаторами динамиков как полями, а i-векторы как значения

enrolledSpeakers = struct;
for s = 1:numSpeakers
    enrolledSpeakers.(string(speakers(s))) = enrolledSpeakersByIdx{s};
end

Верификация

Укажите метод оценки CSS или G-PLDA.

scoringMethod = 'GPLDA';

Ложная частота отклонения (FRR)

Скорость ложного отклонения (FRR) динамика является частотой неправильного отклонения данного динамика. Создайте массив счетов для зарегистрированных векторов динамика и i-векторов того же динамика.

speakersToTest = unique(adsDET.Labels);
numSpeakers = numel(speakersToTest);
scoreFRR = cell(numSpeakers,1);
tic
parfor speakerIdx = 1:numSpeakers
    adsPart = subset(adsDET,adsDET.Labels==speakersToTest(speakerIdx));
    numFiles = numel(adsPart.Files);
    
    ivectorToTest = enrolledSpeakers.(string(speakersToTest(speakerIdx))); %#ok<PFBNS> 
    
    score = zeros(numFiles,1);
    for fileIdx = 1:numFiles
        audioData = read(adsPart);
        
        % Extract features
        Y = helperFeatureExtraction(audioData,afe,normFactors);
        
        % Compute a posteriori log-likelihood
        logLikelihood = helperGMMLogLikelihood(Y,ubm);
        
        % Compute a posteriori normalized probability
        amax = max(logLikelihood,[],1);
        logLikelihoodSum = amax + log(sum(exp(logLikelihood-amax),1));
        gamma = exp(logLikelihood - logLikelihoodSum)';
        
        % Compute Baum-Welch statistics
        n = sum(gamma,1);
        f = Y * gamma - n.*(ubmMu);
        
        % Extract i-vector
        w = pinv(I + (TS.*repelem(n(:),numFeatures))' * T) * TSi * f(:);
        
        % Intersession Compensation
        w = projectionMatrix*w;
        
        % Score
        if strcmpi(scoringMethod,'CSS')
            score(fileIdx) = dot(ivectorToTest,w)/(norm(w)*norm(ivectorToTest));
        else
            score(fileIdx) = gpldaScore(gpldaModel,w,ivectorToTest);
        end
    end
    scoreFRR{speakerIdx} = score;
end
fprintf('FRR calculated (%0.0f seconds).\n',toc)
FRR calculated (17 seconds).

Ложная частота приемки (FAR)

Коэффициент ложного принятия (FAR) является частотой того, что высказывания, не принадлежащие зарегистрированному оратору, неправильно принимаются как принадлежащие зарегистрированному оратору. Создайте массив счетов для зарегистрированных динамиков и i-векторов различных динамиков.

speakersToTest = unique(adsDET.Labels);
numSpeakers = numel(speakersToTest);
scoreFAR = cell(numSpeakers,1);
tic
parfor speakerIdx = 1:numSpeakers
    adsPart = subset(adsDET,adsDET.Labels~=speakersToTest(speakerIdx));
    numFiles = numel(adsPart.Files);
    
    ivectorToTest = enrolledSpeakers.(string(speakersToTest(speakerIdx))); %#ok<PFBNS> 
    score = zeros(numFiles,1);
    for fileIdx = 1:numFiles
        audioData = read(adsPart);
        
        % Extract features
        Y = helperFeatureExtraction(audioData,afe,normFactors);
        
        % Compute a posteriori log-likelihood
        logLikelihood = helperGMMLogLikelihood(Y,ubm);
        
        % Compute a posteriori normalized probability
        amax = max(logLikelihood,[],1);
        logLikelihoodSum = amax + log(sum(exp(logLikelihood-amax),1));
        gamma = exp(logLikelihood - logLikelihoodSum)';
        
        % Compute Baum-Welch statistics
        n = sum(gamma,1);
        f = Y * gamma - n.*(ubmMu);
        
        % Extract i-vector
        w = pinv(I + (TS.*repelem(n(:),numFeatures))' * T) * TSi * f(:);
        
        % Intersession compensation
        w = projectionMatrix * w;
        
        % Score
        if strcmpi(scoringMethod,'CSS')
            score(fileIdx) = dot(ivectorToTest,w)/(norm(w)*norm(ivectorToTest));
        else
            score(fileIdx) = gpldaScore(gpldaModel,w,ivectorToTest);
        end
    end
    scoreFAR{speakerIdx} = score;
end
fprintf('FAR calculated (%0.0f seconds).\n',toc)
FAR calculated (17 seconds).

Равная вероятность ошибок (EER)

Чтобы сравнить несколько систем, вам нужна одна метрика, которая объединяет эффективность FAR и FRR. Для этого вы определяете равную частоту ошибок (EER), которая является порогом, где совпадают кривые FAR и FRR. На практике порог EER может оказаться не лучшим выбором. Для примера, если верификация динамика используется как часть подхода с множественной аутентификацией для передачи проводов, FAR, скорее всего, будет более взвешенным, чем FRR.

amin = min(cat(1,scoreFRR{:},scoreFAR{:}));
amax = max(cat(1,scoreFRR{:},scoreFAR{:}));

thresholdsToTest = linspace(amin,amax,1000);

% Compute the FRR and FAR for each of the thresholds.
if strcmpi(scoringMethod,'CSS')
    % In CSS, a larger score indicates the enroll and test ivectors are
    % similar.
    FRR = mean(cat(1,scoreFRR{:})<thresholdsToTest);
    FAR = mean(cat(1,scoreFAR{:})>thresholdsToTest);
else
    % In G-PLDA, a smaller score indicates the enroll and test ivectors are
    % similar.
    FRR = mean(cat(1,scoreFRR{:})>thresholdsToTest);
    FAR = mean(cat(1,scoreFAR{:})<thresholdsToTest);
end

[~,EERThresholdIdx] = min(abs(FAR - FRR));
EERThreshold = thresholdsToTest(EERThresholdIdx);
EER = mean([FAR(EERThresholdIdx),FRR(EERThresholdIdx)]);

figure
plot(thresholdsToTest,FAR,'k', ...
     thresholdsToTest,FRR,'b', ...
     EERThreshold,EER,'ro','MarkerFaceColor','r')
title(sprintf('Equal Error Rate = %0.4f, Threshold = %0.4f',EER,EERThreshold))
xlabel('Threshold')
ylabel('Error Rate')
legend('False Acceptance Rate (FAR)','False Rejection Rate (FRR)','Equal Error Rate (EER)','Location','best')
grid on
axis tight

Вспомогательные функции

Редукция данных и нормализация

function [features,numFrames] = helperFeatureExtraction(audioData,afe,normFactors)
    % Input:
    % audioData   - column vector of audio data
    % afe         - audioFeatureExtractor object
    % normFactors - mean and standard deviation of the features used for normalization. 
    %               If normFactors is empty, no normalization is applied.
    %
    % Output
    % features    - matrix of features extracted
    % numFrames   - number of frames (feature vectors) returned
    
    % Normalize
    audioData = audioData/max(abs(audioData(:)));
    
    % Protect against NaNs
    audioData(isnan(audioData)) = 0;
    
    % Isolate speech segment
    idx = detectSpeech(audioData,afe.SampleRate);
    features = [];
    for ii = 1:size(idx,1)
        f = extract(afe,audioData(idx(ii,1):idx(ii,2)));
        features = [features;f]; %#ok<AGROW> 
    end

    % 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 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(xMinusMu.*(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 = temp + log(gmm.ComponentProportion)';
end

Счет G-PLDA

function score = gpldaScore(gpldaModel,w1,wt)
% Center the data
w1 = w1 - gpldaModel.mu;
wt = wt - gpldaModel.mu;

% Whiten the data
w1 = gpldaModel.WhiteningMatrix*w1;
wt = gpldaModel.WhiteningMatrix*wt;

% Length-normalize the data
w1 = w1./vecnorm(w1);
wt = wt./vecnorm(wt);

% Score the similarity of the i-vectors based on the log-likelihood.
VVt = gpldaModel.EigenVoices * gpldaModel.EigenVoices';
SVVt = gpldaModel.Sigma + VVt;

term1 = pinv([SVVt VVt;VVt SVVt]);
term2 = pinv(SVVt);

w1wt = [w1;wt];
score = w1wt'*term1*w1wt - w1'*term2*w1 - wt'*term2*wt;
end

Ссылки

[1] Reynolds, Douglas A., et al. «Динамик Верификации использования адаптированных Смешанных гауссовских моделей». Цифровая обработка сигналов, том 10, № 1-3, январь 2000, стр. 19-41. DOI.org (Crossref), doi: 10.1006/dspr.1999.0361.

[2] Kenny, Patrick, et al. Joint Factor Analysis Verson Eigenchannels in Speaker Recognition (неопр.) (недоступная ссылка). Транзакции IEEE по обработке звука, речи и языка, том 15, № 4, май 2007, стр. 1435-47. DOI.org (Crossref), doi:10.1109/TASL.2006.881693.

[3] Kenny, P., et al. «Исследование изменчивости интерспикера при верификации динамика». Транзакции IEEE по обработке звука, речи и языка, том 16, № 5, июль 2008, стр. 980-88. DOI.org (Crossref), doi:10.1109/TASL.2008.925147.

[4] Dehak, Najim, et al. «Анализ фактора переднего плана для верификации динамика». Транзакции IEEE по обработке звука, речи и языка, том 19, № 4, май 2011, стр. 788-98. DOI.org (Crossref), doi:10.1109/TASL.2010.2064307.

[5] Матейка, Павел, Ондрей Глембек, Фабио Кастальдо, M.j. Алам, Олдрич Плхот, Патрик Кенни, Лукас Бургет и Ян Серноки. «Full-Covariation UBM and Heavy-Tailed PLDA in i-Vector Speaker Верификации». 2011 Международная конференция IEEE по акустике, речи и обработке сигналов (ICASSP), 2011. https://doi.org/10.1109/icassp.2011.5947436.

[6] Snyder, David, et al. «X-Vectors: Robust DNN Embeddings for Speaker Recognition». 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), IEEE, 2018, pp. 5329-33. DOI.org (Crossref), doi:10.1109/ICASSP.2018.8461375.

[7] Лаборатория обработки сигналов и речевой связи. Доступ к 12 декабря 2019 года. https://www.spsc.tugraz.at/databases-and-tools/ptdb-tug-pitch-tracking-database-from-graz-university-of-technology.html.

[8] Variani, Ehsan, et al. Глубокие нейронные сети для Text-Dependent Speaker Верификации. 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), IEEE, 2014, pp. 4052-56. DOI.org (Crossref), doi:10.1109/ICASSP.2014.6854363.

[9] Дехак, Наджим, Реда Дехак, Джеймс Р. Гласс, Дуглас А. Рейнольдс и Патрик Кенни. Оценка косинусоидальности без методов нормализации счета. Одиссея (2010).

[10] Verma, Pulkit и Pradip K. Das. «I-векторы в приложениях обработки речи: опрос». International Journal of Speech Technology, vol. 18, № 4, Dec. 2015, pp. 529-46. DOI.org (Crossref), doi: 10.1007/s10772-015-9295-3.

[11] D. Garcia-Romero and C. Espy-Wilson, «Analysis of I-vector Length Normalization in Speaker Recognition Systems». Interspeech, 2011, pp. 249-252.

[12] Кенни, Патрик. Bayesian Speaker Verification with Heavy-Tailed Priors (неопр.) (недоступная ссылка). Odyssey 2010 - Семинар по распознаванию и распознаванию языков, Брно, Чешская Республика, 2010.

[13] Sizov, Aleksandr, Kong Aik Lee, and Tomi Kinnunen. «Унификация вероятностных вариантов линейного дискриминантного анализа в биометрической аутентификации». Lecture Notes in Computer Science Structural, Syntactic, and Statistical Pattern Recognition, 2014, 464-75. https://doi.org/10.1007/978-3-662-44415-3_47.

[14] Раян, Падманабхан, Антон Афанасьев, Вилле Хаутамяки, и Томи Киннунен. 2014. «От одного до нескольких I-векторов регистрации: практические варианты оценки PLDA для верификации спикера». Цифровая обработка сигналов 31 (август): 93-101. https://doi.org/10.1016/j.dsp.2014.05.001.