% SpatialIRsMap_sps.m Spatial acoustic static mapping of an IR recording
% This script implements acoustic source mapping of an IR using SPS format 
% got from a spherical microphone array. It shows Matlab figures and saves 
% png files containing pseudocolor SPL maps.
%
% INPUT
% - sps_format.mat -> file containing the SPS format for IR, generated by GenerateSpsFormat.m
% - panoramic_360_image.jpg
%
% OUTPUT
% - acousticMap.png
%
% Authors: Daniel Pinardi, Marco Binelli.
% University of Parma, 2017.
% Last review: 05/2023, Daniel Pinardi

clear variables
close all
clc

addpath( './lib' )                                                          % add to path the library with needed common files
a_configureFileNamesAndDirectories                                          % configure filenames, directories, resources

set( 0, 'DefaultTextInterpreter', 'none' );                                 % for plot title visualization
warning off                                                                 % to avoid warnings on image size


%% AUDIO PARAMETERS
% set fullscale value
% EM32_to_SPS32_2018    = 119
% EM32_to_SPS122_2018   = 122
fullScaleLevel          = 130;                                              % dB

% set EM32 gain to adjust FS level. Write -1 to avoid.
EM32gain                = -1;                                               % integer from 0 to 30

% Trimming
trimEnable              = 1;                                                % enable trimming if 1

% Signal processing
plotSignalEnable        = 0;
aWeightingEnable        = 1;                                                % enable to perform A filtering
nfft                    = 16384;										    % resolution of FFT
overlap                 = 0.5;												% Default is 0.5

% LOW PASS FILTER PARAMETERS -> 0 to disable
lowPassCutFreq          = 0;                                                % Cutoff Frequency [Hz]

% BAND PASS FILTER PARAMETERS -> 0 (one of the two) to disable
bandpassLowCutFreq      = 0;                                                % Low Cutoff Frequency [Hz]
bandpassHighCutFreq     = 400;                                              % High Cutoff Frequency [Hz]

% OCTAVE BAND FILTER PARAMETERS -> 0 (one of the two) to disable
% Rounded center freq are:
% [ 31 63 125 250 500 1000 2000 4000 8000 16000 ]
% [  1  2   3   4   5    6    7    8    9    10 ]
firstOctaveBand         = 0;                                                % select the first band to analyze
lastOctaveBand          = 8;                                                % select the last band to analyze


%% IMAGE PARAMETERS
% pointing correction: + -> to the right / - -> to the left
horizontalTilt          = 0; % degrees

% correction of the offset between optical and acoustic centers
imageOffsetCorrection   = 0;                                                % 1 -> enable
radiusCorrection        = 1.5;      % [m]
verticalOffset          = 0.2;      % [m] 
outputWindowWidth       = 2048;                                             % desired width of background image

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

% Surf mapping
useFixedSplPlotBounds   = 0;                                                % if 0 each surface plot has color normalization 
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; 0 to disable. Show isolevel lines.
contourSpan             = 1; 												% [dB] Range of isolevel lines.   
valueTickEnable         = 0;												% 1 to enable; 0 to disable. Show values of isolevel lines.
tickSize                = 16;												% Font dimension.
tickBold                = 1;												% 1 to enable; 0 to disable. Font appearance.
colorbarEnable          = 1;												% 1 to enable; 0 to disable. Show colorbar.


%% LOAD DATA FROM FILES
% define filename to load and folder output
mapsFolder              = mapsFolderIrs;
filename                = encodedIrsMatFileName;
fprintf('\nIrs elaboration started\n');
 
% SPS-noise
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 );


%% IMAGE PROCESSING
% constant and other parameters
surfpoint               = 360;
correctionIndex         = round( horizontalTilt * surfpoint / 360 );

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

% correction of offset between acoustic center and image center
if imageOffsetCorrection == 1
    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

if outputBlackWhite == 1
    background          = flipud( backgroundImage( :, :, [ 1 1 1 ] ) ); % flip and conversion to black white
else
    background          = flipud( backgroundImage );                    % flip
end


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


%% FULL SCALE LEVEL
if EM32gain >= 0 && EM32gain <= 30
	fullScaleLevel		= 0.0004*( EM32gain^3 ) - 0.0078*( EM32gain^2 ) - 0.8821*( EM32gain ) + fullScaleLevel; % dB
end


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

