individual_masking_thresholds.m

音频的水印加入算法

function  [LTt, LTn] = Individual_masking_thresholds(X, Tonal_list, ...
												  Non_tonal_list) 
%
%   Authors: Fabien A.P. Petitcolas ([e-mail])
%            Teddy Furon ([e-mail])
%
%   Corrections:
%            Michael Arnold ([e-mail])
%            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 [phone], Telefax 
%    [phone]. 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 Map

if 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