2.

4 Psychoacoustic Model 35

LM=20
100
LM=40
LM=60
LM=80
Sound pressure level (dB)

80
LM=100

60

40

20

0
−4 0 4 8 12
Critical band rate (Bark)

Fig. 2.17 Spreading function in ISO/IEC Psychoacoustic Model 1

Figure 2.17 shows spreading functions in Model 1 for different levels of the masker.
It is seen that the higher SPL the masker has, the more asymmetric the curve looks.
Specifically, higher frequencies exhibit more masking than lower frequencies when
the level of masker is high. This two-piece linear spreading function is a good
approximation to the masking thresholds of TMT in Fig. 2.16.
In addition, four models described above for spreading functions, i.e., two-slope
SF, Schroeder SF, Psychoacoustic Model 1 SF, and Model 2 SF, are compared
at a level of 80 dB in Fig. 2.18. Among these four models, two-slope spreading
function is the most conservative one, and Model 1 spreading function allows for
more upward spreading of masking than others [51].

2.4.1.2 Implementation of Psychoacoustic Model 1

In different application scenarios, psychoacoustic model can be implemented in

different ways to satisfy the criteria required. ISO/IEC MPEG-1 Standard [77]
utilizes two informative psychoacoustic models, Psychoacoustic Model 1 and 2, to
determine the MMT for inaudibility. Typically, Model 1 is applied to MPEG Layers
I and II and Model 2 to MPEG Layer III. Both models are commonly in use and
well performed. Psychoacoustic Model 1 proposed a low-complication method to
analyze spectral data and output SMR, whereas Psychoacoustic Model 2 performs
a more detailed analysis at the expense of greater computational complexity
36 2 Principles of Psychoacoustics

Two−slope SF
80 Schroeder SF
Model 1 SF
Model 2 SF
Sound pressure level (dB)

60

40

20

0
−8 −4 0 4 8 12
Critical band rate (Bark)

Fig. 2.18 Comparison of four spreading functions relative to an 80 dB masker

[31, 78, 79]. Hence, Psychoacoustic Model 1 for Layer I is later employed in our
audio watermarking scheme in consideration of its higher efficiency.
In our case, the input to Psychoacoustic Model 1 is one frame of audio signal
and the corresponding output is its MMT. The whole procedure of implementation
consists of six steps [72, 73, 77, 80]:
1. FFT analysis and SPL normalization
2. Identification of tonal and nontonal maskers
3. Decimation of invalid tonal and nontonal maskers
4. Calculation of individual masking thresholds
5. Calculation of global masking threshold
6. Determination of the MMT

The details of each step are expounded as follows:

• STEP 1: FFT analysis and SPL normalization

For an accurate analysis of frequency components, fast Fourier transform (FFT) is
performed to obtain a high-resolution spectral estimate of incoming frame x .n/.
In Psychoacoustic Model 1, the input frame has a size of N D 512 points. To
minimize the leakage effect, x .n/ is multiplied with a modified Hanning window
w .n/ defined by
r r   
8 8 1 2 n
w .n/ D hann .N / D  1  cos 0  n  N  1; (2.9)
3 3 2 N
2.4 Psychoacoustic Model 37



where hann .N / D 12 1  cos 2Nn is the N -point Hanning window. Factor
q
8
3
is a gain to compensate the average power of w .n/, so that hw .n/2 i 
X
N 1 h i
2
1
N
w .n/ D 1. Then, power spectral density (PSD) of x .n/ is computed as
nD0

ˇ "N 1  #ˇ2
ˇ1 X 2 nk ˇˇ
ˇ N
PSD .k/ =dB D 10 log10 ˇ x .n/ w .n/ exp j ˇ 0k< :
ˇN N ˇ 2
nD0
(2.10)
After that, PSD estimate PSD .k/ is normalized to a SPL level of 96 dB, i.e., the
maximal is limited to 96 dB.

P .k/ =dB D 96  max fPSD .k/g C PSD .k/

(2.11)
D 4P C PSD .k/ ;

where 4P D 96  max fPSD .k/g. It is because we have no prior knowledge

regarding actual playback levels, the absolute pressure level of a sound can only
be specified by comparing to a reference. To this
end, a sinusoid with amplitude
equal to half of PCM quantizer spacing A0 D 4 2
is defined as having a SPL of
0 dB, i.e., 20 log10 .A0 =A0 / D 0 dB. Consequently, for 16-bit
 PCM data, a  sinusoid
. 2 /
16 1 4
with amplitude equal to the overload level of quantizer Amax D 2
would
16

have a SPL of about 96 dB, i.e., 20 log10 .Amax =A0 / D 20 log10 2  1  96 dB

