オーディオの特徴に関する逐次特徴選択

この例では次を使用します。

この例では、数字の音声認識タスクに適用される特徴選択の標準的なワークフローを説明します。

逐次特徴選択では、与えらえた特徴セットでネットワークに学習させた後、精度が最大になるまで特徴を段階的に追加または削除します [1]。この例では、Free Spoken Digit データセット [2] を使用して、逐次特徴選択を数字音声認識タスクに適用します。

ストリーミング数字音声の認識

この例を進めるために、まず、事前学習済みのネットワーク、ネットワークの学習に使用する audioFeatureExtractor (Audio Toolbox) オブジェクト、および特徴の正規化係数を読み込みます。

load("network_Audio_SequentialFeatureSelection.mat","bestNet","afe","normalizers");

マイクからオーディオを読み取るために audioDeviceReader (Audio Toolbox) を作成します。3 つの dsp.AsyncBuffer (DSP System Toolbox) オブジェクトを作成します。1 つは、マイクから読み取ったオーディオをバッファーするため、1 つは、音声検出用に入力オーディオの短時間エネルギーをバッファーするため、もう 1 つは予測をバッファーするために使用します。

fs = afe.SampleRate;

deviceReader = audioDeviceReader(SampleRate=fs,SamplesPerFrame=256);

audioBuffer = dsp.AsyncBuffer(fs*3);
steBuffer = dsp.AsyncBuffer(1000);
predictionBuffer = dsp.AsyncBuffer(5);

ストリーミングオーディオ、推論中のネットワーク出力の確率、および予測を表示するためのプロットを作成します。

fig = figure;

streamAxes = subplot(3,1,1);
streamPlot = plot(zeros(fs,1));
ylabel("Amplitude")
xlabel("Time (s)")
title("Audio Stream")
streamAxes.XTick = [0,fs];
streamAxes.XTickLabel = [0,1];
streamAxes.YLim = [-1,1];

analyzedAxes = subplot(3,1,2);
analyzedPlot = plot(zeros(fs/2,1));
title("Analyzed Segment")
ylabel("Amplitude")
xlabel("Time (s)")
set(gca,XTickLabel=[])
analyzedAxes.XTick = [0,fs/2];
analyzedAxes.XTickLabel = [0,0.5];
analyzedAxes.YLim = [-1,1];

probabilityAxes = subplot(3,1,3);
probabilityPlot = bar(0:9,0.1*ones(1,10));
axis([-1,10,0,1])
ylabel("Probability")
xlabel("Class")

ストリーミング数字音声 (0 ～ 9 の数字) の認識を 20 秒間実行します。ループを実行している間、数字をどれか読み上げて精度をテストします。

最初に、短時間エネルギーのしきい値を定義します。このしきい値より小さい信号には音声が含まれていないとみなします。

