function [h] = FxRoomImpulse(fs, mic, n, r, rm, src);
%RIR   Room Impulse Response.
%   [h] = RoomImpulse(FS, MIC, N, R, RM, SRC) performs a simple room impulse
%         response calculation.
%
%      FS =  sample rate.
%      MIC = row vector giving the x,y,z coordinates of
%            the microphone.  
%      N =   The program will account for (2*N+1)^3 virtual sources 
%      R =   reflection coefficient for the walls, in general -1<R<1.
%      RM =  row vector giving the dimensions of the room.  
%      SRC = row vector giving the x,y,z coordinates of 
%            the sound source.
%
%   EXAMPLE:
%
%      >>fs=44100;
%      >>rm=[20 19 21];
%      >>mic=[19 18 1.6];
%      >>src=[5 2 1];
%      >>r=0.91;
%      >>n=12;
%      >>h=RoomImpulse(fs, mic, n, r, rm, src);
%
%   NOTES:
%
%   1) To implement this filter, you will need to do a fast 
%      convolution.  The program FCONV.m will do this. 
%   2) All distances are in meters.
%      
%
%  Author: Andreas Floros
%  Digital Audio Processing Applications Group
%  Dept. of Audiovisual Arts, Ionian University
%  v.1.0 - Last Modified: 12/05/2009.

nn=[-n:1:n];                          % Index for the sequence
rms= nn+0.5-0.5*(-1).^nn;             
srcs=(-1).^(nn);                      
xi=[srcs*src(1)+rms*rm(1)-mic(1)];    
yj=[srcs*src(2)+rms*rm(2)-mic(2)];     
zk=[srcs*src(3)+rms*rm(3)-mic(3)];     

[i,j,k]=meshgrid(xi,yj,zk);           % convert vectors to 3D matrices
d=sqrt(i.^2+j.^2+k.^2);               
time=d./343;                          
time=round(time*fs);                  

b=1./(4*pi*((d).^(2)));               
[e,f,g] = meshgrid(nn, nn, nn);       % convert vectors to 3D matrices
c=r.^(abs(e)+abs(f)+abs(g));          
e= b.*c;                              

h=full(sparse(time(:),1,e(:)));       
h=h/max(abs(h));                      % Normalize h

