% SpatialNoiseVideoMap_sps_video.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. Coherence elaboration is supported.
% Background is a video. It saves video overlay with pseudocolor SPL maps.
%
% INPUT
% - sps_format.mat -> file containing the SPS format generated by b_generateEncodedFormat
% - panoramic_360_video.mp4
%
% OUTPUT
% - acusticNoiseVideo.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';
audioCoherence          = 0;        % 1 to enable coherent 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        = 0;        % 1 to enable: plot time signal under video maps. Disable for VR
bufferSize              = 16384;    % [samples]
nfft                    = 16384;    % [samples]

% coherence analysis
totalTFenable           = 1;        % if 1, use only one transfer function calculated on whole signal. Good for time invariant system!
plotTFenable            = 0;

% reference signals to be used for coherence analysis
ref_01                  = '';
ref_02                  = 'static';
filtRef                 = 0;        % if 0 -> filt SPS, if 1 -> filt REF
% TIP: mimo + totTF -> filtRef = 0 ( -> filtSPS )
% TIP: mimo + movTF -> filtRef = 1

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

% video
videoOffsetCorrection   = 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 
outputWindowWidth       = 1200;     % 0 to autoselect maximum possible, depending on screen resoution. Aspect ratio will be 2:1
keepVideoNativeRes      = 1;        % 1 -> video keep resolution (scaled map), 0 -> video scaled to map

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

% Surf mapping
useFixedSplPlotBounds   = 0;        % if 1 each surface plot has color normalization with fixedLowerLimit and fixedUpperLimit
fixedLowerLimit         = 80;       % user defined lower limit
fixedUpperLimit         = 90;       % user defined upper limit
noiseGate               = 8;       % 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 backgroundImageFile backgroundImageFileName invSweepFile invSweepFileName irsFolder irsMatFileName ...
irsWavFileName mapsFolderIrs recordedSweepFile recordedSweepFolder spsNoiseWavFileName spsNoiseFolder ...
spsFilterMatrixFile spsFilterMatrixFileName spsIrsFolder spsIrsMatFileName spsIrsWavFileName outputFolder

fprintf( '\nCoherent Noise elaboration started\n' );
fprintf('\nLoading reference file\n' );
[ refRec, fsReference ] = audioread( recordedReferenceSignalFile );         % read reference recorded files

% extract desired reference signal, depending on previous selection
[ refMicNumber, referenceName, referenceMic ] = extractRef( refRec, ref_01, ref_02 );
refMicAudio             = referenceMic;
% enable MIMO analysis
if refMicNumber > 1 
    referenceName       = 'multipleRef';
    multipleRef         = 1;
    fprintf( '\nMIMO coherence enabled\n' )
else
    multipleRef         = 0;
end
clear refRec
   
% 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° video
disp( ['Loading background video file: ' backgroundVideoFile ] );
backgroundVideo         = VideoReader( backgroundVideoFile );               % load panoramic video
frameRate               = backgroundVideo.FrameRate;                        % set colormap framerate equal to pano video
backgroundHeightOrig    = backgroundVideo.Height;                           % get pano video height
backgroundWidthOrig     = backgroundVideo.Width;                            % get pano video width

% check monitor dimension and resize if necessary
[ outputWindowWidth ]   = screenSize( outputWindowWidth, timeSignalEnable );
resizeFactor            = outputWindowWidth / backgroundWidthOrig;          % define resizing factor for pano video
backgroundHeight        = round( backgroundHeightOrig * resizeFactor );     % new pano video value
backgroundWidth         = round( backgroundWidthOrig * resizeFactor );      % new pano video value

