% SpatialNoiseVideoMap_ambisonics_image.m Spatial acoustic mapping of a recording
% This script implements source acoustic mapping using recordings of noise signal 
% got from a microphone array encoded into Ambisonics.
% Background is a static image. It saves video overlay with pseudocolor SPL maps.
%
% INPUT
% - encoded_format.mat -> file containing the ambisonics format generated by b_generateEncodedFormat
% - panoramic_360_image.jpg
%
% OUTPUT
% - acusticVideo.mp4
%
% Author: Daniel Pinardi.
% University of Aalto, Helsinki - Finlan, 2019.
% Last review: 01/2020, Daniel Pinardi

clear all
close all
clc

addpath( './lib' )                                                          % add to path the library with needed common files
a_configureFileNamesAndDirectories                                            % configure filenames, directories, resources
opengl( 'save', 'hardware' )												% set parameter for video processing
set( 0, 'DefaultTextInterpreter', 'none' );                                 % setting for map titles
warning off                                                                 % to avoid warnings on image size


%% AUDIO PARAMETERS
% set fullscale value
% EIGENMIKE -> GAIN = x  -  FS_SPS = 115    ; FS_ambix = 115
% CYLINDER  -> GAIN = x  -  FS_SPS =        ; FS_ambix = 
% RICOH     -> GAIN = 0     FS_SPS =        ; FS_ambix = 
% AMBEO     -> GAIN = 0     FS_SPS =        ; FS_ambix = 
% HEAD      -> GAIN = 0     FS_SPS =        ; FS_ambix = 
EM32gain                = 20;                                               % integer from 0 to 30
fullScaleLevel          = 115;                                              % dB

% Signal processing
trimEnable              = 1;        % 1 to enable signal trimming
aWeightingEnable        = 1;        % 1 to enable A filtering

% Bandpass filtering -> 0 (one of the two) to disable
bandpassLowCutFreq      = 0;        % Low Cutoff Frequency [Hz]
bandpassHighCutFreq     = 6000;     % High Cutoff Frequency [Hz]

% Audio processing
audioMono               = 0;        % 1 to enable: add omni audio to .mp4
ambix_o1                = 0;        % 1 to enable: add ambix 1st order to map video
ambix_o3                = 1;        % 1 to enable: add ambix 3rd order to map video
outBitsPerSample        = 32;       % 16, 24 or 32

timeSignalEnable        = 1;        % 1 to enable: plot time signal under video maps. Disable for VR
bufferSize              = 16384;    % [samples]
autoselectOverlap       = 1;        % enable for VR video
overlap                 = 0.975;    % user defined overlap                            
frameRate               = 30;       % Default 30 - min value for VLC is 6. 30 for VR

%% AMBISONICS PARAMETER
fileAlreadyBformat      = 0; % 1 -> enable: use if the file is not A-format recording but already encoded outside
                             % file must be placed in ArrayRecording and filename must be set in
                             % configureFileNamesAndDirectories -> recordedFileName
resolutionMap           = 5; % degree
maxOrderToUse           = 3;

% Frequency range order limitation
o2_lowF                 = 500;      % [Hz]
o3_lowF                 = 1000;     % [Hz]
o4_lowF                 = 2000;     % [Hz]
o4_highF                = 8000; 	% [Hz]
o3_highF                = 9000;     % [Hz]
o2_highF                = 10000;    % [Hz]

% Map type enable: ONLY ONE ENABLED
pwd                     = 1; % plane wave decomposition
mvdr                    = 0;
intVect                 = 0; % only o1
music                   = 0;
musicSourceNum          = 2; % suggested 1 for o1 and 2 for o3

%% IMAGE/VIDEO PARAMETER
% pointing correction: + -> turn map to the right / - -> turn map to the left
horizontalTilt          = 0;        % degrees

