b@4: function result = getloudness( vector, Fs, timescale, plotBOOL )
b@2: %GETLOUDNESS returns loudness levels according to EBU R 128 specification
b@2: % Function getLoudness accepts an audio file vector (vector), a sampling
b@2: % rate (Fs), a timescale ('M' - Momentary, 'S' - Short, 'I' - Integrated)
b@2: % and a boolean of whether or not to plot the results, and returns either a
b@2: % scalar for Integrated timescales, or a vector of loudnesses with a
b@2: % refresh rate of 10 Hz.
b@2: %
b@2: % Calling sequence:
b@2: %   result = getLoudness(vector, Fs, timescale, plotBOOL)
b@2: 
b@2: % 	% Defining filter coefficients for the pre-filter, according to ITU-R BS.1770-1 
b@2: %     % Modelling a spherical head
b@2: %     prefilterBcoeffs = [1.53512485958697, -2.69169618940638, 1.19839281085285];
b@2: %     prefilterAcoeffs = [1.0, -1.69065929318241, 0.73248077421585];
b@2: %     
b@2: %     % Defining filter coefficients for the RLB, according to ITU-R BS.1770-1
b@2: %     rlbBcoeffs = [1.0, -2.0, 1.0];
b@2: %     rlbAcoeffs = [1.0, -1.99004745483398, 0.99007225036621];
b@2:     
b@2:  % MY 44.1kHZ Coefficients   
b@2:     
b@2: %     prefilterBcoeffs = [1.53093728623786, -2.65120001232265, 1.16732946528280];
b@2: %     prefilterAcoeffs = [1.0, -1.66410930225081, 0.711176041448821];
b@2: 
b@2: % by Pedro Pestana; modified by Brecht De Man
b@2: 
b@2:     rlbBcoeffs = [[0.994975383507587;], [-1.98995076701517;], [0.994975383507587;]];
b@2:     rlbAcoeffs = [1, -1.98992552008493, 0.989976013945421];
b@2: 
b@2:     % loudness prefilter function below this one (see Pestana 2013)
b@4:     [prefilterBcoeffs, prefilterAcoeffs] = getloudnessprefilter(10, Fs);
b@2: 
b@2:    
b@2:     % Weighting coefficients with channel format [L, R, C, Ls, Rs]
b@2:     nrChannels = size(vector,2);
b@2:     switch nrChannels
b@2:         case 1,
b@2:             G = [1.0];
b@2:         case 2,
b@2:             G = [1.0 1.0];
b@2:         case 5,
b@2:             G = [1.0 1.0 1.0 1.41 1.41];
b@2:         case 6,
b@2:             G = [1.0 1.0 1.0 1.41 1.41]; % double check... AES is L,R,C,LFE,Ls,Rs
b@2:         otherwise,
b@2:             fprintf('Unsuported number of channels\n');
b@2:     end
b@2:     
b@2:     pfVector = filter(prefilterBcoeffs, prefilterAcoeffs, vector); %Apply pre-filter
b@2: 	rlbVector = filter(rlbBcoeffs, rlbAcoeffs, pfVector);% Apply RLB filter
b@2:     vectorSquared = rlbVector.^2;
b@2:     
b@2:     switch timescale
b@2:         case 'M'
b@2:             % MODE: MOMENTARY
b@2:             winDur = 400; % window time in miliseconds
b@2:             winSize = Fs*winDur/1000; % window duration in samples
b@2:             olapSize = Fs*100/1000; % not specified
b@2:             maxI = floor( length(vector)/olapSize - winSize/olapSize + 1); % number of iterations
b@2:             for i = 1:maxI
b@2:                 blockEnergy(i,:) = sum( vectorSquared( (i-1)*olapSize+1 : (i-1)*olapSize+winSize ,:) )./winSize;
b@2:                 time(i) = (i*olapSize)/Fs;
b@2:             end
b@2:             gTimesZ = bsxfun(@times, blockEnergy, G);
b@2:             loudness = -0.691 + 10 .* log10(sum(gTimesZ,2));
b@2:             
b@2:             if(plotBOOL)
b@2:                 plot(time, loudness, 'b-', 'LineWidth', 2);
b@2:                 axis([0 length(vector)/Fs -70 0]);
b@2:                 xlabel('Time (s)');
b@2:                 ylabel('Momentary Loudness (dB)')
b@2:                 grid on;
b@2:             end
b@2:             
b@2:             result = loudness;
b@2:             
b@2:         case 'S'
b@2:             % MODE: SHORT
b@2:             winDur = 3000; % window time in miliseconds
b@2:             winSize = Fs*winDur/1000; % window duration in samples
b@2:             olapSize = Fs*100/1000; % 10 Hz refresh rate as per spec
b@2:             maxI = floor( length(vector)/olapSize - winSize/olapSize + 1); % number of iterations
b@2:             for i = 1:maxI
b@2:                 blockEnergy(i,:) = sum( vectorSquared( (i-1)*olapSize+1 : (i-1)*olapSize+winSize ,:) )./winSize;
b@2:                 time(i) = (i*olapSize)/Fs;
b@2:             end
b@2:             gTimesZ = bsxfun(@times, blockEnergy, G);
b@2:             loudness = -0.691 + 10 .* log10(sum(gTimesZ,2));
b@2:             
b@2:             if(plotBOOL)
b@2:                 plot(time, loudness, 'r-.', 'LineWidth', 2);
b@2:                 axis([0 length(vector)/Fs -70 0]);
b@2:                 xlabel('Time (s)');
b@2:                 ylabel('Short Term Loudness (dB)')
b@2:                 grid on;
b@2:             end
b@2:             
b@2:             result = loudness;
b@2:         
b@2:         case 'I' 
b@2: 
b@2:             % MODE: INTEGRATED
b@2:             
b@2:             % Gating block interval duration
b@2:             winDur = 400; % in miliseconds
b@2:             winSize = Fs*winDur/1000; % in samples
b@2:             overlap = 0.5; % overlap. Specs: 1/(1-n) must be a natural number exc. 1
b@2:             olapSize = overlap*winSize; % block start interval
b@2:             maxI = floor( length(vector)/olapSize - winSize/olapSize + 1 ); % number of iterations
b@2:     
b@2:             for i = 1:maxI
b@2:                 % equation 3
b@2:                 blockEnergy(i,:) = sum(vectorSquared( (i-1)*olapSize+1 : (i-1)*olapSize+winSize,:))./winSize;
b@2:             end
b@2:             
b@2:             % equation 4 - loudness
b@2:             gTimesZ = bsxfun(@times, blockEnergy, G);  %element-by-elemnt array multiple of G and blockEnergy
b@2:             loudness = -0.691 + 10 .* log10(sum(gTimesZ,2));
b@2:     
b@2:             % applying the absolute gate
b@2:             absGate = loudness >= -70;
b@2:             absgatedLoudness = loudness(absGate);
b@2:             absgatedEnergy = blockEnergy(absGate,:);
b@2:             overallAbsLoudness = -0.691 + 10 .* log10( sum((sum(absgatedEnergy)./size(absgatedEnergy,1)).*G) );
b@2:     
b@2:             % applying the relative gate 8 dB down from overallAbsLoudness
b@2:             relGateLevel = overallAbsLoudness - 8;
b@2:             relGate = loudness >= relGateLevel;
b@2:             relgatedLoudness = loudness(relGate);
b@2:             relgateEnergy = blockEnergy(relGate,:);
b@2:             overallRelLoudness = -0.691 + 10 .* log10( sum((sum(relgateEnergy)./size(relgateEnergy,1)).*G) );
b@2:     
b@2:             result = overallRelLoudness;
b@2:            
b@2:             
b@2:         case 'ITU' 
b@2: 
b@2:             % MODE: ITU (brecht)
b@2:             
b@2:             % Gating block interval duration
b@2:             winDur = 400; % in miliseconds
b@2:             winSize = Fs*winDur/1000; % in samples
b@2:             overlap = 0.75; % overlap. Specs: 1/(1-n) must be a natural number exc. 1
b@2:             olapSize = overlap*winSize; % block start interval
b@2:             maxI = floor( length(vector)/olapSize - winSize/olapSize + 1 ); % number of iterations
b@2:     
b@2:             for i = 1:maxI
b@2:                 % equation 3
b@2:                 blockEnergy(i,:) = sum(vectorSquared( (i-1)*olapSize+1 : (i-1)*olapSize+winSize,:))./winSize;
b@2:             end
b@2:             
b@2:             % equation 4 - loudness
b@2:             gTimesZ = bsxfun(@times, blockEnergy, G);  %element-by-elemnt array multiple of G and blockEnergy
b@2:             loudness = -0.691 + 10 .* log10(sum(gTimesZ,2));
b@2:     
b@2:             % applying the absolute gate
b@2:             absGate = loudness >= -70;
b@2:             absgatedLoudness = loudness(absGate);
b@2:             absgatedEnergy = blockEnergy(absGate,:);
b@2:             overallAbsLoudness = -0.691 + 10 .* log10( sum((sum(absgatedEnergy)./size(absgatedEnergy,1)).*G) );
b@2:     
b@2:             % applying the relative gate 10 dB down from overallAbsLoudness
b@2:             relGateLevel = overallAbsLoudness - 10;
b@2:             relGate = loudness >= relGateLevel;
b@2:             relgatedLoudness = loudness(relGate);
b@2:             relgateEnergy = blockEnergy(relGate,:);
b@2:             overallRelLoudness = -0.691 + 10 .* log10( sum((sum(relgateEnergy)./size(relgateEnergy,1)).*G) );
b@2:     
b@2:             result = overallRelLoudness;
b@2:             
b@2:         case 'B' 
b@2: 
b@2:             % MODE: INTEGRATED - modified by Brecht De Man to follow ITU-R BS.1770-3
b@2:             % and Pestana 2013
b@2:             
b@2:             % GAIN???
b@2:             
b@2:             % Gating block interval duration
b@2:             winDur = 280; % in miliseconds (see Pestana 2013)
b@2:             winSize = Fs*winDur/1000; % in samples
b@2:             overlap = 0.75; %! % overlap. Specs: 1/(1-n) must be a natural number exc. 1
b@2:             olapSize = overlap*winSize; % block start interval
b@2:             maxI = floor( length(vector)/olapSize - winSize/olapSize + 1 ); % number of iterations
b@2:     
b@2:             for i = 1:maxI
b@2:                 % equation 3
b@2:                 blockEnergy(i,:) = sum(vectorSquared( (i-1)*olapSize+1 : (i-1)*olapSize+winSize,:))./winSize;
b@2:             end
b@2:             
b@2:             % equation 4 - loudness
b@2:             gTimesZ = bsxfun(@times, blockEnergy, G);  %element-by-elemnt array multiple of G and blockEnergy
b@2:             loudness = -0.691 + 10 .* log10(sum(gTimesZ,2));
b@2:     
b@2:             % applying the absolute gate
b@2:             absGate = loudness >= -70;
b@2:             absgatedLoudness = loudness(absGate);
b@2:             absgatedEnergy = blockEnergy(absGate,:);
b@2:             overallAbsLoudness = -0.691 + 10 .* log10( sum((sum(absgatedEnergy)./size(absgatedEnergy,1)).*G) );
b@2:     
b@2:             % applying the relative gate 10 dB (!) down from overallAbsLoudness
b@2:             relGateLevel = overallAbsLoudness - 10; % !
b@2:             relGate = loudness >= relGateLevel;
b@2:             relgatedLoudness = loudness(relGate);
b@2:             relgateEnergy = blockEnergy(relGate,:);
b@2:             overallRelLoudness = -0.691 + 10 .* log10( sum((sum(relgateEnergy)./size(relgateEnergy,1)).*G) );
b@2:     
b@2:             result = overallRelLoudness;
b@2:             
b@2:         otherwise  
b@2:            
b@2:             % MODE: My short-time measure
b@2:             
b@2:             winDur = 3000; % window time in miliseconds
b@2:             winSize = Fs*winDur/1000; % window duration in samples
b@2:             blockEnergy = sum(vectorSquared(:))./winSize;
b@2: 
b@2:             % equation 4 - loudness
b@2:             gTimesZ = bsxfun(@times, blockEnergy, G);  %element-by-elemnt array multiple of G and blockEnergy
b@2:             loudness = -0.691 + 10 .* log10(sum(gTimesZ,2));
b@2:             
b@2:             winSize=0;
b@2:             olapSize=0;
b@2:             
b@2:             result = loudness; 
b@2:             
b@2:     end
b@2: end
b@2: 
b@2: 
b@2: 
b@4: function [prefilterBcoeffs, prefilterAcoeffs] = getloudnessprefilter(gainIndB, fs)
b@2: % GET LOUDNESS calculates coefficients for prefilter with arbitrary gain as 
b@2: % specified in EBU R 128 and the modification proposed in Pestana 2013. 
b@2: % 
b@2: % by Brecht De Man at Centre for Digital Music on 24 April 2014
b@2: 
b@2: if nargin < 2
b@2:     fs = 44100; % standard 44.1kHz
b@2:     if nargin < 1
b@2:         gainIndB = 4.0; % as recommended in EBU R 128 / ITU-R BS.1770-3
b@2:     end
b@2: end
b@2: 
b@2: f0      = 1681.9;
b@2: %Q = ?
b@2: %BW = ?
b@2: S       = 1;
b@2:  
b@2: A       = sqrt(10^(gainIndB/20));
b@2: w0      = 2*pi*f0/fs;
b@2: alpha   = sin(w0/2) * sqrt( (A+1/A) * (1/S-1) + 2);
b@2:  
b@2: b0      = A*((A+1) + (A-1)*cos(w0) + 2*sqrt(A)*alpha);
b@2: b1      = -2*A*((A-1) + (A+1)*cos(w0));
b@2: b2      = A*((A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha);
b@2: a0      = (A+1) - (A-1)*cos(w0) + 2*sqrt(A)*alpha;
b@2: a1      = 2*((A-1) - (A+1)*cos(w0));
b@2: a2      = (A+1) - (A-1)*cos(w0) - 2*sqrt(A)*alpha;
b@2:  
b@2: % b0 =  (1 - cos(w0))/2;
b@2: % b1 =   1 - cos(w0);
b@2: % b2 =  (1 - cos(w0))/2;
b@2: % a0 =   1 + alpha;
b@2: % a1 =  -2*cos(w0);
b@2: % a2 =   1 - alpha;
b@2:  
b@2: % b0      = b0/a0;
b@2: % b1      = b1/a0;
b@2: % b2      = b2/a0;
b@2: % a0      = a0/a0;
b@2: % a1      = a1/a0;
b@2: % a2      = a2/a0;
b@2:  
b@2: prefilterBcoeffs    = [b0, b1, b2];
b@2: prefilterAcoeffs    = [a0, a1, a2];
b@2: 
b@2: 
b@2: % % DEBUG: comparison with ITU/EBU standard hard coded coefficients
b@2: % [Hpf,Wpf] = freqz(prefilterBcoeffs, prefilterAcoeffs, 20000, fs); % response in radians per sample
b@2: %  
b@2: % subplot(2,1,1)
b@2: % semilogx(20*log10(abs(Hpf)));
b@2: % axis([0 20000 -10 15]);
b@2: % grid on
b@2: % set(gca, 'XTickMode','manual');
b@2: % set(gca, 'XTickLabel',num2str(get(gca,'XTick')'));
b@2: %  
b@2: % 
b@2: % % fs=48k, 4dB gain
b@2: % prefilterBcoeffs = [1.53512485958697, -2.69169618940638, 1.19839281085285];
b@2: % prefilterAcoeffs = [1.0, -1.69065929318241, 0.73248077421585];
b@2: %  
b@2: % [Hpf,Wpf] = freqz(prefilterBcoeffs, prefilterAcoeffs, 20000,48000); % response in radians per sample
b@2: %  
b@2: % subplot(2,1,2)
b@2: % semilogx(20*log10(abs(Hpf)));
b@2: % axis([0 20000 -10 10]);
b@2: % grid on
b@2: % set(gca, 'XTickMode','manual');
b@2: % set(gca, 'XTickLabel',num2str(get(gca,'XTick')'));
b@2: 
b@2: 
b@2: % see 
b@2: % S. Mansbridge, S. Finn, and J. D. Reiss, "Implementation and Evaluation 
b@2: % of Autonomous Multi-track Fader Control," 132nd Convention of the Audio 
b@2: % Engineering Society, 2012. 
b@2: % f_c = 1681.97;
b@2: % Omega = tan(pi*f_c/fs);
b@2: % V_H = 1.58; % probably gain (?+4dB)
b@2: % V_L = 1.26; % switch with V_B?
b@2: % V_B = 1; 
b@2: % Q   = 0.71;
b@2: % 
b@2: % b(1) = V_L * Omega^2 + V_B*Omega/Q + V_H;
b@2: % b(2) = 2*(V_L*Omega^2 - V_H);
b@2: % b(3) = V_L*Omega^2 - V_B*Omega/Q + V_H;
b@2: % a(1) = Omega^2 + Omega/Q + 1;
b@2: % a(2) = 2*(Omega^2 - 1);
b@2: % a(3) = Omega^2 - Omega/Q + 1;
b@2: % 
b@2: % b = b /a(1);
b@2: % a = a /a(1);
b@2: 
b@2: % f  = 1682; % Hz
b@2: % S  = 1; % steep without noticeable notches/peaks
b@2: % A  = 10^(gain/40);
b@2: % w0 =2*pi*f/fs;
b@2: % alpha = sin(w0)/2 * sqrt( (A+ 1/A)*(1/S - 1) + 2 );
b@2: % 
b@2: % % filter coefficients
b@2: % b(1)=    A*( (A+1) + (A-1)*cos(w0) + 2*sqrt(A)*alpha );
b@2: % b(2) = -2*A*( (A-1) + (A+1)*cos(w0)                   );
b@2: % b(3) =    A*( (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha );
b@2: % a(1) =        (A+1) - (A-1)*cos(w0) + 2*sqrt(A)*alpha;
b@2: % a(2) =    2*( (A-1) - (A+1)*cos(w0)                   );
b@2: % a(3) =        (A+1) - (A-1)*cos(w0) - 2*sqrt(A)*alpha;
b@2: % 
b@2: 
b@2: 
b@2: % DEBUG
b@2: %fvtool(b,a,'Fs',fs);
b@2: 
b@2: end