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

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

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

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

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

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

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

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

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

Этот пример использует Базу данных Отслеживания Тангажа из Технологического университета Граца (PTDB-TUG) [7]. Набор данных состоит из 20 английских носителей языка, читающих 2 342 фонетически богатых предложения из корпуса 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 MFCCs, 20 дельт-MFCCs и 20 дельт дельты MFCCs. Используйте длину окна дельты 9. Извлеките функции из окон Hann на 25 мс с транзитным участком на 10 мс.

numCoeffs = 20;
deltaWindowLength = 9;

windowDuration = 0.025;
hopDuration = 0.01;

windowSamples = вокруг (windowDuration*fs);
hopSamples = вокруг (hopDuration*fs);
overlapSamples = windowSamples - hopSamples;

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

Извлеките функции из аудио, считанного из учебного datastore. Функции возвращены как numHops- 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).

Покажите коэффициенты нормализации

Используйте функцию помощника, 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 обычно, содержат приблизительно 2 048 компонентов.

numComponents = 64;
if speedUpExample
    numComponents = 32;
end
альфа = единицы (1, numComponents)/numComponents;
mu = randn (numFeatures, numComponents);
vari = rand (numFeatures, numComponents) + eps;
ubm = struct ('ComponentProportion'\alpha,'mu'\mu,'sigma', vari);

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

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

            % Compute a posteriori normalized probability
            amax = макс. (логарифмическая правдоподобность, [], 1);
            logLikelihoodSum = amax + журнал (сумма (exp (логарифмическая-правдоподобность-amax), 1));
            гамма = 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 + сумма (logLikelihoodSum);
        end
    end
    
    % Print current log-likelihood
    fprintf'Training UBM: %d/%d complete (%0.0f seconds), Log-likelihood = %0.0f\n', проход, maxIter, toc, L)
    
    % MAXIMIZATION
    N = макс. (N, eps);
    ubm.ComponentProportion = макс. (N/sum (N), eps);
    ubm.ComponentProportion = ubm.ComponentProportion/sum (ubm.ComponentProportion);
    ubm.mu = F./N;
    ubm.sigma = макс. (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

Вычислите Baum-валлийскую статистику

Baum-валлийскими статистическими данными является 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.

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

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) isa CF×CF диагональная матрица, блоки которой Nc(s)I(c=1,...C).

  • F(s) isa 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) является Baum-валлийской статистической величиной второго порядка в центре. Однако обновление Σ часто пропускается на практике, когда это оказывает мало влияния. Этот пример не обновляется Σ.

