室内音響関係 / Room Acoustics Related Stuff

残響曲線 / Reverberation Envelope

function y = reverb_envelope(x)
% REVERB_ENVELOPE   Calculates reverberation envelope
%    y=reverb_envelope(x) calculates reverberation envelope y from
%    an impulse response x.  This function can be used to see how
%    the reverberation decays in the corresponding room.
%
%    Reference: Ooga, Yamazaki, and Kaneda "Acoustic Systems and
%               Digital Processing for Them" (p.163)
%
% 2007-05-28 MARUI Atsushi

if ndims(x)>2
  return;
end

if size(x,1) < size(x,2)
  x = x';
end

y = flipud(cumsum(flipud(x .^ 2)));

残響時間 / Reverberation Time

ISO 3382に準拠したものではなく、デモ的に用意しています。

function rt = reverb_time(x, fs, draw_curve, T)
%REVERB_TIME   Compute Reverberation Time (RT60)
%   rt = reverb_time(x, fs);
%   Compute RT60 from T30 (least square fit on reverberation envelope
%   at -5dB and -35dB).
%
%   rt = reverb_time(x, fs, true); draws figure of reverb time
%
%   rt = reverb_time(x, fs, true, [t1 t2 t3]); draws figure of reverb
%   time with least square fit on reverb envelope from t1 dB to t2 dB
%   and extrapolate the line to t3 dB (by default, t1=-5, t2=-35, t3=-60).
%
%   2008-11-17 MARUI Atsushi

if nargin>3
  T1 = T(1);
  T2 = T(2);
  T3 = T(3);
else
  T1 = -5;
  T2 = -35;
  T3 = -60;
end

if nargin>2
  re = reverb_envelope(x, fs, draw_curve);
else
  re = reverb_envelope(x, fs);
end
mm = max(10*log10(abs(re)));
t = ([0:length(re)-1]/fs)';
z = 10*log10(abs(re))-mm;
ind = [sum(z >= T1)  sum(z >= T2)];

% linear fit (-5dB ~ -35dB) and find the time of being -60dB
p = polyfit(z(ind(1):ind(2)), t(ind(1):ind(2)), 1);
rt = polyval(p, T3);

if nargin>2 && draw_curve==true
  p = polyfit(t(ind(1):ind(2)), z(ind(1):ind(2)), 1);
  hold on;
  h = plot(t, polyval(p, t), 'k--');
  set(h, 'Color', [.5 .5 .5]);
  title(sprintf('Reverberation Envelope (RT60=%.3fsec)', rt));
  hold off;
end

C80

ISO 3382にしたがったものではなく、デモ的に用意しています。ISO 3382では、(1)直接音を見つける、(2)50msや80msで分割、(3)それぞれをオクターブ帯域に分割、(4)C80の計算、というステップで行います。下記プログラムはオクターブ帯域に分割していませんし、直接音を見つける部分も若干あやしいです。いずれ直します。

This program does not conform to ISO 3382, but a simpler implementation. The real ISO 3382 C-value computation has following steps: (1) find direct sound, (2) divide impulse response to “early” and “late” at 50 ms, 80 ms, or other time, (3) separate into octave bands, (4) calculate ratio of early and late. The program below does not divide into octave bands and algorithm to find direct sound is not working as what ISO 3382 specifies. I'll correct them when I have a chance.

function C = clarity(x, ms, fs)
%CLARITY   Compute V-value (clarity, definition, ...)
%   clarity(x,ms,fs) compute C-value for specific separation time
%   (eg. 80ms for C80) for audio signal x with sampling frequency fs.
%
%   2010-04-16 MARUI Atsushi

% find the sample number of the direct sound
[z, zind] = max(log10(abs(diff(x(1:fs)))) > -4);   % roughly look for the first peak
z = diff(x(1:fs));
for n=zind:length(x)
  if z(n-1) * z(n) < 0, break;   % find zero crossing of the difference of the impulse response
  end
end
zind = n;   % zero crossing is where the direct sound is
xx = x(zind:end);


te = round(ms/1000 * fs);
numer = sum(xx(1:te) .^ 2);
denom = sum(xx(te+1:end) .^ 2);
C = 10 * log10 (numer / denom);