if trimEnable == 1
    
    fig                 = figure( 'windowstate', 'maximized' );
    hax                 = axes;
    plot( spsMean( :, 1 ) )
    xlim( [ 0 size( spsMean, 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

	spsFormat			= spsFormat( 1 + noiseStart : noiseEnd, : );      % trim noise
    
else % no trimming
    
    fprintf( 'Continuing without trimming.\n' );
	figure( 'windowstate', 'maximized' ) 
    plot( spsMean );
    xlim( [ 0 size( spsMean, 1 ) ] )
    xlabel( 'Samples' )
    grid on
    title( 'Untrimmed signal' );
    fprintf( 'Check signal...then press any key to continue or Ctrl-C to abort\n\n' );
    pause
    
end
close all


%% FILTERING PROCESS
freqAxis                = linspace( 0, fs/2, nfft/2+1 );

% Complex sum
complexSum              = computeComplexSum( spsFormat, nfft, overlap );
                       
% Total unfiltered SPL
[ splTotal, limitFreq ] = computeSPLunfiltered( complexSum, fullScaleLevel, spsMean, fs, nfft, overlap, plotSignalEnable);
     
% Lowpass filtered SPLG
if lowPassCutFreq ~= 0
    [ splLowPassFilt ]  = computeSPLlowpassfilt( complexSum, fullScaleLevel, freqAxis, lowPassCutFreq, spsMean, ...
                                                 fs, nfft, overlap, plotSignalEnable, limitFreq );
end

% Bandpass filtered SPL
if bandpassLowCutFreq ~= 0 && bandpassHighCutFreq ~= 0
    [ splBandPassFilt ] = computeSPLbandpassfilt( complexSum, fullScaleLevel, freqAxis, bandpassLowCutFreq, bandpassHighCutFreq, ...
                            spsMean, fs, nfft, overlap, plotSignalEnable, limitFreq );
end

% Octave band filtered SPL
if firstOctaveBand ~= 0 && lastOctaveBand ~= 0
    [ splOctaveFilt, centerFreq ] = computeSPLoctavefilt( fs, firstOctaveBand, lastOctaveBand, plotSignalEnable, ...
                                freqAxis, complexSum, spsMean, nfft, overlap, limitFreq, fullScaleLevel );
else
    centerFreq          = 0;
end

% Check plot if enabled
if plotSignalEnable == 1
    fprintf( 'Check filtering...then press any key to continue or Ctrl-C to abort.\n\n' );
    pause
    close all
end


%% GENERATION OF SPL MATRIX
% first row             -> total value of SPL (without filtering)
% second row            -> total value of SPL (with low pass filtering)
% third row             -> total value of SPL (with band pass filtering)
% fourth to 'x' row     -> total value of SPL for each octave band computed
spl                     = splTotal;
if exist( 'splLowPassFilt', 'var' )
    spl                 = [ spl; splLowPassFilt ];
end
if exist( 'splBandPassFilt', 'var' )
    spl                 = [ spl; splBandPassFilt ];
end
if exist( 'splOctaveFilt', 'var' )
    spl                 = [ spl; splOctaveFilt ];
end

% extend spl value to duplicated points at top and bottom. Only uniform
% division 32 and 122 can be extended, since they have virtual mics at poles.
if virtualMicNumber == 122 || ( virtualMicNumber == 32 && sps32type == 1 )
    [ spl ]             = extendSPLdir( spl );
    virtualMicNumberExtended = size( spl, 2 );
elseif virtualMicNumber == 64
    [ spl ]             = extendSPLdir64( spl );
    virtualMicNumberExtended = size( spl, 2 );
else 
    virtualMicNumberExtended = 0;
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 );
elseif virtualMicNumber == 64
	fprintf( 'Building 64 virtual mics grid\n\n' );
	[ extX, extY, extXI, extYI, XI, YI, X, Y ] = GridBuilder64( backgroundWidth, backgroundHeight, surfpoint, virtualMicNumber );
elseif virtualMicNumber == 122
    fprintf( 'Building 122 virtual mics grid\n\n' );
    [ extX, extY, extXI, extYI, XI, YI, X, Y ] = GridBuilder122( backgroundWidth, backgroundHeight, surfpoint, virtualMicNumber );
end


%% PLOT
numberPlot              = size(spl, 1);                                     % number of maps to be computed
refName                 = '';

% personalized colorscale
cmap                = jet( 64 );
cmap( 1, : )        = [ 0 0 0 ]; % black full scale

% show grid over background image and save image
figure( numberPlot + 1 );
gridOverImage( backgroundImage, virtualMicNumber, virtualMicNumberExtended, mapsFolder, X, Y, sps32type );

% plot acoustic map
% 1                 -> total unfiltered SPL
% 2                 -> total SPL low pass filtered
% 3                 -> total SPL band pass filtered
% from 4 to end     -> SPL with octave filtering

for y = 1 : numberPlot
    splPlot             = spl( y, : );                                      % every cycle takes one row

    [ ZI ]              = splInterpolation( double( splPlot ), virtualMicNumber, extX, extY, extXI, extYI, surfpoint, sps32type );
    
    if correctionIndex > 0  
        ZI              = ZI( :, [ correctionIndex+1 : surfpoint, 1:correctionIndex ] );
    elseif correctionIndex < 0
        ZI              = ZI( :, [ surfpoint+correctionIndex+1 : surfpoint, 1:surfpoint+correctionIndex ] );
    end
    
    upperLimit          = ceil( max( max ( ZI ) ) );
    lowerLimit          = upperLimit - noiseGate;
    
    figure( y );
    
    % compute color map
    staticColorMap( background, backgroundWidth, backgroundHeight, XI, YI, ZI, lowerLimit, upperLimit, contourSpan, ...
                    alpha, contourEnable, useFixedSplPlotBounds, fixedLowerLimit, fixedUpperLimit, valueTickEnable, cmap, ...
                    tickSize, tickBold, colorbarEnable );
        
    % print title for each map and save in png format, ordering by position
    saveColorMap( recordedFileName, encodingFilterMatrixFileName, mapsFolder, lowPassCutFreq, ...
                  bandpassLowCutFreq, bandpassHighCutFreq, centerFreq, y, refName, aWeightingEnable )
              
              
end

rmpath( './lib' )
fprintf( 'Ended correctly.\n' )
 