3

Does someone have a matlab script to get the impulse response from the recorded measurements using log sweep? I have the input file (sweep) and the output and I need to get the impulse response and then calculate reverberation times for octave bands.

MBaz
  • 15,314
  • 9
  • 30
  • 44
rfbrgs
  • 51
  • 4
  • 2
    don't have a script. fft both the input file and the output file (perhaps zero-padding both). divide the fft of the output with the fft of the input. then ifft the result. – robert bristow-johnson Nov 19 '16 at 20:59
  • thanks! but and then how can I get the reverberation times for octave bands? my output file has lots of zero values, and then I got NaN when I use dB...how to solve this? – rfbrgs Nov 19 '16 at 22:09
  • This sounds like a more complicated problem. Could you please rephrase your question so it is applicable to any language? We are trying to be language-agnostic on DSP SE. Nonetheless, please have a look here first. What you have to do is to convolve your recording with time-reversed and scaled signal. When it comes to calculation of reverberation time in octaves, it's a different question I believe. – jojeck Nov 19 '16 at 22:21
  • both division by zero and the logarithm of zero is a problem. usually what i do is add a very small positive number (like $\epsilon = 10^{-12}$) to the positive value in the denominator or the $\log(\cdot)$ argument. – robert bristow-johnson Nov 20 '16 at 02:42
  • @robertbristow-johnson can you turn your comment into an answer here or here with a few math details? Would be super cool, thanks in advance! – Basj Dec 28 '16 at 22:16
  • that's a little work. the comment above is pretty straight forward, if you don't mind FFT-ing an entire sweep and the synchronously recorded response to the sweep. this is actually what they call in audioland a "2-channel FFT" spectrum analyzer. and the other two that i know anything about is the swept-frequency sinusoid and maximum-length sequences. there is specific math for the linear-swept complex sinusoid in continuous time that draws on the fact that either a gaussian pulse or a linear-freq sweep (or combination) is self-similar with the Fourier Transform. – robert bristow-johnson Dec 29 '16 at 00:11
  • this is on research gate somewhere. but i can't find the address. – robert bristow-johnson Dec 29 '16 at 00:13
  • here is the research gate copy. – robert bristow-johnson Dec 29 '16 at 02:01

3 Answers3

1
referenceFilename = 'reference.wav'; 
responseFilename = 'response.wav'; 

[reference, referenceFs] = audioread(referenceFilename);
[nReferenceSamples, nReferenceChannels] = size(reference);

[response, responseFs] = audioread(responseFilename);
[nResponseSamples, nResponseChannels] = size(response);

if referenceFs ~= responseFs
  error('reference audio and response audio files must have same sample rate');
end


nFFT = 2^( ceil( log2(nReferenceSamples+nResponseSamples) ) );

x = zeros(nFFT, 1);
y = zeros(nFFT, 1);

x(1:nReferenceSamples) = reference(:,1);    % use only left channel
y(1:nResponseSamples) = response(:,1);      % use only left channel


h = ifft( fft(y) ./ fft(x) );

h = h(1:nFFT/2);                            % truncate latter half

plot(h);

sound(h, responseFs);
robert bristow-johnson
  • 20,661
  • 4
  • 38
  • 76
0

Once you obtain the IR you must square it, integrate it backwards (Schroeder integral), use octave filters and calculate the slope which is the RT.

J Almagro
  • 1
0

A deconvolution algorithm would work. In the simplest form it is the inverse Fourier transform of the result of dividing the the Fourier transform of the output and the Fourier transform of the input. This may not work if the FT of the input has zeros in it but there are algorithms referenced in the link below which address this issue.

https://en.m.wikipedia.org/wiki/Deconvolution

zoulzubazz
  • 88
  • 4