% SpatialNoiseVideoMap_sps_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 SPS.
% Background is a static image. It saves video overlay with pseudocolor SPL maps.
%
% INPUT
% - sps_format.mat -> file containing the SPS format generated by b_generateEncodedFormat
% - panoramic_360_image.jpg
%
% OUTPUT
% - acusticVideo.mp4
%
% Author: Daniel Pinardi.
% University of Parma, 2017-18
% 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
ambixMatrixFileName     = 'em32_ambix-o3_2018.wav';
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


%% IMAGE/VIDEO PARAMETERS
% pointing correction: + -> turn map to the right / - -> turn map to the left
horizontalTilt          = 0;        % degrees
surfpoint               = 120;      % number of horiz. interpolation points: increase up to 360 for horiz. correction

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
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

fprintf( '\nNoise elaboration started\n' );

% SPS-noise
filename                = encodedNoiseMatFileName;
disp( [ 'Loading SPS format file: ' filename ] );
load( filename );
spsFormat               = encodedFormat;
clear encodedFormat
virtualMicNumber        = size( spsFormat, 2 );                             % number of virtual microphones
if virtualMicNumber == 32
    prompt              = 'SPS type (0 -> capsule directions ; 1 -> balloon): ';
    sps32type           = input( prompt );
else
    sps32type           = 2;
end

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

% Check audio and load filtering matrix if necessary
[ arrayRec, ~ ]         = audioread( recordedNoiseFile );
inChNum                 = size( arrayRec, 2 );
if audioMono == 0 && ambix_o1 == 0 && ambix_o3 == 0 % audio is not to be calculated: output will be mute.
    fprintf ( 'Output video will be mute.\n' )
elseif audioMono || ambix_o1 || ambix_o3
    % audio has to be calculated: load filtering matrix
    ambixMatrix         = [ resourcesFolder ambixMatrixFileName ];
    disp( [ 'Loading Ambisonics matrix: ' ambixMatrix ] );
    [ ambixIR, ~ ]      = audioread( ambixMatrix );
    outChNum            = size( ambixIR, 2 );
    clear ambixIR
    [ ambixIR, fs ]     = load_ir( ambixMatrix, inChNum, outChNum ); % ambix 3rd order is [ m-IN x n-OUT x Nsamples ]
end


%% 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
    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


%% FULL SCALE LEVEL
fullScaleLevel          = 0.0004*( EM32gain^3 ) - 0.0078*( EM32gain^2 ) - 0.8821*( EM32gain ) + fullScaleLevel; % dB
clear EM32gain


%% A-weighting
if aWeightingEnable
    fprintf( 'Performing A-weighting.\n' );
    AWeighting          = weightingFilter( 'A-weighting', fs );
    for i = 1 : virtualMicNumber
        spsFormat( :, i ) = AWeighting( spsFormat( :, i ) );
    end 
    fprintf( 'A-weighting 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
spsMean                 = mean( spsFormat, 2 );                             % mean signal for fft plot

if trimEnable
    fprintf( 'Perform trimming\n\n' );
    fig             = figure;
    hax             = axes;
    plot( spsMean )
    xlim( [ 0 size( spsMean, 1 ) ] )
    xlabel( 'Samples' )
    grid on
    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 );
    spsFormat       = spsFormat( 1 + noiseStart : noiseEnd, : );        % trim noise
    if size( arrayRec, 2 ) ~= 0
        arrayRec     = arrayRec( 1 + noiseStart : noiseEnd, : );        % trim audio
    end
    spsMean             = mean( spsFormat, 2 );                             % update value

else % no trimming
    
    fprintf( 'Continuing without trimming.\n' );
    plot( spsMean );
    xlim( [ 0 length( spsMean ) ] )
    grid on
    xlabel( 'Samples' )
    title( 'Untrimmed signal' );
    fprintf( 'Check signal...then press any key to continue or Ctrl-C to abort\n\n' );
    pause
    
end
close all

spsLength               = size( spsFormat, 1 ); 
numbersFrame            = floor( spsLength / stepSize );                    % integer number of frames
noiseEnd                = stepSize * numbersFrame;                          % update value
spsFormat               = spsFormat( 1 : noiseEnd, : );                     % trim noise
if size( arrayRec, 2 ) ~= 0
    arrayRec             = arrayRec( 1 : noiseEnd, : );                       % trim audio
