annotate aux/getLoudness.m @ 4:b28ffd29e6e1

Audio file preparation for listening test

function result = getloudness( vector, Fs, timescale, plotBOOL )
%GETLOUDNESS returns loudness levels according to EBU R 128 specification
% Function getLoudness accepts an audio file vector (vector), a sampling
% rate (Fs), a timescale ('M' - Momentary, 'S' - Short, 'I' - Integrated)
% and a boolean of whether or not to plot the results, and returns either a
% scalar for Integrated timescales, or a vector of loudnesses with a
% refresh rate of 10 Hz.
%
% Calling sequence:
%   result = getLoudness(vector, Fs, timescale, plotBOOL)

% 	% Defining filter coefficients for the pre-filter, according to ITU-R BS.1770-1 
%     % Modelling a spherical head
%     prefilterBcoeffs = [1.53512485958697, -2.69169618940638, 1.19839281085285];
%     prefilterAcoeffs = [1.0, -1.69065929318241, 0.73248077421585];
%     
%     % Defining filter coefficients for the RLB, according to ITU-R BS.1770-1
%     rlbBcoeffs = [1.0, -2.0, 1.0];
%     rlbAcoeffs = [1.0, -1.99004745483398, 0.99007225036621];
    
 % MY 44.1kHZ Coefficients   
    
%     prefilterBcoeffs = [1.53093728623786, -2.65120001232265, 1.16732946528280];
%     prefilterAcoeffs = [1.0, -1.66410930225081, 0.711176041448821];

