percentile_level.m

N%-percentile level of a wavefile using F,S or I temporal constants and A, C, Z weighting networks

function Nptauw = percentile_level(x,p,tau,w,N,Fs)
% This function computes the A- C- or non-weighted percentil levels
% See ISO 1996-1.
%
% To calibrate results just add the sound presure level Lpcal, corresponding to 0 dBFS.
% (this value depends on your electroacoustic system). Consider a 1 kHz signal at full
% scale in Matlab correspond to Lpcal in your acoustic measurement point
% (and the frequency response of your system is plain).
%
% Notice that your measuremnt point could be before recording the analysed wavefile
% or after reproducing it with loudspeakers or headphones.
%
% Usage:
% Nptauw = percentile_level(x,p,N,tauf,w,Fs)
%
% where
% Nptauw is a matrix of p*100 percentile levels where columns are
% audio channels
% x: signal use column for each signal channe or a wave filename
% including path
% p: percentile. Use 10 for N10 or 90 for N90
% tau: time constant for the average (movig average filter)
% 'F': Fast (default); 'S': Slow; 'I': Impulse; or numeric in
% seconds
% w: weighting network
% 'A' (default): A-weighted; 'C': C-weighted; 'Z': non-weighted
% N: # of sampes by frame (2^16 default)*
% Fs: sampling rate (44100 default)*
% *: Be carefull, these parameter affects the precision of A- and C-weighting
% networks
% Note: no provision is given as regards calibration. The user
% should use a calibrated system in order to get meaningful results.
%% Example
% x = [.03*randn(1e6,1).*abs(sin((1:1e6)'/1e6*2*pi)) ...
% .02*randn(1e6,1).*abs(cos((1:1e6)'/1e6*2*pi/5))+.02];
% Lpcal = 103.2;
% NpfasZ = percentile_level(x) + Lpcal;
% p=[1 5 10 50 90 95 99];
% plot(p,NpfasZ);
%
% -------------------------
% Author: Ernesto Accolti
% Date: 2013-09-23
% Revision: 2013-09-26, 2014-05-25
%

...

if ischar(x)
[siz, Fs] = wavread(x,'size');
% if ch, siz(2) = 1; end
else
siz = size(x);
if isempty(Fs),
Fs =44100;
disp('Fs = 44100 is arbitrary assigned');
end
end
samples = siz(1);
channels = siz(2);

...

%% Variables
SC = floor(samples/N); % frame number
b = (-150:.1:0)'; % histogram bins
Nptauw = zeros(length(p),channels); % initializate
zz = zeros(1,channels); % initialize moving average filter
auxfZ=zeros(length(b),channels); % initialize histogram bufer

% 1st order butterworth filter coefficients for moving average
Af = [1, -tau / (1/Fs + tau)]; Bf = [1/Fs / (1/Fs + tau)];
%% Constants
ff = 20*2.^(0:.01:9.97);
C = 12194^2 * ff.^2 ./ (ff.^2 + 20.6^2) ./ (ff.^2 + 12194^2);
A = ff.^2 .* C ./ sqrt(ff.^2 + 107.7^2) ./ sqrt(ff.^2 + 737.9^2);
A = 1.25766693643638 * A;
ff = [0,15,ff, Fs/2]/(Fs/2);
NNff = 2.^12;

%% weighting network design
switch lower(w)
case 'a'
A = [0, 0,A , 0];
wd = fir2(NNff,ff,A);
case 'c'
A = [0, 0,C , 0];
wd = fir2(NNff,ff,A);
case or('z','l')
wd=[];
otherwise
disp('Unknown weighting filter, using A weighting')
A = [0, 0,A , 0];
wd = fir2(NNff,ff,A);
end
zw = zeros(length(wd)-1,channels); % initialize weighting filter
%% Main
for i=1:SC
% read a frame
if ischar(x)
xx=wavread(x,[(i-1)*N+1 (i-1)*N+N]); % Reads each frame i
else
xx = x((i-1)*N+1 : (i-1)*N+N,:); % Reads each frame i
end
% filter weight
if not(isempty(wd))
[xx,zw] = filter(wd, 1, xx,zw); % weighting
end
% Filter squared signal (energy average)
[en,zz] = filter(Bf, Af, xx.^2,zz); % energy moving avernage
% Take histogram
if i*N > Fs*tau,
if (i-1)*N +1 < Fs*tau, ini = ceil(Fs*tau-(i-1)*N);
else ini =1; end
Lp_frame_i = 10*log10(en(ini:end,:));
auxfZ = hist(Lp_frame_i,b) + auxfZ;
end
end

% Last frame
if samples> N*SC;
% read a frame
if ischar(x)
xx=wavread(x,[SC*N+1 samples]); % Reads last frame
else
xx = x(SC*N+1:samples,:); % Reads last frame
end
% filter weighting network
if not(isempty(wd))
[xx,zw] = filter(wd, 1, xx,zw); % energy moving avernage
end
% Filter squared signal (energy average)
en = filter(Bf, Af, xx.^2,zz); % energy moving avernage
% Take histogram
if samples > Fs*tau,
if SC*N +1 < Fs*tau, ini = ceil(Fs*tau);
else ini =1; end
auxfZ = hist(10*log10(en(ini:end,:)),b) + auxfZ;
end
end

auxfZ = cumsum(auxfZ /(samples-ceil(Fs*tau)+1));
%% Extract N-percentiles from histogram
for i=1:length(p)
for j=1:channels
Nptauw(i,j) = b(find((auxfZ(:,j))>=1-p(i)/100,1));
end
end