%
% FUNCTION 6.5 : "cp0602_thr_dB_vectors"
%
% Evaluates the bandwidth of the input 'PSD' with values in
% the frequency range given by the frequency vector 'f'
% Bandwidth is evaluated according to the given 'threshold'
% (in dB)
%
% The function returns: a truncated frequency vector 'freq_th_dB'
% corresponding to the bandwidth at -'threshold' dB of the input
% PSD, a truncated PSD vector 'PSD_th_dB' containing the 
% corresponding components of the input PSD 'BW', and the
% -'threshold' dB bandwidth
% 'f_high' the higher limit of the bandwidth
% 'f_low' the lower limit of the bandwidth
%
% 
% Programmed by Guerino Giancola / Luca De Nardis

function [f_th_dB, PSD_th_dB] = cp0602_thr_dB_vectors(f,...
   PSD, threshold)

N = length(f);  
df = f(length(f))/length(f);

% ------------------------------------------------
% Step One - Evaluation of the frequency bandwidth
% ------------------------------------------------

[Ppeak,index] = max(PSD);  % Ppeak is the peak value of the
                           % PSD
f_peak = index * df;       % peak frequency   

Pth = Ppeak*10^(threshold/10); % Pth is the value of the
                               % PSD corresponding to the
                               % given threshold

% iterative algorithm for evaluating high and low
% frequencies

imax = index;
P0h = PSD(index);

while (P0h>Pth)&(imax<=N)
    imax = imax + 1;
    P0h = PSD(imax);
end 

f_high = (imax-1) * df;             % high frequency

imin = index;
P0l = PSD(index);

while (P0l>Pth)&(imin>1)&(index>1)
    imin = imin - 1;
    P0l = PSD(imin);
end 

f_low = (min(index,imin)-1) * df;   % low frequency

% end of iterative algorithm

BW = f_high - f_low;           % signal frequency bandwidth

% -----------------------------------------------------
% Step Two - Cutting the vectors in the -'threshold' dB
%            bandwidth
% -----------------------------------------------------

f_th_dB = f(imin:imax);
PSD_th_dB = PSD(imin:imax);