
% MultiChConvolve.m Generation of IRs from array recorded sweep.
% This script takes sweep recorded files, deconvolve with inverse filters, trims IRs 
% and store the results in a .mat or .wav file 
%
% INPUT
% multi ch sweep recorded .wav file from arrays
% inverse sweep .wav
%
% OUTPUT
% IRs trimmed
%
% Authors: Daniel Pinardi, Marco Binelli.
% University of Parma, 2017
% Last review: 01/2020, Daniel Pinardi

clear variables
close all
clc

%% TRIMMING PARAMETERS FOR IRs
irTargetLength              = 144000; % IR number of samples after cut
irPeakAnticipation          = 10000; % number of samples before IR peak (approximately)


%% SAVE ALSO WAV FILE: POSSIBLE VALUES '0 -> disable' , '1 -> enable'
saveWavEnable               = 1;


%% LOAD DATA FROM FILES
addpath( './lib' )
a_configureFileNamesAndDirectories

% load array sweep recorded
disp( [ 'Loading multi-channels recorded sweep files: ' recordedSweepFile ] );
[ recordedSweep, fs ]       = audioread( recordedSweepFile );
numberChannels              = size( recordedSweep, 2 );  
fprintf( [ 'Number of channels of array recording: ' num2str( numberChannels ) '\n' ] )

% load inverse sweep 
disp( [ 'Loading inverse sweep file: ' invSweepFile ] );
[ inverseSweep, ~ ]         = audioread( invSweepFile );


%% DECONVOLUTION
% Recorded sweep are convolved by inverse sweep giving IRs                       
inverseSweepLength          = length( inverseSweep );
recordedSweepLength         = length( recordedSweep );
irLength                    = recordedSweepLength + inverseSweepLength - 1;
irs                         = zeros( irLength, numberChannels );            % pre-allocate matrix for IRs

% Initialize string to show completed convolution percentage 
fprintf( '\nDeconvolving audio: ' )
percentageString            = '0%%';
percentageStringLenght      = length( percentageString );
fprintf( '%s', percentageString )
oldCompletedPercentage      = -1; 

for channelIndex = 1 : numberChannels
    
    % Updating string to show completed convolution percentage 
    CompletedPercentage     = round( channelIndex / numberChannels * 100 );
    if oldCompletedPercentage ~= CompletedPercentage
        for i = 1:percentageStringLenght
            fprintf( '\b' )
        end
        percentageString    = sprintf( '%d%%', CompletedPercentage );
        percentageStringLenght  = length( percentageString );
        fprintf( '%s', percentageString )
    end
    oldCompletedPercentage  = CompletedPercentage;
    
    % perform deconvolution
    irs( :, channelIndex )  = fd_conv( recordedSweep( :, channelIndex ), inverseSweep );
    
end

fprintf( '\n' )

%% IR trim
% plot before trimming
ok                          = 0;
i                           = 1;
while ok == 0
    
    if i > 1
        prompt              = 'Number of samples to keep: ';
        irTargetLength      = input( prompt );
        prompt              = 'Number of samples before peak: ';
        irPeakAnticipation  = input( prompt );
    end
    
    figure( 'windowstate', 'maximized' )
    subplot( 2, 1, 1 )
    plot( irs( :, 1 ) );
    xlim( [ 0 irLength ] )
    xlabel( 'Samples' )
    grid on
    title( 'Irs before trimming (only channel 1)' );

    [ ~ , peakIndex ]       = max( abs( irs( :, : ) ) );                    % index of peak value of the IRs
    peakIndex               = round( mean( peakIndex ) );                   % mean value of peak index
    iStart                  = peakIndex - irPeakAnticipation;               % start index for IR trim
    iStop                   = iStart + irTargetLength - 1;                  % stop index for IR trim
    irsShort                = zeros( irTargetLength, numberChannels );      % pre-allocate matrix for trimmed IRs

    for i = 1 : numberChannels
        irsShort( :, i )    = irs( iStart : iStop, i );                     % perform trimming
        irsShort( :, i )    = irsShort( :, i ) - mean( irsShort( :, i ) );  % remove DC
    end

    % plot after trimming
    subplot( 2, 1, 2 )
    plot( irsShort );
    xlim( [ 0 irTargetLength ] )
    xlabel( 'Samples' )
    grid on
    title( [ 'Irs after trimming - ' num2str( irTargetLength ) ' samples' ] );

    prompt                  = 'Check trimmed signal. Is ok? (1 -> accept 0 -> redo trimming): ';
    ok                      = input( prompt );    
    
    if ok == 0
        i                   = i + 1;
        close all
    end
    
end

close all
clear irs
irs                         = irsShort;                                     % rename
clear irsShort


%% CONVOLUTION GAIN
if isempty( sweepFileName ) ~= 1
	% load array sweep recorded
	disp( [ 'Loading played sweep: ' sweepFile ] );
	[ playedSweep, fs ]     = audioread( sweepFile );
	delta                   = fd_conv( playedSweep, inverseSweep );
	convGain                = 20 * log10( max( abs( delta ) ) );
	fprintf( [ 'Convolution gain: ' num2str( convGain ) ' dB.\n' ] )
	convFactor              = 10 ^ ( convGain / 20 );
	irs                     = irs ./ convFactor;
else
	convFactor				= max( max( abs( irs ) ) );
	convGain				= 20 * log10( convFactor );
	fprintf( [ 'Convolution gain is: ' num2str( convGain ) ' dB\n' ] )
	prompt					= 'Rescaling IRs to 0 dB ? [1 to accept, 0 to avoid]: ';
	doRescaling				= input( prompt );
	if doRescaling == 1
		irs					= irs ./ convFactor;
	end
end

%% Exporting results in .mat file
if ~isfolder( irsFolder )
    mkdir( irsFolder )
end
disp( [ 'Saving multi-channels IRs MAT files to: ' irsMatFileName ] );
save( irsMatFileName, 'irs', 'fs' );


%% Normalize and save multich .wav file
if saveWavEnable == 1
    maxAbsValue             = max( max( abs( irs ) ) );                     % search for absolute max value
    irs                     = irs./ maxAbsValue;                            % normalize IRs between -1 and +1
    disp( [ 'Saving multi-channels IRs WAV files to: ' irsWavFileName ] );
    audiowrite( irsWavFileName, irs, fs, 'BitsPerSample', 32 );
end

rmpath( './lib' )
fprintf( '\nEnded correctly.\n\n' );
