
% 1/3 octave filter bank


function Db = calc_densities(Xb, Yb, fmb, fs)

    l = min(size(Xb,2), size(Yb,2));

    for n = 1:length(fmb)
        Db(n,1) = Xb(n,[1:l]) * Yb(n,[1:l])' / ( 0.2316 * fmb(n) * 2/fs );    
    end
    
end


% bandwidth calc:
% halfband = (2^(1/NbandsPerOctave))^(1/2)
% full = halfband - 1/halfband 
% = 0.7071 for NbandsPerOctave = 1
% = 0.3483 for NbandsPerOctave = 2
% = 0.2316 for NbandsPerOctave = 3
% = 0.1156 for NbandsPerOctave = 6
% = 0.0693 for NbandsPerOctave = 10
% 1.013 is an overlap factor for filter order Nf = 10