% by Pedro Pestana; modified by Brecht De Man

    rlbBcoeffs = [[0.994975383507587;], [-1.98995076701517;], [0.994975383507587;]];
    rlbAcoeffs = [1, -1.98992552008493, 0.989976013945421];

    % loudness prefilter function below this one (see Pestana 2013)
    [prefilterBcoeffs, prefilterAcoeffs] = getloudnessprefilter(10, Fs);

   
    % Weighting coefficients with channel format [L, R, C, Ls, Rs]
    nrChannels = size(vector,2);
    switch nrChannels
        case 1,
            G = [1.0];
        case 2,
            G = [1.0 1.0];
        case 5,
            G = [1.0 1.0 1.0 1.41 1.41];
        case 6,
            G = [1.0 1.0 1.0 1.41 1.41]; % double check... AES is L,R,C,LFE,Ls,Rs
        otherwise,
            fprintf('Unsuported number of channels\n');
    end
    
    pfVector = filter(prefilterBcoeffs, prefilterAcoeffs, vector); %Apply pre-filter
	rlbVector = filter(rlbBcoeffs, rlbAcoeffs, pfVector);% Apply RLB filter
    vectorSquared = rlbVector.^2;
    
    switch timescale
        case 'M'
            % MODE: MOMENTARY
            winDur = 400; % window time in miliseconds
            winSize = Fs*winDur/1000; % window duration in samples
            olapSize = Fs*100/1000; % not specified
            maxI = floor( length(vector)/olapSize - winSize/olapSize + 1); % number of iterations
            for i = 1:maxI
                blockEnergy(i,:) = sum( vectorSquared( (i-1)*olapSize+1 : (i-1)*olapSize+winSize ,:) )./winSize;
                time(i) = (i*olapSize)/Fs;
            end
            gTimesZ = bsxfun(@times, blockEnergy, G);
            loudness = -0.691 + 10 .* log10(sum(gTimesZ,2));
            
            if(plotBOOL)
                plot(time, loudness, 'b-', 'LineWidth', 2);
                axis([0 length(vector)/Fs -70 0]);
                xlabel('Time (s)');
                ylabel('Momentary Loudness (dB)')
                grid on;
            end
            
            result = loudness;
            
        case 'S'
            % MODE: SHORT
            winDur = 3000; % window time in miliseconds
            winSize = Fs*winDur/1000; % window duration in samples
            olapSize = Fs*100/1000; % 10 Hz refresh rate as per spec
            maxI = floor( length(vector)/olapSize - winSize/olapSize + 1); % number of iterations
            for i = 1:maxI
                blockEnergy(i,:) = sum( vectorSquared( (i-1)*olapSize+1 : (i-1)*olapSize+winSize ,:) )./winSize;
                time(i) = (i*olapSize)/Fs;
            end
            gTimesZ = bsxfun(@times, blockEnergy, G);
            loudness = -0.691 + 10 .* log10(sum(gTimesZ,2));
            
            if(plotBOOL)
                plot(time, loudness, 'r-.', 'LineWidth', 2);
                axis([0 length(vector)/Fs -70 0]);
                xlabel('Time (s)');
                ylabel('Short Term Loudness (dB)')
                grid on;
            end
            
            result = loudness;
        
        case 'I' 

            % MODE: INTEGRATED
            
            % Gating block interval duration
            winDur = 400; % in miliseconds
            winSize = Fs*winDur/1000; % in samples
            overlap = 0.5; % overlap. Specs: 1/(1-n) must be a natural number exc. 1
            olapSize = overlap*winSize; % block start interval
            maxI = floor( length(vector)/olapSize - winSize/olapSize + 1 ); % number of iterations
    
            for i = 1:maxI
                % equation 3
                blockEnergy(i,:) = sum(vectorSquared( (i-1)*olapSize+1 : (i-1)*olapSize+winSize,:))./winSize;
            end
            
            % equation 4 - loudness
            gTimesZ = bsxfun(@times, blockEnergy, G);  %element-by-elemnt array multiple of G and blockEnergy
            loudness = -0.691 + 10 .* log10(sum(gTimesZ,2));
    
            % applying the absolute gate
            absGate = loudness >= -70;
            absgatedLoudness = loudness(absGate);
            absgatedEnergy = blockEnergy(absGate,:);
            overallAbsLoudness = -0.691 + 10 .* log10( sum((sum(absgatedEnergy)./size(absgatedEnergy,1)).*G) );
    
            % applying the relative gate 8 dB down from overallAbsLoudness
            relGateLevel = overallAbsLoudness - 8;
            relGate = loudness >= relGateLevel;
            relgatedLoudness = loudness(relGate);
            relgateEnergy = blockEnergy(relGate,:);
            overallRelLoudness = -0.691 + 10 .* log10( sum((sum(relgateEnergy)./size(relgateEnergy,1)).*G) );
    
            result = overallRelLoudness;
           
            
        case 'ITU' 

            % MODE: ITU (brecht)
            
            % Gating block interval duration
            winDur = 400; % in miliseconds
            winSize = Fs*winDur/1000; % in samples
            overlap = 0.75; % overlap. Specs: 1/(1-n) must be a natural number exc. 1
            olapSize = overlap*winSize; % block start interval
            maxI = floor( length(vector)/olapSize - winSize/olapSize + 1 ); % number of iterations
    
            for i = 1:maxI
                % equation 3
                blockEnergy(i,:) = sum(vectorSquared( (i-1)*olapSize+1 : (i-1)*olapSize+winSize,:))./winSize;
            end
            
            % equation 4 - loudness
            gTimesZ = bsxfun(@times, blockEnergy, G);  %element-by-elemnt array multiple of G and blockEnergy
            loudness = -0.691 + 10 .* log10(sum(gTimesZ,2));
    
            % applying the absolute gate
            absGate = loudness >= -70;
            absgatedLoudness = loudness(absGate);
            absgatedEnergy = blockEnergy(absGate,:);
            overallAbsLoudness = -0.691 + 10 .* log10( sum((sum(absgatedEnergy)./size(absgatedEnergy,1)).*G) );
    
            % applying the relative gate 10 dB down from overallAbsLoudness
            relGateLevel = overallAbsLoudness - 10;
            relGate = loudness >= relGateLevel;
            relgatedLoudness = loudness(relGate);
            relgateEnergy = blockEnergy(relGate,:);
            overallRelLoudness = -0.691 + 10 .* log10( sum((sum(relgateEnergy)./size(relgateEnergy,1)).*G) );
    
            result = overallRelLoudness;
            
        case 'B' 

            % MODE: INTEGRATED - modified by Brecht De Man to follow ITU-R BS.1770-3
            % and Pestana 2013
            
            % GAIN???
            
            % Gating block interval duration
            winDur = 280; % in miliseconds (see Pestana 2013)
            winSize = Fs*winDur/1000; % in samples
            overlap = 0.75; %! % overlap. Specs: 1/(1-n) must be a natural number exc. 1
            olapSize = overlap*winSize; % block start interval
            maxI = floor( length(vector)/olapSize - winSize/olapSize + 1 ); % number of iterations
    
            for i = 1:maxI
                % equation 3
                blockEnergy(i,:) = sum(vectorSquared( (i-1)*olapSize+1 : (i-1)*olapSize+winSize,:))./winSize;
            end
            
            % equation 4 - loudness
            gTimesZ = bsxfun(@times, blockEnergy, G);  %element-by-elemnt array multiple of G and blockEnergy
            loudness = -0.691 + 10 .* log10(sum(gTimesZ,2));
    
            % applying the absolute gate
            absGate = loudness >= -70;
            absgatedLoudness = loudness(absGate);
            absgatedEnergy = blockEnergy(absGate,:);
            overallAbsLoudness = -0.691 + 10 .* log10( sum((sum(absgatedEnergy)./size(absgatedEnergy,1)).*G) );
    
            % applying the relative gate 10 dB (!) down from overallAbsLoudness
            relGateLevel = overallAbsLoudness - 10; % !
            relGate = loudness >= relGateLevel;
            relgatedLoudness = loudness(relGate);
            relgateEnergy = blockEnergy(relGate,:);
            overallRelLoudness = -0.691 + 10 .* log10( sum((sum(relgateEnergy)./size(relgateEnergy,1)).*G) );
    
            result = overallRelLoudness;
            
        otherwise  
           
            % MODE: My short-time measure
            
            winDur = 3000; % window time in miliseconds
            winSize = Fs*winDur/1000; % window duration in samples
            blockEnergy = sum(vectorSquared(:))./winSize;

            % equation 4 - loudness
            gTimesZ = bsxfun(@times, blockEnergy, G);  %element-by-elemnt array multiple of G and blockEnergy
            loudness = -0.691 + 10 .* log10(sum(gTimesZ,2));
            
            winSize=0;
            olapSize=0;
            
            result = loudness; 
            
    end