end


%% ZERO PADDING
zeroPadMatrix           = zeros( bufferSize/2, virtualMicNumber );          % center of buffer must be at 0
spsFormat               = [ zeroPadMatrix; spsFormat; zeroPadMatrix ];      % apply zero padding
spsMean                 = mean( spsFormat, 2 );                             % update
if size( arrayRec, 2 ) ~= 0
    zeroPadMatrix       = zeros( bufferSize/2, inChNum );                   % center of buffer must be at 0
    arrayRec             = [ zeroPadMatrix; arrayRec; zeroPadMatrix ];        % apply zero padding
end


%% Build 122 or 32 points grid
if virtualMicNumber == 32 && sps32type == 0
    fprintf( 'Building 32 virtual mics grid - EM directions\n\n' );
    [ extX, extY, extXI, extYI, XI, YI, X, Y ] = GridBuilderEM32( backgroundWidth, backgroundHeight, surfpoint, virtualMicNumber );
elseif virtualMicNumber == 32 && sps32type == 1
    fprintf( 'Building 32 virtual mics grid\n\n' );
    [ extX, extY, extXI, extYI, XI, YI, X, Y ] = GridBuilder32( backgroundWidth, backgroundHeight, surfpoint, virtualMicNumber );
else
    fprintf( 'Building 122 virtual mics grid\n\n' );
    [ extX, extY, extXI, extYI, XI, YI, X, Y ] = GridBuilder122( backgroundWidth, backgroundHeight, surfpoint, virtualMicNumber );
end


%% DATA PROCESSING
fprintf( [ 'NFFT value setted equal to bufferSize: ' num2str( bufferSize ) '\n'] )
nfft                    = bufferSize;

freq                    = linspace( 0, fs/2, nfft/2 + 1 );
tmp                     = abs( freq - bandpassLowCutFreq );
[ ~, idxLowBP ]         = min( tmp );                                           % index of closest value
tmp                     = abs( freq - bandpassHighCutFreq );
[ ~, idxHighBP ]        = min( tmp );
freqWindowing           = [ zeros( virtualMicNumber, idxLowBP ) ...
                            ones( virtualMicNumber, idxHighBP - idxLowBP + 1 ) ...
                            zeros( virtualMicNumber, ( nfft / 2 + 1 ) - idxHighBP - 1 ) ]';

if bandpassHighCutFreq == 0 || bandpassLowCutFreq == 0
    freqWindowing( freqWindowing == 0 ) = 1;
end

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

% 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;
    
    [ ZI ]              = frameProcessing( blockIndex, stepSize, bufferSize, spsFormat, nfft, freqWindowing, ...
                                           virtualMicNumber, fullScaleLevel, extX, extY, extXI, extYI, surfpoint, sps32type );
            
    splMatrix( :, :, blockIndex ) = ZI;

end

splMatrix               = single( splMatrix );

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

upperLimit              = max( max( max( splMatrix ) ) );                       % find absolute max level
lowerLimit              = upperLimit - noiseGate;                               % range to plot

if useFixedSplPlotBounds == 1
    v                   = fixedLowerLimit : contourSpan : fixedUpperLimit;      % iso-line levels
    splMatrix( splMatrix < fixedLowerLimit ) = 0;                           % set to 0 SPL values < lowerLimit
else                         
    v                   = lowerLimit+1 : contourSpan : upperLimit;              % iso-line levels
    splMatrix( splMatrix < lowerLimit ) = 0;                                % set to 0 SPL values < lowerLimit
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

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

