find_tonal_components.m

音频的水印加入算法

function [Flags, Tonal_list, Non_tonal_list] = Find_tonal_components(X)
%
%   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. 
%----------------------------------------------------------------------------
%
% [Flags, Tonal_list, Non_tonal_list] = Find_tonal_components(X) 
% Identifies and lists both tonal and non - tonal components of the audio
% signal. 
% 
% -- INPUT --
% X: Power density spectrum vector in dB scale of one block (with FFT_SIZE
% samples).
%
% -- OUTPUT --
% Flags: Vector of size FFT_SIZE/2 containing the markers TONAL or NON_TONAL. 
%
% Tonal_list: Matrix with two columns for the TONAL components. The first
% column for the frequency indexes (INDEX) in linear scale and the second one
% for the sound pressure level (SPL) [dB].
%
% Non_tonal_list: Matrix with two columns for the NON_TONAL components. The
% first column for the frequency indexes (INDEX) in linear scale and the second
% one for the sound pressure level (SPL) [dB].
% ------------
%  
%   See also Decimation
%---------------------------------------------------------------------------- 

global FFT_SIZE MIN_POWER NOT_EXAMINED IRRELEVANT INDEX SPL TONAL BARK
global NON_TONAL TH Map CB

% List tonal components and compute sound pressure level 
% Index array without RIGHT margin
jnd = 1:254; 
% Index array without LEFT margin
ind = 2:256; 

% a) Labelling of local maxima
% Check X(k-1) < X(k): IF X(k-1) < X(k) 
% THEN delay_1(k) = 1 ELSE delay_1(k) = 0 
delay_1 = ceil(0.5 + sign(X(ind) - X(ind - 1)) / 2);
                                                    
% Check X(k-1) < X(k) AND X(k) >= X(k+1): 
% IF X(k) = local maxima THEN Flags(k) = 1 ELSE Flags(k) = 0 (exclude
% margins) 
Flags = [0 delay_1(jnd) &~ delay_1(jnd + 1) 0];

% b) Labelling of tonal components and calculation of the sound pressure level
% IF X(k) == lokal maxima && X(k) - X(k+j) >= 7dB THEN Flags(k) = 1 ELSE
% Flags(k) = 0
 
% X(k) - X(k+j) >= 7dB
% The multiplication of Flags with the array of "0" and "1" implements the
% AND operation. Every additional multiplication AND the Flags from previous
% calculation.

% j = -2 , +2
ind = 3:254;
Flags = Flags .* [0 0 ceil(0.5 + sign(X(ind) - X(ind - 2) - 7) / 2) 0 0];
Flags = Flags .* [0 0 ceil(0.5 + sign(X(ind) - X(ind + 2) - 7) / 2) 0 0];

% j = -3, +3
ind = 63:253; B = zeros(1, 3);
%  one(n,m) is for saveing the results form previous calculations 
tmp = [ones(1, 253 - size(ind, 2)) ...
	   ceil(0.5 + sign(X(ind) - X(ind - 3) - 7) / 2) ...
	   B];
Flags = Flags .* tmp;
tmp= [ones(1, 253 - size(ind, 2)) ...
	  2.* ceil(0.5 + sign(X(ind) - X(ind + 3) - 7) / 2) ...
	  B];
Flags = Flags .* tmp;

% j = -4, +4
ind = 127:252; B = zeros(1, 4);
tmp = [ones(1, 252 - size(ind, 2)) ...
	   ceil(0.5 + sign(X(ind) - X(ind - 4) - 7) / 2) ...
	   B];
Flags = Flags.*tmp;
tmp = [ones(1, 252 - size(ind, 2)) ...
	   ceil(0.5 + sign(X(ind) - X(ind + 4) - 7) / 2) ...
	   B];
Flags = Flags.*tmp;

%j = -5,+5
ind = 127:251; B = zeros(1, 5);
tmp = [ones(1, 251 - size(ind, 2)) ...
	   ceil(0.5 + sign(X(ind) - X(ind - 5) - 7) / 2) ...
	   B];
Flags = Flags.*tmp;
tmp = [ones(1, 251 - size(ind, 2)) ...
	   ceil(0.5 + sign(X(ind) - X(ind + 5) - 7) / 2) ...
	   B];
Flags = Flags.*tmp;

%j = -6,+6
ind = 127:250; B = zeros(1, 6);
tmp = [ones(1, 250 - size(ind, 2)) ...
	   ceil(0.5 + sign(X(ind) - X(ind - 6) - 7) / 2) ...
	   B];
Flags = Flags.*tmp;
tmp = [ones(1, 250 - size(ind, 2)) ...
	   3.*ceil(0.5 + sign(X(ind) - X(ind + 6) - 7) / 2) ...
	   B];
Flags = Flags.*tmp;

% K holds the indexes of the non-zero elements of flags
K = find(Flags)';