[51].
An example of the initial and normalized PSD estimates as well as the threshold
in quiet are shown in Fig. 2.19, where the frequencies of two graphs are plotted on
linear and Bark scales, respectively. Note that in psychoacoustic models, an offset
depending on the overall bit rate is employed for the threshold in quiet. It is equal to
12 dB for bit rates no less than 96 kbits=s and 0 dB for bit rates less than 96 kbits=s
per channel [77]. Sound tracks used in our experiments are of CD quality, whose
bit rates are normally greater than 96 kbits=s. Therefore, by comparing Fig. 2.19 to
Figs. 2.7 and 2.8, the threshold in quiet in Fig. 2.19 is shifted downward by 12 dB.
• STEP 2: Identification of tonal and nontonal maskers
On account of “asymmetry of masking,” it is required to discern frequency
components as tonal (i.e., sinusoidal) and nontonal (i.e., noise-like) maskers. Tonal
maskers are selected from local maxima of normalized PSD estimate, P .k/. A local
maxima refers to the maximum PSD within its two neighbors:

P .k/  P .k C 1/ and P .k/  P .k  1/ (2.12)

38 2 Principles of Psychoacoustics

a
Initial PSD
100
Normalized PSD
Threshold in Quiet

50
Sound pressure level (dB)

−50

−100

−150
0 4000 8000 12000 16000 20000
Frequency (Hz)

b
Initial PSD
100
Normalized PSD
Threshold in Quiet

50
Sound pressure level (dB)

−50

−100

−150
1 3 5 7 9 11 13 15 17 19 21 23 25
Critical band rate (Bark)

Fig. 2.19 Initial and normalized PSD estimates (a) Frequency on linear scale (b) Frequency on
Bark scale
2.4 Psychoacoustic Model 39

If the value of a local maxima is at least 7 dB greater than that of its neighboring
components within a certain Bark range Dk , such a maxima will be marked as a
tonal masker. All the tonal components comprise the “tonal” set, STM :

STM D fP .k/ j ŒP .k/  P .k ˙ Dk /  7 dBg ; (2.13)

where Dk varies with different frequency indices.11

8
< f˙2g ; 2 < k < 63 $ 2Fs
N
 63F
N
s
kHz
Dk 2 f˙2; ˙3g ; 63  k < 127 $ 63Fs
 N kHz
127Fs
:
: N
f˙2; ˙3; : : : ; ˙6g ; 127  k  250 $ 127Fs
N
 250F
N
s
kHz

One point to note is that [77] did not specify the value of Dk for 251  k  256,
because the maskers within this range are already dominated by the threshold in
quiet (as seen in Fig. 2.19) and have no contribution to masking threshold. Actually,
it is the first criterion for decimation in Step 3.
As the effect of masking is additive in the logarithmic domain, the SPL of each
tonal component is calculated by
h P .k1/ P .k/ P .kC1/
i
PTM .k/ =dB D 10 log10 10 10 C 10 10 C 10 10 : (2.14)

In addition, the remaining components within each critical band12 are treated to be
nontonal. So we sum up their intensities as the SPL of a single nontonal masker for
each critical band, PNM :
  X  P .j / 
PNM k =dB D 10 log10 10 10 8P .j / … STM ; (2.15)
j

where that k is the frequency index nearest to the geometric mean13 of each critical
band. Correspondingly, all the nontonal components are put into the “nontonal”
set, SNM .

11
The frequency edges are calculated based on the sampling frequency Fs .
12
Critical band boundaries vary with the Layer and sampling frequency. ISO/IEC IS 11172-3 [77]
has tabulated such parameters in Table D.2a–f. In our case, Table D.2b for Layer I at a sampling
frequency of 44.1 kHz is adopted.
!1=M
Y
M
13
The geometric mean of a data set Œa1 ; a2 ; : : : ; aM  is defined as am . It is sometimes
mD1
!1=M " #
Y
M X
M
called the log-average, i.e., am D 10ˆ 1
M
log10 .am / .
mD1 mD1
40 2 Principles of Psychoacoustics

a 120
Normalized power spectrum
Local maxima
100 Tonal components
Non−tonal components
Threshold in Quiet
80 Removed components
Sound pressure level (dB)

60

40

20

−20
0 4000 8000 12000 16000 20000
Frequency (Hz)
b 120
Normalized power spectrum
Local maxima
100 Tonal components
Non−tonal components TM
Threshold in Quiet TM 2
TM3
1
80
Sound pressure level (dB)

Removed components

60
16.55 17.08 17.57 18.01 Barks

40 0.49 Bark

