MATLAB Answers

Impulse response acoustic information calculator

16 views (last 30 days)
Bill Tate
Bill Tate on 6 Jul 2020
Answered: Chidvi Modala on 14 Jul 2020
Hey I am having trouble getting this code to work, I can not get the first line to work in irstats. it seems that it can not read my audio file and gives a error of. Also does anyone know what Varargin means ? :
function [rt,drr,cte,cfs,edt] = irStats(filename,varargin)
%IRSTATS Calculate RT, DRR, Cte, and EDT for impulse response file
%
% RT = IOSR.ACOUSTICS.IRSTATS(FILENAME) returns the reverberation time
% (to -60 dB) using a method based on ISO 3382-1:2009. The function uses
% reverse cumulative trapezoidal integration to estimate the decay curve,
% and a linear least-square fit to estimate the slope between 0 dB and
% -60 dB. Estimates are taken in octave bands and the overall figure is
% an average of the 500 Hz and 1 kHz bands.
%
% FILENAME should be the full path to an audio file or the name of an
% audio file on the Matlab search path. The file can be of any format
% supported by the AUDIOREAD function, and have any number of channels;
% estimates (and plots) will be returned for each channel.
%
% The function returns a 1xN vector of RTs, where N is the number of
% channels in the audio file.
%
% The function determines the direct sound as the peak of the squared
% impulse response.
%
% [RT,DRR] = IOSR.ACOUSTICS.IRSTATS(FILENAME) returns the
% direct-to-reverberant-ratio DRR for the impulse; DRR is a 1xN vector.
% This is calculated in the following way:
%
% DRR = 10 * log10( X(T0-C:T0+C)^2 / X(T0+C+1:end)^2 )
%
% where X is the approximated integral of the impulse, T0 is the time of
% the direct impulse, and C=2.5ms [1].
%
% [RT,DRR,CTE] = IOSR.ACOUSTICS.IRSTATS(FILENAME) returns the
% early-to-late index CTE for the impulse; CTE is a 1xN vector.
% This is calculated in the following way:
%
% CTE = 10 * log10( X(T0-C:T0+TE)^2 / X(T0+TE+1:end)^2 )
%
% where TE is 50 ms.
%
% [RT,DRR,CTE,CFS] = IOSR.ACOUSTICS.IRSTATS(FILENAME) returns the
% octave-band centre frequencies CFS used in the calculation of RT.
%
% [RT,DRR,CTE,CFS,EDT] = IOSR.ACOUSTICS.IRSTATS(FILENAME) returns the
% early decay time EDT, which is the same size as RT. The slope of the
% decay curve is determined from the fit between 0 and -10 dB. The decay
% time is calculated from the slope as the time required for a 60 dB
% decay.
%
% ... = IOSR.ACOUSTICS.IRSTATS(...,'PARAMETER',VALUE) allows numerous
% parameters to be specified. These parameters are:
%
% 'graph' : {false} | true
% Controls whether decay curves are plotted. Specifically, graphs
% are plotted of the impulse response, decay curves, and linear
% least-square fit for each octave band and channel of the audio
% file. If the EDT output is specified, the EDT fit will also be
% plotted.
% 'te' : {0.05} | scalar
% Specifies the early time limit (in seconds).
% 'spec' : {'mean'} | 'full'
% Determines the nature of RT and EDT outputs. With spec='mean'
% (default) the reported RT and EDT are the mean of the 500 Hz
% and 1 kHz bands. With spec='full', the function returns the
% RT and EDT as calculated for each octave band returned in
% CFS; RT and EDT have size [M N] where M=length(CFS).
% 'y_fit' : {[0 60]} | two-element vector
% Specifies the decibel range over which the decay curve should
% be evaluated. For example, 'y_fit' may be [-5 -25] or [-5 -35]
% corresponding to the RT20 and RT30 respectively.
% 'correction' : {0.0025} | scalar
% Specifies the correction parameter C (in seconds) given above
% for DRR and CTE calculations. Values of up to 10 ms have been
% suggested in the literature.
%
% Octave-band filters are calculated according to ANSI S1.1-1986 and IEC
% standards. Note that the OCTDSGN function recommends centre frequencies
% fc in the range fs/200 < fc < fs/5.
%
% The author would like to thank Feifei Xiong for his input on the
% correction parameter.
%
% References
%
% [1] Zahorik, P., 2002: 'Direct-to-reverberant energy ratio
% sensitivity', The Journal of the Acoustical Society of America,
% 112, 2110-2117.
%
% See also AUDIOREAD, OCTDSGN.
% Copyright 2016 University of Surrey.
%% validate inputs and set options
% check dependency
if exist('octdsgn','file')~=2
web('http://uk.mathworks.com/matlabcentral/fileexchange/69-octave','-new','-browser')
error('iosr:irStats:OctaveToolbox','Please download and install the OCTAVE toolbox from: http://uk.mathworks.com/matlabcentral/fileexchange/69-octave')
end
% check file exists
assert(exist(filename,'file')==2, 'iosr:irStats:invalidFile', ['iosr.acoustics.irStats: ' filename ' does not exist'])
% set defaults
options = struct(...
'graph',false,...
'te',0.05,...
'spec','mean',...
'y_fit',[0 -60],...
'correction',0.0025);
% read parameter/value inputs
if nargin>1 % if parameters are specified
% read the acceptable names
optionNames = fieldnames(options);
% count arguments
nArgs = length(varargin);
if round(nArgs/2)~=nArgs/2
error('iosr:irStats:nameValuePair','IRSTATS needs propertyName/propertyValue pairs')
end
% overwrite defults
for pair = reshape(varargin,2,[]) % pair is {propName;propValue}
IX = strcmpi(pair{1},optionNames); % find match parameter names
if any(IX)
% do the overwrite
options.(optionNames{IX}) = pair{2};
else
error('iosr:irStats:unknownOption','%s is not a recognized parameter name',pair{1})
end
end
end
% check options size and type
assert(isvector(options.y_fit) && numel(options.y_fit)==2, 'iosr:irStats:invalidYfit', '''y_fit'' must be a two-element vector.')
assert(isscalar(options.correction), 'iosr:irStats:invalidCorrection', '''correction'' must be a scalar.')
assert(isscalar(options.te), 'iosr:irStats:invalidTe', '''te'' must be a scalar.')
assert(ischar(options.spec), 'iosr:irStats:invalidSpec', '''spec'' must be a char array.')
assert(islogical(options.graph) && numel(options.graph)==1, 'iosr:irStats:invalidGraph', '''graph'' must be logical.')
% check for reasonable values
assert(all(options.y_fit<=0), 'iosr:irStats:invalidYfit', '''y_fit'' values must be less than or equal to 0.')
assert(options.te>=0, 'iosr:irStats:invalidTe', '''te'' must be greater than or equal to 0.')
assert(options.correction>=0, 'iosr:irStats:invalidCorrection', '''correction'' must be greater than or equal to 0.')
%% read in audio file
% read in impulse
[x,fs] = audioread(filename);
assert(fs>=5000, 'iosr:irStats:invalidFs', 'Sampling frequency is too low. FS must be at least 5000 Hz.')
assert(sum(abs(x(:)))>0, 'iosr:irStats:silence', 'The signal appears to be silent.')
% set te in samples
te = round(options.te*fs);
% Check sanity of te
assert(te<length(x), 'iosr:irStats:invalidTe', 'The specified early time limit te is longer than the duration of the impulse!')
% get number of channels
numchans = size(x,2);
%% set up octave-band filters
% octave-band center frequencies
cfs = [31.25 62.5 125 250 500 1000 2000 4000 8000 16000];
% octave-band filter order
N = 3;
% limit centre frequencies so filter coefficients are stable
cfs = cfs(cfs>fs/200 & cfs<fs/5);
cfs = cfs(:);
% calculate filter coefficients
a = zeros(length(cfs),(2*N)+1);
b = zeros(length(cfs),(2*N)+1);
for f = 1:length(cfs)
[b(f,:),a(f,:)] = octdsgn(cfs(f),fs,N);
end
%% perform calculations
% empty matrices to fill
z = zeros([length(cfs) size(x)]);
rt_temp = zeros([length(cfs) numchans]);
edt = zeros([length(cfs) numchans]);
t0 = zeros(1,numchans);
drr = zeros(1,numchans);
cte = zeros(1,numchans);
correction = round(options.correction*fs);
% filter and integrate
for n = 1:numchans
% find direct impulse
peak = find(x(:,n).^2==max(x(:,n).^2));
if numel(peak)>1
warning('iosr:irStats:multiplePeaks','More than one peak found. Choosing first peak. Are you sure this is an impulse response?')
end
t0(n) = peak(1);
% draw figure
if options.graph
scrsz = get(0,'ScreenSize');
figpos = [((n-1)/numchans)*scrsz(3) scrsz(4) scrsz(3)/2 scrsz(4)];
figure('Name',['Channel ' num2str(n)],'OuterPosition',figpos);
end
% evaluate impulse in each octave band
for f = 1:length(cfs)
y = filter(b(f,:),a(f,:),x(:,n)); % octave-band filter
temp = cumtrapz(y(end:-1:1).^2); % decay curve
z(f,:,n) = temp(end:-1:1);
[rt_temp(f,n),E_rt,fit_rt] = calc_decay(z(f,t0:end,n),options.y_fit,60,fs,cfs(f)); % estimate RT
[edt(f,n),E_edt,fit_edt] = calc_decay(z(f,t0:end,n),[0,-10],60,fs,cfs(f)); % estimate EDT
if options.graph % plot
% time axes for different vectors
ty = ((0:length(y)-1)-t0(n))./fs;
tE_rt = (0:length(E_rt)-1)./fs;
tE_edt = (0:length(E_edt)-1)./fs;
% plot
subplot(length(cfs),2,(2*f)-1)
plot(ty,y,'k') % octave-band impulse
if f==1
title({'Impulse response'; ''; [num2str(cfs(f)) ' Hz octave band']})
else
title([num2str(cfs(f)) ' Hz octave band'])
end
if f==length(cfs)
xlabel('Time [s]')
else
set(gca,'xticklabel',[]);
end
ylabel('Amplitude')
set(gca,'position',[1 1 1 1.05].*get(gca,'position'),'xlim',[min(ty) max(ty)]);
subplot(length(cfs),2,2*f)
% energy decay and linear least-square fit
if nargout==5
% plot EDT fit if EDT wanted
plot(tE_rt,E_rt,'-k',tE_rt,fit_rt,'--r',tE_edt,fit_edt,':b')
else
plot(tE_rt,E_rt,'-k',tE_rt,fit_rt,'--r')
end
% title for top row
if f==1
title({'Decay curve'; ''; [num2str(cfs(f)) ' Hz octave band']})
else
title([num2str(cfs(f)) ' Hz octave band'])
end
% x label for bottom row
if f==length(cfs)
xlabel('Time [s]')
else
set(gca,'xticklabel',[]);
end
ylabel('Energy [dB]')
set(gca,'position',[1 1 1 1.05].*get(gca,'position'),'ylim',[-70 0],'xlim',[0 max(tE_rt)]);
% choose legend according to EDT request
fitstr = num2str(abs(diff(options.y_fit)));
if nargout==5
legend('Energy decay curve',['Linear fit (RT' fitstr ')'],'Linear fit (EDT)','location','northeast')
else
legend('Energy decay curve',['Linear fit (RT' fitstr ')'],'location','northeast')
end
end
end
% DRR
if nargout>=2
drr(n) = 10.*log10(...
trapz(x(max(1,t0(n)-correction):t0(n)+correction,n).^2)/...
trapz(x(t0(n)+correction+1:end,n).^2)...
);
end
% Cte
if nargout>=3
if t0(n)+te+1>size(x,1)
warning('iosr:irStats:teOutOfRange',['Early time limit (te) out of range in channel ' num2str(n) '. Try lowering te.'])
cte(n) = NaN;
else
cte(n) = 10.*log10(...
trapz(x(max(1,t0(n)-correction):t0(n)+te).^2)/...
trapz(x(t0(n)+te+1:end,n).^2)...
);
end
end
end
%% write output
switch lower(options.spec)
case 'full'
rt = rt_temp;
case 'mean'
rt = mean(rt_temp(cfs==500 | cfs==1000,:)); % overall RT
edt = mean(edt(cfs==500 | cfs==1000,:)); % overall EDT
otherwise
error('iosr:irStats:unknownSpec','Unknown ''spec'': must be ''full'' or ''mean''.')
end
end
function [t,E,fit] = calc_decay(z,y_fit,y_dec,fs,f)
% CALC_DECAY calculate decay time from decay curve
% Returns the time for a specified decay y_dec calculated
% from the fit over the range y_fit. The input is the
% integral of the impulse sample at fs Hz. The function also
% returns the energy decay curve in dB and the corresponding
% fit.
E = 10.*log10(z); % put into dB
E = E-max(E); % normalise to max 0
if any(isinf(E))
E = E(1:find(isinf(E),1,'first')-1); % remove trailing infinite values
end
% find yfit x-range
IX1 = findDB(E,max(y_fit),1,f);
IX2 = findDB(E,min(y_fit),length(E),f);
IX = IX1:IX2;
% calculate fit over yfit
diff_y = abs(diff(y_fit)); % dB range diff
x = reshape(IX,1,length(IX));
y = reshape(E(IX),1,length(IX));
p = polyfit(x,y,1);
fitLength = max(length(E),(1.1*diff_y/abs(length(E)*p(1)))*length(E)); % evaluate fit over sufficient dynamic range
fit = polyval(p,1:fitLength); % actual fit
fit0 = fit-max(fit); % fit anchored to 0dB
t = (y_dec/diff_y)*findDB(fit0,-diff_y,[],f)/fs; % estimate decay time
fit = fit(1:length(E));
end
function IX = findDB(E,dB,default,f)
% FINDDB find dB value in energy decay
IX = find(E<=dB,1,'first');
if isempty(IX)
if isempty(default)
error('iosr:irStats:dynamicRange','Impulse response has insufficient dynamic range at %i Hz to evaluate at %i dB.',f,dB)
else
warning('iosr:irStats:dynamicRange''Impulse response has insufficient dynamic range at %i Hz to evaluate at %i dB. Evaluating at %.1f dB instead.',f,dB,E(default))
IX = default;
end
end
end

  0 Comments

Sign in to comment.

Answers (1)

Chidvi Modala
Chidvi Modala on 14 Jul 2020
varargin is an input parameter given to a function definition statement in order to accept any number of input arguments.
For example:
function XYZ(varargin)
% Function implementation
end
If the above function is called with 2 input arguments, say XYZ(2,3)
then varargin{1} takes the value of 2 and varargin{2} takes the value of 3 inside XYZ function implementation
For more information on varargin, you can refer to this link
Can you provide us what exactly is the error that you are getting while running irStats function?

  0 Comments

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!