b@13: function integratedloudness = loudness(signal, fs)
brecht@14: %LOUDNESS returns loudness level (scalar) according to EBU R 128 
b@13: % specification and accepts an audio file vector (vector) and a sampling
b@13: % rate (fs). 
b@13: % 
b@13: % by Brecht De Man on Monday 25 April 2016 
b@13: 
b@13:     nrSamples  = size(signal,1);
b@13:     nrChannels = size(signal,2);
b@13: 
b@13:     % loudness prefilters (functions below this one)
b@13:     [b1, a1] = loudnessprefilter1(fs);
b@13:     [b2, a2] = loudnessprefilter2(fs);
b@13:    
b@13:     % Weighting coefficients with channel format [L, R, C, Ls, Rs]
b@13:     G = [1.0 1.0 1.0 1.41 1.41 0.0 0.0]; % temporary
b@13:     
b@13:     % Apply pre-filter
b@13:     signal_intermediate = filter(b1, a1, signal); 
b@13:     signal_filtered = filter(b2, a2, signal_intermediate); 
b@13: 
b@13:     % Gating block interval duration
b@13:     T_g = .4; % seconds
b@13:     Gamma_a = -70.0; % absolute threshold: -70 LKFS
b@13:     overlap = 0.75; % relative overlap (0.0 - 1.0)
b@13:     step = 1 - overlap;
b@13:     
b@13:     T = nrSamples/fs; % length of measurement interval in seconds
b@13:     j_range = 1:floor((T-T_g)/(T_g*step)+1);
b@13:     z = zeros(nrChannels, max(size(j_range)));
b@13:     for i = 1:nrChannels % for each channel i
b@13:         for j=j_range % for each window j
b@13:             lbound = floor(fs*T_g*(j-1)*step+1); 
b@13:             hbound = floor(fs*T_g*((j-1)*step+1)); 
b@13:             z(i,j) = (1/(T_g*fs))*sum(signal_filtered(lbound:hbound, i).^2);
b@13:         end
b@13:     end  
b@13: 
b@13:     G_current = G(1:nrChannels); % discard weighting coefficients G_i unused channels
b@13:     l = zeros(size(j_range,1), 1);
b@13:     for j=j_range % for each window j
b@13:         l(j) = -.691 + 10.0*log10(sum(G_current*z(:,j)));
b@13:     end
b@13:         
b@13:     % throw out anything below absolute threshold:
b@13:     z_avg = mean(z(:,l>Gamma_a),2);
b@13:     Gamma_r = -.691 + 10.0*log10(G_current*z_avg) - 10.0;
b@13:     % throw out anything below relative threshold:
b@13:     z_avg = mean(z(:,l>Gamma_r),2);
b@13:     integratedloudness = -.691 + 10.0*log10(G_current*z_avg);
b@13:             
b@13: end
b@13: 
b@13: function [b, a] = loudnessprefilter1(fs)
b@13: % LOUDNESSPREFILTER1 calculates coefficients for the first prefilter 
b@13: % specified in EBU R 128. 
b@13: 
b@13:     % pre-filter 1 coefficients
b@13:     f0 = 1681.9744509555319;
brecht@14:     G  = 3.99984385397334;
b@13:     Q  = 0.7071752369554193;
b@13:     K  = tan(pi * f0 / fs);
b@13:     Vh = 10.0^(G / 20.0); % TODO: precompute
b@13:     Vb = Vh^0.499666774155; % TODO: precompute
b@13:     a0_ = 1.0 + K / Q + K * K;
b@13:     b0 = (Vh + Vb * K / Q + K * K) / a0_;
b@13:     b1 = 2.0 * (K * K -  Vh) / a0_;
b@13:     b2 = (Vh - Vb * K / Q + K * K) / a0_;
b@13:     a0 = 1.0;
b@13:     a1 = 2.0 * (K * K - 1.0) / a0_;
b@13:     a2 = (1.0 - K / Q + K * K) / a0_;
b@13:     
b@13:     b = [b0, b1, b2];
b@13:     a = [a0, a1, a2];
b@13:     
b@13:     % DEBUG
b@13:     %fvtool(b,a,'fs',fs);
b@13: 
b@13: end
b@13: 
b@13: function [b, a] = loudnessprefilter2(fs)
b@13: % LOUDNESSPREFILTER2 calculates coefficients for the second prefilter 
b@13: % specified in EBU R 128. 
b@13: 
b@13:     % pre-filter 2 coefficients
b@13:     f0 = 38.13547087613982;
b@13:     Q  = 0.5003270373253953;
b@13:     K  = tan(pi * f0 / fs);
b@13:     a0 = 1.0;
b@13:     a1 = 2.0 * (K * K - 1.0) / (1.0 + K / Q + K * K);
b@13:     a2 = (1.0 - K / Q + K * K) / (1.0 + K / Q + K * K);
b@13:     b0 = 1.0;
b@13:     b1 = -2.0;
b@13:     b2 = 1.0;
b@13:     
b@13:     b = [b0, b1, b2];
b@13:     a = [a0, a1, a2];
b@13:     
b@13:     % DEBUG
b@13:     %fvtool(b,a,'fs',fs);
b@13: 
b@13: end