function [FilterCoeffs,SpectrumVec] = MakeIMResp(Fs,beta,X_src,X_rcv,room,cc,LimAtten_dB,measT60,silentflag)
%MakeIMResp  impulse response based on image method calculations
%
% [h,H] = MakeIMRes(Fs,beta,source,sink,room,c,LimAtten,measT60)
% 
% This function generates the impulse response between a sound source and
% an acoustic receiver, based on various environmental parameters such as
% source and sensor positions, enclosure's dimension and reflection
% coefficients, etc. Input parameters are defined as follows:
%
%       Fs: scalar, sampling frequency in hertz. Eg: 8000.
%     beta: vector of dimension 6, corresponding to each wall's reflection
%           coefficient: [x1 x2 y1 y2 z1 z2]. Index 1 indicates wall closest
%           to origin. Set to [0 0 0 0 0 0] to obtain anechoic response (only 
%           direct path). Eg: [0.75 0.75 0.85 0.25 0.3 0.9].
%   source: vector dimension 3, indicating the location of the source in
%           space: [x y z]. Eg: [1 1 1.5].
%     sink: vector of dimension 3, indicating the location of the microphone
%           in space: [x y z]. Eg: [2 2 1.5].
%     room: vector of dimension 3, indicating the rectangular room dimensions:
%           [x_length y_length z_length]. Eg: [4 4 3].
%        c: scalar, speed of acoustic waves. Eg: 342.
% LimAtten: scalar, parameter in dB determining how much of the resulting
%           transfer function is cropped: the impulse response is computed
%           until the time index where its overall energy content has
%           decreased by 'LimAtten' dB, after which the computations stop.
%           Not relevant if beta=zeros(1,6). Eg: 50.
%  measT60: scalar, measured reverberation time in seconds. The above beta
%           parameters typically generate a reverberation time in the TF
%           that is slightly different from the Sabine or Eyring formula
%           using the given beta values. If the T60 value has been
%           practically measured (e.g. using IMRevTimeAnalysis.m) for the
%           current setup, set 'measT60' to this value for a more accurate
%           computation of the TF. Otherwise leave empty, i.e. define as [].
% 
% This function returns the time coefficients of the filter (transfer
% function) in 'h' and its corresponding frequency response 'H'. The filter
% coefficients are real and non-normalised. The first value in vector 'h',
% h(1), corresponds to time t=0. The number of coefficients returned is
% variable and results from the value of 'LimAtten' defined by the user:
% the filter length will be as large as necessary to capture all the
% relevant highest-order reflections. 
% 
% Implementation based on Allen and Berkley's Image Method for Efficiently
% Simulating Small-room Acoustics. See: J.B. Allen and D.A. Berkley, "Image
% method for efficiently simulating small-room acoustics", J. Acoust. Soc.
% Am., issue 65, number 4, April 1979. The computations are carried out in
% the frequency domain, which allows fractional delays for all reflections.

% Release date: September 2006
% Author: Eric A. Lehmann, WATRI, Perth, Australia (www.watri.org.au)
%         Eric.Lehmann@watri.org.au -- Tel. +61 (0)8 6488 4642
%
% Copyright 2006 Eric A. Lehmann
%
% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2 of the License, or (at your
% option) any later version.
% 
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
% Public License for more details.
% 
% You should have received a copy of the GNU General Public License along
% with this program; if not, write to the Free Software Foundation, Inc.,
% 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA


% Explanations for the following code -------------------------------------
% The following implementation of the image method principle has been
% speficically optimised for execution speed. The following code is based
% on the observation that for a particular dimension, the delays from the
% image sources to the receiver increases monotonically as the absolute
% value of the image index (m, n, or l) increases. Hence, all image sources
% for whose indices are above or equal to a specific limit index (for which
% the received delay is above the relevant cut-off value) can be discarded.
% The following code checks, for each dimension, the delay of each received
% path and automatically determines when to stop, thus avoiding unnecessary
% computations.
% This approach is more effective than a similar approach based on the
% amplitude level of the received image source signal. The latter approach
% will typically cut off image source that are below a given amplitude
% threshold, but it is possible that two such sources would still produce a
% combined amplitude larger than the threshold in the resulting TF (this
% might typically leads to significant errors in the tail of the TF). The
% technique used in the code below effectively crops the end of the TF,
% leaving the uncropped part of the TF unaltered (the amount of TF cropped
% depends on the 'LimAtten' parameter).
% The resulting number of considered image sources hence automatically
% results from environmental factors, such as the room dimensions, the
% source and sensor positions, and the walls' reflection coefficients. As a
% result, the length of the computed transfer function has an optimally
% minimum length (no extra padding with negligibly small values).
%--------------------------------------------------------------------------

global SpectrumVec FreqPoints      % not too pretty, but this avoids passing potentially large vectors to frequently called functions...

%---- Check user input:
if X_rcv(1)>=room(1) | X_rcv(2)>=room(2) | X_rcv(3)>=room(3) | X_rcv(1)<=0 | X_rcv(2)<=0 | X_rcv(3)<=0,
    error('Receiver must be within the room boundaries.');
elseif X_src(1)>=room(1) | X_src(2)>=room(2) | X_src(3)>=room(3) | X_src(1)<=0 | X_src(2)<=0 | X_src(3)<=0,
    error('Source must be within the room boundaries.');
elseif ~isempty(find(beta>=1)) | ~isempty(find(beta<0)),
    error('Parameter ''beta'' must be in range [0...1[.');
end

if nargin<9,
    silentflag = 0;     % set to 1 to disable on-screen messages
end

X_src = X_src(:);       % Source location
X_rcv = X_rcv(:);       % Receiver location
beta = beta(:);
Rr = 2*room(:);         % Room dimensions
DPdel = norm(X_rcv - X_src)/cc;         % direct path delay in [s]