% image
imageOffsetCorrection   = 0;    	% 1 to enable. Correct vertical offset between acoustical and optical center
radiusCorrection        = 1.5;      % [m]. The radius of the sphere where correction is performed
verticalOffset          = 0.15;     % [m]. The amount of offset to be compensated 
mapWindowWidth          = 1200;     % 0 to autoselect maximum possible, depending on screen resoution. Aspect ratio will be 2:1
outputWidth             = 0;        % 0 -> keep the same of colormap, xxxx -> rescale image and colormap

outputBlackWhite        = 1;        % 1 -> blackwhite, 0 -> color
alpha                   = 0.3;      % parameter for transparency

% Surf mapping
useAdimensialPlotBounds = 0;        % colormap is in the range 0 - 1
adimensionalLowerLimit  = 0.05;     % upperlimit is always 1
useFixedSplPlotBounds   = 0;        % if 0 each surface plot has color normalization with fixedLowerLimit and fixedUpperLimit
fixedLowerLimit         = 80;       % user defined lower limit
fixedUpperLimit         = 90;       % user defined upper limit
noiseGate               = 5;        % max difference from maxlevel to plot
contourEnable           = 1;        % 1 to enable isolevel lines
contourSpan             = 1;        % [dB] define isoline difference
valueTickEnable         = 0;        % 1 to enable: show SPL value of isolevel lines
tickSize                = 14;
tickBold                = 1;
colorbarEnable          = 0;
injectMetadata          = 0;        % 1 -> enable, 0 -> disable


%% CHECK IF FFMPEG IS PRESENT IN C:\WINDOWS
if exist( 'C:\Windows\ffmpeg.exe', 'file' ) == 0
    error( 'FFMPEG not found. Please download from " https://ffmpeg.zeranoe.com/builds/ " and copy ffmpeg.exe in C:\Windows' );
end


%% LOAD DATA FROM FILES
% clear unused variables
clear backgroundVideoFile backgroundVideoFileName invSweepFile invSweepFileName irsFolder irsMatFileName irsWavFileName
clear mapsFolderIrs mapsFolderNoise outputFolder recordedNoiseFolder recordedReferenceSignalFileName recordedReferenceSignalFolder
clear recordedSweepFile recordedSweepFolder spsFilterMatrixFile spsFilterMatrixFileName spsIrsFolder spsIrsMatFileName
clear spsIrsWavFileName spsNoiseFolder spsNoiseWavFileName

% define filename to load and folder output
filename                = encodedNoiseMatFileName;
fprintf( '\nNoise elaboration started\n' );

% Ambisonics format
if fileAlreadyBformat
    [ ambixFormat, fs ] = audioread( [ './ArrayRecording/' recordedFileName ] );
else
    disp( [ 'Loading Ambix format file: ' encodedNoiseMatFileName] );
    load( encodedNoiseMatFileName );
    ambixFormat         = encodedFormat;
    clear encodedFormat
end

ambixVmic               = size( ambixFormat, 2 );                           % number of virtual microphones
ambixOrder              = sqrt( ambixVmic ) - 1;                            % ambisonics order available
if maxOrderToUse < ambixOrder
    ambixVmic           = ( maxOrderToUse + 1 ) ^ 2;
    ambixFormat         = ambixFormat( :, 1 : ambixVmic );
    ambixOrder          = maxOrderToUse;
end

if ambix_o3 == 1 && ambixOrder < 3
    fprintf( 'Output soundtrack downgraded to ambix-o1.\n' )
    ambix_o3            = 0;
    ambix_o1            = 1;
end

% Background 360° image
disp( [ 'Loading background image file: ' backgroundImageFile ] );
backgroundImage         = imread( backgroundImageFile );


%% IMAGE PROCESSING
% check monitor dimension and resize if necessary
[ mapWindowWidth ]      = screenSize( mapWindowWidth, timeSignalEnable );                                                         % clear unused variables

% check image aspect ratio and perform correction if needed ( output format is 2x1 ) and resize image
[ backgroundImage, backgroundHeight, backgroundWidth, ~ ] = image2x1resize( backgroundImage, mapWindowWidth );