20

−20
1 3 5 7 9 11 13 15 17 19 21 23 25
Critical band rate (Bark)

Fig. 2.20 Tonal and nontonal maskers (a) Frequency on a linear scale (b) Frequency on Bark scale

Tonal and nontonal maskers are denoted by pentagram and asterisk symbols in
Fig. 2.20, respectively. Particularly, the associated critical band for each masker is
indicated in the graph on Bark scale.
• STEP 3: Decimation of invalid tonal and nontonal maskers
On considering their possible contributions to masking threshold, the sets of tonal
and nontonal maskers are examined according to two criteria as follows:
2.4 Psychoacoustic Model 41

One rule is that any tonal and nontonal maskers below the threshold in quiet
are removed. That is, only the maskers that satisfy Eq. (2.16) are retained, where
ATH .k/ is the SPL of threshold in quiet at frequency index k:

PTM; NM .k/  ATH .k/ : (2.16)

For example, one of each tonal and nontonal maskers between 24 and 25 Barks is
discarded, as shown in Fig. 2.20b.
The other rule is to simplify any group of maskers occurring within a distance
of 0.5 Bark: only the masker with the highest SPL is preserved and the rest are
eliminated.

PTM; NM .k/ D arg max PTM; NM .k C k0 / : (2.17)

k0 2Œ0:5;0:5

For example, two pairs of tonal maskers between 17 and 19 Barks, fTM1 ; TM2 g
and fTM2 ; TM3 g, are inspected. As shown in an enlarged drawing on the right of
Fig. 2.20b, the distance between fTM1 ; TM2 g is 0.49 Bark, and TM1 has a lower
SPL than TM2 . Therefore, TM2 is preserved, whereas TM1 is removed. Similarly,
we dispose of TM3 but retain TM2 for fTM2 ; TM3 g.
ˇ
ˇ Distance W17:57  17:08 D 0:49 Bark
TM2  fTM1 ; TM2 g ˇˇ
SPL W PTM1 < PTM2
ˇ
ˇ Distance W18:01  17:57 D 0:44 Bark
TM2  fTM2 ; TM3 g ˇˇ
SPL W PTM2 > PTM3

In Fig. 2.20, the invalid tonal and nontonal maskers being decimated are denoted by
a circle.
• STEP 4: Calculation of individual masking thresholds
After eliminating invalid maskers, individual masking threshold is computed for
each tonal and nontonal masker. An individual masking threshold L .j; i / refers to
the masker at frequency index j contributing to masking effect on the maskee at
frequency index i . It corresponds to L Œz .j / ; z .i /, where z .j / and z .i / are the
masker and maskee’s frequencies in Bark scale. In MPEG psychoacoustic models,
only a subset of samples over the whole spectrum are considered to be maskees and
involved in the calculation of global masking threshold. The number and frequencies
of maskees also depend on the Layer and sampling frequency, as tabulated from
Table D.1a–f in [77]. In our case, Table D.1b for Layer I at a sampling frequency of
44.1 kHz is adopted, where 106 maskees are taken into account.
The individual masking thresholds for tonal and nontonal maskers, LTM
Œz .j / ; z .i / and LNM Œz .j / ; z .i /, are calculated by

LTM Œz .j / ; z .i / =dB D PTM Œz .j / C 4TM Œz .j / C SF Œz .j / ; z .i / (2.18)

42 2 Principles of Psychoacoustics

LNM Œz .j / ; z .i / =dB D PNM Œz .j / C 4NM Œz .j / C SF Œz .j / ; z .i / ; (2.19)

where PTM Œz .j / and PNM Œz .j / are the SPLs of tonal and nontonal maskers
at a Bark scale of z .j /, respectively. The term 4X is called masking index, an
offset between the excitation pattern and actual masking threshold. As mentioned
in Sect. 2.3.1, the excitation pattern needs to be shifted by an appropriate amount in
order to obtain the masking curve relative to the masker. Because tonal and nontonal
maskers have different masking capability, i.e., the noise is a better masker than pure
tone, the masking indices of tonal and nontonal maskers are defined separately as
follows [77]:

4TM Œz .j / D 6:025  0:275z .j / (2.20)

4NM Œz .j / D 2:025  0:175z .j / : (2.21)

The term SF Œz .j / ; z .i / is the spreading function discussed already in Sect. 2.4.1.1.