steThreshold = 0.015;
idxVec = 1:fs;
tic
while toc < 20
    
    % Read in a frame of audio from your device.
    audioIn = deviceReader();
    
    % Write the audio into a the buffer.
    write(audioBuffer,audioIn);
    
    % While 200 ms of data is unused, continue this loop.
    while audioBuffer.NumUnreadSamples > 0.2*fs
        
        % Read 1 second from the audio buffer. Of that 1 second, 800 ms
        % is rereading old data and 200 ms is new data.
        audioToAnalyze = read(audioBuffer,fs,0.8*fs);
        
        % Update the figure to plot the current audio data.
        streamPlot.YData = audioToAnalyze;

        ste = mean(abs(audioToAnalyze));
        write(steBuffer,ste);
        if steBuffer.NumUnreadSamples > 5
            abc = sort(peek(steBuffer));
            steThreshold = abc(round(0.4*numel(abc)));
        end
        if ste > steThreshold
            
            % Use the detectSpeeech function to determine if a region of speech
            % is present.
            idx = detectSpeech(audioToAnalyze,fs);
            
            % If a region of speech is present, perform the following.
            if ~isempty(idx)
                % Zero out all parts of the signal except the speech
                % region, and trim to 0.5 seconds.
                audioToAnalyze = trimOrPad(audioToAnalyze(idx(1,1):idx(1,2)),fs/2);
                
                % Normalize the audio.
                audioToAnalyze = audioToAnalyze/max(abs(audioToAnalyze));
                
                % Update the analyzed segment plot
                analyzedPlot.YData = audioToAnalyze;

                % Extract the features and transpose them so that time is
                % across columns.
                features = (extract(afe,audioToAnalyze))';

                % Normalize the features.
                features = (features - normalizers.Mean) ./ normalizers.StandardDeviation;
                
                % Call classify to determine the probabilities and the
                % winning label.
                features(isnan(features)) = 0;
                [label,probs] = classify(bestNet,features);
                
                % Update the plot with the probabilities and the winning
                % label.
                probabilityPlot.YData = probs;
                write(predictionBuffer,probs);

                if predictionBuffer.NumUnreadSamples == predictionBuffer.Capacity
                    lastTen = peek(predictionBuffer);
                    [~,decision] = max(mean(lastTen.*hann(size(lastTen,1)),1));
                    probabilityAxes.Title.String = num2str(decision-1);
                end
            end
        else
            % If the signal energy is below the threshold, assume no speech
            % detected.
             probabilityAxes.Title.String = "";
             probabilityPlot.YData = 0.1*ones(10,1);
             analyzedPlot.YData = zeros(fs/2,1);
             reset(predictionBuffer)
        end
        
        drawnow limitrate
    end
end

この例の残りでは、ストリーミング検出に使用したネットワークがどのように学習されたか、およびネットワークに渡される特徴がどのように選択されたかを説明します。

学習データセットと検証データセットの作成

Free Spoken Digit データセット (FSDD) をダウンロードします [2]。FSDD には、数字音声 (0 ～ 9) の短いオーディオファイルが含まれています。

downloadFolder = matlab.internal.examples.downloadSupportFile("audio","FSDD.zip");
dataFolder = tempdir;
unzip(downloadFolder,dataFolder)
dataset = fullfile(dataFolder,"FSDD");

録音を指す audioDatastore (Audio Toolbox) を作成します。データセットのサンプルレートを指定します。

ads = audioDatastore(dataset,IncludeSubfolders=true);
[~,adsInfo] = read(ads);
fs = adsInfo.SampleRate;

ファイル名の最初の要素は、ファイルに含まれる数字音声です。ファイル名の最初の要素を取得して categorical に変換し、audioDatastore の Labels プロパティを設定します。

[~,filenames] = cellfun(@(x)fileparts(x),ads.Files,UniformOutput=false);
ads.Labels = categorical(string(cellfun(@(x)x(1),filenames)));

データストアを開発セットと検証セットに分割するには、splitEachLabel (Audio Toolbox) を使用します。データの 80% を開発に割り当て、残りの 20% を検証に割り当てます。

[adsTrain,adsValidation] = splitEachLabel(ads,0.8);

オーディオ特徴抽出器のセットアップ

ウィンドウ幅 30 ms、更新頻度 10 ms でオーディオ特徴を抽出する audioFeatureExtractor (Audio Toolbox) オブジェクトを作成します。この例でテストしたい特徴をすべて true に設定します。

win = hamming(round(0.03*fs),"periodic");
overlapLength = round(0.02*fs);

afe = audioFeatureExtractor( ...
    Window=win, ...
    OverlapLength=overlapLength, ...
    SampleRate=fs, ...
    ...
    linearSpectrum=false, ...
    melSpectrum=false, ...
    barkSpectrum=false, ...
    erbSpectrum=false, ...
    ...
    mfcc=true, ...
    mfccDelta=true, ...
    mfccDeltaDelta=true, ...
    gtcc=true, ...
    gtccDelta=true, ...
    gtccDeltaDelta=true, ...
    ...
    spectralCentroid=true, ...
    spectralCrest=true, ...
    spectralDecrease=true, ...
    spectralEntropy=true, ...
    spectralFlatness=true, ...
    spectralFlux=true, ...
    spectralKurtosis=true, ...
    spectralRolloffPoint=true, ...
    spectralSkewness=true, ...
    spectralSlope=true, ...
    spectralSpread=true, ...
    ...
    pitch=false, ...
    harmonicRatio=false, ...
    zerocrossrate=false, ...
    shortTimeEnergy=false);

