Image Source Method RIR and Energy Decay Time

Hi, I'm trying to implement an room acoustics software working with image source method for my thesis work. I managed to finish up to some level but I'm stuck with the energy decay curve. Upto now I created image sources for a rectangular room by traditional plane formulas and distances. I'm 100% sure about the locations and distances of the images and I calculated the pressure of the each image source. Then I determined the unique distances of the image sources and got the rms values of the so called shells. When I calculated the energy decay curve and calculated the -60dB drop from the first reflection the value I got is much lower then the value calculated by the Sabines formula. The part of the code is as below. I do not know what I am doing wrong. I will be glad if someone can show me the way or point what I'm doing wrong. Any help will be appreciated.

Initial variables dist:distances of the image sources pressure:pressures of the image sources calculated with the formula sqrt(1-alpha)^i/(4*pi*dist) where i is the order of the reflection and alpha is the absorption coefficient of the surfaces.

time=(dist)./340; t=unique(time);

pres=[];

   p=[];  
   d=0;
   e=1;
   while e<=length(length(pval(:,1)))
      for i=1:length(t)
         a=find(time==t(i));
         rms=sqrt(sum(pressure(a).^2));
         p=[p;rms];
         c=length(a);
         d=c+d;
         e=d+1;
      end
   end

After that code I'm assuming that I got t vs p graph which is also called echogram or impulse response diagram without the direct sound and without convolution with a signal.

Then I tried to calculate the energy decay curve described by Schroeders backwards integral as:

P=flipud(p); Pcum=cumsum(P.^2,1); %cumulative backwards summation of the p^2 P=flipud(Pcum);

and the rest of the code is:

    E(:,i)=10*log10(((P(:)./P(1))));

A=max(find(t<t(1)+0.02)); %find 20ms after the first reflection pf20ms = polyfit(t(1:A),E((1:A)),1); %curve fitting for the 0-20ms time interval plate=polyfit(t(A:length(t)),E((A:length(t))),1); % curve fitting for the rest of the impulse response pfull=polyfit(t,E,1); % curve fitting for the complete impulse response

RTcalcf20ms=(-60-pf20ms(1,2))/pf20ms(1,1); RTcalcplate=(-60-plate(1,2))/plate(1,1); RTcalcpfull=(-60-pfull(1,2))/pfull(1,1);

The values that I got for reflection order 5 and absorption coeffiecent of 0.1 for a room of 17x28x12 are between very about 10 secs for first 20 ms and 0.9 secs for late part where the calculated Reverberation Time from Sabine's formula give 4.58secs.. Thanks

Answer 1

[Daniel Shub]

It is not clear to me what your question is, but that doesn't generally stop me from guessing. The description of what you are doing sounds fine.

The FEX has a number of submissions dealing with reverberation: http://www.mathworks.com/matlabcentral/fileexchange/?term=reverberation

You could compare your implementation of the room image model with others, to make sure you are 110% sure that yours is correct. You could then compare your implementation of Schroeder's backwards integral with others. You could also run your implementation of the backwards integral against real RIRs and compare to Sabine's predictions.