% computeSPLoctavefilt.m This function compute octave band filtered spl from complex sum

function [ splOctaveFilt, centerFreq ] = computeSPLoctavefilt( fs, firstOctaveBand, lastOctaveBand, plotSignalEnable, ...
                             freqTot, complexSum, spsMean, nfft, overlap, limitFreq, fullScaleLevel )
                                
    % Octave band filters generation
    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' );
        halfFrf             = avgFft( spsMean, nfft, overlap );
        amplitudeFrfdB      = 20 * log10( mean( abs( halfFrf ), 3 ) );
        freq                = linspace( 0, fs/2, nfft/2+1 );
    end

    for i = 1 : octaveBandsNumber
        
        lowCutFreq          = centerFreq( i ) / sqrt( 2 );
        highCutFreq         = lowCutFreq * 2;

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

        rmsOctaveFilt( i, : ) = ( sqrt( sum( complexSum( :, idxLowOctave:idxHighOctave), 2 ) ) / sqrt( 2 ) )';

        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
            xlim( [ 20 20000 ] )
            xlabel( 'Frequency [Hz]' )
            ylim( limitFreq )
            ylabel( '[dB]' )
            hold off

        end     

    end

    splOctaveFilt           = fullScaleLevel + 20 * log10( rmsOctaveFilt );     % convert rms to SPL
    
end
