percentile_level.m

N%-percentile level of a wavefile using F,S or I temporal constants and A, C, Z weighting networks
230 Downloads
Updated 28 May 2014

View License

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

Cite As

Ernesto (2024). percentile_level.m (https://www.mathworks.com/matlabcentral/fileexchange/46745-percentile_level-m), MATLAB Central File Exchange. Retrieved .

MATLAB Release Compatibility
Created with R2012a
Compatible with any release
Platform Compatibility
Windows macOS Linux
Categories
Find more on Audio I/O and Waveform Generation in Help Center and MATLAB Answers

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
Version Published Release Notes
1.1.0.0

Updated comments

1.0.0.0