if videoOffsetCorrection
    [ xImage, yImage ]  = meshgrid( 1 : backgroundWidth, 1 : backgroundHeight );
    middleIndexY        = round( backgroundWidth / 2 );                     % Find index of middle element
    elevation           = yImage / middleIndexY * 180 - 90;                 % elevation from -90 to +90
    yCorrected          = ( atand( ( radiusCorrection * sind(  elevation ) - verticalOffset ) ...
                          ./ ( radiusCorrection * cosd(  elevation ) ) ) + 90) / 180 * middleIndexY;
    correctionMatrix    = cat( 3, xImage, yCorrected ); 
    resamp              = makeresampler( 'linear', 'fill' );
end   
clear yCorrected elevation middleIndexY xImage yImage radiusCorrection verticalOffset
 
% Audio track of video
[ audioVideo, fsVideo ] = audioread( backgroundVideoFile );                 % load audio track from recorded pano video
audioVideo              = audioVideo( :, 1 );
clear backgroundVideoFile

% array recorded
[ arrayRec, ~ ]         = audioread( recordedNoiseFile );
inChNum                 = size( arrayRec, 2 );

% Check audio and load filtering matrix if necessary
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 ]
    if audioCoherence == 1 && totalTFenable == 0
        audioCoherence  = 0;
        fprintf( 'Total TF disabled --> coherent audio disabled\n' ) 
    end
end


%% FULL SCALE LEVEL: EMIB gain is not linear. Equation has been calculated previously.
fullScaleLevel          = 0.0004*( EM32gain^3 ) - 0.0078*( EM32gain^2 ) - 0.8821*( EM32gain ) + fullScaleLevel; % [dB]
clear EM32gain


%% Evaluate stepSize and overlap
stepSize                = floor( fs / frameRate );                          % normally 1600 or 1601 samples
overlapVideo            = ( 1 - stepSize / bufferSize ) * 100;
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( overlapVideo ) '. \n' ] );


%% SYNC REF
spsLength               = length( spsFormat );
% check temporal alignment
if refMicNumber == 1
    
    delay               = syncRef( spsFormat( :, 1 ), referenceMic );
    if delay > 0
        referenceMic    = referenceMic( delay : end );
        refMicAudio     = refMicAudio( delay : end );
    else
        referenceMic    = [ zeros( abs( delay ), 1 ) ; referenceMic ];
        refMicAudio     = [ zeros( abs( delay ), 1 ) ; refMicAudio ];
    end
    % adjust signals length
    refMicLength        = length( referenceMic );
    if refMicLength > spsLength
        referenceMic    = referenceMic( 1 : length( spsFormat ), : );
        refMicAudio     = refMicAudio( 1 : length( spsFormat ), : );
    else
        referenceMic    = [ referenceMic ; zeros( spsLength - refMicLength, refMicNumber ) ];
        refMicAudio     = [ refMicAudio ; zeros( spsLength - refMicLength, refMicNumber ) ];
    end

else
    
    referenceMicTemp            = zeros( spsLength, 2 );
    refMicAudioTemp             = zeros( spsLength, 2 );
    for i = 1 : 2
        
        referenceMicTemp_i      = referenceMic( :, i );
        refMicAudioTemp_i       = refMicAudio( :, i );
        delay                   = syncRef( spsFormat( :, 1 ), referenceMicTemp_i );
        if delay > 0
            referenceMicTemp_i  = referenceMicTemp_i( delay : end );
            refMicAudioTemp_i   = refMicAudioTemp_i( delay : end );
        else
            referenceMicTemp_i  = [ zeros( abs( delay ), 1 ) ; referenceMicTemp_i ];
            refMicAudioTemp_i   = [ zeros( abs( delay ), 1 ) ; refMicAudioTemp_i ];
        end
        % adjust signals length
        refMicLength            = length( referenceMicTemp_i );
        if refMicLength > spsLength
            referenceMicTemp_i  = referenceMicTemp_i( 1 : length( spsFormat ), : );
            refMicAudioTemp_i   = refMicAudioTemp_i( 1 : length( spsFormat ), : );
        else
            referenceMicTemp_i  = [ referenceMicTemp_i ; zeros( spsLength - refMicLength, 1 ) ];
            refMicAudioTemp_i   = [ refMicAudioTemp_i ; zeros( spsLength - refMicLength, 1 ) ];
        end
        referenceMicTemp( :, i ) = referenceMicTemp_i;
        refMicAudioTemp( :, i )  = refMicAudioTemp_i;
    end
    
    referenceMic        = referenceMicTemp;
    refMicAudio         = refMicAudioTemp;
    clear referenceMicTemp referenceMicTemp_i refMicAudioTemp_i refMicAudioTemp
    
