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);

获得 IR 后，您必须对其求平方，将其反向积分（施罗德积分），使用倍频程滤波器并计算斜率即 RT。

反卷积算法会起作用。在最简单的形式中，它是将输出的傅里叶变换和输入的傅里叶变换相除的结果的傅里叶逆变换。如果输入的 FT 中有零，这可能不起作用，但下面的链接中引用了解决此问题的算法。

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