%---- Define enough frequency points for resulting time impulse response:
if ~isequal(beta,zeros(6,1)),
    if isempty(measT60),    % if no practical T60 measurement available, use Sabine estimate
        V = prod(room);
        aa_sab = (2-beta(1)^2-beta(2)^2)*room(2)*room(3) + (2-beta(3)^2-beta(4)^2)*room(1)*room(3) + (2-beta(5)^2-beta(6)^2)*room(1)*room(2);
        T60val = 0.161*V/aa_sab;      % Sabine's reverberation time in [s]
    else
        T60val = measT60;   % practical T60 measurement determines real energy decay in TF!
    end
    foo = LimAtten_dB * T60val / 60;    % desired length of TF (TF decays by 60dB for T60 seconds after direct path delay)
    MaxDelay = DPdel + foo;             % maximum delay in TF: direct path plus TForder
else
    MaxDelay = 2*DPdel;         % anechoic case: allow for 2 times direct path in TF 
end
TForder = ceil(MaxDelay*Fs);    % total TF length

FreqPoints = linspace(0,Fs/2,TForder).';
SpectrumVec = zeros(TForder,1);

%---- summation over room dimensions:
if ~silentflag, fprintf('   [MakeIMResp] Computing transfer function '); end;
for a = 0:1
    for b = 0:1
        for d = 0:1
            if ~silentflag, fprintf('.'); end;
            
            m = 1;              % Check delay values for m=1 and above
            FoundLValBelowLim = Check_lDim(a,b,d,m,X_rcv,X_src,Rr,cc,MaxDelay,beta);
            while FoundLValBelowLim==1,
                m = m+1;
                FoundLValBelowLim = Check_lDim(a,b,d,m,X_rcv,X_src,Rr,cc,MaxDelay,beta);
            end
            
            m = 0;              % Check delay values for m=0 and below
            FoundLValBelowLim = Check_lDim(a,b,d,m,X_rcv,X_src,Rr,cc,MaxDelay,beta);
            while FoundLValBelowLim==1,
                m = m-1;
                FoundLValBelowLim = Check_lDim(a,b,d,m,X_rcv,X_src,Rr,cc,MaxDelay,beta);
            end

        end
    end
end
if ~silentflag, fprintf('\n'); end;

%---- Inverse Fourier transform:
SpectrumVec(1) = SpectrumVec(1)/2;      % remove DC component in resulting time coefficients.
freqvec = -i*2*pi*linspace(0,1/2,TForder);
FilterCoeffs = zeros(TForder,1);
for ii=1:TForder,
    freq = exp((ii-1)*freqvec);
    FilterCoeffs(ii) = real(freq*SpectrumVec);
end
FilterCoeffs = FilterCoeffs/TForder;


%============
function [FoundLValBelowLim] = Check_lDim(a,b,d,m,X_rcv,X_src,Rr,cc,MaxDelay,beta)

FoundLValBelowLim = 0;

l = 1;              % Check delay values for l=1 and above
FoundNValBelowLim = Check_nDim(a,b,d,l,m,X_rcv,X_src,Rr,cc,MaxDelay,beta);
while FoundNValBelowLim==1,
    l = l+1;
    FoundNValBelowLim = Check_nDim(a,b,d,l,m,X_rcv,X_src,Rr,cc,MaxDelay,beta);
end
if l~=1, FoundLValBelowLim = 1; end;

l = 0;              % Check delay values for l=0 and below
FoundNValBelowLim = Check_nDim(a,b,d,l,m,X_rcv,X_src,Rr,cc,MaxDelay,beta);
while FoundNValBelowLim==1,
    l = l-1;
    FoundNValBelowLim = Check_nDim(a,b,d,l,m,X_rcv,X_src,Rr,cc,MaxDelay,beta);
end
if l~=0, FoundLValBelowLim = 1; end;


%============
function [FoundNValBelowLim] = Check_nDim(a,b,d,l,m,X_rcv,X_src,Rr,cc,MaxDelay,beta)

global SpectrumVec FreqPoints

FoundNValBelowLim = 0;

n = 1;          % Check delay values for n=1 and above
dist = norm( [2*a-1; 2*b-1; 2*d-1].*X_src + X_rcv - Rr.*[n;l;m] );
foo_time = dist/cc;
while foo_time<=MaxDelay,    % if delay is below TF length limit for n=1, check n=2,3,4...
    foo_amplitude = prod(beta.^abs([n-a; n; l-b; l; m-d; m])) / (4*pi*dist);
    SpectrumVec = SpectrumVec + foo_amplitude * exp(i*2*pi*foo_time*FreqPoints);
    n = n+1;
    dist = norm( [2*a-1; 2*b-1; 2*d-1].*X_src + X_rcv - Rr.*[n;l;m] );
    foo_time = dist/cc;
end
if n~=1, FoundNValBelowLim = 1; end;

n = 0;          % Check delay values for n=0 and below
dist = norm( [2*a-1; 2*b-1; 2*d-1].*X_src + X_rcv - Rr.*[n;l;m] );
foo_time = dist/cc;
while foo_time<=MaxDelay,    % if delay is below TF length for n=0, check n=-1,-2,-3...
    foo_amplitude = prod(beta.^abs([n-a; n; l-b; l; m-d; m])) / (4*pi*dist);
    SpectrumVec = SpectrumVec + foo_amplitude * exp(i*2*pi*foo_time*FreqPoints);
    n = n-1;
    dist = norm( [2*a-1; 2*b-1; 2*d-1].*X_src + X_rcv - Rr.*[n;l;m] );
    foo_time = dist/cc;
end
if n~=0, FoundNValBelowLim = 1; end;