end

    
%% SYNC AUDIO / VIDEO
fprintf( 'Sync audio - video.\n' );
[ delay ]               = syncRef( arrayRec( :, 1 ), audioVideo );
fprintf( [ 'Delay is: ' num2str( delay ) ' samples.\n' ] );
% apply delay
if delay >= 0
    timeToSkip          = delay / fsVideo;
else
    timeToSkip          = 0;
    spsFormat           = spsFormat( -delay+1 : end, : );                   % apply delay
    arrayRec            = arrayRec( -delay+1 : end, : );                    % apply delay
    referenceMic        = referenceMic( -delay+1 : end, : );                % apply delay   
    refMicAudio         = refMicAudio( -delay+1 : end, : );                 % apply delay 
    spsLength           = length( spsFormat );                              % update value
    audioVideoLength    = length( audioVideo );                             % update value
    if spsLength > audioVideoLength
        spsFormat       = spsFormat( 1 : audioVideoLength, : );
        arrayRec        = arrayRec( 1 : audioVideoLength, : );
        referenceMic    = referenceMic( 1 : audioVideoLength, : );
        refMicAudio     = refMicAudio( 1 : audioVideoLength, : );
    end
end 
clear audioVideo


%% A-weighting
if aWeightingEnable
    fprintf( 'Performing A-weighting.\n' );
    AWeighting          = weightingFilter( 'A-weighting', fs );
    for i = 1 : virtualMicNumber
        spsFormat( :, i ) = AWeighting( spsFormat( :, i ) );
    end  
    for i = 1 : refMicNumber
        referenceMic( :, i ) = AWeighting( referenceMic( :, i ) );
    end
    fprintf( 'A-weighting ended.\n' );
end

spsMean                 = mean( spsFormat, 2 );                             % mean signal for fft plot


%% TRIM
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
    referenceMic        = referenceMic( 1 + noiseStart : noiseEnd, : );         % trim references
    refMicAudio         = refMicAudio( 1 + noiseStart : noiseEnd, : );          % trim references
    if any( arrayRec, 'all' ) ~= 0
        arrayRec        = arrayRec( 1 + noiseStart : noiseEnd, : );            % trim audio
    end

    spsMean             = mean( spsFormat, 2 );                             % update value

else % no trimming
    
    fprintf( 'Continuing without trimming.\n' );   
    noiseStart          = 0;
    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
referenceMic            = referenceMic( 1 : noiseEnd, : );                  % trim references
refMicAudio             = refMicAudio( 1 : noiseEnd, : );
if any( arrayRec, 'all' ) ~= 0
    arrayRec            = arrayRec( 1 : noiseEnd, : );                      % trim audio
end
timeToSkip              = timeToSkip + noiseStart / fsVideo;                % update value
fprintf( [ 'Skipped video: ' num2str( timeToSkip ) ' s.\n' ] );


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

zeroPadMatrix           = zeros( bufferSize/2, refMicNumber );              % center of buffer must be at 0
referenceMic            = [ zeroPadMatrix; referenceMic; zeroPadMatrix ];   % apply zero padding
refMicAudio             = [ zeroPadMatrix ; refMicAudio ; zeroPadMatrix ];
if any( arrayRec, 'all' ) ~= 0
    zeroPadMatrix       = zeros( bufferSize/2, inChNum );                   % center of buffer must be at 0
    arrayRec            = [ zeroPadMatrix; arrayRec; zeroPadMatrix ];       % apply zero padding
