% HTHRES Hearing threshold
%
% H. S. Malvar - Nov'00, (c) H. S. Malvar
%
% Syntax: Ht = hthres(X,fs)
%
% X = vector of frequency magnitude components
% fs = sampling frequency, in Hz
% Ht = vector of thresholds, in rms
function Ht = hthres(X,fs)
Dabs = 60; % in dB, how much to bring down Fletcher-Munson curves;
% it depends on the assumed playback level
Thmin = -60; % in dB, how far down can any threshold go
Thmax = 60; % in dB, how far up can any threshold go
Rfac = 8; % Power reduction factor in dB for masking within the
% same bark frequency band
% Bark subbands upper limits
Bh = [100 200 300 400 510 630 770 920 1080 1270 1480 1720 2000 2320 2700 3150 3700 4400 5300 6400 7700 9500 12000 15500 22050];
P = X.*X; % power spectrum
TdB = 0*P; % Threshold in dB
Nbands = length(X);
f=[0:Nbands-1]'*fs/(2*(Nbands-1)); % subband center frequencies, in Hz
% Loop to generate Bark power spectrum and Bark center thresholds
dindx = Nbands*2/fs; % # of coefficients per Hz
stop = 0;
i1 = 1;
il=[]; iu=[];
for i = 1:25 % scan all 25 Bark bands
% i-th Bark subband covers original subbands i1 to i2
i2 = round(Bh(i)*dindx);
if i2 > Nbands
i2 = Nbands;
stop = 1; % stop loop if signal bandwidth is reached
end
il = [il;i1];
iu = [iu;i2];
% Average RMS signal amplitude over Bark band [i1,i2]
Arms = sqrt(mean(P(i1:i2)));
% average signal amplitude over Bark band [i1,i2], in dB
Adb = 20*log10(Arms+eps);
% Tr = relative threshold given power level within i-th band
Tr = Adb - Rfac; % center threshold
Sbu(i) = Tr;
i1 = i2+1;
if stop break; end;
end
Lsb = i; % number of Barks covered for the given fs
% Cross-Bark spreading
Sbu = Sbu';
Sb = Sbu;
pl = -100;
pr = -100;
pli = Sbu(1);
for i = 1:Lsb
if i < Lsb; pr = Sbu(i+1)-25; else pr = -100; end
if pli > pl; pl = pli; end
pl = pl - 10;
if pl > Sb(i); Sb(i) = pl; end
if pr > Sb(i); Sb(i) = pr; end
pli = Sbu(i);
end
% Adjust thresholds considering absolute masking
for i = 1:Lsb
i1 = il(i);
i2 = iu(i);
Tr = Sb(i);
% Ta = Fletcher-Munson absolute threshold
Ta = mean(thrabs(f(i1:i2)))-Dabs;
if Tr < Ta; Tr = Ta; end
% Clip at minimum and maximum levels
if Tr < Thmin; Tr = Thmin; end
if Tr > Thmax; Tr = Thmax; end
Sb(i) = Tr; % save dB thresholds in vector Sb
TdB(i1:i2) = Tr; % same threshold for all bands within Bark band
end
Ht = 10.^(TdB/20); % convert from power to rms levels
e-REdING. Biblioteca de la Escuela Superior de Ingenieros de Sevilla.
