| [rta, Lya, rtpa]=reverb_time(SP, Fs, fc, t_high, t_low, num_x_filter, pbs, pab, method)
|
function [rta, Lya, rtpa]=reverb_time(SP, Fs, fc, t_high, t_low, num_x_filter, pbs, pab, method)
% % Syntax:
% % [rta, Lya, rtpa]=reverb_time(SP, Fs, fc, t_high, t_low, num_x_filter,pbs, pab, method);
% %
% % This program supports two methods of calculating the reverberation
% % time: the speaker on speaker off method and the balloon pop method.
% % The speaker on speaker off method is considered to be much more
% % accurate. It is highly recommended to use the speaker on speaker off
% % method if possible.
% %
% % This program calculates reverberation times for 1/3 octave broad
% % band noise for multiple microphones. See the description of the
% % variable fc.
% %
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% % Description of the Input Variables
% %
% % SP is the time record in units Pa and has size (num_mics,
% % num_samples). SP is automatically rearranged so that the number
% % of elements in the second dimension are greater than the first
% % dimension.
% %
% % Fs (Hz) is the samping rate.
% %
% % fc (Hz) is the third octave band center frequency.
% %
% % Filters have anomalies and program returns meaningless data
% % for fc > Fs/3.
% %
% % It is difficult to excite source below 31 Hz with enough energy
% % for a meaningful reverberation time.
% %
% % For the balloon pop method, the impulsive reverberant decay must
% % reach the background level.
% %
% % For the Balloon pop method measurement times of 10 seconds or
% % longer are recommended for fc < 200 Hz.
% %
% %
% % t_high (s) is used in the speaker on spekaer off method.
% % t_high is the duration of signal on at the beginning of the time
% % record. This is used to calculate the average signal level.
% %
% % t_low (s) is used by both methods. t_low is the duration of quiet
% % time at the end of the time record. t_low is used to calculate the
% % average background level.
% %
% % num_x_filter is the number of times to filter the time record.
% % The minimum value of num_x_filter is 1.
% % Typically a value of 2 is sufficient to remove the contributions
% % from other third octave bands.
% %
% % If there is spurious background noise, espcially in the low
% % frequencies fc < 100 Hz, then additional filtering may be
% % beneficial. num_x_filter=10 can remove up to 120 dB of noise.
% % Unfortunately the third octave band filter changes the calculated
% % reverberation times, so be careful.
% %
% % Increasing num_x_filter does not necessarily improve accuracy.
% % Increasing the input signal speaker level or balloon pop level
% % may improve the accuracy for both methods.
% %
% % pbs is an acronym for the percentage of the total range in dB below
% % the average signal level. pbs is used to discard spurious
% % high level sound data at the speaker cutoff. pbs=20 usually works.
% %
% % pab is an acronym for the percentage of the range above the average
% % background level. pab is used to discard spurious background
% % noise where the reveberant decay approaches the background level.
% % This is more of a problem in environments with very low background
% % For example, measuring reverberation time in a hemi-anechoic
% % chamber.
% %
% % method selects the speaker on speaker off method or the balloon pop
% % method. 1 for speaker, otherwise balloon pop.
% %
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% % Description of the Output Variables
% %
% % rta % The array of the reverberation times
% %
% % Lya % The array of processed time records
% %
% % rtpa % The array of time records of the solution of the
% % % reverberation time calculation
% %
% % SP2 % Resampled and filtered time record
% %
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Example=''; % Simulated reverberation data
% % For the Speaker on Speaker off method
%
% % This example is for a four microphone with each microphone recording at
% % a different level. The reverberation time is set to a different for each
% % microphone to illustrate the decay line fitting.
%
% tau=[2, 2.3, 3.0, 2.4]; % time constants for exponential decay
%
% tb=5; % seconds, duration of decay signal
%
% Fs=50000; % sampling frequency rate
%
% t=1/Fs*(1:(tb*Fs)); % seconds time array for calculating
% % decay signal
%
% A=[100, 110, 115, 130];
% decay1=A(1).*exp(-1./tau(1).*t).*randn(1,length(t));
% decay2=A(2).*exp(-1./tau(2).*t).*randn(1,length(t));
% decay3=A(3).*exp(-1./tau(3).*t).*randn(1,length(t));
% decay4=A(4).*exp(-1./tau(4).*t).*randn(1,length(t));
%
% % The test signal has
% % 1) nominal one second of signal high
% % 2) reverberant decay
% % 3) nominal one second of background noise
% %
% SP=zeros(4, 2*Fs+length(t));
% SP(1, :)=[A(1)*randn(1,Fs) decay1 A(1).*exp(-1./tau(1).*tb)*randn(1, Fs);];
% SP(2, :)=[A(2)*randn(1,1.5*Fs) decay2 A(2).*exp(-1./tau(2).*tb)*randn(1, 0.5*Fs);];
% SP(3, :)=[A(3)*randn(1,1.2*Fs) decay2 A(3).*exp(-1./tau(3).*tb)*randn(1, 0.8*Fs);];
% SP(4, :)=[A(4)*randn(1,0.8*Fs) decay2 A(4).*exp(-1./tau(4).*tb)*randn(1, 1.2*Fs);];
%
% figure(1); for e1=1:4; subplot(4, 1, e1); plot(SP(e1, :)); hold on; end;
%
% Fs=50000; % (Hz) Sampling rate
%
% fc=400; % (Hz) Center frequency of the third-octave band
%
% t_high=0.5; % (s) duration of loud time at the beginning of the time record
% % used to calculate the average signal level,
% % which is usead to estimate the speaker cutoff time
% % A longer t_high time is better when there is filter phase delay
% % especially at low frequencies (Fc < 100) or
% % if the speaker turns onslowly
%
% t_low=0.5; % (s) duration of quiet time at the end of the time record
% % used to calculate the average background level for
% % calculating the time when the background noise level is reached
% % a longer time gets the background level above
% % the dips and peaks in the background level from
% % sudden quiet and background dog barks etc...
%
% num_x_filter=2; % Number of times to filter the data. Minimum value is 1
% % typically a value of 2 is sufficient to remove noise.
% % num_x_filter=10 can remove 120 dB of noise.
% % Increasing num_x_filter changes the cacualted
% % reverberation time.
%
% pbs=15; % Percentage of the total range in dB below the average signal
% % level to discard from fitting to a line. This avoids
% % spurious high level sound data at the speaker cutoff.
%
% pab=40; % Percentage of the range above the average background level
% % to discard. Increases the background cutoff level above
% % spurious background noise where the reveberation decay
% % approaches the background level. A value of 40 worked in
% % a laboratory hemi-anechoic chamber and may work for
% % acoustical environments with very low background noise levels.
%
% pab=15; % A value of 15 worked for an animal shelter with dogs barking
% % in the surrounding rooms.
%
% method=1; % 1 is the speaker on speaker off method.
% % otherwise is the balloon pop method
%
% [rta, Lya, rtpa]=reverb_time(SP, Fs, fc, t_high, t_low,num_x_filter,pbs,pab, method);
%
% % Exact reverberation times are
% t_60=log(1000).*tau;
%
% % The percent error is
% Percent_error=100.*(1-(t_60./rta'));
% %
% %
%
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Example=''; % Simulated reverberation data
% % For the Balloon Pop Method
%
% % This example is for a four microphone with each microphone recording at
% % a different level. The reverberation time is set to a different for each
% % microphone to illustrate the decay line fitting.
%
% tau=[2, 2.3, 3.0, 2.4]; % time constants for exponential decay
%
% tb=5; % seconds, duration of decay signal
%
% Fs=50000; % sampling frequency rate
%
% t=1/Fs*(1:(tb*Fs)); % seconds time array for calculating
% % decay signal
%
% A=[100, 110, 115, 130];
% decay1=A(1).*exp(-1./tau(1).*t).*randn(1,length(t));
% decay2=A(2).*exp(-1./tau(2).*t).*randn(1,length(t));
% decay3=A(3).*exp(-1./tau(3).*t).*randn(1,length(t));
% decay4=A(4).*exp(-1./tau(4).*t).*randn(1,length(t));
%
% % The test signal has
% % 1) nominal one-tenth of a second of background noise
% % 2) reverberant decay
% % 3) nominal one second of background noise
%
% SP=zeros(4, 2*Fs+length(t));
% SP(1, :)=[0.01*A(1)*randn(1,0.1*Fs) decay1 A(1).*exp(-1./tau(1).*tb)*randn(1, Fs);];
% SP(2, :)=[0.01*A(2)*randn(1,0.1*Fs) decay2 A(2).*exp(-1./tau(2).*tb)*randn(1, Fs);];
% SP(3, :)=[0.01*A(3)*randn(1,0.1*Fs) decay2 A(3).*exp(-1./tau(3).*tb)*randn(1, Fs);];
% SP(4, :)=[0.01*A(4)*randn(1,0.1*Fs) decay2 A(4).*exp(-1./tau(4).*tb)*randn(1, Fs);];
%
% Fs=50000;
% fc=400;
% t_high=0.5;
% t_low=0.5;
% num_x_filter=2;
% pbs=15;
% pab=15;
% method=2;
%
% [rta, Lya, rtpa]=reverb_time(SP, Fs, fc, t_high, t_low,num_x_filter,pbs,pab, method);
%
% % Exact reverberation times are
% t_60=log(1000).*tau;
%
% % The percent error is
% Percent_error=100.*(1-(t_60./rta'));
%
% %
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% % Uses the Signal Processing Toolbox
% % Uses the following progrmas which can be found on matlab central
% % Uses the octave program set from Christophe Couvreur
% % Uses programs from Alexandros Leontitsis
% % Uses the modified programs from Alexandros Leontitsis
% %
% % List of Subprograms
% % findextrema From Schuberth Schuberth found on matlab central
% % oct3dsgn From Christophe Couvreur found on matlab central
% % sub_mean
% % LMTSreg A modified program original from Alexandros Leontitsis
% % found on matlab central
% % LMSloc From Alexandros Leontitsis found on matlab central
% % LMS_trim
% % rand_int
% % LMTSregor A modified program original from Alexandros Leontitsis
% % found on matlab central
% %
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% % Program was written by Edward L. Zechmann
% %
% % date 16 October 2007
% %
% % modified 13 January 2008 updated comments
% % changed from filter.m to filtfilt.m
% % Added impulsive nosie method
% %
% % modified 14 January 2008 Reduced number of points allowed for
% % the robust line fit
% % updated comments
% %
% % modified 13 February 2008 Added the balloon pop method and
% % example
% % Made the fit adaptive for both methods
% % Added a leverage point to the
% % Speaker on speaker off method when the
% % decay line fit does not properly
% % intersect the signal high time record
% %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% % Please feel free to modify this code.
% %
if nargin < 9
method=1;
end
fc=round(fc);
[SP]=convert_double(SP);
[m1, n1]=size(SP);
if m1 > n1
SP=SP';
[m1 n1]=size(SP);
end
if nargin < 6
num_x_filter=1;
else
if num_x_filter < 1
num_x_filter=1;
end
end
dt=1/Fs;
%Resample timeseries so that filter has better performance
p=1;
q=1;
flag1=0;
Fsn=Fs;
if fc < Fs/50
Fsn=40*fc;
flag1=1;
elseif fc > Fs/10
Fsn=10*fc;
flag1=1;
end
if flag1 == 1
dtn=1/Fsn;
[p,q] = rat(dt/dtn);
buf=resample(SP(1, :), p, q);
n2=length(buf);
SP2=zeros(m1,n2);
SP2(1, :)=buf;
clear('buf');
if m1 > 1
for e1=2:m1;
SP2(e1, :)=resample(SP(e1, :), p, q);
end
end
Fs2=Fs*p/q;
else
Fs2=Fs;
SP2=SP;
end
% Calculate the filter coefficients for the third octave filter
[Bc,Ac] = oct3dsgn(fc,Fs2,3);
clear('SP');
for e1=1:m1;
for e2=1:num_x_filter;
SP2(e1, :)=filtfilt(Bc,Ac,SP2(e1, :)); %Process data with a 1/3 octave bandpass filter
end
end
if flag1 == 1
buf=resample(SP2(1, :), q, p);
SP=zeros(m1,length(buf));
SP(1, :)=buf;
clear('buf');
if m1 > 1
for e1=2:m1;
SP(e1, :)=resample(SP2(e1, :), q, p);
end
end
else
SP=SP2;
Fs=Fs2;
end
clear('SP2');
% return SP to size [m1, n1]
% check size of SP zero pad if necessary
[m2 n2]=size(SP);
if n1 ~= n2
if n2 > n1
SP=SP(1:m1, 1:n1);
else
a=min(min(abs(SP)));
for e1=1:m1;
for e2=1:n1;
if e2 <= n2
SP(e1, e2)=SP(e1, e2);
else
SP(e1, e2)=a;
end
end
end
end
end
[m1 n1]=size(SP);
SP=SP-mean(SP, 2)*ones(1, n1);
a=1:m1;
SP = abs(hilbert(SP')');
I_high2=floor(t_high/dt);
I_low2=floor(t_low/dt);
rta=zeros(m1,1);
Lya=zeros(m1,n1);
rtpa=zeros(m1,n1);
for e1=1:m1;
if isequal(method, 1)
% Process data for the speaker on speaker off method
[buf1, y22 ]=sub_mean(SP(e1, 1:n1), Fs, 25);
else
% Process data for the balloon pop method
[buf1, y22 ]=sub_mean(SP(e1, 1:n1), Fs, 100);
end
Ly=20*log10(y22(1, :)./(20*10^(-6)));
% reset the values for I_high and I_low adn the other quantities
I_high=I_high2;
I_low=I_low2;
IX_est=0;
IX3=n1;
count=0;
t1=0;
t2=n1-I_low;
IX=1;
flag=0;
while (logical(abs(IX3-IX_est)/IX3 > 0.1) || logical(abs(I_high-t1)/I_high > 0.1) || logical(abs(I_low-(n1-t2))/I_low > 0.1)) && logical(count < 3)
if isequal(count, 0)
IX3=n1-I_low;
end
count=count+1;
I_high=max([I_high, IX]);
buf1=[n1-t2, n1-IX3];
buf1=buf1(buf1>0);
I_low=min(buf1);
if I_low > n1
I_low=n1-1;
end
if I_low < 1
I_low=1;
end
if I_high > n1
I_high=n1-1;
end
if I_high < 1
I_high=1;
end
if isequal(method, 1)
% Level at beginning
Lb=sqrt(sum(y22(1, 1:I_high)'.^2)./I_high);
else
% Level of the balloon pop
Lb=max(y22(1, :));
end
SPLb=20*log10(Lb./(20*10^(-6)));
% Level at the end
Le=sqrt(sum(y22(1, (n1-I_low):n1)'.^2)./I_low);
SPLe=20*log10(Le./(20*10^(-6)));
% Find the beginning and ending of the decay
% Set cutoff level pbs percent of dB range below the signal level
% Set cutoff level pab percent of dB range above the background level
% Find last data points where the level is below the signal level or
% below the balloon pop level
if isequal(method, 1)
IX=find( Ly > SPLb-pbs/100*(SPLb-SPLe), 1, 'last' );
% find data points where Ly is below the background level
IX2=find( Ly < SPLe+pab/100*(SPLb-SPLe) );
% find the first point below the balloon pop level, but above the
% background level
IX3=IX2(find(IX2 > IX, 1 ));
else
[max_LY IX]=max(Ly);
IX2=find( Ly < SPLe+pab/100*(SPLb-SPLe) );
IX3=IX2(find(IX2 > IX, 1 ));
end
% Make sure there are plenty of data points between IX and IX3
if IX3 - IX < 1000;
if IX3+1000 > n1
IX=IX-1000;
else
IX3=IX3+1000;
end
end
if length(IX:IX3) > 5
if isequal(method, 1)
% SPeaker On Speaker Off Method
IX3=min([max([IX3, IX_est]), n1]);
if IX3 > n1
IX3=n1-1;
end
if IX3 < 1
IX3=1;
end
rtd=Ly(1, IX:IX3);
tt=dt*(IX:IX3);
[LMSout,blms,Rsq]=LMTSreg(rtd',tt', 1000, 100000);
p=[blms(2) blms(1)];
t2=floor(Fs*(SPLe-p(2))/p(1));
IX_est=t2;
else
% Balloon Pop Method
IX3=max(IX3, IX_est);
if IX3 > n1
IX3=n1-1;
end
if IX3 < 1
IX3=1;
end
rtd=Ly(1, IX:IX3);
tt=dt*(IX:IX3);
num_pts=length(rtd);
num_bins=min([num_pts/5, 1000]);
rbp=floor(num_pts/num_bins);
bins=(1:rbp:num_pts);
epts=zeros((length(bins)-1), 1);
for e2=1:(length(bins)-1);
[buf, eptsix]=max(rtd(bins(e2):bins(e2+1)));
epts(e2)=bins(e2)-1+eptsix;
end
[ma,mi]=findextrema(rtd);
rtd_n1=length(rtd);
ma=unique(round(ma));
ma=intersect(ma, epts);
if length(ma) > 10
if ma(1) < 1
ma(1)=1;
end
if ma(end) > rtd_n1
ma(end)=rtd_n1;
end
ma=unique(ma);
[LMSout,blms,Rsq]=LMTSregor(-max_LY+rtd(ma)',-tt(1)+tt(ma)', 1000, 100000);
p=[blms(1) max_LY];
else
[LMSout,blms,Rsq]=LMTSregor(-max_LY+rtd',-tt(1)+tt', 1000, 100000);
p=[blms(1) max_LY];
end
t2=floor(Fs*(SPLe-p(2))/p(1));
IX_est=IX+t2;
end
else
count=1000000;
p(1)=60;
p(2)=SPLb;
rtp=SPLb*ones(1, n1);
end
t1=ceil(Fs*(SPLb-p(2))/p(1)); % first intersection point
t2=floor(Fs*(SPLe-p(2))/p(1)); % second intersection point
if ((-t1+IX)/IX > 0.3 || logical(p(1) > 0)) && isequal(method, 1)
% Method failed
% a leverage point is necessary to ensure a good fit
% The intesection point of the signal high cutoff limit and
% the instantaneous sound level data will be the leverage point
flag=1;
count=4;
[LMSout,blms,Rsq]=LMTSregor(-Ly(IX)+rtd',-dt*(IX)+tt', 1000, 100000);
p=[blms(1) Ly(IX)];
t1=ceil(Fs*(SPLb-p(2))/p(1)); % first intersection point
t2=floor(Fs*(SPLe-p(2))/p(1)); % second intersection point
end
rty = polyval(p,dt*(t1:t2));
if isequal(method, 1) && ~isequal(flag, 1)
rtp=[ SPLb*ones(1, t1-1) rty SPLe*ones(1, (n1-t2))];
else
rtp=[ SPLb*ones(1, IX+t1-1) rty SPLe*ones(1, (n1-IX-t2))];
end
end
rt=-60/p(1); % Estimated Reveberation Time
rta(e1)=rt; % Array of estimated reverberation time measurements
Lya(e1, :)=Ly; % Time record of processed time record
rtpa(e1, :)=rtp(1, 1:n1); % Time record of reverberation time solution
end
|
|
