exponenta event banner

Проверка динамика с использованием i-векторов

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

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

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

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

Разработка системы i-Vector

В этом примере разрабатывается стандартная 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  

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

[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 МФКК, 20 дельта-МФКК и 20 дельта-дельта МФКК. Используйте длину дельта-окна 9. Извлеките элементы из окон 25 мс Hann с 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',fs, ...
    'Window',hann(windowSamples,'periodic'), ...
    'OverlapLength',overlapSamples, ...
    ...
    'mfcc',true, ...
    'mfccDelta',true, ...
    'mfccDeltaDelta',true);
setExtractorParams(afe,'mfcc','DeltaWindowLength',deltaWindowLength,'NumCoeffs',numCoeffs)

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

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

Обучение

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

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, обычно содержат около 2048 компонентов.

numComponents = 64;
if speedUpExample
    numComponents = 32;
end
alpha = ones(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
tic
for iter = 1:maxIter
    tic
    % EXPECTATION
    N = zeros(1,numComponents);
    F = zeros(numFeatures,numComponents);
    S = zeros(numFeatures,numComponents);
    L = 0;
    parfor ii = 1:numPar
        adsPart = partition(adsTrain,numPar,ii);
        while hasdata(adsPart)
            audioData = read(adsPart);
            
            % Extract features
            Y = helperFeatureExtraction(audioData,afe,normFactors);
 
            % Compute a posteriori log-liklihood
            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;
            s = (Y.*Y) * 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
    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.ComponentProportion = 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

Расчет статистики Baum-Welch

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

(ы) =∑tγt (с)

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

  • Yt - вектор признаков в момент времени t.

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

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

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

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 (ы), как описано в [3], так что

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

  • F (s) - C F × 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 (ы) lT-1 (ы)

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

Tc = Ac-1 × Start

T=[T1T2⋮TC]

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

Σ = (∑sN (s))-1 ((∑sS (s)) - диагональ (Κ×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′Σ-1 NT) T′Σ-1 F

На данном этапе все еще рассматривается каждый учебный файл в качестве отдельного докладчика. Однако на следующем шаге, когда вы обучаете матрицу проекции уменьшать размерность и увеличивать межговорящие различия, 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=1S1ns∑i=1ns (wis-ws ) (wis-ws ‾) ′

где

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

  • w‾=1N∑s=1S∑i=1nswis - это средний i-вектор для всех говорящих.

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

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

Sbv = λ Swv

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

performLDA = true;
if performLDA
    tic
    numEigenvectors = 16;

    Sw = zeros(size(projectionMatrix,1));
    Sb = zeros(size(projectionMatrix,1));
    wbar = mean(cat(2,w{:}),2);
    for ii = 1:numel(w)
        ws = w{ii};
        wsbar = mean(ws,2);
        Sb = Sb + (wsbar - wbar)*(wsbar - wbar)';
        Sw = Sw + cov(ws',1);
    end
    
    [A,~] = eigs(Sb,Sw,numEigenvectors);
    A = (A./vecnorm(A))';

    ivectorsTrain = A * ivectorsTrain;
    
    w = mat2cell(ivectorsTrain,size(ivectorsTrain,1),utterancePerSpeaker);
    
    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=1S∑s=1S1ns∑i=1ns (wis-ws ) (wis-ws ‾) ′

где

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

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

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

W-1 = BB ′

performWCCN = true;
if performWCCN
    tic
    alpha = 0.9;
    
    W = zeros(size(projectionMatrix,1));
    for ii = 1:numel(w)
        W = W + cov(w{ii}',1);
    end
    W = W/numel(w);
    
    W = (1 - alpha)*W + alpha*eye(size(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 - матрица факторной нагрузки, также известная как собственные функции.

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

numEigenVoices = 16;

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

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

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

N=∑i=1Kni

μ=1N∑i,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 = eye(size(ivectorsMatrix,1));
end

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

ivectorsMatrix = ivectorsMatrix./vecnorm(ivectorsMatrix);

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

S=∑ijφijφijT

S = ivectorsMatrix*ivectorsMatrix';

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

ivectors = mat2cell(ivectorsMatrix,D,utterancePerSpeaker);

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

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-вектора, которые были центрированы, отбелены и нормализованы по длине, оценка рассчитывается как:

оценка (w1, wt) = [w1TwtT] [Λ + VVTVVVTΛ + VVT] [w1wt] -w1T [Λ + VVT] -1w1-wtT [Λ + VVT] -1wt + C

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

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

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

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

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

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

В pratice тестовые i-векторы и, в зависимости от вашей системы, векторы регистрации не используются при обучении модели G-PLDA. В следующем разделе анализа используются ранее невидимые данные для регистрации и проверки. Поддерживающая функция gpleyScore инкапсулирует указанные выше шаги оценки и дополнительно выполняет центрирование, отбеливание и нормализацию. Сохранение обученной модели G-PLDA в качестве структуры для использования с функцией поддержки 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-векторов говорящих и 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.199.0361.

[2] Кенни, Патрик и др. «Совместный факторный анализ в сравнении с собственными каналами при распознавании говорящих». IEEE Transactions on Audio, Speech and Language Processing, vol. 15, no. 4, May 2007, pp. 1435-47. DOI.org (Crossref), doi:10.1109/TASL.2006.881693.

[3] Кенни, П., и др. «Исследование межговорящей изменчивости в проверке говорящих». Транзакции 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] Матейка, Павел, Ондрей Глембек, Фабио Кастальдо, М.Дж. Алам, Олдрич Плчот, Патрик Кенни, Лукас Бургет и Ян Черноки. 2011 Международная конференция IEEE по акустике, речи и обработке сигналов (ICASSSP), 2011. https://doi.org/10.1109/icassp.2011.5947436.

[6] Снайдер, Дэвид, и др. Международная конференция IEEE 2018 по акустике, обработке речи и сигналов (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, Эхсан, и др. Международная конференция IEEE 2014 по акустике, речи и обработке сигналов (ICASSP), IEEE, 2014, стр. 4052-56. DOI.org (Crossref), doi:10.1109/ICASSP.2014.6854363.

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

[10] Верма, Пулкит и Прадип К. Дас. «I-векторы в приложениях обработки речи: опрос». Международный журнал речевых технологий, том 18, № 4, декабрь 2015 года, стр. 529-46. DOI.org (Crossref), doi: 10.1007/s10772-015-9295-3.

[11] Д. Гарсия-Ромеро и К. Эспи-Уилсон, «Анализ нормализации длины I-вектора в системах распознавания говорящих». Interspeech, 2011, стр. 249-252.

[12] Кенни, Патрик. «Байесовская проверка спикера с тяжелохвостыми приорами». Одиссея 2010 - Семинар по распознаванию спикеров и языков, Брно, Чехия, 2010.

[13] Сизов, Александр, Конг Айк Ли и Томи Киннунен. «Унификация вероятностных вариантов линейного дискриминантного анализа при биометрической проверке подлинности». Лекции по компьютерной науке Структурное, синтаксическое и статистическое распознавание образов, 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.