% Check if there are tonal components. 
if K, 
  Flags(1, K + 1) = IRRELEVANT;
  Flags(1, K - 1) = IRRELEVANT;
  Flags(1, K - 2) = IRRELEVANT;
  Flags(1, K + 2) = IRRELEVANT; 
  % Set the index and calculate sound pressure level of the tonal list.
  Tonal_list(:, INDEX) = K;
  Tonal_list(:, SPL) = 10.*log( 10.^(X(K - 1) / 10) + ...
								10.^(X(K) / 10) + ... 
								10.^(X(K + 1) / 10) )'./log(10);
else
  % Set empty array to enable check in decimation function.
  Tonal_list = [];
end

K = find(Flags == 2)';
if K, 
  Flags(1, K - 3) = IRRELEVANT;
  Flags(1, K + 3) = IRRELEVANT;
  % ?Reset value of flag back to TONAL (=1).
  Flags(1, K) = TONAL;
end

K = find(Flags == 6)';
if K, 
  Flags(1, K - 2) = IRRELEVANT;
  Flags(1, K + 2) = IRRELEVANT;
  Flags(1, K - 3) = IRRELEVANT;
  Flags(1, K + 3) = IRRELEVANT;
  
  Flags(1, K - 4) = IRRELEVANT;
  Flags(1, K + 4) = IRRELEVANT;
  Flags(1, K - 5) = IRRELEVANT;
  Flags(1, K + 5) = IRRELEVANT;
  Flags(1, K - 6) = IRRELEVANT;
  Flags(1, K + 6) = IRRELEVANT;
  Flags(1, K) = TONAL;
end
F = Flags;

% List the non tonal components and compute power.  All the spectral lines that
% have not been examined during the previous search are summed together to form
% the non - tonal component.
Non_tonal_list = []; 

% For each critical band, compute the power in non-tonal components. 
% CB : Critical Band Boundaries indexes for 44.1kHz see pp. 123 D.2b 
% TH : Frequency, critical band rates and absolute threshold see pp. 117
% D.1b  

% Loop the critical bands for a specific sampling rate.
for i = 1:length(CB(:, 1)) - 1 
  % CB is a 25 element column vector. It contains the indexes for the top
  % end of each critical band expressed in the non-linear indexes. There are 25
  % critical bands. 
  %
  % TH is a 106x3 matrix. TH(:, 1) are frequency INDEXES in LINEAR scale at
  % the top of each critical band. 
  %
  % K contains all frequency indexes in linear scale in the critical band.
  % INDEX = 1
  % lowIdx = TH(CB(i), INDEX)
  % highIdx = TH(CB(i + 1), INDEX)
  K = TH(CB(i), INDEX):TH(CB(i + 1), INDEX) - 1;

  % Kf: list of indices in the critical band which are = 0 (Non-tonal noise
  % components)
  Kf = find([ones(1, K(1) - 1) ...
			 Flags(K) ...
			 ones(1, length(Flags) - K(length(K)))] == 0)';

  if Kf,      
	% To avoid calculation of log10(0) if X(Kf) == 0, add MIN_POWER to the sum!
	power = 10*log( sum([10.^(X(Kf)/10) 10^(MIN_POWER/10)]) )./log(10); 

	% Calculate the distance between each non-tonal component Kf and the lower
    % limit of the current critical band to obtain the weights.
	% Calculate the power of the non-tonal components in the critical band.

	% OLD calculation
	% weight = sum( 10.^(X(Kf)/10).*( TH(Map(Kf), BARK)' - i ) );

	% NEW calculation
	weight = sum( 10.^(X(Kf)/10).*( TH(Map(Kf), BARK)' - TH(CB(i), BARK)' ) );

    % K(length(K)) + 1 - K(1) is the length of critical band in LINEAR scale.
	% Calculate index nearest to the geometrical mean of the critical band
	index = K(1) + round( weight / 10^(power/10)*( K(length(K)) + 1 - K(1) ) );
	Flags(1, Kf) = IRRELEVANT; 	
  else 
	% There are no non-tonal components in the current critical band.
	% Take the mean as index for the non-tonal component. round() rounds 
	index  = round( mean([K(1) K(length(K)) + 1]) );
	power = MIN_POWER;
  end 
  
  % If index < length(Flags) AND index < 1 set index = 1.
  % If index > length(Flags) set index = length(Flags).
  index = max(min(index, length(Flags)), 1); 

  % If the calculated index for the non-tonal component is already flagged as
  % TONAL use the next one.
  if Flags(index) == TONAL
 	index = index + 1;
  end
  % Set index and sound pressure level of the non-tonal list and Flag array.
  Non_tonal_list(i, INDEX) = index; 
  Non_tonal_list(i, SPL) = power; 
  Flags(index) = NON_TONAL; 
end