Создайте переменную сигмы.

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
    тик
    numEigenvectors = 16;

    Коротковолновый = нули (размер (projectionMatrix, 1));
    Сб = нули (размер (projectionMatrix, 1));
    wbar = среднее значение (кошка (2, w {:}), 2);
    for ii = 1:numel (w)
        ws = w {ii};
        wsbar = среднее значение (ws, 2);
        Сб = Сб + (wsbar - wbar) * (wsbar - wbar)';
        Коротковолновый = коротковолновый + cov (ws',1);
    end
    
    [A, ~] = eigs (Сб, коротковолновый, numEigenvectors);
    A = (A./vecnorm (A))';

    ivectorsTrain = * ivectorsTrain;
    
    w = mat2cell (ivectorsTrain, размер (ivectorsTrain, 1), utterancePerSpeaker);
    
    projectionMatrix = * 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
    тик
    альфа = 0.9;
    
    W = нули (размер (projectionMatrix, 1));
    for ii = 1:numel (w)
        W = W + cov (w {ii} ', 1);
    end
    W = W/numel (w);
    
    W = (1 - альфа) *W + alpha*eye (размер (W, 1));
    
    B = chol (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

Примените матрицу проекции к составу.

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

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

ϕij=μ+Vyi+εij

yiΝ(0,Ι)

εijΝ(0,Λ-1)

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

Задайте количество eigenvoices, чтобы использовать. Обычно числа между 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 (whiteningType,'ZCA')
    S = cov (ivectorsMatrix');
    [~, sD, sV] = svd (S);
    W = diag (1. / (sqrt (diag (sD)) + eps)) *sV';
    ivectorsMatrix = W * ivectorsMatrix;
elseif strcmpi (whiteningType,'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 как

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

Инициализируйте eigenvoices матрицу, 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 eigenvoice матрица. C термин является учтенными константами и может быть пропущен на практике.

speakerIdx = 2;
utteranceIdx = 1;
w1 = ivectors {speakerIdx} (: utteranceIdx);

speakerIdx = 1;
utteranceIdx = 10;
вес = ivectors {speakerIdx} (: utteranceIdx);

VVt = V*V';
SigmaPlusVVt = pinv (Lambda) + VVt;

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

w1wt = [w1; вес];
выиграйте = w1wt '*term1*w1wt - w1 '*term2*w1 - вес '*term2*wt
score = 52.4507

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

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

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

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

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

  1. Извлечение признаков

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

  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-векторам динамика и 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] Рейнольдс, Дуглас А., и др. “Верификация динамика Используя Адаптированные смешанные гауссовские модели”. Цифровая обработка сигналов, издание 10, № 1-3, январь 2000, стр 19–41. DOI.org (Crossref), doi:10.1006/dspr.1999.0361.

[2] Кенни, Патрик, и др. “Объединенный Факторный анализ По сравнению с Eigenchannels в Распознавании Динамика”. Транзакции IEEE на Аудио, Речи и Обработке Языка, издании 15, № 4, май 2007, стр 1435–47. DOI.org (Crossref), doi:10.1109/TASL.2006.881693.

[3] Кенни, P., и др. “Исследование Изменчивости Междинамика в Верификации Динамика”. Транзакции IEEE на Аудио, Речи, и Обработке Языка, издании 16, № 5, июль 2008, стр 980–88. DOI.org (Crossref), doi:10.1109/TASL.2008.925147.

[4] Dehak, Najim, и др. “Выйдите напрямую Факторный анализ для Верификации Динамика”. Транзакции IEEE на Аудио, Речи, и Обработке Языка, издании 19, № 4, май 2011, стр 788–98. DOI.org (Crossref), doi:10.1109/TASL.2010.2064307.

[5] Matejka, Павел, Ondrej Glembek, Фабио Кастальдо, M.j. Алам, Oldrich Plchot, Патрик Кенни, Лукаш Бургет и Ян Серноки. “Full-Covariance UBM и PLDA с тяжелым хвостом в Верификации Динамика i-вектора”. 2 011 Международных конференций IEEE по вопросам Акустики, Речи и Обработки сигналов (ICASSP), 2011. https://doi.org/10.1109/icassp.2011.5947436.

[6] Снайдер, Дэвид, и др. “X-векторы: Устойчивые Вложения DNN для Распознавания Динамика”. 2 018 Международных конференций IEEE по вопросам Акустики, Речи и Обработки сигналов (ICASSP), IEEE, 2018, стр 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, Эхсан, и др. “Глубокие нейронные сети для Маленькой Верификации Динамика Текстового Зависимого Места”. 2 014 Международных конференций IEEE по вопросам Акустики, Речи и Обработки сигналов (ICASSP), IEEE, 2014, стр 4052–56. DOI.org (Crossref), doi:10.1109/ICASSP.2014.6854363.

[9] Dehak, Najim, Réda Dehak, Джеймс Р. Стекло, Дуглас А. Рейнольдс и Патрик Кенни. “Подобие косинуса, выигрывающее без методов нормализации счета”. Одиссея (2010).

[10] Verma, Pulkit и Прэдип К. Дас. “I-векторы в Речевых Приложениях обработки: Обзор”. Международный журнал Речевой Технологии, издания 18, № 4, декабрь 2015, стр 529–46. DOI.org (Crossref), doi:10.1007/s10772-015-9295-3.

[11] Д. Гарсия-Ромеро и К. замечать-Вильсон, “Анализ Нормализации I-длины-вектора в Системах Распознавания Динамика”. Межречь, 2011, стр 249–252.

[12] Кенни, Патрик. "Байесова верификация динамика с уголовным прошлым с тяжелым хвостом". Одиссея 2010 - семинар распознавания динамика и языка, Брно, Чешская Республика, 2010.

[13] Сизов, Александр, Кун Айк Ли и Томи Киннунен. “Объединяя Вероятностные Варианты Линейного дискриминантного анализа в Биометрической Аутентификации”. Читайте лекции Примечаниям в Информатике Структурное, Синтаксическое, и Статистическое Распознавание образов, 2014, 464–75. https://doi.org/10.1007/978-3-662-44415-3_47.

[14] Раджан, Padmanabhan, Антон Афанасьев, Ville Hautamäki и Tomi Kinnunen. 2014. “От Одного до Нескольких I-векторов Приема: Практический PLDA Выигрыш Вариантов для Верификации Динамика”. Цифровая обработка сигналов 31 (август): 93–101. https://doi.org/10.1016/j.dsp.2014.05.001.