Code covered by the BSD License  

Highlights from
Impulse response acoustic information calculator

image thumbnail

Impulse response acoustic information calculator

by

 

10 Jul 2013 (Updated )

Calculate RT, DRR, Cte, and EDT for impulse response file

IR_stats(filename,varargin)
function [rt,drr,cte,cfs,edt] = IR_stats(filename,varargin)
% Calculate RT, DRR, Cte, and EDT for impulse response file
% 
%   RT = IR_STATS(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] = IR_STATS(FILENAME) returns the direct-to-reverberant-ratio
%   DRR for the impulse; DRR is the same size as RT. 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] = IR_STATS(FILENAME) returns the early-to-late index CTE
%   for the impulse; CTE is the same size as RT. 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] = IR_STATS(FILENAME) returns the octave-band centre
%   frequencies CFS used in the calculation of RT.
% 
%   [RT,DRR,CTE,CFS,EDT] = IR_STATS(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.
% 
%   ... = IR_STATS(...,'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.

% =========================================================================
% Last changed:     $Date: 2015-06-19 12:05:24 +0100 (Fri, 19 Jun 2015) $
% Last committed:   $Revision: 387 $
% Last changed by:  $Author: ch0022 $
% =========================================================================

    %% validate inputs and set options

    % check file exists
    assert(exist(filename,'file')==2,['IR_stats: ' 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('IR_STATS 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('%s is not a recognized parameter name',pair{1})
           end
        end
    end
    
    %% read in audio file

    % read in impulse
    [x,fs] = audioread(filename);
    assert(fs>=5000,'Sampling frequency is too low. FS must be at least 5000 Hz.')

    % set te in samples
    te = round(options.te*fs);

    % Check sanity of te
    assert(te<length(x),'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
        t0(n) = find(x(:,n).^2==max(x(:,n).^2)); % find direct impulse
        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
        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); % estimate RT
            [edt(f,n),E_edt,fit_edt] = calc_decay(z(f,t0:end,n),[0,-10],60,fs); % 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.*log(...
                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(['Early time limit (te) out of range in channel ' num2str(n) '. Try lowering te.'])
                cte(n) = NaN;
            else
                cte(n) = 10.*log(...
                    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('Unknown ''spec'': must be ''full'' or ''mean''.')
    end

end

function [t,E,fit] = calc_decay(z,y_fit,y_dec,fs)
% 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
    E = E(1:find(isinf(E),1,'first')-1); % remove trailing infinite values
    IX = find(E<=max(y_fit),1,'first'):find(E<=min(y_fit),1,'first'); % find yfit x-range
    if isempty(IX)
        error('Impulse response has insufficient dynamic range to evaluate to %i dB',min(y_fit))
    end

    % calculate fit over yfit
    x = reshape(IX,1,length(IX));
    y = reshape(E(IX),1,length(IX));
    p = polyfit(x,y,1);
    fit = polyval(p,1:2*length(E)); % actual fit
    fit2 = fit-max(fit); % fit anchored to 0dB

    diff_y = abs(diff(y_fit)); % dB range diff
    t = (y_dec/diff_y)*find(fit2<=-diff_y,1,'first')/fs; % estimate decay time

    fit = fit(1:length(E));

end

Contact us