end

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


%% FILTERING WINDOW
if nfft < bufferSize || totalTFenable == 0
    nfft                = bufferSize;
    fprintf( [ 'NFFT value setted equal to bufferSize: ' num2str( bufferSize ) '\n'] )
end

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

if multipleRef
    freqWindowing       = repmat( freqWindowing, [ 1, 1, 2 ] );
end


%% TRANSFER FUNCTION PROCESSING
if totalTFenable
    
    transferFunction    = tfestimate( referenceMic, spsFormat, rectwin( nfft ), [], nfft );
    
    if plotTFenable
        figure
        if multipleRef
            for i = 1 : 2
                tf_i    = transferFunction( :, :, i );
                tfIR    = [ zeros( 1, virtualMicNumber ) ; tf_i( 2 : end - 1, : ); zeros( 1, virtualMicNumber ) ; ...
                            flipud( conj( tf_i( 2 : end - 1, : ) ) ) ];
                tfIR    = fftshift( ifft( tfIR ) );      
                subplot( 2, 1, i )
                plot( tfIR )
                grid on
                xlim( [ 0 size( tfIR, 1 ) ] ) 
                xlabel( 'Samples' )
                if i == 1
                    title( 'TF plot - ref 1' )
                else
                    title( 'TF plot - ref 2' )
                end
            end
        else
            tfIR        = [ zeros( 1, virtualMicNumber ) ; transferFunction( 2 : end - 1, : ); zeros( 1, virtualMicNumber ) ; ...
                            flipud( conj( transferFunction( 2 : end - 1, : ) ) ) ];
            tfIR        = fftshift( ifft( tfIR ) );      
            plot( tfIR )
            grid on
            xlim( [ 0 size( tfIR, 1 ) ] ) 
            xlabel( 'Samples' )
            title( 'TF plot' )
        end
    end
end


%% FRAME PROCESSING
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 )

if totalTFenable
    
    transferFunction            = transferFunction .* freqWindowing;
    
    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_totalCoh( blockIndex, stepSize, bufferSize, spsFormat, referenceMic, filtRef, nfft, ...
                                                            virtualMicNumber, fullScaleLevel, extX, extY, extXI, extYI, surfpoint, ...
                                                            transferFunction, sps32type );

        splMatrix( :, :, blockIndex )   = ZI;

    end