end



function [prefilterBcoeffs, prefilterAcoeffs] = getloudnessprefilter(gainIndB, fs)
% GET LOUDNESS calculates coefficients for prefilter with arbitrary gain as 
% specified in EBU R 128 and the modification proposed in Pestana 2013. 
% 
% by Brecht De Man at Centre for Digital Music on 24 April 2014

if nargin < 2
    fs = 44100; % standard 44.1kHz
    if nargin < 1
        gainIndB = 4.0; % as recommended in EBU R 128 / ITU-R BS.1770-3
    end
end

f0      = 1681.9;
%Q = ?
%BW = ?
S       = 1;
 
A       = sqrt(10^(gainIndB/20));
w0      = 2*pi*f0/fs;
alpha   = sin(w0/2) * sqrt( (A+1/A) * (1/S-1) + 2);
 
b0      = A*((A+1) + (A-1)*cos(w0) + 2*sqrt(A)*alpha);
b1      = -2*A*((A-1) + (A+1)*cos(w0));
b2      = A*((A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha);
a0      = (A+1) - (A-1)*cos(w0) + 2*sqrt(A)*alpha;
a1      = 2*((A-1) - (A+1)*cos(w0));
a2      = (A+1) - (A-1)*cos(w0) - 2*sqrt(A)*alpha;
 
% b0 =  (1 - cos(w0))/2;
% b1 =   1 - cos(w0);
% b2 =  (1 - cos(w0))/2;
% a0 =   1 + alpha;
% a1 =  -2*cos(w0);
% a2 =   1 - alpha;
 
% b0      = b0/a0;
% b1      = b1/a0;
% b2      = b2/a0;
% a0      = a0/a0;
% a1      = a1/a0;
% a2      = a2/a0;
 
prefilterBcoeffs    = [b0, b1, b2];
prefilterAcoeffs    = [a0, a1, a2];


% % DEBUG: comparison with ITU/EBU standard hard coded coefficients
% [Hpf,Wpf] = freqz(prefilterBcoeffs, prefilterAcoeffs, 20000, fs); % response in radians per sample
%  
% subplot(2,1,1)
% semilogx(20*log10(abs(Hpf)));
% axis([0 20000 -10 15]);
% grid on
% set(gca, 'XTickMode','manual');
% set(gca, 'XTickLabel',num2str(get(gca,'XTick')'));
%  
% 
% % fs=48k, 4dB gain
% prefilterBcoeffs = [1.53512485958697, -2.69169618940638, 1.19839281085285];
% prefilterAcoeffs = [1.0, -1.69065929318241, 0.73248077421585];
%  
% [Hpf,Wpf] = freqz(prefilterBcoeffs, prefilterAcoeffs, 20000,48000); % response in radians per sample
%  
% subplot(2,1,2)
% semilogx(20*log10(abs(Hpf)));
% axis([0 20000 -10 10]);
% grid on
% set(gca, 'XTickMode','manual');
% set(gca, 'XTickLabel',num2str(get(gca,'XTick')'));


% see 
% S. Mansbridge, S. Finn, and J. D. Reiss, "Implementation and Evaluation 
% of Autonomous Multi-track Fader Control," 132nd Convention of the Audio 
% Engineering Society, 2012. 
% f_c = 1681.97;
% Omega = tan(pi*f_c/fs);
% V_H = 1.58; % probably gain (?+4dB)
% V_L = 1.26; % switch with V_B?
% V_B = 1; 
% Q   = 0.71;
% 
% b(1) = V_L * Omega^2 + V_B*Omega/Q + V_H;
% b(2) = 2*(V_L*Omega^2 - V_H);
% b(3) = V_L*Omega^2 - V_B*Omega/Q + V_H;
% a(1) = Omega^2 + Omega/Q + 1;
% a(2) = 2*(Omega^2 - 1);
% a(3) = Omega^2 - Omega/Q + 1;
% 
% b = b /a(1);
% a = a /a(1);

% f  = 1682; % Hz
% S  = 1; % steep without noticeable notches/peaks
% A  = 10^(gain/40);
% w0 =2*pi*f/fs;
% alpha = sin(w0)/2 * sqrt( (A+ 1/A)*(1/S - 1) + 2 );
% 
% % filter coefficients
% b(1)=    A*( (A+1) + (A-1)*cos(w0) + 2*sqrt(A)*alpha );
% b(2) = -2*A*( (A-1) + (A+1)*cos(w0)                   );
% b(3) =    A*( (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha );
% a(1) =        (A+1) - (A-1)*cos(w0) + 2*sqrt(A)*alpha;
% a(2) =    2*( (A-1) - (A+1)*cos(w0)                   );
% a(3) =        (A+1) - (A-1)*cos(w0) - 2*sqrt(A)*alpha;
% 


% DEBUG
%fvtool(b,a,'Fs',fs);

end