annotate aux/loudness.m @ 15:24be5e9ce25b tip

Update README

function integratedloudness = loudness(signal, fs)
%LOUDNESS returns loudness level (scalar) according to EBU R 128 
% specification and accepts an audio file vector (vector) and a sampling
% rate (fs). 
% 
% by Brecht De Man on Monday 25 April 2016 

    nrSamples  = size(signal,1);
    nrChannels = size(signal,2);

    % loudness prefilters (functions below this one)
    [b1, a1] = loudnessprefilter1(fs);
    [b2, a2] = loudnessprefilter2(fs);
   
    % Weighting coefficients with channel format [L, R, C, Ls, Rs]
    G = [1.0 1.0 1.0 1.41 1.41 0.0 0.0]; % temporary
    
    % Apply pre-filter
    signal_intermediate = filter(b1, a1, signal); 
    signal_filtered = filter(b2, a2, signal_intermediate); 

    % Gating block interval duration
    T_g = .4; % seconds
    Gamma_a = -70.0; % absolute threshold: -70 LKFS
    overlap = 0.75; % relative overlap (0.0 - 1.0)
    step = 1 - overlap;
    
    T = nrSamples/fs; % length of measurement interval in seconds
    j_range = 1:floor((T-T_g)/(T_g*step)+1);
    z = zeros(nrChannels, max(size(j_range)));
    for i = 1:nrChannels % for each channel i
        for j=j_range % for each window j
            lbound = floor(fs*T_g*(j-1)*step+1); 
            hbound = floor(fs*T_g*((j-1)*step+1)); 
            z(i,j) = (1/(T_g*fs))*sum(signal_filtered(lbound:hbound, i).^2);
        end
    end  

    G_current = G(1:nrChannels); % discard weighting coefficients G_i unused channels
    l = zeros(size(j_range,1), 1);
    for j=j_range % for each window j
        l(j) = -.691 + 10.0*log10(sum(G_current*z(:,j)));
    end
        
    % throw out anything below absolute threshold:
    z_avg = mean(z(:,l>Gamma_a),2);
    Gamma_r = -.691 + 10.0*log10(G_current*z_avg) - 10.0;
    % throw out anything below relative threshold:
    z_avg = mean(z(:,l>Gamma_r),2);
    integratedloudness = -.691 + 10.0*log10(G_current*z_avg);
            
end

function [b, a] = loudnessprefilter1(fs)
% LOUDNESSPREFILTER1 calculates coefficients for the first prefilter 
% specified in EBU R 128. 

    % pre-filter 1 coefficients
    f0 = 1681.9744509555319;
    G  = 3.99984385397334;
    Q  = 0.7071752369554193;
    K  = tan(pi * f0 / fs);
    Vh = 10.0^(G / 20.0); % TODO: precompute
    Vb = Vh^0.499666774155; % TODO: precompute
    a0_ = 1.0 + K / Q + K * K;
    b0 = (Vh + Vb * K / Q + K * K) / a0_;
    b1 = 2.0 * (K * K -  Vh) / a0_;
    b2 = (Vh - Vb * K / Q + K * K) / a0_;
    a0 = 1.0;
    a1 = 2.0 * (K * K - 1.0) / a0_;
    a2 = (1.0 - K / Q + K * K) / a0_;
    
    b = [b0, b1, b2];
    a = [a0, a1, a2];
    
    % DEBUG
    %fvtool(b,a,'fs',fs);

end

function [b, a] = loudnessprefilter2(fs)
% LOUDNESSPREFILTER2 calculates coefficients for the second prefilter 
% specified in EBU R 128. 

    % pre-filter 2 coefficients
    f0 = 38.13547087613982;
    Q  = 0.5003270373253953;
    K  = tan(pi * f0 / fs);
    a0 = 1.0;
    a1 = 2.0 * (K * K - 1.0) / (1.0 + K / Q + K * K);
    a2 = (1.0 - K / Q + K * K) / (1.0 + K / Q + K * K);
    b0 = 1.0;
    b1 = -2.0;
    b2 = 1.0;
    
    b = [b0, b1, b2];
    a = [a0, a1, a2];
    
    % DEBUG
    %fvtool(b,a,'fs',fs);

end