else
    
    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_movingCoh( blockIndex, referenceMic, spsFormat, stepSize, filtRef, bufferSize, nfft, ...
                                                             freqWindowing, virtualMicNumber, fullScaleLevel, extX, extY, extXI, ...
                                                             extYI, surfpoint, sps32type );

        splMatrix( :, :, blockIndex )   = ZI;

    end
    
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' );
referenceName           = deblank( referenceName );
videoTitle              = [ recordedFileName( 1:end-4 ) '_' referenceName ];    % name of file video
videoFile               = VideoWriter( videoTitle, 'MPEG-4' );                  % initializes video
videoFile.FrameRate     = frameRate;                                            % set frame rate
videoFile.Quality       = 100;                                                  % set quality
open( videoFile );
backgroundVideo.CurrentTime = timeToSkip;                                   % index of starting point
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         = outputWindowWidth / ( 2 * 0.8 );                      % define video window
    videoWindow         = figure( 'units', 'pixel', 'innerposition', [ 0 0 outputWindowWidth 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

        % take each video frame as background image
        plotBG          = subplot( 'Position', [ 0 0.2 1 0.8 ] );
        background      = readFrame( backgroundVideo );
        background      = imresize( background, resizeFactor );             % resize

        if outputBlackWhite
            background  = (1+alpha)*flipud( background( :, :, [ 1 1 1 ] ) ); % flip and conversion to black white
        else
            background  = (1+alpha)*flipud( background );                   % flip 
        end
        
        if videoOffsetCorrection
            background  = uint8( tformarray( background, [], resamp, [2 1], [1 2], [], correctionMatrix, .3 ) );
        end
        
        % subplot colormap
        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
    
    close( videoFile );                                                     % finalizes video
    close all
    
else
    
    videoHeigth         = outputWindowWidth / 2;                            % set video window dimension to 2:1 for VR
    videoWindow         = figure( 'units', 'pixel', 'innerposition', [ 0 0 outputWindowWidth videoHeigth  ], 'menubar', 'none' );
    
    if videoOffsetCorrection
        
        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

            % take each video frame as background image
            plotBG          = subplot( 'Position', [ 0 0 1 1 ] );
            background      = readFrame( backgroundVideo );
            background      = flipud( background( :, :, [ 1 1 1 ] ) );      % conversion to black white
            background      = imresize( background, resizeFactor );         % resize
            background      = uint8( tformarray( background, [], resamp, [2 1], [1 2], [], correctionMatrix, .3 ) );
            
            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
        close( videoFile );                                                 % finalizes video 
        close all
        
    else
        
        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 );                                                 % finalizes video 
        close all
        
        
        %% GENERATE MUTE OVERLAY
        backgroundVideo         = [ videoTitle '.mp4' ];
        backgroundVideoOutput   = [ 'overlay_' videoTitle '.mp4' ];
        duration                = round( numbersFrame / frameRate );
        video                   = [ './resources/' backgroundVideoFileName ];
        copyfile( video, './', 'f' )
          
        % conversion of background from color to B/W
        if outputBlackWhite
            system( [ 'ffmpeg -i ' backgroundVideoFileName ' -ss ' num2str( timeToSkip ) ' -t ' num2str( duration ) ' -vf hue=s=0 -an -crf 22 -preset ultrafast tempBackground.mp4' ] );
        else 
            system( [ 'ffmpeg -i ' backgroundVideoFileName ' -ss ' num2str( timeToSkip ) '  -t ' num2str( duration ) ' -an -crf 22 -preset ultrafast tempBackground.mp4' ] );
        end
        
        % perform overlaying
        if keepVideoNativeRes
            system( [ 'ffmpeg -i ' backgroundVideo ' -i tempBackground.mp4 -filter_complex " '...
            '[0:v]setpts=PTS-STARTPTS, scale=' num2str( backgroundWidthOrig ) 'x' num2str( backgroundHeightOrig ) '[top]; '...
            '[1:v]setpts=PTS-STARTPTS, scale=' num2str( backgroundWidthOrig ) 'x' num2str( backgroundHeightOrig ) '[bottom]; '...
            '[top]split[m][a]; [m][a]alphamerge[keyed]; [bottom][keyed]overlay=eof_action=endall" -crf 22 -preset ultrafast ' backgroundVideoOutput ] );
        else
            
            if backgroundHeightOrig ~= backgroundWidthOrig / 2 

                fprintf( [ 'Background video aspect ratio is not 2:1. Processing rescaling...\n\n' ] )
                y_target                = backgroundWidthOrig / 2;
                Y_video                 = ( backgroundHeightOrig - y_target ) / 2;
                copyfile 'tempBackground.mp4' 'tempBackgroundTemp.mp4'
                delete( 'tempBackground.mp4' )
                system( [ 'ffmpeg -i tempBackgroundTemp.mp4 -filter:v "crop=' num2str( backgroundWidthOrig ) ':' num2str( y_target ) ':0:' num2str( Y_video ) ...
                          '" -c:a copy -crf 22 -preset ultrafast tempBackground.mp4' ] );
                delete( 'tempBackgroundTemp.mp4' )

            end
            
            % check videomap dimensions
            videomap                = VideoReader( backgroundVideo );       % load video map
            videoMapH               = videomap.Height;                      % get pano video map height
            videoMapW               = videomap.Width;                       % get pano video map width
            
            system( [ 'ffmpeg -i ' backgroundVideo ' -i tempBackground.mp4 -filter_complex " '...
            '[0:v]setpts=PTS-STARTPTS, scale=' num2str( videoMapW ) 'x' num2str( videoMapH ) '[top]; '...
            '[1:v]setpts=PTS-STARTPTS, scale=' num2str( videoMapW ) 'x' num2str( videoMapH ) '[bottom]; '...
            '[top]split[m][a]; [m][a]alphamerge[keyed]; [bottom][keyed]overlay=eof_action=endall" -crf 22 -preset ultrafast ' backgroundVideoOutput ] );
        
        end
        delete( backgroundVideoFileName )
        delete( 'tempBackground.mp4' )
        delete( backgroundVideo )
        
    end
    
end


%% AUDIO PROCESSING
if audioMono == 1 || ambix_o1 == 1 || ambix_o3 == 1
    
    if audioCoherence
        
        transferFunction        = tfestimate( refMicAudio, arrayRec, rectwin( nfft ), [], nfft );
        fade                    = nfft / 8;              % length of fadeIn & fadeOut
        winType                 = 'hann';
        eval( sprintf( 'fadein = %s(%d);', winType, 2 * fade ) )
        fadein                  = fadein( 1 : fade );
        eval( sprintf( 'fadeout = %s(%d);', winType, 2 * fade ) )
        fadeout                 = fadeout( fade + 1 : end );
        win                     = [ fadein ; ones( nfft - 2 * fade, 1 ) ; fadeout ];
        winMatrix               = repmat( win, [ 1 inChNum ] );
        if multipleRef
            fprintf( '\nFiltering array recording with static MIMO Transfer Function\n' )
            for i = 1 : refMicNumber
                tf_i            = transferFunction( :, :, i );
                tf_i            = [ zeros( 1, inChNum ) ; tf_i( 2 : end - 1, : ); zeros( 1, inChNum ) ; flipud( conj( tf_i( 2 : end - 1, : ) ) ) ];
                tfIR            = fftshift( ifft( tf_i ) ); 
                tfIR            = tfIR .* winMatrix;
                refMicAudioMat  = repmat( refMicAudio( :, i ), [ 1, inChNum ] );
                arrayRec_i( :, :, i ) = fftfilt( tfIR, refMicAudioMat );
            end
            arrayRec            = sum( arrayRec_i, 3 );
        else
            fprintf( '\nFiltering array recording with static Transfer Function\n' )
            transferFunction    = [ zeros( 1, inChNum ) ; transferFunction( 2 : end - 1, : ); zeros( 1, inChNum ) ; ...
                                    flipud( conj( transferFunction( 2 : end - 1, : ) ) ) ];
            tfIR                = fftshift( ifft( transferFunction ) ); 
            tfIR                = tfIR .* winMatrix;
            refMicAudioMat      = repmat( refMicAudio, [ 1, inChNum ] );
            arrayRec            = fftfilt( tfIR, refMicAudioMat );
        end
            
    end
    
    % encoding
    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 ), : );

    % bandpass filtering    
    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
    
    if audioCoherence

        if ambix_o1 || ambix_o3
            encodedFormat( :, 1 )  = -encodedFormat( :, 1 );
            encodedFormat( :, 3 )  = -encodedFormat( :, 3 );
        end
        
        if ambix_o3
            encodedFormat( :, 5 )  = -encodedFormat( :, 5 );
            encodedFormat( :, 7 )  = -encodedFormat( :, 7 );
            encodedFormat( :, 9 )  = -encodedFormat( :, 9 );
            encodedFormat( :, 11 ) = -encodedFormat( :, 11 );
            encodedFormat( :, 13 ) = -encodedFormat( :, 13 );
            encodedFormat( :, 15 ) = -encodedFormat( :, 15 );  
        end

    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
if timeSignalEnable || videoOffsetCorrection
    videoFileName       = [ videoTitle '.mp4' ];
else
    videoFileName       = backgroundVideoOutput;
end

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