%% Extract Features Using VGGish
% The VGGish network requires you to preprocess and extract features from audio signals by converting them to the sample rate the network was trained on, and then extracting log mel spectrograms. This example walks through the required preprocessing and feature extraction to match the preprocessing and feature extraction used to train VGGish. The vggishFeatures function performs these steps for you.
%
% Read in an audio signal to classify. Resample the audio signal to 16 kHz and then convert it to single precision.

[audioIn,fs0] = audioread('Ambiance-16-44p1-mono-12secs.wav');

fs = 16e3;
audioIn = resample(audioIn,fs,fs0);

audioIn = single(audioIn);

%
% Define mel spectrogram parameters and then extract features using the melSpectrogram function.

FFTLength = 512;
numBands = 64;
frequencyRange = [125 7500];
windowLength = 0.025*fs;
overlapLength = 0.015*fs;

melSpect = melSpectrogram(audioIn,fs, ...
    'Window',hann(windowLength,'periodic'), ...
    'OverlapLength',overlapLength, ...
    'FFTLength',FFTLength, ...
    'FrequencyRange',frequencyRange, ...
    'NumBands',numBands, ...
    'FilterBankNormalization','none', ...
    'WindowNormalization',false, ...
    'SpectrumType','magnitude', ...
    'FilterBankDesignDomain','warped');

%
% Convert the mel spectrogram to the log scale.

melSpect = log(melSpect + single(0.001));

%
% Reorient the mel spectrogram so that time is along the first dimension as rows.

melSpect = melSpect.';
[numSTFTWindows,numBands] = size(melSpect)

%
% Partition the spectrogram into frames of length 96 with an overlap of 48. Place the frames along the fourth dimension.

frameWindowLength = 96;
frameOverlapLength = 48;

hopLength = frameWindowLength - frameOverlapLength;
numHops = floor((numSTFTWindows - frameWindowLength)/hopLength) + 1;

frames = zeros(frameWindowLength,numBands,1,numHops,'like',melSpect);
for hop = 1:numHops
    range = 1 + hopLength*(hop-1):hopLength*(hop - 1) + frameWindowLength;
    frames(:,:,1,hop) = melSpect(range,:);
end

%
% Create a VGGish network.

net = vggish;

%
% Call predict to extract feature embeddings from the spectrogram images. The feature embeddings are returned as a numFrames-by-128 matrix, where numFrames is the number of individual spectrograms, and 128 is the number of elements in each feature vector.

features = predict(net,frames);

[numFrames,numFeatures] = size(features)

%
% Compare visualizations of the mel spectrogram and the VGGish feature embeddings.

melSpectrogram(audioIn,fs, ...
    'Window',hann(windowLength,'periodic'), ...
    'OverlapLength',overlapLength, ...
    'FFTLength',FFTLength, ...
    'FrequencyRange',frequencyRange, ...
    'NumBands',numBands, ...
    'FilterBankNormalization','none', ...
    'WindowNormalization',false, ...
    'SpectrumType','magnitude', ...
    'FilterBankDesignDomain','warped');
surf(features,'EdgeColor','none')
view([90,-90])
axis([1 numFeatures 1 numFrames])
xlabel('Feature')
ylabel('Frame')
title('VGGish Feature Embeddings')

% Copyright 2020 The MathWorks, Inc.