層と学習オプションの定義

この例で使用する深層学習層の一覧と trainingOptions を定義します。最初の層である sequenceInputLayer は、単なるプレースホルダーです。最初の層は、逐次特徴選択でどの特徴をテストするかに応じて、適切なサイズの sequenceInputLayer に置き換えられます。

numUnits = 100;
layers = [ ...
    sequenceInputLayer(1)
    bilstmLayer(numUnits,OutputMode="last")
    fullyConnectedLayer(numel(categories(adsTrain.Labels)))
    softmaxLayer
    classificationLayer];

options = trainingOptions("adam", ...
    LearnRateSchedule="piecewise", ...
    Shuffle="every-epoch", ...
    Verbose=false, ...
    MaxEpochs=20);

逐次特徴選択

逐次特徴選択の基本形式では、与えらえた特徴セットでネットワークに学習させた後、精度が向上しなくなるまで特徴を段階的に追加または削除します [1]。

前方選択

4 つの特徴のセットに前方選択を適用する単純なケースについて考えてみましょう。最初の前方選択ループでは、ネットワークを学習し検証精度を比較することによって、4 つの特徴がそれぞれ個別にテストされます。検証精度が最も高い特徴が記録されます。2 番目の前方選択ループでは、最初のループで得られた最適な特徴が残りの特徴とそれぞれ組み合わされます。次に、特徴の各ペアを使用して学習が実行されます。2 番目のループの精度が最初のループの精度より向上しなかった場合、選択プロセスは終了します。そうでない場合、新しい最適な特徴セットが選択されます。前方選択ループは、精度が向上しなくなるまで続けられます。

後方選択

後方特徴選択では、最初にすべての特徴が含まれる特徴セットで学習させた後、特徴を削除して精度が向上するかどうかをテストします。

逐次特徴選択の実行

補助関数 (sequentialFeatureSelection、trainAndValidateNetwork、および trimOrPad) は、前方逐次特徴選択または後方逐次特徴選択を実装します。学習データストア、検証データストア、オーディオ特徴抽出器、ネットワーク層、ネットワークオプション、および方向を指定します。原則として、小さい特徴セットが予想される場合は前方を選択し、大きい特徴セットが予想される場合は後方を選択します。

direction = 'forward';
[logbook,bestFeatures,bestNet,normalizers] = sequentialFeatureSelection(adsTrain,adsValidation,afe,layers,options,direction);

HelperFeatureExtractor からの出力 logbook は、テストされたすべての特徴の構成および対応する検証精度が格納されたテーブルです。

logbook

logbook=62×2 table
                                  Features                                   Accuracy
    _____________________________________________________________________    ________

    "mfccDelta, spectralKurtosis, spectralRolloffPoint"                       98.25  
    "mfccDelta, spectralRolloffPoint"                                         97.75  
    "mfccDelta, spectralEntropy, spectralRolloffPoint"                        97.75  
    "mfccDelta, spectralDecrease, spectralKurtosis, spectralRolloffPoint"     97.25  
    "mfccDelta, mfccDeltaDelta"                                                  97  
    "mfccDelta, gtccDeltaDelta, spectralRolloffPoint"                            97  
    "mfcc, mfccDelta, spectralKurtosis, spectralRolloffPoint"                    97  
    "mfcc, mfccDelta"                                                         96.75  
    "mfccDelta, gtccDeltaDelta, spectralKurtosis, spectralRolloffPoint"       96.75  
    "mfccDelta, spectralRolloffPoint, spectralSlope"                           96.5  
    "mfccDelta"                                                               96.25  
    "mfccDelta, spectralKurtosis"                                             96.25  
    "mfccDelta, spectralSpread"                                               96.25  
    "mfccDelta, spectralDecrease, spectralRolloffPoint"                       96.25  
    "mfccDelta, spectralFlatness, spectralKurtosis, spectralRolloffPoint"     96.25  
    "mfccDelta, gtccDeltaDelta"                                                  96  
      ⋮