if timeSignalEnable

    verticalLimit       = [ min( spsMean ) max( spsMean ) ];
    videoHeigth         = mapWindowWidth / ( 2 * 0.8 );                      % define video window
    videoWindow         = figure( 'units', 'pixel', 'innerposition', [ 0 0 mapWindowWidth videoHeigth  ], 'menubar', 'none' );

    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 ] );
        videoColorMap( background, backgroundWidth, backgroundHeight, XI, YI, ZI, lowerLimit, upperLimit, plotBG, ...
                       v, alpha, contourEnable, useFixedSplPlotBounds, fixedLowerLimit, fixedUpperLimit, valueTickEnable, ...
                       cmap, tickSize, tickBold, colorbarEnable );
        
        % subplot time signal
        plotTime                = subplot( 'Position', [ 0 0 1 0.2 ] );                 % subplot with time signal
        timeSignalPlot( spsMean, 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' );

    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' );
        videoColorMapFFMPEG( XI, YI, ZI, lowerLimit, upperLimit, v, contourEnable, useFixedSplPlotBounds, fixedLowerLimit, ...
                             fixedUpperLimit, valueTickEnable, cmap, 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 ];
    duration                = stepSize * numbersFrame / fs;
    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' );

    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 ] );
        videoColorMap( background, backgroundWidth, backgroundHeight, XI, YI, ZI, lowerLimit, upperLimit, plotBG, ...
                       v, alpha, contourEnable, useFixedSplPlotBounds, fixedLowerLimit, fixedUpperLimit, valueTickEnable, ...
                       cmap, 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
    
    fprintf( '\nConvolving array rec with Ambix FIR matrix.\n' );
    limitLength                 = 200000;
    fileToProcessLength         = length( arrayRec );

    if  fileToProcessLength    <= limitLength
        % Reshaping matrix to use convolve library function: matrix_conv
        fileToProcessMatrix     = reshape( arrayRec', 1, inChNum, fileToProcessLength );   
        encodedFormat           = matrix_conv( fileToProcessMatrix, ambixIR );   % convolving with SPS matrix
    else
        % Reshaping matrix to use convolve library function: oa_multichannel_conv
        fileToProcessMatrix     = reshape( arrayRec', inChNum, fileToProcessLength, 1 );   
        encodedFormat           = oa_multichannel_conv( fileToProcessMatrix, ambixIR );% convolving with SPS matrix
    end
    encodedFormat               = squeeze( encodedFormat )';                    % eliminate exceeding matrix dimension
    
    % readjust length
    for in = 1 : inChNum
        for out = 1 : outChNum
            [ ~, delayInOut ]   = max( abs( ambixIR( in, out, : ) ) );
            delay( in, out )    = delayInOut;
        end
    end
    delay                       = round( mean( delay, 'all' ) );
    encodedFormat               = encodedFormat( delay+1 : delay + length( arrayRec ), : );

    if bandpassLowCutFreq ~= 0 && bandpassHighCutFreq ~= 0
    
        fprintf( 'Performing ambix audio filtering.\n' );
        [ b, a ]                = butter( 10, [ bandpassLowCutFreq bandpassHighCutFreq ] / ( fs/2 ), 'bandpass' );
        encodedFormat           = filter( b, a, encodedFormat );                % perform filtering

    end
    
    % NORMALIZE AND SAVE
    if audioMono % extract, normalize and save
        omniAudio           = encodedFormat( :, 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' );
    end

    if ambix_o1
        outputWavAmbix_o1   = [ recordedFileName( 1:end-4 ) '_ambix_o1.wav' ];  % define filename
        ambixAudio_o1       = encodedFormat( :, 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' );
    end
    
    if ambix_o3
        outputWavAmbix_o3   = [ recordedFileName( 1:end-4 ) '_ambix_o3.wav' ];  % define filename
        maxAbsValue         = max( max( abs( encodedFormat ) ) );               % search for absolute max value
        encodedFormat       = encodedFormat ./ maxAbsValue;                     % normalize between -1 and +1
        fprintf( [ '\nSaving 3rd order Ambix format in WAV file: ' outputWavAmbix_o3 ] );
        audiowrite( outputWavAmbix_o3, encodedFormat, fs, 'BitsPerSample', outBitsPerSample );   % save file
        fprintf( '\nAmbix format saved correctly.\n' );
    end
       
end

%% ADD AUDIO TO OVERLAY
videoFileName           = [ videoTitle '.mp4' ];

if audioMono
	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
    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
    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

if audioMono || ambix_o1 || ambix_o3
    delete( videoFileName )
end

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

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