% SpatialIrsMap_ambisonics.m Spatial acoustic static mapping of a NOISE recording
% This script implements acoustic source mapping of a noise recording
% using Ambisonics format of any order got from microphone array.
% It shows Matlab figures and saves png files containing pseudocolor SPL maps.
%
% INPUT
% - ambisonics_format.mat -> file containing the Ambisonics format generated by b_generateEncodedFormat.m
% - panoramic_360_image.jpg
%
% OUTPUT
% - acousticMap.png
%
% Author: Daniel Pinardi
% University of Aalto, Helsinki - Finland, 2019.
% 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
fullScaleLevel          = 130;                                              % dB

% Trimming of the recorded file
trimEnable              = 1;                                                % enable trimming if 1
clickToTrim             = 0;                                                % 1 to enable graphic trim function, else digit number of samples

% Signal processing
plotSignalEnable        = 0;												% show additional plots
aWeightingEnable        = 1;                                                % enable to perform A filtering
nfft                    = 16384;											% resolution of the 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          = 9;                                                % select the last band to analyze


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

maxOrderToUse           = 7; % limit the maximum Ambisonics order to use in the processing


%% IMAGE PARAMETER
resolutionPlot          = 2; % degree

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

% correction of the vertical offset between audio and video
imageOffsetCorrection   = 0;                                                % 1 -> enable
radiusCorrection        = 1.5;                                              % [m]
verticalOffset          = 0.2;                                              % [m] 
outputWindowWidth       = 1024;                                             % desired width of background image

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

% Surf mapping
useAdimensialPlotBounds = 0;                                                % if 1, the normalization is between 0 and 1
adimensionalLowerLimit  = 0;
useFixedSplPlotBounds   = 0;                                                % if 0 each surface plot has color normalization 
fixedUpperLimit         = 72;                                               % user defined upper limit
fixedLowerLimit         = 67;                                               % user defined lower 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


%%  OTHER PARAMETERS
% Frequency range order limitation
o2_lowF                 = 20;		% [Hz]
o2_highF                = 20000;    % [Hz]
o3_lowF                 = 20;		% [Hz]
o3_highF                = 20000;    % [Hz]
o4_lowF                 = 20;		% [Hz]
o4_highF                = 20000; 	% [Hz]
o5_lowF                 = 20;		% [Hz]
o5_highF                = 20000;    % [Hz]
o6_lowF                 = 20;		% [Hz]
o6_highF                = 20000;    % [Hz]
o7_lowF                 = 20;		% [Hz]
o7_highF                = 20000; 	% [Hz]


%% LOAD DATA FROM FILE
% define filename to load and folder output
mapsFolder              = mapsFolderIrs;
fprintf( '\nNoise elaboration started\n' );

% Ambisonics format
if fileAlreadyBformat == 1
    [ ambixFormat, fs ] = audioread( [ './ArrayRecording/' recordedFileName ] );
else
    disp( [ 'Loading Ambix format file: ' encodedIrsMatFileName] );
    load( encodedIrsMatFileName );
    ambixFormat         = encodedFormat;
    clear encodedFormat
end
virtualMicNumber        = size( ambixFormat, 2 );                           % number of virtual microphones
ambixOrder              = sqrt( virtualMicNumber ) - 1;                     % ambisonics order

if maxOrderToUse < ambixOrder
    virtualMicNumber    = ( maxOrderToUse + 1 ) ^ 2;
    ambixFormat         = ambixFormat( :, 1 : virtualMicNumber );
    ambixOrder          = maxOrderToUse;
end

fprintf( [ 'Analysis limited to order ' num2str( ambixOrder ) '.\n\n' ] )

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


%% IMAGE PROCESSING
% 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 A-weighting.\n' );
    AWeighting          = weightingFilter( 'A-weighting', fs );
    for i = 1 : virtualMicNumber
        ambixFormat( :, i ) = AWeighting( ambixFormat( :, i ) );
    end 
    fprintf( 'A-weighting ended.\n' );
end