sequentialFeatureSelection からの出力 bestFeatures には、最適な特徴が true に設定された struct が含まれています。

bestFeatures

bestFeatures = struct with fields:
                    mfcc: 0
               mfccDelta: 1
          mfccDeltaDelta: 0
                    gtcc: 0
               gtccDelta: 0
          gtccDeltaDelta: 0
        spectralCentroid: 0
           spectralCrest: 0
        spectralDecrease: 0
         spectralEntropy: 0
        spectralFlatness: 0
            spectralFlux: 0
        spectralKurtosis: 1
    spectralRolloffPoint: 1
        spectralSkewness: 0
           spectralSlope: 0
          spectralSpread: 0

この struct を使用して audioFeatureExtractor を設定できます。

set(afe,bestFeatures)
afe

afe = 
  audioFeatureExtractor with properties:

   Properties
                     Window: [240×1 double]
              OverlapLength: 160
                 SampleRate: 8000
                  FFTLength: []
    SpectralDescriptorInput: 'linearSpectrum'
        FeatureVectorLength: 15

   Enabled Features
     mfccDelta, spectralKurtosis, spectralRolloffPoint

   Disabled Features
     linearSpectrum, melSpectrum, barkSpectrum, erbSpectrum, mfcc, mfccDeltaDelta
     gtcc, gtccDelta, gtccDeltaDelta, spectralCentroid, spectralCrest, spectralDecrease
     spectralEntropy, spectralFlatness, spectralFlux, spectralSkewness, spectralSlope, spectralSpread
     pitch, harmonicRatio, zerocrossrate, shortTimeEnergy


   To extract a feature, set the corresponding property to true.
   For example, obj.mfcc = true, adds mfcc to the list of enabled features.

sequentialFeatureSelection は、最も性能が良いネットワークと、選択された特徴に対応する正規化係数も出力します。ネットワーク、構成済みの audioFeatureExtractor、および正規化係数を保存するには、次の行のコメントを解除します。

% save('network_Audio_SequentialFeatureSelection.mat','bestNet','afe','normalizers')

まとめ

この例では、再帰型ニューラルネットワーク (LSTM または BiLSTM) の逐次特徴選択のワークフローについて説明しました。これは、CNN および RNN-CNN のワークフローに簡単に適用することができます。

サポート関数

学習ネットワークと検証ネットワーク

function [trueLabels,predictedLabels,net,normalizers] = trainAndValidateNetwork(adsTrain,adsValidation,afe,layers,options)
% Train and validate a network.
%
%   INPUTS:
%   adsTrain      - audioDatastore object that points to training set
%   adsValidation - audioDatastore object that points to validation set
%   afe           - audioFeatureExtractor object.
%   layers        - Layers of LSTM or BiLSTM network
%   options       - trainingOptions object
%
%   OUTPUTS:
%   trueLabels      - true labels of validation set
%   predictedLabels - predicted labels of validation set
%   net             - trained network
%   normalizers     - normalization factors for features under test

% Copyright 2019 The MathWorks, Inc.

% Convert the data to tall arrays.
tallTrain = tall(adsTrain);
tallValidation = tall(adsValidation);

