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 50ms, 80ms, 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);