fft_analysis.m

音频的水印加入算法

function [X, Delta] = FFT_Analysis(Input)
%
%   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. 
%----------------------------------------------------------------------------
% [X, Delta] = FFT_Analysis(Input)
%
% A Hanning window is applied before computing the FFT.  Finaly, a
% normalisation of the sound pressure level (SPL) is done such that the maximum
% value of the spectrum is 96dB; the dB added are stored in Delta output.
%
% One should take care that the Input is not zero at all samples.
% Otherwise W will be -INF for all samples.
%
% -- INPUT --
% Input: Vector or matrix of float values in the range [-1, 1] representing the
% audio samples. If Input is a matrix with Blocksize (= FFT_SIZE = 512) rows
% the columns of the matrix represent one block of audio samples (with
% Blocksize samples). The number of columns corresponds to the number of blocks
% analyzed simultaenously. 
%
% -- OUTPUT --
% X: Vector or matrix of the power density spectrum normalized to max(X) =
% +96dB. According to Input X is a row, a column vector or a matrix. In the
% latter case the columns of the matrix correspond to different blocks analysed
% simultaneously.
  
% Delta: Delta = 96dB - max(X). Corresponding to the input X, Delta is a scalar
% or a matrix with dimension size(X, 1) x size(X, 2) (= Blocksize x Number Of
% Blocks) with identical rows. Each column has identical entries corresponding
% to the difference to 96dB for each block.
%--------------

global FFT_SIZE MIN_POWER;
% Initialising output 
Delta = []; 
X = [];

% Extract dimension of input array. Expand Hanning window to matrix of
% identical rows according to the input array.
[rows, cols] = size(Input);

% Check if input consists of only one row.
if rows == 1 

  % Prepare the Hanning window according to the number of columns. Adjust
  % Hanning Window of Matlab according to MPEG Standard.
  h = sqrt(8 / 3) * hann(cols, 'periodic')'; 

else

  % If more than one column, take columns as input vectors. Expand window
  % to matrix of identical cols according to the input array. 
  h = repmat( sqrt(8 / 3) * hann(rows, 'periodic'), 1, cols);

end

% Power density spectrum. If the power density is smaller -200 (MIN_POWER) ->
% set to -200. If X is a matrix the calculation will be performed columnwise.

Input_FFT = fft( (Input + eps).*(h + eps) )/FFT_SIZE;
X = max( 20*log( abs(Input_FFT))./log(10), MIN_POWER );

% Normalization to the reference sound pressure level of 96 dB.
% X is vector -> Delta is a scalar.
% X is matrix -> Delta is a row vector containing the difference of
% each column of X to +96 dB. 
 
Delta = 96 - max(X); 

% Expand Delta row vector to matrix with dimension: rows x size(X, 2) (=
% Blocksize x Number Of Blocks) for further calculations.
if cols > 1
  % The rows of array Delta will be identical.
  Delta = repmat(Delta, rows, 1);
end
  
% Normalize X array to SPL +96dB.
X = X + Delta;