b@0: function lkfs = loudness_itu(x, fs)
b@0: %LOUDNESS_ITU   compute loudness (LKFS) based on ITU-R BS.1770-2
b@0: %   LOUDNESS_ITU(x, fs, mode) compute loudness based on ITU-R BS.1770-2
b@0: %   specification.  x is an input signal (mono/stereo/5.1ch) with sampling
b@0: %   frequency fs.
b@0: %
b@0: %   Input signal has to be either in mono (M), stereo (L and R), or
b@0: %   5.1ch(L, R, C, LFE, Ls, and Rs) in respective channel order.
b@0: %
b@0: %   Example: loudness calculation
b@0: %   [x,fs] = wavread('soundfile.wav');
b@0: %   lkfs = LOUDNESS_ITU(x, fs);
b@0: %
b@0: %   2012-01-06 MARUI Atsushi
b@0: 
b@0: % constants
b@0: blockSize = 400;   % in ms
b@0: overlapSize = 0.75;   % in percentage
b@0: channelWeights = [1 1 1 0 sqrt(2) sqrt(2)];   % in L, R, C, LFE, Ls, Rs order
b@0: absoluteThreshold = -70;   % in dB
b@0: relativeThreshold = -10;   % in dB
b@0: 
b@0: % preparation
b@0: if size(x,1)~=length(x)
b@0:     x = x';
b@0: end
b@0: numch = size(x,2);
b@0: 
b@0: if fs~=48000
b@0:     x = resample(x, 48000, fs);
b@0:     fs = 48000;
b@0: end
b@0: 
b@0: switch(numch)
b@0:     case 5
b@0:         chwat = channelWeights([1 2 3 5 6]);
b@0:     otherwise
b@0:         chwat = channelWeights(1:numch);
b@0: end
b@0: 
b@0: % K-filter
b@0: B1 = [
b@0:     1.53512485958697
b@0:     -2.69169618940638
b@0:     1.19839281085285
b@0:     ];
b@0: A1 = [
b@0:     1.0
b@0:     -1.69065929318241
b@0:     0.73248077421585
b@0:     ];
b@0: B2 = [
b@0:     1.0
b@0:     -2.0
b@0:     1.0
b@0:     ];
b@0: A2 = [
b@0:     1.0
b@0:     -1.99004745483398
b@0:     0.99007225036621
b@0:     ];
b@0: y = filter(B2,A2,filter(B1,A1,x));
b@0: 
b@0: % Mean square
b@0: numBlock = ceil(length(y)/ blockSize) + 1;
b@0: yy = zeros(numBlock * blockSize, numch);
b@0: yy(1:length(y),:) = y;
b@0: j = 0:(length(y) - blockSize)/(blockSize * (1-overlapSize));
b@0: z = zeros(length(j), numch);
b@0: for n1 = 1:length(j)
b@0:     yyy = yy(blockSize*j(n1)*(1-overlapSize)+1:blockSize*(j(n1)*(1-overlapSize)+1)+1,:);
b@0:     z(n1,:) = sum(yyy .^ 2) / blockSize;
b@0: end
b@0: l = -0.691 + 10*log10(sum(repmat(chwat,length(z),1) .* z, 2));
b@0: 
b@0: % Gating (absolute)
b@0: Jg = l > absoluteThreshold;
b@0: zz = repmat(channelWeights(1:numch),length(z),1) .* z;
b@0: L_KG = -0.691 + 10*log10(chwat * sum(zz(Jg,:), 1)' / sum(Jg));
b@0: 
b@0: % Gating (relative)
b@0: Gamma_r = L_KG + relativeThreshold;
b@0: Jg = l > Gamma_r;
b@0: L_KG = -0.691 + 10*log10(chwat * sum(zz(Jg,:), 1)' / sum(Jg));
b@0: 
b@0: % OPTIONAL: draw pretty figure
b@0: if false
b@0:     t = (blockSize*(j*(1-overlapSize)+1)+1)/fs;
b@0:     h = plot(t, l, 'color', [.5 .5 .5]);
b@0:     hold on;
b@0:     h = line([t(1) t(end)], [L_KG L_KG]);
b@0:     set(h, 'Color', 'b');
b@0:     set(h, 'LineStyle', '--');
b@0:     set(h, 'LineWidth', 1);
b@0:     h = line([t(1) t(end)], [Gamma_r Gamma_r]);
b@0:     set(h, 'Color', 'b');
b@0:     set(h, 'LineStyle', ':');
b@0:     set(gca, 'Xlim', [t(1) t(end)]);
b@0:     hold off;
b@0:     legend({'momentary', 'integrated', 'relative gate'});
b@0: end
b@0: 
b@0: % output
b@0: lkfs = L_KG;