% correction of offset between acoustic center and image center
if imageOffsetCorrection
    fprintf( 'Perform image correction\n\n' );
    [ backgroundImage ] = imageVertOffsetCorr( backgroundImage, radiusCorrection, verticalOffset );
end

if outputBlackWhite == 1
    background          = (1+alpha)*flipud( backgroundImage( :, :, [ 1 1 1 ] ) ); % flip and conversion to black white
else
    background          = (1+alpha)*flipud( backgroundImage );              % flip   
end
clear radiusCorrection verticalOffset


%% A-weight filtering
ambixFormatAudio        = ambixFormat;
if aWeightingEnable
    fprintf( 'Generating A-weight filter.\n' );
    AWeighting          = weightingFilter( 'A-weighting', fs );             % use Matlab built-in A_weigth filter
    fprintf( 'A-weight filtering started.\n' );
    for i = 1 : ambixVmic
        ambixFormat( :, i ) = AWeighting( ambixFormat( :, i ) );
    end
    fprintf( 'A-weighting ended.\n' );
end


%% Bandpass filtering
if bandpassLowCutFreq ~= 0 && bandpassHighCutFreq ~= 0
    fprintf( 'Generating band-pass filter.\n' );
    hBandPass           = fdesign.bandpass( 'N,F3dB1,F3dB2', 10, bandpassLowCutFreq, bandpassHighCutFreq, fs );
    HBandPass           = design( hBandPass, 'butter' );
    fprintf( 'Band-pass filtering started.\n' );
    ambixFormat         = filter( HBandPass, ambixFormat );
    ambixFormatAudio    = filter( HBandPass, ambixFormatAudio );
    fprintf( 'Band-pass filtering ended.\n' );
end


%% Evaluate stepSize and overlap
if autoselectOverlap
    stepSize            = floor( fs / frameRate );
    overlap             = ( 1 - stepSize / bufferSize ) * 100;    
else
    stepSize            = floor( bufferSize * ( 1 - overlap ) );                % calculate stepSize
end
fprintf( [ '\nFrame rate is ', num2str( frameRate ) ' fps, fs is ' num2str( fs ) ' Hz and stepsize is ' num2str( stepSize ) ' samples. \n' ] );
fprintf( [ '\nBuffer size is ' num2str( bufferSize ) ' samples and video overlap is ', num2str( overlap ) '. \n' ] );


%% TRIM
if trimEnable
    
    fig                 = figure;
    hax                 = axes;
    plot( ambixFormat( :, 1 ) )
    xlim( [ 0 size( ambixFormat, 1 ) ] )
    grid on
    xlabel( 'Samples' )
    title( 'Audio recorded' )

    prompt              = 'AUDIO samples to skip ( 0 is accepted ): ';
    noiseStart          = input( prompt );
    line( [ noiseStart noiseStart ], get( hax, 'YLim' ), 'Color', [ 1 0 0 ] )
    prompt              = 'last AUDIO sample: ';
    noiseEnd            = input( prompt );

    close all
    ambixFormat         = ambixFormat( 1 + noiseStart : noiseEnd, : );      % trim noise
    ambixFormatAudio    = ambixFormatAudio( 1 + noiseStart : noiseEnd, : ); % trim noise
end

ambixLength             = size( ambixFormat, 1 );
numbersFrame            = floor( ambixLength / stepSize );                  % integer number of frames
noiseEnd                = stepSize * numbersFrame;                          % update value
ambixFormat             = ambixFormat( 1 : noiseEnd, : );                   % trim noise
ambixFormatAudio        = ambixFormatAudio( 1 : noiseEnd, : );              % trim noise


%% FULL SCALE LEVEL
fullScaleLevel          = 0.0004*( EM32gain^3 ) - 0.0078*( EM32gain^2 ) - 0.8821*( EM32gain ) + fullScaleLevel; % dB
clear EM32gain
fullScale               = fullScaleLevel + 20 * log10( rms( ambixFormat( :, 1 ) ) );