Psychoacoustic Model 1 employs spreading function in Eq. (2.8), rewritten in the
following expression:
8
ˆ
ˆ 17d z  0:4PX Œz .j / C 11; 3  d z < 1
ˆ
ˆ
<.0:4P Œz .j / C 6/ d z; 1  d z < 0
X
10 log10 SF .d z/ =dB D ;
ˆ
ˆ17d z;
ˆ 0  dz < 1
:̂
17d z C 0:15PX Œz .j / .d z  1/ ; 1  dz < 8
(2.22)
where d z is the distance from the maskee to masker, d z D z .i /  z .j /,
as defined in Sect. 2.4.1.1. PX Œz .j / refers to PTM Œz .j / in the case of tonal
masker, otherwise PNM Œz .j / for nontonal masker. Notice that for reasons of
implementation complexity, the masking is no longer considered if d z < 3 Bark or
d z  8 Bark and thereby LTM Œz .j / ; z .i / and LNM Œz .j / ; z .i / are set to 1 dB
outside the above ranges [77].
Figure 2.21 shows the individual masking thresholds for both tonal and nontonal
maskers survived from the decimation.
• STEP 5: Calculation of global masking threshold
The global masking threshold is the combination of individual masking thresholds
and the threshold in quiet. Since the mixture of masking is additive, the global
masking threshold at frequency index i is calculated according to
2 3
ATH.i / X
NTM
LTM Œz.j /;z.i / X
N NM
LNM Œz.j /;z.i /
LG .i / =dB D 10 log10 410 10 C 10 10 C 10 10 5;
j D1 j D1
(2.23)
2.4 Psychoacoustic Model 43

a
120
Tonal components
Masking thresholds for tonal components
100 Non−tonal components
Masking thresholds for non−tonal components

80
Sound pressure level (dB)

60

40

20

−20
0 4000 8000 12000 16000 20000
Frequency (Hz)

b 120
Tonal components
Masking thresholds for tonal components
100 Non−tonal components
Masking thresholds for non−tonal components

80
Sound pressure level (dB)

60

40

20

−20
1 3 5 7 9 11 13 15 17 19 21 23 25
Critical band rate (Bark)

Fig. 2.21 Individual masking thresholds (a) Frequency on linear scale (b) Frequency on Bark
scale

where ATH .i / is the SPL of threshold in quiet at frequency index i , NTM and NNM
are the number of tonal and nontonal maskers, and LTM Œz .j / ; W and LNM Œz .j / ; W
are their corresponding individual masking thresholds.
44 2 Principles of Psychoacoustics

The global masking threshold is denoted by a bold dashed black line in

Fig. 2.22.
• STEP 6: Determination of the MMT
The MMT is derived from the global masking threshold. As mentioned in Step 4,
the global masking threshold LG is computed on only a subset of samples (here
106 samples) over the frequency spectrum, i.e., 1  i  106. Then these spectral
subsamples are mapped onto 32 uniform subbands, as shown in Fig. 2.23. Each
subband contains N=232
D 512=2
32
D 8 samples. Therefore, the minimum masking
level in the nth subband .1  n  32/ is determined by the following expression:

LMin .n/ =dB D min LG .i / ; (2.24)

fid .i/2subband n

where fid .i / is the frequency index corresponding to the i th subsample. After

spreading every LMin .n/ .1  n  32/ over its subband with 8 samples, we get
the MMT LMMT :

LMMT .m/ D LMin .n/ m D Œ8 .n  1/ C 1 W 8n: (2.25)

2.4.1.3 Comparison Between Psychoacoustic Model 1 and Model 2

The general idea of implementation on Psychoacoustic Model 2 is similar to Model

1. However, the concrete operations of calculating MMT in Psychoacoustic Model 2
are quite different from that of Model 1, as depicted in the following steps [78,81]:
• STEP 1: FFT analysis and calculation of complex spectrum
The input to Model 2 is a set of 1,024 samples, twice longer than 512-point frame
in Model 1. Before performing FFT, a Hanning window is applied as well.
• STEP 2: Definition of threshold calculation partitions and spreading function
The notion of “threshold calculation partitions” is a significant difference in Model
2. Instead of identifying the tonal and nontonal maskers in each critical band in
Model 1, Model 2 groups the frequency lines into so-called threshold calculation
partitions. Such partitions are also of nonlinear widths, but with finer frequency
resolution than critical band. Each partition has a width of either one FFT line
(at low frequencies) or 1/3 critical band (at high frequencies), whichever is wider
[78]. According to this criterion, there are 57 partitions at a sampling frequency of
44.1 kHz by calculation. The result complies with Table D.3b in [77].
The spreading function in Model 2 is described by Eq. (2.7) and one specific
spreading function is defined for each partition. Note that 10 log10 SF .d z/ in
Eq. (2.7) is level-independent and thereby suitable for alleviating the computational
burden of convolution in Step 4.