function  [LTt, LTn] = Individual_masking_thresholds(X, Tonal_list, ...												  Non_tonal_list) %%   Authors: Fabien A.P. Petitcolas (fapp2@cl.cam.ac.uk)%            Teddy Furon (furont@thmulti.com)%%   Corrections:%            Michael Arnold (arnold@igd.fhg.de)%            Fraunhofer Institute for Computer Graphics (IGD)     % %   References: %    [1] Information technology -- Coding of moving pictures and associated %        audio for digital storage media at up to 1,5 Mbits/s -- Part3: audio. %        British standard. BSI, London. October 1993.Implementation of %        ISO/IEC 11172-3:1993. BSI, London. First edition 1993-08-01. %%   Legal notice: %    This computer program is based on ISO/IEC 11172-3:1993, Information %    technology -- Coding of moving pictures and associated audio for digital %    storage media at up to about 1,5 Mbit/s -- Part 3: Audio, with the %    permission of ISO. Copies of this standards can be purchased from the %    British Standards Institution, 389 Chiswick High Road, GB-London W4 4AL,  %    Telephone:+ 44 181 996 90 00, Telefax:+ 44 181 996 74 00 or from ISO, %    postal box 56, CH-1211 Geneva 20, Telephone +41 22 749 0111, Telefax %    +4122 734 1079. Copyright remains with ISO. %----------------------------------------------------------------------------%%   [LTt, LTn] = Individual_masking_thresholds(X, Tonal_list, Non_tonal_list) %   Individual masking thresholds for both tonal and non-tonal  %   components are set to -infinity since the masking function has an %   infinite attenuation beyond -3 and +8 barks, that is the component %   has no masking effect on frequencies beyond thos ranges [1, pp. 113--114] %% -- INPUT --% X: Power density spectrum in dB scale of one block with FFT_SIZE samples.%% Tonal_list: Matrix with two columns for the TONAL components. The first% column for the indexes (INDEX) and the second one for the sound pressure% level (SPL) of the corresponding components.%% Non_tonal_list: Matrix with two columns for the NON_TONAL components. The% first column for the indexes (INDEX) and the second one for the sound% pressure level (SPL) of the corresponding components.%% -- OUTPUT -- % LTt: Matrix of size (length Tonal_list)x(length(TH(:, 1) = 106)% LTn: Matrix of size (length Non_tonal_list)x(length(TH(:, 1) = 106)% ------------global MIN_POWER INDEX BARK TH Mapif isempty(Tonal_list)    LTt = []; else   % Create (length Tonal_list)x(length(TH(:, 1) = 106) matrix with MIN_POWER  % entries.   LTt = zeros(length(Tonal_list(:, 1)), length(TH(:, 1))) + MIN_POWER; end % Create (length Non_tonal_list)x(length(TH(:, 1) = 106) matrix with MIN_POWER% entries.LTn = zeros(length(Non_tonal_list(:, 1)), length(TH(:, 1))) + MIN_POWER;  % Only a subset of the samples are considered for the calculation of the% global masking threshold. The number of these samples depends on the% sampling rate and the encoding layer. All the information needed is in TH% which contains the frequencies, critical band rates and absolute% threshold.% Calculation of the loudness threshold for tonal maskers LTt.if not(isempty(Tonal_list))  % Create zi (Length of Tonal_list)x106 matrix with equal rows containing all  % Bark-Values  zi = ones(length(Tonal_list(:, 1)), 1) * TH(:, BARK)';     % Create zj (length of Tonal List)x106 matrix, each column hold the Bark  % value of the tonal masker. This is done by mapping the frequency  % index of the masker onto the corresponding Bark value.  zj = TH(Map(Tonal_list(:, INDEX)), BARK) * ones(1, length(TH(:, BARK)));    % Generate matrix describing the distance to the masker in Bark values!   dz = zi - zj;     % Generate zero array of size dz.  vf = zeros(size(dz));       % Masking function    % Generate (length Tonal_list)x106 matrix containing equal columns with  % the SPLs of the tonal maskers.  XT = X(Tonal_list(:, INDEX)) * ones(1, size(dz, 2));     vf = vf + (17 * (dz + 1) - (0.4 * XT + 6)) .* (( -3 <= dz) & (dz <- 1));    vf = vf + (0.4 * XT + 6) .* dz .* (( -1 <= dz) & (dz < 0));	  vf = vf -(17 * dz .* ((0 <= dz) & (dz < 1))); 	  vf = vf + (-(dz - 1) .* (17 - 0.15 * XT) - 17) .* ((1 <= dz) & (dz < 8));    % Select relevant component of the masking function according to the bark  % range.  select = (dz >= -3) & (dz < 8);  % Calculate masking index for tonal maskers avtm.   avtm = (-1.525 - 0.275 .* zj - 4.5);   % Calculate individual masking threshold. Add MIN_POWER to the components,  % which are out of the BARK range -3 <= dz < 8.  LTt = (XT + avtm + vf) .* select + MIN_POWER .* (~select); end    % Calculation of the loudness threshold for non-tonal maskers LTn.if not(isempty(Non_tonal_list))    % Create zi (Length of Non_tonal_list)x106 matrix with equal rows containing  % all Bark-Values.  zi = ones(length(Non_tonal_list(:, 1)), 1) * TH(:, BARK)';    % Create zj (length of Non_tonal List)x106 matrix, each column hold the Bark  % value of the non-tonal masker. This is done by mapping the frequency  % index of the masker onto the corresponding Bark value.  zj = TH(Map(Non_tonal_list(:, INDEX)), BARK) * ones(1, length(TH(:, BARK)));   % Generate matrix describing the distance to the masker in Bark values!  dz = zi - zj;     % Generate zero array of size dz.  vf = zeros(size(dz));       % Masking function     % Generate (length Non_tonal_list)x106 matrix containing equal columns with  % the SPLs of the non-tonal maskers.  XT = X(Non_tonal_list(:, INDEX)) * ones(1, size(dz, 2));	  vf = vf + (17 * (dz + 1) - (0.4 * XT + 6)) .* ((-3 <= dz) & (dz < -1));	  vf = vf + (0.4 * XT + 6) .* dz .* ((-1 <= dz) & (dz < 0));	  vf = vf - (17 * dz .* ((0 <= dz) & (dz < 1))); 	  vf = vf + (-(dz - 1) .* (17 - 0.15 * XT) - 17) .* ((1 <= dz) & (dz < 8));	  % Select relevant component of the masking function according to the bark  % range.  select = (dz >= -3) & (dz < 8);  % Calculate masking index for non-tonal maskers avnm.   avnm = (-1.525 - 0.175 .* zj - 0.5);   % Calculate individual masking threshold. Add MIN_POWER to the components,  % which are out of the BARK range -3 <= dz < 8.  LTn = (XT + avnm + vf) .* select + MIN_POWER .* (~select); end