%% ZERO PADDING
zeroPadMatrix           = zeros( bufferSize/2, ambixVmic );                 % center of buffer must be at 0
ambixFormat             = [ zeroPadMatrix; ambixFormat; zeroPadMatrix ];    % apply zero padding
ambixFormatAudio        = [ zeroPadMatrix; ambixFormatAudio; zeroPadMatrix ];% apply zero padding


%% Adjust ch num
if intVect == 1
    fprintf( 'Order set to 1 for intensity vector analysis.' )
    useAdimensialPlotBounds = 1;   
    valueTickEnable     = 0;
    useFixedSplPlotBounds = 0;
    M_sh2pv             = beamWeightsPressureVelocity( 'real' );
    freqDomain          = 0;
    ambixFormat         = ambixFormat( :, 1 : 4 );
    ambixFormatAudio    = ambixFormatAudio( :, 1 : 4 );
    if ambix_o3 == 1
        ambix_o3        = 0;
        ambix_o1        = 1;
    end
else
    freqDomain          = 1;
    M_sh2pv             = 0;
end
ambixVmic               = size( ambixFormat, 2 );


%% GRID CREATION
aziElev2aziPolar        = @(dirs) [ dirs( :, 1 ) pi/2-dirs( :, 2 ) ];       
grid_dirs               = grid2dirs( resolutionMap, resolutionMap, 0, 0 ); % Grid of directions to evaluate DoA estimation
grid_dirs( :, 1 )       = - grid_dirs( :, 1 );
gridMap                 = Fdirs2grid( grid_dirs( :, 1 ), resolutionMap, resolutionMap, 1 );
gridMapEl               = size( gridMap, 1 );
gridMapAz               = size( gridMap, 2 );


%% DATA PROCESSING
ambixOmni               = ambixFormat( :, 1 );
nfft                    = bufferSize;

fprintf( [ 'Processing ' num2str( numbersFrame ) ' frames.\n' ] );
splMatrix               = zeros ( gridMapEl, gridMapAz, numbersFrame );     % preallocation

% Define window
hannWin                 = hann( bufferSize );
energyWinFactor         = ( rms( hannWin ) ) ^-1;
hannWin                 = repmat( hannWin, [ 1 ambixVmic ] );

% set frequency limit for orders > 1
freq                    = linspace( 0, fs/2, nfft/2 + 1 );
[ ~, index_o2_lowF ]    = min( abs( freq - o2_lowF ) );
[ ~, index_o3_lowF ]    = min( abs( freq - o3_lowF ) );
[ ~, index_o4_lowF ]    = min( abs( freq - o4_lowF ) );
[ ~, index_o2_highF ]   = min( abs( freq - o2_highF ) );
[ ~, index_o3_highF ]   = min( abs( freq - o3_highF ) );
[ ~, index_o4_highF ]   = min( abs( freq - o4_highF ) );

% Percentage initializing
oldCompletedPercentage  = -1;
fprintf( '\nFrames processing: ' )
varStr                  = '0%%';
varStrLen               = length( varStr );
fprintf( '%s', varStr )

for blockIndex = 1 : numbersFrame
    
    %%%%%%%%%%% PERCENTAGE PROCESSING %%%%%%%%%%%  
    completedPercentage = round( blockIndex/numbersFrame*100 );
    if oldCompletedPercentage ~= completedPercentage
        for i = 1 : varStrLen
            fprintf( '\b' )
        end
        varStr          = sprintf( '%d%%', completedPercentage );
        varStrLen       = length( varStr );
        fprintf( '%s', varStr )
    end
    oldCompletedPercentage = completedPercentage;
    
    if freqDomain == 1
        [ DoA ]         = frameProcessing_ambisonics_fd( ambixFormat, blockIndex, stepSize, bufferSize, pwd, mvdr, ...
                                                         music, grid_dirs, musicSourceNum, hannWin, energyWinFactor, ...
                                                         index_o2_lowF, index_o3_lowF, index_o4_lowF, index_o2_highF, ...
                                                         index_o3_highF, index_o4_highF, ambixOrder );
    else % intensity vector analysis
        
        [ DoA ]         = frameProcessing_ambisonics_td( ambixFormat, blockIndex, stepSize, bufferSize, grid_dirs, ...
                                                         M_sh2pv);
    end
    
    splMatrix( :, :, blockIndex ) = Fdirs2grid( DoA, resolutionMap, resolutionMap, 1 );

