% MODELLATE.m - estimate late field parameters.

% (c) Copyright 1999 Kind of Loud Technologies, LLC.  All rights reserved.
%
% Created: 28-Jul-1999.
% Revised: 28-Jul-1999, v1.0, JSA.
% Revised: 13-May-2003, v1.1, JSA - windowed power spectrum, T60 plot added.
% Revised: 12-May-2005, v1.2, JSA - edited for Music 424 demo.
% Version: v1.2.


%% initialization
%%

%% analysis variables
nbins = 2048;	%% analysis window length, samples
nskip = 44;	%% interval between successive analysis frames, samples

fbc = 125*2.^[0:0.5:7]';	%% analysis bin band centers, Hz
abins = length(fbc);		%% analysis bandwidth, bins

rmin = 10^(-30/20);	%% reverberation minimum amplitude, fraction

%% impulse response
ir;		%% impulse response, sample array
fs;		%% sampling rate, Hz
nsamp = length(ir);		%% impulse response length, samples
nframes = ceil((nsamp-nbins)/nskip);	%% impulse response length, frames


%% load impulse response
%%
ir = ir - mean(ir);
ir = ir / max(abs(ir));

figure(1); plot([0:nsamp-1]*(1000/fs), ir); grid;
title(['impulse response']);
xlabel('time - milliseconds'); ylabel('response amplitude');
pause;


%% form psd, smoothed psd
%%
index = kron(ones(1,nframes), [1:nbins]) + kron(nskip*[0:nframes-1], ones(1,nbins));
irperm = reshape(ir(index,1), nbins, nframes);
psd = abs(fft(irperm .* (hanning(nbins)*ones(1,nframes)), 2*nbins)).^2;

wbc = fbc*2/fs;	%% analysis bin band centers, cycles/sample
wbe = [0; sqrt(fbc(1:abins-1).*fbc(2:abins))*2/fs; 1];	%% analysis bin band edges, cycles/sample
psds = zeros(abins,nframes);
for i = [1:abins],
	smooth = ([0:nbins-1]/nbins >= wbe(i)) & ([0:nbins-1]/nbins < wbe(i+1));
	psds(i,:) = smooth * psd(1:nbins,:) / sum(smooth);
end;

%% plot smoothed psd
figure(1); image([0:nframes-1]*nskip*(1000/fs), [0.5 abins-0.5]/abins, 10*log10(psds/max(max(psds)))+64); grid; axis('xy');
title(['response spectra, ', int2str(round(1000*nbins/fs)), '-msec. frames.']);
xlabel('time - msec.'); ylabel('frequency - Hz');
set(gca, 'YTick', [0.5:2:abins-0.5]/abins);
ystring = int2str(fix(fbc(1))); for i = [3:2:abins], ystring = [ystring, '|', int2str(fix(fbc(i)))]; end;
set(gca, 'YTickLabel', ystring);
cbar = colorbar('vert');
set(get(cbar, 'Title'), 'String', 'power - dB');
pause;


%% find q, tau
%% note: ([0:nframes-1]'-(nframes-1)/2) in tempTQ diagonalizes inverse
%%
q = zeros(abins,1);
itau = zeros(abins,1);
for b = [1:abins],
	temp = find(sqrt(psds(b,:)/max(psds(b,:))) > rmin);
	[junk, ndx] = max(psds(b,:));
	index = [ndx+1:max(temp)]';
	basis = [-60*index*nskip/fs, ones(size(index))];
	tempTQ = basis \ (10*log10(psds(b,index))');
	itau(b) = tempTQ(1);		%% tau60 inverse
	q(b) = 10.^(tempTQ(2)/10);	%% equalization filter power

	figure(2);
	plot([0:nframes-1]*nskip*(1000/fs), 10*log10(psds(b,:)), 'og', (index-1)*nskip*1000/fs, basis*tempTQ, '-b'); grid;
	title(['measured and modeled power spectrum, ', int2str(fix(fbc(b))), ' Hz.']);
	xlabel('time - milliseconds'); ylabel('power - dB');
	pause;
end;

%% plot q
figure(2); semilogx((2*1000\fs)*wbc, 10*log10(q), '-r'); grid;
title(['q - resonance spectrum.']);
xlabel('frequency - kHz'); ylabel('power - dB');
set(gca, 'XLim', [0.1 20]);
pause;

%% plot tau
figure(2); semilogx((2*1000\fs)*wbc, 1./itau); grid;
title(['tau - 60-dB decay time.']);
xlabel('frequency - kHz'); ylabel('decay time - seconds');
set(gca, 'XLim', [0.1 20]);