function [output,output2]=rt(x,fs,bands,dBstart,dBend,plotparameter)
%----------------------------------
%function rt
%output=rt(x,fs,bands,dBstart,dBend,plotparameter)
%x=signal
%fs=samplig frequency
%bands=1 (third octave bands) or 2 (octaves)   
%dBstart and dBend: example -5 -15->EDT
%plotparameter=1 ->the function will plot de decay curves
%----------------------------------

if nargin<6
   disp('----------------------------------');
   disp('function rt');
   disp('output=rt(x,fs,bands,dBstart,dBend,plotparameter)');
   disp('x=signal');
   disp('fs=samplig frequency');
   disp('bands=1 (third octave bands) or 2 (octaves)');   
   disp('dBstart and dBend: example -5 -15->EDT');
   disp('plotparameter=1 ->the function will plot the decay curves');
   disp('alternatively [output,output2]=rt... gives a second output with the central frequency list');
   disp('----------------------------------');
   return
end
warning off

np1=size(x);
if min(np1)>1
   disp('ERROR')
   disp('this function only works with mono signals')
   return
end

if np1(2)>np1(1)
   x=x';
end

if bands==1
   fcent=1000*(2.^((-10:7)/3));
else
   fcent=1000*(2.^((-3:2)));
end

output2=fcent;

disp('calling funtion myfilter')

y=filtra(x,fs,bands);
yy=cumsum(y.^2);

for ii=1:length(fcent);
   yy(:,ii)=10*log10(1-yy(:,ii)/max(yy(:,ii)));
end



for ii=1:length(fcent)
   [useless,posi_dBend]=min(abs(yy(:,ii)-dBend));
   [useless,posi_dBstart]=min(abs(yy(:,ii)-dBstart));
   

   output(ii)=60/(-dBend+dBstart)*(posi_dBend-posi_dBstart)/fs;
end


if plotparameter==1
   plot((0:length(yy')-1)/fs,yy')
   
   if bands==1
      legend('100','125','160','200','250','315','400','500','630','800','1000','1250','1600','2000','2500','3150','4000','5000');
   else
      legend('125','250','500','1000','2000','4000');
   end
axis( [0  (length(yy')-1)/fs  -60 0])   
end
warning on
return