%% TRIMMING
ambixOmni               = ambixFormat( :, 1 );                              % omni signal
if trimEnable == 1
    fprintf( 'Perform trimming\n\n' );
    if clickToTrim == 1
        ambixFormat     = noiseTrim( ambixFormat, 0 );
    else
        fig             = figure( 'windowstate', 'maximized' );
        hax             = axes;
        plot( ambixOmni )
        xlim( [ 0 size( ambixOmni, 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 );
        ambixFormat     = ambixFormat( 1 + noiseStart : noiseEnd, : );       % trim noise
    end
    
    ambixOmni           = ambixFormat( :, 1 );                              % update value

else % no trimming
    
    fprintf( 'Continuing without trimming.\n' );
    
    if plotSignalEnable == 1
		figure( 'windowstate', 'maximized' )
        plot( ambixOmni );
        xlim( [ 0 size( ambixOmni, 1 ) ] )
        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
    
end
close all
        

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


%% FILTER PROCESSING
freq                    = linspace( 0, fs/2, nfft/2 + 1 );

% FFT
ambixFrf                = avgFft( ambixFormat, nfft, overlap );
ambixFrfComplex         = squeeze( mean( ambixFrf, 3 ) );
complexSum              = real( ambixFrf ).^2 + imag( ambixFrf ).^2;
ambixFrfReal            = mean( complexSum, 3 );

% set frequency limit for orders > 1
[ ~, index_o2_lowF ]    = min( abs( freq - o2_lowF ) );
[ ~, index_o2_highF ]   = min( abs( freq - o2_highF ) );

[ ~, index_o3_lowF ]    = min( abs( freq - o3_lowF ) );
[ ~, index_o3_highF ]   = min( abs( freq - o3_highF ) );

[ ~, index_o4_lowF ]    = min( abs( freq - o4_lowF ) );
[ ~, index_o4_highF ]   = min( abs( freq - o4_highF ) );

[ ~, index_o5_lowF ]    = min( abs( freq - o5_lowF ) );
[ ~, index_o5_highF ]   = min( abs( freq - o5_highF ) );

[ ~, index_o6_lowF ]    = min( abs( freq - o6_lowF ) );
[ ~, index_o6_highF ]   = min( abs( freq - o6_highF ) );

[ ~, index_o7_lowF ]    = min( abs( freq - o7_lowF ) );
[ ~, index_o7_highF ]   = min( abs( freq - o7_highF ) );

if ambixOrder == 2
    ambixFrfComplex( 1 : index_o2_lowF, 5 : 9 )         = 0;
    ambixFrfComplex( index_o2_highF : end, 5 : 9 )      = 0;
elseif ambixOrder == 3
    ambixFrfComplex( 1 : index_o2_lowF, 5 : 9 )         = 0;
    ambixFrfComplex( index_o2_highF : end, 5 : 9 )      = 0;
    ambixFrfComplex( 1 : index_o3_lowF, 10 : 16 )       = 0;
    ambixFrfComplex( index_o3_highF : end, 10 : 16 )    = 0;
elseif ambixOrder == 4
    ambixFrfComplex( 1 : index_o2_lowF, 5 : 9 )         = 0;
    ambixFrfComplex( index_o2_highF : end, 5 : 9 )      = 0;
    ambixFrfComplex( 1 : index_o3_lowF, 10 : 16 )       = 0;
    ambixFrfComplex( index_o3_highF : end, 10 : 16 )    = 0;
    ambixFrfComplex( 1 : index_o4_lowF, 17 : 25 )       = 0;
    ambixFrfComplex( index_o4_highF : end, 17 : 25 )    = 0;
elseif ambixOrder == 5
    ambixFrfComplex( 1 : index_o2_lowF, 5 : 9 )         = 0;
    ambixFrfComplex( index_o2_highF : end, 5 : 9 )      = 0;
    ambixFrfComplex( 1 : index_o3_lowF, 10 : 16 )       = 0;
    ambixFrfComplex( index_o3_highF : end, 10 : 16 )    = 0;
    ambixFrfComplex( 1 : index_o4_lowF, 17 : 25 )       = 0;
    ambixFrfComplex( index_o4_highF : end, 17 : 25 )    = 0;
	ambixFrfComplex( 1 : index_o5_lowF, 26 : 36 )       = 0;
    ambixFrfComplex( index_o5_highF : end, 26 : 36 )    = 0;
elseif ambixOrder == 6
    ambixFrfComplex( 1 : index_o2_lowF, 5 : 9 )         = 0;
    ambixFrfComplex( index_o2_highF : end, 5 : 9 )      = 0;
    ambixFrfComplex( 1 : index_o3_lowF, 10 : 16 )       = 0;
    ambixFrfComplex( index_o3_highF : end, 10 : 16 )    = 0;
    ambixFrfComplex( 1 : index_o4_lowF, 17 : 25 )       = 0;
    ambixFrfComplex( index_o4_highF : end, 17 : 25 )    = 0;
	ambixFrfComplex( 1 : index_o5_lowF, 26 : 36 )       = 0;
    ambixFrfComplex( index_o5_highF : end, 26 : 36 )    = 0;
	ambixFrfComplex( 1 : index_o6_lowF, 37 : 49 )       = 0;
    ambixFrfComplex( index_o6_highF : end, 37 : 49 )    = 0;
elseif ambixOrder == 7
    ambixFrfComplex( 1 : index_o2_lowF, 5 : 9 )         = 0;
    ambixFrfComplex( index_o2_highF : end, 5 : 9 )      = 0;
    ambixFrfComplex( 1 : index_o3_lowF, 10 : 16 )       = 0;
    ambixFrfComplex( index_o3_highF : end, 10 : 16 )    = 0;
    ambixFrfComplex( 1 : index_o4_lowF, 17 : 25 )       = 0;
    ambixFrfComplex( index_o4_highF : end, 17 : 25 )    = 0;
	ambixFrfComplex( 1 : index_o5_lowF, 26 : 36 )       = 0;
    ambixFrfComplex( index_o5_highF : end, 26 : 36 )    = 0;
	ambixFrfComplex( 1 : index_o6_lowF, 37 : 49 )       = 0;
    ambixFrfComplex( index_o6_highF : end, 37 : 49 )    = 0;
	ambixFrfComplex( 1 : index_o7_lowF, 50 : 64 )       = 0;
    ambixFrfComplex( index_o7_highF : end, 50 : 64 )    = 0;
end

% unfiltered signal
sphCovUnfiltered			= cov( ambixFrfComplex );
rmsOmniUnfiltered			= ( sum( ambixFrfReal( :, 1 ), 1 )' ) ./ 2;
[ DoA_unfiltered, ~ ]		= sphPWDmap( sphCovUnfiltered, grid_dirs, 0 );
algorithm					= 'pwd';


if plotSignalEnable == 1
    halfFrf                 = avgFft( ambixOmni, nfft, overlap );
    amplitudeFrfdB          = 20 * log10( mean( abs( halfFrf ), 3 ) );
    freq                    = linspace( 0, fs/2, nfft/2+1 ); 
    k                       = round( 20 / ( fs / nfft ) );
    j                       = round( 20000 / ( fs / nfft ) );
    limitFreq               = [ min( amplitudeFrfdB( k:j ) ) max( amplitudeFrfdB( k:j ) ) ];
    figure( 'windowstate', 'maximized' )
    semilogx( freq, amplitudeFrfdB )
    grid on
    xlabel( 'Frequency [Hz]' )
    title( 'Unfiltered SPL' )
    xlim( [ 20 20000 ] )
    ylim( limitFreq )
    ylabel( '[dB]' )
end

% Low pass filtering
if lowPassCutFreq ~= 0
  
    tmp                 = abs( freq - lowPassCutFreq );
    [ ~, idxLP ]        = min( tmp );                                   % index of closest value
    ambixFrfLpReal      = ambixFrfReal( 1 : idxLP, : );
    ambixFrfLpComplex   = ambixFrfComplex( 1 : idxLP, : );
    sphCovLpFilt        = cov( ambixFrfLpComplex );
    rmsOmniLowPassFilt  = ( sum( ambixFrfLpReal( :, 1 ), 1 )' ) ./ 2;

    [ DoA_LpFilt, ~ ]   = sphPWDmap( sphCovLpFilt, grid_dirs, 0 );

    if plotSignalEnable == 1
        figure( 'windowstate', 'maximized' )
        semilogx( freq, amplitudeFrfdB )
        hold on
        line( [ lowPassCutFreq lowPassCutFreq ], limitFreq, 'Color', 'red' );
        title( 'Low pass' )
        grid on
        xlabel( 'Frequency [Hz]' )
        ylabel( '[dB]' )
        xlim( [ 20 20000 ] )
        ylim( limitFreq )
        hold off
    end
end

% Band pass filtering
if bandpassLowCutFreq ~= 0 && bandpassHighCutFreq ~= 0
    
    tmp                 = abs( freq - bandpassLowCutFreq );
    [~, idxLowBP]       = min( tmp );                                       % index of closest value
    tmp                 = abs( freq - bandpassHighCutFreq );
    [~, idxHighBP]      = min( tmp );                                       % index of closest value
    ambixFrfBpReal      = ambixFrfReal( idxLowBP : idxHighBP, : );
    ambixFrfBpComplex   = ambixFrfComplex( idxLowBP : idxHighBP, : );
    sphCovBpFilt        = cov( ambixFrfBpComplex );
    rmsOmniBandPassFilt = ( sum( ambixFrfBpReal( :, 1 ), 1 )' ) ./ 2;
   
    [ DoA_BpFilt, ~ ]   = sphPWDmap( sphCovBpFilt, grid_dirs, 0 );

    if plotSignalEnable == 1
        figure( 'windowstate', 'maximized' )
        semilogx( freq, amplitudeFrfdB )
        hold on
        line( [ bandpassLowCutFreq bandpassLowCutFreq ], limitFreq, 'Color', 'red' );
        line( [ bandpassHighCutFreq bandpassHighCutFreq ], limitFreq, 'Color', 'red' );
        title( 'Band pass ')
        grid on
        xlabel( 'Frequency [Hz]' )
        ylabel( '[dB]' )
        xlim( [ 20 20000 ] )
        ylim( limitFreq )
        hold off
    end
end

centerFreq                  = 0;
% Octave band filters generation
if firstOctaveBand ~= 0 && lastOctaveBand ~= 0
    bandsPerOctave          = 1;
    octaveBandFilterOrder   = 10;                                               % Filter Order
    centerFreq              = 1000;                                             % Center Frequency (Hz)
    f                       = fdesign.octave( bandsPerOctave, 'Class 1', 'N,F0', octaveBandFilterOrder, centerFreq, fs );
    centerFreq              = validfrequencies( f );                            % calc. approx center frequencies
    centerFreq              = centerFreq( firstOctaveBand : lastOctaveBand );   % selected freq range
    octaveBandsNumber       = length( centerFreq );                             % number of bands selected

    % Octave band filters computation
    if plotSignalEnable == 1
        m                   = ceil( octaveBandsNumber / 2 );
        figure( 'windowstate', 'maximized' )
    end

    DoA_OctaveFilt			= zeros( size( grid_dirs, 1 ), octaveBandsNumber );
    
	rmsOmniOctaveFilt		= zeros( 1, octaveBandsNumber );
    for i = 1 : octaveBandsNumber
        lowCutFreq          = centerFreq( i ) / sqrt( 2 );
        highCutFreq         = lowCutFreq * 2;

        tmp                 = abs( freq - lowCutFreq );
        [~, idxLowOctave]   = min( tmp );                                       % index of closest value
        tmp                 = abs( freq - highCutFreq );
        [~, idxHighOctave]  = min( tmp );                                       % index of closest value

        sphCovOctaveFilt            = cov( ambixFrfComplex( idxLowOctave : idxHighOctave, : ) );
        rmsOmniOctaveFilt( i )      = ( sum( ambixFrfReal( idxLowOctave : idxHighOctave, 1 ), 1 )' ) ./ 2;
		[ DoA_OctaveBandFilt, ~ ]   = sphPWDmap( sphCovOctaveFilt, grid_dirs, 0 );
   
        DoA_OctaveFilt( :, i )  = DoA_OctaveBandFilt;
        
        if plotSignalEnable == 1
            % plot FFT
            subplot( m, 2, i);
            tmp                 = abs( freq - lowCutFreq );
            [ ~, idxLowOctave ] = min( tmp );                                   % index of closest value
            tmp                 = abs( freq - highCutFreq );
            [ ~, idxHighOctave ] = min( tmp );                                  % index of closest value

            semilogx( freq( idxLowOctave : idxHighOctave ), amplitudeFrfdB( idxLowOctave : idxHighOctave ) );	 	 
            hold on
            line( [ centerFreq(i) centerFreq(i) ], limitFreq, 'Color', 'red' );
            title( [ 'Filtering octave band: ' num2str( centerFreq(i) ) ] );
            grid on
            xlabel( 'Frequency [Hz]' )
            ylabel( '[dB]' )
            xlim( [ 20 20000 ] )
            ylim( limitFreq )
            hold off
        end     
    end
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
DoA_unfiltered          = abs( DoA_unfiltered );
maxDoA_unfiltered       = max( DoA_unfiltered );
if useAdimensialPlotBounds == 1
    spl                 = DoA_unfiltered ./ maxDoA_unfiltered;
else
    spl					= DoA_unfiltered ./ maxDoA_unfiltered * rmsOmniUnfiltered;
    spl					= 10 * log10( spl ) + fullScaleLevel;
end

if exist( 'DoA_LpFilt', 'var' )
    DoA_LpFilt          = abs( DoA_LpFilt );
    maxDoA_LpFilt       = max( DoA_LpFilt );
    if useAdimensialPlotBounds == 1
        DoA_LpFilt      = DoA_LpFilt ./ maxDoA_LpFilt;
    else
        DoA_LpFilt		= DoA_LpFilt ./ maxDoA_LpFilt * rmsOmniLowPassFilt;
        DoA_LpFilt		= 10 * log10( DoA_LpFilt ) + fullScaleLevel;
    end
    DoA_LpFilt( isnan( DoA_LpFilt ) ) = 0;
    spl                 = [ spl, DoA_LpFilt ];
end
if exist( 'DoA_BpFilt', 'var' )
    DoA_BpFilt          = abs( DoA_BpFilt );
    maxDoA_BpFilt       = max( DoA_BpFilt );
    if useAdimensialPlotBounds == 1
        DoA_BpFilt      = DoA_BpFilt ./ maxDoA_BpFilt;
    else
        DoA_BpFilt		= DoA_BpFilt ./ maxDoA_BpFilt * rmsOmniBandPassFilt;
        DoA_BpFilt		= 10 * log10( DoA_BpFilt ) + fullScaleLevel;
    end
    DoA_BpFilt( isnan( DoA_BpFilt ) ) = 0;
    spl                 = [ spl, DoA_BpFilt ];
end
if exist( 'DoA_OctaveFilt', 'var' )
    DoA_OctaveFilt      = abs( DoA_OctaveFilt );
    maxDoA_OctaveFilt   = max( DoA_OctaveFilt );
    if useAdimensialPlotBounds == 1
        DoA_OctaveFilt  = DoA_OctaveFilt ./ maxDoA_OctaveFilt;
    else
        DoA_OctaveFilt	= DoA_OctaveFilt ./ maxDoA_OctaveFilt .* rmsOmniOctaveFilt;
        DoA_OctaveFilt	= 10 * log10( DoA_OctaveFilt ) + fullScaleLevel;
    end
    DoA_OctaveFilt( isnan( DoA_OctaveFilt ) ) = 0;
    spl                 = [ spl, DoA_OctaveFilt ];
end


%% GENERATION OF COLORMAPS
correctionIndex         = round( horizontalTilt * ( 360 / resolutionPlot + 1 ) / 360 );
% personalized colorscale
cmap                    = jet( 64 );
cmap( 1, : )            = [ 0 0 0 ]; % black full scale

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

for plot_i = 1 : numberPlot
    
    splPlot             = spl( :, plot_i );                                 % every cycle takes one row
    if useAdimensialPlotBounds == 1
        upperLimit      = 1;
        lowerLimit      = adimensionalLowerLimit;
    else
        upperLimit          = ceil( max( max ( splPlot ) ) );
        lowerLimit          = upperLimit - noiseGate;
        if lowerLimit < 0
            lowerLimit  = 0;
        end
    end
    
    ambisonicsStaticColorMap( background, backgroundWidth, backgroundHeight, contourSpan, lowerLimit, upperLimit, ...
                              alpha, contourEnable, useFixedSplPlotBounds, fixedLowerLimit, fixedUpperLimit, ...
                              valueTickEnable, cmap, splPlot, resolutionPlot, resolutionPlot, tickSize, tickBold, ... 
                              correctionIndex, colorbarEnable );
                     
    % print title for each map and save in png format, ordering by position
    saveColorMap( recordedFileName, encodingFilterMatrixFileName, mapsFolder, lowPassCutFreq, bandpassLowCutFreq, ...
                  bandpassHighCutFreq, centerFreq, plot_i, refName, aWeightingEnable, algorithm )

end

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