Get from Ico-github-logo

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video

Highlights from
Impulse response acoustic information calculator

  • functionTemplate(input_args)
    FUNCTIONTEMPLATE Summary of this function goes on this H1 line
  • iosr.acoustics.irStats(fi...
    IRSTATS Calculate RT, DRR, Cte, and EDT for impulse response file
  • iosr.acoustics.rtEst(abs_...
    RTEST Estimate reverberation time based on room size and absorption
  • iosr.auditory.azimuth2itd...
    AZIMUTH2ITD Convert azimuth in degrees to ITD
  • iosr.auditory.binSearch(S...
    Conduct a binary search
  • iosr.auditory.calcIld(L,R...
    CALCILD Calculate normalised interaural level difference
  • iosr.auditory.chXcorr(hc_...
    CHXCORR Calculate cross-correlograms with a wide range of options.
  • iosr.auditory.chXcorr2(hc...
    CHXCORR2 Calculate cross-correlograms with a range of options.
  • iosr.auditory.createWindo...
    Create a Hann or exp. window with specified onsets/offsets
  • iosr.auditory.dupWeight(f)
    DUP_WEIGHT Calculate duplex weighting coefficients for ITD and ILD
  • iosr.auditory.erbRate2hz(x)
    ERBRATE2HZ Convert ERB rate to Hz.
  • iosr.auditory.freqMulti(f)
    FREQMULTI Calculate frequency coefficient for ITD-azimuth warping
  • iosr.auditory.gammatoneFa...
    GAMMATONEFAST Produce an array of responses from gammatone filters via FFT
  • iosr.auditory.hz2erbRate(x)
    HZ2ERBRATE Convert Hz to ERB rate
  • iosr.auditory.instItd(l,r...
    INSTITD Calculate instantaneous ITD
  • iosr.auditory.iso226(phon...
    ISO226 ISO 226:2003 Normal equal-loudness-level contours
  • iosr.auditory.itd2azimuth...
    ITD2AZIMUTH Convert ITD to azimuth
  • iosr.auditory.lindemannIn...
    LINDEMANNINH Signal pre-processing for Lindemann's cross-correlation
  • iosr.auditory.loudWeight(...
    LOUDWEIGHT Calculate loudness weighting coefficients based on ISO 226
  • iosr.auditory.makeErbCFs(...
    MAKEERBCFS Make a series of center frequencies equally spaced in ERB-rate.
  • iosr.auditory.meddisHairC...
    Calculate Ray Meddis' hair cell model for a number of channels.
  • iosr.auditory.perceptualC...
    PERCEPTUALCENTROID Perceptual spectral centroid
  • iosr.auditory.xcorrLindem...
    XCORRLINDEMANN Cross-correlation based on Lindemann's precedence model
  • iosr.bss.applyIdealMasks(...
    APPLYIDEALMASKS Calculate and apply ideal masks via STFT
  • iosr.bss.applyMask(s,m,nf...
    APPLYMASK Apply a time-frequency mask to an STFT
  • iosr.bss.calcImr(m,im)
    CALCIMR Calculates the Ideal Mask Ratio (IMR)
  • iosr.bss.calcSnr(output,t...
    CALCSNR Calculate the separation SNR
  • iosr.bss.cfs2fcs(cfs,fs)
    CFS2FCS Calculate gammatone crossover frequencies.
  • iosr.bss.generateMixtures...
    GENERATEMIXTURES Generate arrays of mixtures from targets and interferers.
  • iosr.bss.getFullMask(m,fr...
    GETFULLMASK Convert frame rate mask to a sample-by-sample mask
  • iosr.bss.idealMasks(st,si...
    IDEALMASKS Calculate ideal time-frequency masks from STFTs
  • iosr.bss.resynthesise(x,f...
    RESYNTHESISE Resynthesise a target from a time-frequency mask
  • iosr.dsp.autocorr(x,Q,dim)
    AUTOCORR Perform autocorrelation via FFT
  • iosr.dsp.convFft(a,b,shape)
    CONVFFT Convolve two vectors using FFT multiplication
  • iosr.dsp.istft(s,nfft,hop...
    ISTFT Calculate the Inverse Short-Time Fourier Transform
  • iosr.dsp.lapwin(L,b)
    LAPWIN Laplace window.
  • iosr.dsp.localpeaks(x,mode)
    LOCALPEAKS Find local peaks and troughs in a vector
  • iosr.dsp.ltas(x,fs,varargin)
    LTAS calculate the long-term average spectrum of a signal
  • iosr.dsp.matchEQ(x,fs,mag...
    MATCHEQ Match the LTAS of a signal to an arbitrary spectral magnitude
  • iosr.dsp.rcoswin(N)
    RCOSWIN Raised cosine window.
  • iosr.dsp.rms(x,dim)
    RMS Calculate the rms of a vector or matrix
  • iosr.dsp.sincFilter(x,Wn,...
    SINCFILTER Apply a near-ideal low-pass or band-pass brickwall filter
  • iosr.dsp.smoothSpectrum(X...
    SMOOTHSPECTRUM Apply 1/N-octave smoothing to a frequency spectrum
  • iosr.dsp.stft(x,nfft,hop,fs)
    STFT Calculate the short-time Fourier transform of a signal
  • iosr.dsp.vsmooth(x,frame,...
    VSMOOTH Smooth a vector using mathematical functions
  • iosr.figures.chMap(M)
    CHMAP Create a monochrome-compatible colour map
  • iosr.figures.cmrMap(M)
    CMRMAP Create a monochrome-compatible colour map
  • iosr.figures.multiwaveplo...
    MULTIWAVEPLOT Stacked line plots from a matrix or vectors
  • iosr.figures.subfigrid(nr...
    SUBFIGRID Create axis positions for subfigures
  • iosr.general.cell2csv(C,f...
    CELL2CSV Output a cell array to a CSV file
  • iosr.general.checkMexComp...
    CHECKMEXCOMPILED Check if mex file is compiled for system
  • iosr.general.getContents(...
    GETCONTENTS Get the contents of a specified directory
  • iosr.general.updateConten...
    UPDATECONTENTS Create a Contents.m file including subdirectories
  • iosr.general.urn(varargin)
    URN Generate random number sequence without duplicates
  • iosr.install
    INSTALL Set search paths, and download and install dependencies.
  • iosr.statistics.getRmse(X...
    GETRMSE Calculate the root-mean-square error between input data
  • iosr.statistics.laprnd(va...
    LAPRND Pseudorandom numbers drawn from the Laplace distribution
  • iosr.statistics.qqPlot(va...
    QQPLOT Quantile-quantile plot with patch option
  • iosr.statistics.quantile(...
    QUANTILE Quantiles of a sample via various methods.
  • iosr.statistics.tab2box(X...
    TAB2BOX Prepare tabular data for boxPlot function
  • iosr.statistics.trirnd(va...
    TRIRND Pseudorandom numbers drawn from the triangular distribution
  • iosr.svn.buildSvnProfile(...
    BUILDSVNPROFILE Read data from files tagged with SVN keywords
  • iosr.svn.headRev(folders,...
    HEADREV Retrieve the head revision for specified files
  • iosr.svn.readSvnKeyword(f...
    READSVNKEYWORD Read data from a file tagged with an SVN keyword
  • boxPlot.m
  • classTemplate
    CLASSTEMPLATE Summary of this class goes on this H1 line
  • iosr.bss.mixture
    MIXTURE Class of sound source separation mixture.
  • iosr.bss.source
    SOURCE Class of sound source separation source.
    AUDIO Abstract superclass providing audio-related properties and methods.
  • iosr.statistics.functiona...
    FUNCTIONALSPREADPLOT Draw a functional plot showing data spread
  • iosr.statistics.statsPlot
    STATSPLOT An abstract superclass for classes that plot statistics
  • example.m
    determine STFT parameters
  • Contents.m
  • View all files
5.0 | 4 ratings Rate this file 44 Downloads (last 30 days) File Size: 639 KB File ID: #42566 Version: 1.5.4
image thumbnail

Impulse response acoustic information calculator



10 Jul 2013 (Updated )

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

| Watch this File

File Information

NOTE: this function is now available from the IoSR Matlab Toolbox as iosr.acoustics.irStats.
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
      '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.


  [1] Zahorik, P., 2002: 'Direct-to-reverberant energy ratio
      sensitivity', The Journal of the Acoustical Society of America,
      112, 2110-2117.



Octave inspired this file.

Required Products Signal Processing Toolbox
MATLAB release MATLAB 8.1 (R2013a)
Other requirements This function requires the Octave toolbox available on the FX:
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (14)
05 Dec 2016 jing zhang

17 Sep 2016 Christopher Hummersone

@Mohsin can you post a link to the impulse response file that causes the error (or send the file to me directly).

Comment only
17 Sep 2016 Mohsin Ahmed

I am having this error:

>> RT = iosr.acoustics.irStats('K:\EmoDB_distance\Release_RT_estimation_MatlabFileExchange\speech_file\TSP1min48.wav', 'graph', false, 'te', 0.05 , 'y_fit', [-5 -25], 'correction', 0.0025, 'spec', 'mean')
In an assignment A(I) = B, the number of elements in B and I must be the same.

Error in iosr.acoustics.irStats (line 191)
t0(n) = find(x(:,n).^2==max(x(:,n).^2)); % find direct impulse

Comment only
16 Sep 2016 Mohsin Ahmed

I tried the following command

RT = irStats('K:\EmoDB_distance\Release_RT_estimation_MatlabFileExchange\speech_file\merged.wav', 'graph', true, 'te', 0.05 , 'spec', 'full', 'y_fit', [0 60], 'correction', 0.0025)

but cannot make it work. it says
'Undefined function 'irStats' for input arguments of type 'char'.' Could you please help me fix this issue?

Comment only
15 Sep 2016 Christopher Hummersone

@Razvan can you post a minimal working example so that I can try to reproduce and fix the bug? Please also post a link to the impulse response file that causes the error (or send the file to me directly).

Comment only
15 Sep 2016 Razvan Todea

I can not find calc_decay function

>> rt=irStats(filename);
Undefined function 'calc_decay' for input arguments of type 'double'.

Error in irStats (line 201)
[rt_temp(f,n),E_rt,fit_rt] = calc_decay(z(f,t0:end,n),options.y_fit,60,fs,cfs(f)); % estimate

Comment only
26 May 2016 Earl Vickers

Thanks for the bug fix! Seems to work well now even with impulse responses having a low dynamic range.

25 May 2016 Christopher Hummersone

@Earl, can you post a minimal working example so that I can try to reproduce and fix the bug?

Comment only
25 May 2016 Earl Vickers

I often get:

Assignment has more non-singleton rhs dimensions than non-singleton subscripts

Error in IR_stats (line 186)
[rt_temp(f,n),E_rt,fit_rt] = calc_decay(z(f,t0:end,n),options.y_fit,60,fs);

20 Mar 2016 Christopher Hummersone

Hi Augusto. Yes you can. Try:

Rt60 = ir_stats(filename,'spec','full');


Comment only
16 Mar 2016 Augusto Carvalho de Sousa

Hi there,

Is it possible with this .m file to determine de RT60 for the separate octave bands? I was trying to do so but I couldn't figure it out. Tried to pass a filtered IR generated with the audiowrite function but it didn't work out.

Comment only
14 Nov 2014 Christopher Hummersone

Hi Clay. I need a bit more information in order to help. Firstly, did you download the Octave toolbox per the requirements above (and add the files to your MATLAB search path)?

Comment only
13 Nov 2014 Clay Pipkin

An error comes up:

Error in IR_stats (line 164)
[b(f,:),a(f,:)] = octdsgn(cfs(f),fs,N);

Any thoughts?

Comment only
27 Jul 2013 Jose Ercolino

15 Jul 2013 1.1

Minor to tweak to handling out-of-range te values.

19 Jul 2013 1.2

Minor tweak to file and updated screenshot.

10 Dec 2013 1.3

Made a number of tweaks: added more comments, added EDT output, and made slight modification to RT calculation so that y range can be specified and is [0 -60] by default.

11 Dec 2013 1.4

Corrected t0 estimation and H1 line.

03 Sep 2014 1.5

Changed inputs to accept parameter/value pairs, including new 'correction' parameter. Changed default correction.

19 Jun 2015 1.5.1

Added error trap for dynamic range in IR.

19 Jun 2015 1.5.2

A number of improvements and clarifications to the help. Some minor refactoring to improve the clarity of code.

19 Jun 2015 1.5.2

Updated description.

25 May 2016 1.5.3

Added more input checking.

26 May 2016 1.5.4

Made decay time calculation more robust to input files with a low dynamic range (thanks to Earl Vickers for the bug report).

21 Feb 2017 1.5.4

Migrated to GitHub.

Contact us