% Extract features from the training set. Reorient the features so that
% time is along rows to be compatible with sequenceInputLayer.
fs = afe.SampleRate;
tallTrain = cellfun(@(x)trimOrPad(x,fs/2),tallTrain,UniformOutput=false);
tallTrain = cellfun(@(x)x/max(abs(x),[],"all"),tallTrain,UniformOutput=false);
tallFeaturesTrain = cellfun(@(x)extract(afe,x),tallTrain,UniformOutput=false);
tallFeaturesTrain = cellfun(@(x)x',tallFeaturesTrain,UniformOutput=false);  %#ok<NASGU>
[~,featuresTrain] = evalc('gather(tallFeaturesTrain)'); % Use evalc to suppress command-line output.

tallValidation = cellfun(@(x)trimOrPad(x,fs/2),tallValidation,UniformOutput=false);
tallValidation = cellfun(@(x)x/max(abs(x),[],'all'),tallValidation,UniformOutput=false);
tallFeaturesValidation = cellfun(@(x)extract(afe,x),tallValidation,UniformOutput=false);
tallFeaturesValidation = cellfun(@(x)x',tallFeaturesValidation,UniformOutput=false); %#ok<NASGU>
[~,featuresValidation] = evalc('gather(tallFeaturesValidation)'); % Use evalc to suppress command-line output.

% Use the training set to determine the mean and standard deviation of each
% feature. Normalize the training and validation sets.
allFeatures = cat(2,featuresTrain{:});
M = mean(allFeatures,2,"omitnan");
S = std(allFeatures,0,2,"omitnan");
featuresTrain = cellfun(@(x)(x-M)./S,featuresTrain,UniformOutput=false);
for ii = 1:numel(featuresTrain)
    idx = find(isnan(featuresTrain{ii}));
    if ~isempty(idx)
        featuresTrain{ii}(idx) = 0;
    end
end
featuresValidation = cellfun(@(x)(x-M)./S,featuresValidation,UniformOutput=false);
for ii = 1:numel(featuresValidation)
    idx = find(isnan(featuresValidation{ii}));
    if ~isempty(idx)
        featuresValidation{ii}(idx) = 0;
    end
end

% Replicate the labels of the train and validation sets so that they are in
% one-to-one correspondence with the sequences.
labelsTrain = adsTrain.Labels;

% Update input layer for the number of features under test.
layers(1) = sequenceInputLayer(size(featuresTrain{1},1));

% Train the network.
net = trainNetwork(featuresTrain,labelsTrain,layers,options);

% Evaluate the network. Call classify to get the predicted labels for each
% sequence.
predictedLabels = classify(net,featuresValidation);
trueLabels = adsValidation.Labels;

% Save the normalization factors as a struct.
normalizers.Mean = M;
normalizers.StandardDeviation = S;
end

逐次特徴選択

function [logbook,bestFeatures,bestNet,bestNormalizers] = sequentialFeatureSelection(adsTrain,adsValidate,afeThis,layers,options,direction)
%
%   INPUTS:
%   adsTrain    - audioDatastore object that points to training set
%   adsValidate - audioDatastore object that points to validation set
%   afe         - audioFeatureExtractor object. Set all features to test to true
%   layers      - Layers of LSTM or BiLSTM network
%   options     - trainingOptions object
%   direction   - SFS direction, specify as 'forward' or 'backward'
%
%   OUTPUTS:
%   logbook         - table containing feature configurations tested and corresponding validation accuracies
%   bestFeatures    - struct containg best feature configuration
%   bestNet         - Trained network with highest validation accuracy
%   bestNormalizers - Feature normalization factors for best features

% Copyright 2019 The MathWorks, Inc.

afe = copy(afeThis);
featuresToTest = fieldnames(info(afe));
N = numel(featuresToTest);
bestValidationAccuracy = 0;

% Set the initial feature configuration: all on for backward selection
% or all off for forward selection.
featureConfig = info(afe);
for i = 1:N
    if strcmpi(direction,"backward")
        featureConfig.(featuresToTest{i}) = true;
    else
        featureConfig.(featuresToTest{i}) = false;
    end
end

% Initialize logbook to track feature configuration and accuracy.
logbook = table(featureConfig,0,VariableNames=["Feature Configuration","Accuracy"]);

% Perform sequential feature evaluation.
wrapperIdx = 1;
bestAccuracy = 0;
while wrapperIdx <= N
    % Create a cell array containing all feature configurations to test
    % in the current loop.
    featureConfigsToTest = cell(numel(featuresToTest),1);
    for ii = 1:numel(featuresToTest)
        if strcmpi(direction,"backward")
            featureConfig.(featuresToTest{ii}) = false;
        else
            featureConfig.(featuresToTest{ii}) = true;
        end
        featureConfigsToTest{ii} = featureConfig;
        if strcmpi(direction,"backward")
            featureConfig.(featuresToTest{ii}) = true;
        else
            featureConfig.(featuresToTest{ii}) = false;
        end
    end

    % Loop over every feature set.
    for ii = 1:numel(featureConfigsToTest)

        % Determine the current feature configuration to test. Update
        % the feature afe.
        currentConfig = featureConfigsToTest{ii};
        set(afe,currentConfig)

        % Train and get k-fold cross-validation accuracy for current
        % feature configuration.
        [trueLabels,predictedLabels,net,normalizers] = trainAndValidateNetwork(adsTrain,adsValidate,afe,layers,options);
        valAccuracy = mean(trueLabels==predictedLabels)*100;
        if valAccuracy > bestValidationAccuracy
            bestValidationAccuracy = valAccuracy;
            bestNet = net;
            bestNormalizers = normalizers;
        end

        % Update Logbook
        result = table(currentConfig,valAccuracy,VariableNames=["Feature Configuration","Accuracy"]);
        logbook = [logbook;result]; %#ok<AGROW> 

    end

    % Determine and print the setting with the best accuracy. If accuracy
    % did not improve, end the run.
    [a,b] = max(logbook{:,"Accuracy"});
    if a <= bestAccuracy
        wrapperIdx = inf;
    else
        wrapperIdx = wrapperIdx + 1;
    end
    bestAccuracy = a;

    % Update the features-to-test based on the most recent winner.
    winner = logbook{b,"Feature Configuration"};
    fn = fieldnames(winner);
    tf = structfun(@(x)(x),winner);
    if strcmpi(direction,"backward")
        featuresToRemove = fn(~tf);
    else
        featuresToRemove = fn(tf);
    end
    for ii = 1:numel(featuresToRemove)
        loc =  strcmp(featuresToTest,featuresToRemove{ii});
        featuresToTest(loc) = [];
        if strcmpi(direction,"backward")
            featureConfig.(featuresToRemove{ii}) = false;
        else
            featureConfig.(featuresToRemove{ii}) = true;
        end
    end

end

% Sort the logbook and make it more readable.
logbook(1,:) = []; % Delete placeholder first row.
logbook = sortrows(logbook,"Accuracy","descend");
bestFeatures = logbook{1,"Feature Configuration"};
m = logbook{:,"Feature Configuration"};
fn = fieldnames(m);
myString = strings(numel(m),1);
for wrapperIdx = 1:numel(m)
    tf = structfun(@(x)(x),logbook{wrapperIdx,"Feature Configuration"});
    myString(wrapperIdx) = strjoin(fn(tf),", ");
end
logbook = table(myString,logbook{:,"Accuracy"},VariableNames=["Features","Accuracy"]);
end

トリミングまたはパディング

function y = trimOrPad(x,n)
% y = trimOrPad(x,n) trims or pads the input x to n samples. If x is
% trimmed, it is trimmed equally on the front and back. If x is padded, it is
% padded equally on the front and back with zeros. For odd-length trimming or
% padding, the extra sample is trimmed or padded from the back.

% Copyright 2019 The MathWorks, Inc.
a = size(x,1);
if a < n
    frontPad = floor((n-a)/2);
    backPad = n - a - frontPad;
    y = [zeros(frontPad,1);x;zeros(backPad,1)];
elseif a > n
    frontTrim = floor((a-n)/2)+1;
    backTrim = a - n - frontTrim;
    y = x(frontTrim:end-backTrim);
else
    y = x;
end
end

参考文献

[1] Jain, A., and D. Zongker. "Feature Selection: Evaluation, Application, and Small Sample Performance." IEEE Transactions on Pattern Analysis and Machine Intelligence. Vol. 19, Issue 2, 1997, pp. 153-158.

[2] Jakobovski. “Jakobovski/Free-Spoken-Digit-Dataset.” GitHub, May 30, 2019. https://github.com/Jakobovski/free-spoken-digit-dataset.