end

splMatrix               = single( splMatrix );

% set colormap limits
maxDoA                  = max( max( max( abs( splMatrix ) ) ) );

if useAdimensialPlotBounds == 1
    
    splMatrix           = splMatrix ./ maxDoA;
    upperLimit          = 1;
    lowerLimit          = adimensionalLowerLimit;
    v                   = lowerLimit : 0.1 : upperLimit;                    % iso-line levels
    
else
    
    if freqDomain
        ambixFrf        = avgFft( ambixFormat, nfft, 0.5 );
        complexSum      = real( ambixFrf ).^2 + imag( ambixFrf ).^2;
        ambixFrfReal    = mean( complexSum, 3 );
        rmsOmniUnfiltered = ( sum( ambixFrfReal( :, 1 ), 1 )' ) ./ 2;
        splMatrix       = splMatrix ./ maxDoA * rmsOmniUnfiltered;
        splMatrix       = 10 * log10( splMatrix ) + fullScaleLevel;
    else
        splMatrix       = splMatrix ./ maxDoA * ( rms( ambixOmni ) );
        splMatrix       = 20 * log10( splMatrix ) + fullScaleLevel;
    end
    
    if useFixedSplPlotBounds == 1
        upperLimit      = fixedUpperLimit;
        lowerLimit      = fixedLowerLimit;
    else
        upperLimit      = ceil( max( max( max( splMatrix ) ) ) );           % find absolute max level
        lowerLimit      = upperLimit - noiseGate;                           % range to plot
    end
    
    v                   = lowerLimit+1 : contourSpan : upperLimit;          % iso-line levels
    splMatrix( splMatrix < lowerLimit ) = 0;                                % set to 0 SPL values < lowerLimit
        
end

% pointing correction
correctionIndex     = round( horizontalTilt * gridMapAz / 360 );
if correctionIndex > 0  
    splMatrix       = splMatrix( :, [ correctionIndex+1 : gridMapAz, 1:correctionIndex ], : );
elseif correctionIndex < 0
    splMatrix       = splMatrix( :, [ gridMapAz+correctionIndex+1 : gridMapAz, 1:gridMapAz+correctionIndex ], : );
end

fprintf( '\nData processing ended correctly.\n' )


%% VIDEO PROCESSING
fprintf( '\nInitializing video processing.\n' );
videoTitle          = recordedFileName( 1 : end-4 );                        % name of file video
videoFile           = VideoWriter( videoTitle, 'MPEG-4' );                  % initializes video
videoFile.FrameRate = frameRate;                                            % set frame rate
videoFile.Quality   = 100;                                                  % set quality
open( videoFile );
cmap                = jet( 64 );
cmap( 1, : )        = [ 0 0 0 ]; % black full scale
duration            = stepSize * numbersFrame / fs;

% Percentage initializing
oldCompletedPercentage  = -1;
fprintf( '\nVideo processing: ' )
varStr              = '0%%';
varStrLen           = length( varStr );
fprintf( '%s', varStr )

if timeSignalEnable

    verticalLimit   = [ min( ambixOmni ) max( ambixOmni ) ];
    videoHeigth     = mapWindowWidth / ( 2 * 0.8 );                         % define video window
    videoWindow     = figure( 'units', 'pixel', 'innerposition', [ 0 0 mapWindowWidth videoHeigth  ], 'menubar', 'none' );
    imageHeight     = mapWindowWidth / 2;
    
    % reshape diretions for colormap generation
    elev            = -imageHeight/2 : resolutionMap/90*imageHeight/2 : imageHeight/2;
    azi             = -mapWindowWidth/2 : resolutionMap/180*mapWindowWidth/2 : mapWindowWidth/2;
    aziMat          = repmat( azi, size( elev, 2 ), 1 );
    elevMat         = repmat( elev, size( azi, 2 ), 1 )';

    for blockIndex = 1 : numbersFrame
        
        %%%%%%%%%%% PERCENTAGE PROCESSING %%%%%%%%%%%  
        completedPercentage     = round( blockIndex/numbersFrame*100 );
        if oldCompletedPercentage ~= completedPercentage
            for i = 1 : varStrLen
                fprintf( '\b' )
            end
            varStr              = sprintf( '%d%%', completedPercentage );
            varStrLen           = length( varStr );
            fprintf( '%s', varStr )
        end
        oldCompletedPercentage  = completedPercentage;
        
        ZI          = splMatrix( :, :, blockIndex );                        % extract current value to plot

        % subplot colormap
        plotBG      = subplot( 'Position', [ 0 0.2 1 0.8 ] );
        ambisonicsVideoColorMap( background, backgroundWidth, backgroundHeight, azi, elev, ZI, lowerLimit, upperLimit, ...
                                 plotBG, v, alpha, contourEnable, valueTickEnable, cmap, aziMat, elevMat, tickSize, tickBold, colorbarEnable )
        
        % subplot time signal
        plotTime   	= subplot( 'Position', [ 0 0 1 0.2 ] );                 % subplot with time signal
        timeSignalPlot( ambixOmni, verticalLimit, blockIndex, stepSize, bufferSize );

        writeVideo( videoFile, getframe( videoWindow ) );
        
    end
    
elseif outputWidth ~= 0
    
    videoHeigth     = mapWindowWidth / 2;                                % set video window dimension to 2:1 for VR
    videoWindow   	= figure( 'units', 'pixel', 'innerposition', [ 0 0 mapWindowWidth videoHeigth  ], 'menubar', 'none' );

    % reshape diretions for colormap generation
    elev            = -videoHeigth/2 : resolutionMap/90*videoHeigth/2 : videoHeigth/2;
    azi             = -mapWindowWidth/2 : resolutionMap/180*mapWindowWidth/2 : mapWindowWidth/2;
    aziMat          = repmat( azi, size( elev, 2 ), 1 );
    elevMat         = repmat( elev, size( azi, 2 ), 1 )';
    
    for blockIndex = 1 : numbersFrame

        %%%%%%%%%%% PERCENTAGE PROCESSING %%%%%%%%%%%  
        completedPercentage     = round( blockIndex/numbersFrame*100 );
        if oldCompletedPercentage ~= completedPercentage
            for i = 1 : varStrLen
                fprintf( '\b' )
            end
            varStr              = sprintf( '%d%%', completedPercentage );
            varStrLen           = length( varStr );
            fprintf( '%s', varStr )
        end
        oldCompletedPercentage  = completedPercentage;

        ZI              = splMatrix(:,:,blockIndex);                    % extract current value to plot

        plotBG          = subplot( 'Position', [ 0 0 1 1 ] );
        set( plotBG, 'YDir', 'normal' );
        ambisonicsVideoColorMapFFMPEG( azi, elev, ZI, lowerLimit, upperLimit, v, contourEnable, ...
                                       valueTickEnable, cmap, aziMat, elevMat, tickSize, tickBold )

        writeVideo( videoFile, getframe( videoWindow ) );
        
    end
    
    close( videoFile )
    % copy background image to root folder and convert to requested format
    copyfile( backgroundImageFile, './', 'f' )

    if outputBlackWhite
        system( [ 'ffmpeg -y -i ' backgroundImageFileName ' -vf hue=s=0 ' backgroundImageFileName ] )
    end
    
    system( [ 'ffmpeg -i ' backgroundImageFileName ' -vf scale=' num2str( outputWidth ) ':' num2str( outputWidth/2 ) ' temp.png ' ] );
    delete( backgroundImageFileName )
    
    % create background video from image
    mapVideoFile            = [ videoTitle '.mp4' ];
    backgroundVideoOutput   = [ 'overlay_' mapVideoFile ];
    system( [ 'ffmpeg -loop 1 -i temp.png -c:v libx264 -t ' num2str( duration ) ' -pix_fmt yuv420p tempBackground.mp4' ] );
    delete( 'temp.png' )
    clear mapVideo 
    
    % generate overlay
    system( [ 'ffmpeg -y -i ' mapVideoFile ' -i tempBackground.mp4 -filter_complex " '...
    '[0:v]setpts=PTS-STARTPTS, scale=' num2str( outputWidth ) 'x' num2str( outputWidth/2 ) '[top]; '...
    '[1:v]setpts=PTS-STARTPTS, scale=' num2str( outputWidth ) 'x' num2str( outputWidth/2 ) '[bottom]; '...
    '[top]split[m][a]; [m][a]alphamerge[keyed]; [bottom][keyed]overlay=eof_action=endall" ' backgroundVideoOutput ] );        
    delete( 'tempBackground.mp4' )
    delete( mapVideoFile )
    movefile( backgroundVideoOutput, mapVideoFile );
        
else
    
    videoHeigth     = mapWindowWidth / 2;                                % set video window dimension to 2:1 for VR
    videoWindow   	= figure( 'units', 'pixel', 'innerposition', [ 0 0 mapWindowWidth videoHeigth  ], 'menubar', 'none' );
    
    % reshape diretions for colormap generation
    elev            = -videoHeigth/2 : resolutionMap/90*videoHeigth/2 : videoHeigth/2;
    azi             = -mapWindowWidth/2 : resolutionMap/180*mapWindowWidth/2 : mapWindowWidth/2;
    aziMat          = repmat( azi, size( elev, 2 ), 1 );
    elevMat         = repmat( elev, size( azi, 2 ), 1 )';
    
    for blockIndex = 1 : numbersFrame
        
        %%%%%%%%%%% PERCENTAGE PROCESSING %%%%%%%%%%%  
        completedPercentage     = round( blockIndex/numbersFrame*100 );
        if oldCompletedPercentage ~= completedPercentage
            for i = 1 : varStrLen
                fprintf( '\b' )
            end
            varStr              = sprintf( '%d%%', completedPercentage );
            varStrLen           = length( varStr );
            fprintf( '%s', varStr )
        end
        oldCompletedPercentage  = completedPercentage;
        
        ZI          = splMatrix( :, :, blockIndex );                        % extract current value to plot

        % subplot colormap
        plotBG      = subplot( 'Position', [ 0 0 1 1 ] );
        ambisonicsVideoColorMap( background, backgroundWidth, backgroundHeight, azi, elev, ZI, lowerLimit, upperLimit, ...
                                 plotBG, v, alpha, contourEnable, valueTickEnable, cmap, aziMat, elevMat, tickSize, tickBold, colorbarEnable )
                        
        writeVideo( videoFile, getframe( videoWindow ) );
        
    end
    
end

close( videoFile );                                                         % finalizes video 
close all


%% AUDIO PROCESSING
if autoselectOverlap == 0
    fprintf( '\n\nautoselectOverlap = 0 --> Audio disabled\n' )
    audioMono   = 0;
    ambix_o1    = 0;
    ambix_o3    = 0;    
end

if audioMono == 1 || ambix_o1 == 1 || ambix_o3 == 1

    % adjust to video duration
    ambixFormat             = ambixFormatAudio( 1 : duration * fs, : );

    videoFileName           = [ videoTitle '.mp4' ];
    
    % NORMALIZE AND SAVE
    if audioMono 
        
        % extract, normalize and save
        omniAudio           = ambixFormat( :, 1 );                          % keep only ambix omni mic
        maxAbsValue         = max( max( abs( omniAudio ) ) );               % search for absolute max value
        omniAudio           = omniAudio ./ maxAbsValue;                     % normalize between -1 and +1
        outputWavOmni       = [ recordedFileName( 1:end-4 ) '_omni.wav' ];
        fprintf( [ '\nSaving omni audio format in WAV file: ' outputWavOmni ] );
        audiowrite( outputWavOmni, omniAudio, fs, 'BitsPerSample', outBitsPerSample );
        fprintf( '\nOmni audio saved correctly.\n' );

        % add to video
        fprintf( '\nAdding omni audio\n' )
        outputFileName      = [ videoTitle '_omni.mp4' ];
        system( [ 'FFMPEG -i ' videoFileName ' -i ' outputWavOmni ' -map 0:v -map 1:a -c:v copy -shortest ' outputFileName ] );
        delete( outputWavOmni )

    end

    if ambix_o1
        
        % extract, normalize and save
        outputWavAmbix_o1   = [ recordedFileName( 1:end-4 ) '_ambix_o1.wav' ];  % define filename
        ambixAudio_o1       = ambixFormat( :, 1:4 );                        % keep only first four channels
        maxAbsValue         = max( max( abs( ambixAudio_o1 ) ) );           % search for absolute max value
        ambixAudio_o1       = ambixAudio_o1 ./ maxAbsValue;                 % normalize between -1 and +1
        fprintf( [ '\nSaving 1st order Ambix format in WAV file: ' outputWavAmbix_o1 ] );
        audiowrite( outputWavAmbix_o1, ambixAudio_o1, fs, 'BitsPerSample', outBitsPerSample );   % save file
        fprintf( '\nAmbix format saved correctly.\n' );
        
        % add to video
        fprintf( '\nAdding 1st order Ambix audio\n' )
        outputFileName      = [ videoTitle '_ambix_o1.mov' ];
        system( [ 'FFMPEG -i ' videoFileName ' -i ' outputWavAmbix_o1 ' -map 0:v -map 1:a -c:v copy -c:a copy -shortest ' outputFileName ] );
        delete( outputWavAmbix_o1 )
        
    end
    
    if ambix_o3
        
        % extract, normalize and save
        outputWavAmbix_o3   = [ recordedFileName( 1:end-4 ) '_ambix_o3.wav' ];  % define filename
        ambixAudio_o3       = ambixFormat( :, 1:16 );
        maxAbsValue         = max( max( abs( ambixAudio_o3 ) ) );               % search for absolute max value
        ambixAudio_o3       = ambixAudio_o3 ./ maxAbsValue;                     % normalize between -1 and +1
        fprintf( [ '\nSaving 3rd order Ambix format in WAV file: ' outputWavAmbix_o3 ] );
        audiowrite( outputWavAmbix_o3, ambixAudio_o3, fs, 'BitsPerSample', outBitsPerSample );   % save file
        fprintf( '\nAmbix format saved correctly.\n' );
        
        % add to video
        fprintf( '\nAdding 3rd order Ambix audio\n' )
        outputFileName       = [ videoTitle '_ambix_o3.mov' ];
        system( [ 'FFMPEG -i ' videoFileName ' -i ' outputWavAmbix_o3 ' -map 0:v -map 1:a -c:v copy -c:a copy -shortest ' outputFileName ] );
        delete( outputWavAmbix_o3 )
    
    end
       
    delete( videoFileName )
    
end


%% METADATA
if injectMetadata
    source              = './lib/SpatialMediaMetadataInjector.exe';
    copyfile( source, './', 'f' )
    system( 'SpatialMediaMetadataInjector.exe' )
end

delete( 'SpatialMediaMetadataInjector.exe' )
rmpath( './lib' )
fprintf( '\nEnded correctly.\n\n' )
