Romain Hennequin

Bertrand David

Roland Badeau

SCORE INFORMED AUDIO SOURCE SEPARATION USING A PARAMETRIC MODEL OF NON-NEGATIVE SPECTROGRAM

Keywords: audio source separation, music information retrieval, machine learning, non-negative matrix factorization

In this paper we present a new technique for monaural source separation in musical mixtures, which uses the knowledge of the musical score. This information is used to initialize an algorithm which computes a parametric decomposition of the spectrogram based on non-negative matrix factorization (NMF). This algorithm provides time-frequency masks which are used to separate the sources with Wiener filtering.

INTRODUCTION

Underdetermined audio source separation has been a major field of research for the past decades. When applied to musical signals, the goal is usually to obtain tracks corresponding to isolated instruments. Blind source separation was first addressed notably with non-negative matrix factorization (NMF) [START_REF] Virtanen | Monaural sound source separation by nonnegative matrix factorization with temporal continuity[END_REF]. One of the main drawbacks of this technique is the difficulty to cluster the factorized elements and associate them with a source. This is certainly a reason why numerous works introduce additional information to improve separation results. Different kinds of information have been considered: in [START_REF] Smaragdis | Supervised and semisupervised separation of sounds from single-channel mixtures[END_REF], the different spectral shapes of each source are learned on isolated sounds and are then used to decompose the mixture. In [START_REF] Parvaix | A watermarking-based method for informed source separation of audio signals with a single sensor[END_REF], source signals are used as a side information in a coder/decoder scheme.

Recently, the use of an aligned MIDI file to guide source separation was addressed in several works. In [START_REF] Woodruff | Remixing stereo music with score-informed source separation[END_REF], stereo source separation based on spatial cues is improved by the knowledge of the score in order to accurately separate time-frequency bins with overlap. In [START_REF] Raphael | A classifier-based approach to score-guided source separation of musical audio[END_REF], the score of the solo helps to separate it from the accompaniment with a classifier approach. In [START_REF] Smaragdis | Separation by humming: User-guided sound extraction from monophonic mixtures[END_REF], P. Smaragdis proposed a probabilistic latent component analysis (PLCA) for isolating sounds in a mixture from the presentation of a humming query. This query mimics the desired target to be extracted and serves as a prior in the PLCA decomposition of the mixture. This approach is applied in [START_REF] Ganseman | Evaluation of a score-informed source separation system[END_REF], replacing the humming query by an aligned MIDI file and a synthesizer: each track of the MIDI file is synthesized and the provided signals are used as priors in the PLCA decomposition of the mixture. In [START_REF] Every | A spectral-filtering approach to music signal separation[END_REF], harmonic filters are generated from the score: in each analysis frame, the fundamental frequency of each active note in the MIDI file is finely assessed from the peaks in the spectrum.

The research leading to this paper was supported by the CONTINT program of the French National Research Agency (ANR), as a part of the DReaM project (ANR-09-CORD-006-03), and by the Quaero Programme, funded by OSEO, French State agency for innovation.

In this paper, we propose a new approach for score informed source separation, based on a parametric decomposition of the power spectrogram of the mixture. The information extracted from the score is used to initialize the algorithm which provides the decomposition. The algorithm then locally optimizes the parameters (notably the fundamental frequency of each atom at each frame). The obtained decomposition provides a time-frequency mask for each source, which permits to separate the sources from the mixture by Wiener filtering.

The paper is organized as follows. An unusual non-negative framework which uses time-dependent frequency templates for decomposing audio spectrograms is first presented in section 2. We then derive an algorithm for score-informed source separation in section 3. A comparative evaluation of the performance of this algorithm is provided in section 4 where PLCA-based algorithms [START_REF] Ganseman | Evaluation of a score-informed source separation system[END_REF] are employed for benchmarking. Finally, we draw some conclusions and outline future work in section 5.

SPECTROGRAM MODEL

The model of power spectrogram is inspired by NMF but uses timedependent parametric atoms (as presented in detail in [START_REF] Hennequin | Time-dependent parametric and harmonic templates in non-negative matrix factorization[END_REF]). These atoms are harmonic and we only consider the separation of a mixture of harmonic (or quasi-harmonic) instruments thus excluding percussive instruments.

Non-negative Matrix Factorization

Given an F × T non-negative matrix V and an integer R such that F R + RT ≪ F T , NMF approximates V by the product V of an F × R non-negative matrix W and an R × T non-negative matrix H:

∀f ∈ {1, . . . F }, t ∈ {1, . . . T } [V] f t ≈ [ V] f t = R r=1 w f r hrt.
(1) When V is the magnitude or power spectrogram of a musical signal, the templates that are redundant in multiple frames are hopefully most of the time harmonic templates corresponding to musical tones. Thus, each column of W should correspond to a note and each row of H is the time activation associated with that note. However, this property is not guaranteed and generally, further constraints are added [START_REF] Bertin | Enforcing harmonicity and smoothness in Bayesian non-negative matrix factorization applied to polyphonic music transcription[END_REF].

The approximation in (1) is generally quantified with an element-wise divergence to be minimized with respect to H and W:

C(W, H) = D(V|| V) = F f =1 T t=1 d([V] f t , [ V] f t ). (2)
In this paper, the general class of β-divergence (see [START_REF] Hennequin | Time-dependent parametric and harmonic templates in non-negative matrix factorization[END_REF] for its expression) is considered and it is particularized to the Kullback-Leibler divergence (case β = 1) in our experiments.

Model of source spectrograms

As stated previously, the model of the power spectrogram (for each source of the mixture) considered in this paper is the parametric model presented in [START_REF] Hennequin | Time-dependent parametric and harmonic templates in non-negative matrix factorization[END_REF]: the spectrogram of a single instrument (a source) is indexed by an integer k and is decomposed with parametric frequency atoms. It thus writes:

[ Vk ] f t = R r=1 w f krt 0 kf r h krt . (3)
In our model, the parameterization of the atoms is derived for each frame t: the fundamental frequency of a selected frequency template is refined for each frame which is emphasized by the t dependancy of the corresponding symbol f krt 0 . The whole atom thus writes w f krt 0 kf r where r denotes the index of the atom and f the frequency bin.

Atoms are synthesized in the following way:

w f krt 0 kf r = n h p=1 a kp g(f -pf krt 0 ). (4)
Figure 1 represents such an harmonic atom. The function g corresponds to a single spectral component. Thus g is the squared modulus of the Fourier transform of the analysis window (used to compute the spectrogram). a kp is the amplitude of the pth harmonic for every atom of source k: in order to limit the number of parameters in the model, we use the same set of harmonic amplitudes for every atom of a source. n h is the number of harmonics. Each source has its own set of harmonic amplitudes and its own time-varying fundamental frequency. The time-varying fundamental frequency f krt 0 will be finely tuned by the algorithm described in section 3.2.

In this model, the following assumptions are made:

• one assumes that the harmonic part is stationary within a single frame (then the Fourier transform of an harmonic is the Fourier transform of the analysis window centered around the frequency of the harmonic);

• interferences between harmonics are assumed to be weak (this assumption is valid when the fundamental frequency of the atom is not too low);

• one assumes that interferences between peaks of negative frequencies (not taken into account) and positive frequencies are weak;

• one assumes that the frequency aliasing introduced by the signal sampling is weak (which allows using the analytical expression of the continuous Fourier transform of the analysis window).

As the cost function (2) to be minimized is multimodal with respect to the fundamental frequency of each atom, it is impossible to perform a global minimization of this function. It is thus necessary to introduce numerous atoms: in the proposed system, we chose to have an harmonic atom for each semitone of the chromatic scale (i.e. an atom per MIDI note). Atoms and activations can thus be indexed with a MIDI note number and the activations of a source can then be viewed as a piano roll.

Model of the mixture spectrogram

The power spectrogram of the mixture is supposed to be the sum of the spectrograms of the sources (this is a common assumption with non-negative decompositions such as NMF which is generally explained by the approximate statistical independence of the sources). Thus the model of the mixture spectrogram is:

V mix ≈ Vmix = K k=1
Vk , where Vk is the parametric spectrogram of source k (there are K sources) given in equation (3).

SCORE INFORMED SOURCE SEPARATION

The model of parametric spectrogram introduced in the previous sections is used to decompose the mixture spectrogram by means of a multiplicative descent algorithm which aims at minimizing (2), initialized with the information of the score. The score (in this paper a temporally aligned MIDI file) presumably permits to initialize the decomposition in a neighborhood of an optimal decomposition. In this paper, the MIDI file is supposed perfectly aligned with the mixture signal and thus we do not deal with the problem of MIDI alignment which can be done automatically [START_REF] Dannenberg | Polyphonic audio matching for score following and intelligent audio editors[END_REF].

Initialization with the musical score

The MIDI file provides a piano roll for each source. Each of these piano rolls has a direct link with the activations of the atoms of the corresponding sources. Thus, they can be used as an initializing mask for activations: while a note is active in the piano roll, the activation of the corresponding harmonic atom is set to 1. For all other instants, this activation is set to 0. Since we use a multiplicative algorithm (presented in section 3.2), coefficients initialized to 0 will remain 0 over iterations. Thus, it is better to slightly enlarge initialization to 1 before the beginning of a note and after the end of a note: this is necessary to avoid possible alignment errors and to take the possibly slow release of a note into account. An illustration of initializing activation masks is given in figure 2.

Algorithm

The algorithm resembles those used for standard NMF: the decomposition is obtained by minimizing a β-divergence between V mix et Fig. 2. Activation mask for the 3 instruments of a piece: activations of each source are initialized with a binary piano roll extracted from the corresponding MIDI track.

Vmix with respect to the parameter of the model i.e. for each source:

• the fundamental frequency of each atom r at each instant t: f krt 0 , • the set of harmonic amplitudes a kp , • the activation of each note at each instant h krt .

The minimization is done with multiplicative update rules which are successively applied to each of the previous parameters. These rules particularly ensure that the parameters remain non-negative and become constant if the partial derivative of the cost function with respect to the considered parameter goes to zero.

These update rules can be straightforwardly derived from the mono-instrument (a single source) case presented in [START_REF] Hennequin | Time-dependent parametric and harmonic templates in non-negative matrix factorization[END_REF]. Thus we give them without providing the whole derivation. The update rules of parameters for each source k are given by:

f krt 0 ← f krt 0 F krt G krt , h krt ← h krt M krt P krt , a kp ← a kp N kp Q kp ,
with:

G krt = F f =1 n h p=1 a kp pP (f -pf krt 0 )( V mix f t ) β-2 (f V mix f t + pf krt 0 V mix f t ), F krt = F f =1 n h p=1 a kp pP (f -pf krt 0 )( V mix f t ) β-2 (pf krt 0 V mix f t + f V mix f t ), P krt = F f =1 w f krt 0 kf r ( V mix f t ) β-1 , M krt = F f =1 w f krt 0 kf r ( V mix f t ) β-2 V mix f t , Q kp = F f =1 T t=1 R r=1 g(f -pf krt 0 )h krt ( V mix f t ) β-1 , N kp = F f =1 T t=1 R r=1 g(f -pf krt 0 )h krt ( V mix f t ) β-2 V mix f t ,
where P is a non-negative function (see [START_REF] Hennequin | Time-dependent parametric and harmonic templates in non-negative matrix factorization[END_REF]) defined from the derivative of the Fourier transform of the analysis frame g as P (f ) = -g ′ (f ) f .

RESULTS

The performance of the algorithm is assessed by means of experiments done on different synthetic databases, aligned on MIDI files.

It is compared to that of an algorithm based on PLCA [START_REF] Ganseman | Evaluation of a score-informed source separation system[END_REF].

Description of the database

To our knowledge, no publicly available database exists which provides real-recorded mixtures of musical signals, real-recorded separated source signals and the corresponding aligned MIDI files. We thus designed an home-made database aiming at realistically rendering important characteristics of real musical streams, as the possibly important overlap of the sources in the time-frequency domain (in opposition to the randomly built database proposed in [START_REF] Ganseman | Evaluation of a score-informed source separation system[END_REF]), while providing separated tracks. For each piece, the tracks are summed to obtain the whole mix. The database is obtained from 12 MIDI files of royalty-free string quartets of Bach, Beethoven and Boccherini. These MIDI files are processed to synthesize wave files with 2 different methods. The first technique relies on isolated sounds recorded from real musical instruments with three levels of velocity. Their onset is synchronized with the Note On message and the sounds are then faded out at Note Off. The instrument list includes violin, viola and cello. This first technique will be referred to as "M1".

The second technique is based on the TiMidity1 software fed with a common soundfont set (Crisis General Midi 3.02 ). The latter technique will be referred to as "M2". M2 includes a reverberation effect present in the soundfont whereas M1 does not.

It is worth noting that in this paper, we do not test the degradation of the performance of the algorithm with respect to alignment errors. The wave files of the database are downsampled to 11025Hz. The mixture signals are monophonic. The database (MIDI files, separated tracks and mixtures) is available at http://perso.enst.fr/ hennequi/database.zip.

Experiment

The database provides the monophonic mixture and the aligned MIDI file for each piece. Only the 30 first seconds of each piece are processed. Mixture signals are separated using the information of the MIDI file with two algorithms:

• the algorithm that we presented,

• the algorithm based on PLCA presented in [START_REF] Ganseman | Evaluation of a score-informed source separation system[END_REF].

The PLCA-based algorithm requires a training stage where the MIDI tracks of the different files are synthesized and used as priors in the decomposition of the mixture. In order to make the comparison fair between algorithms, the synthesizer used to generate the training separated tracks should differ from the one we use to generate mixture signals. As sounds in our database are synthesized with two different methods, we use one for test and the other for training and vice versa. We also provide the results for the PLCA-based algorithm with the true sources used as priors (same training signals and test signals): they can be thought as an upper bound for the PLCA-algorithm performance, and also as a high-rank reference.

For both algorithms, the mixture signals (and the separated signal tracks for the PLCA-based algorithm) were transformed into spectrograms using a short-time Fourier transform with 92ms-long Hann window, and 75% overlap.

The database provides the original separated sources, which are used to evaluate and compare the performance of both algorithms with the BSS EVAL toolbox [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF]: performance is assessed with signal to interference ratio (SIR), signal to artifact ratio (SAR) and signal to distortion ratio (SDR), all defined in [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF]. The parameters of the PLCA-based algorithm (number of atoms and prior weights) were optimized (on a single piece) to give the best results possible.

Results

Results are given in tables 1 and 2. In both tables, the first row corresponds to the ratios obtained with our algorithm, the second row corresponds to the ratios obtained with the PLCA-based algorithm using different methods to synthesize training signals and test signals, and the third row corresponds to the ratios obtained with the PLCA-based algorithm using the true separated signals as training signals. Table 1 corresponds to the experiment using M1 as synthesizing method for the test signals and table 2 corresponds to the experiment using M2 as synthesizing method for the test signals.

The given ratios are averaged over all pieces and all sources.

In table 1, the values of the ratios show that our algorithm outperforms the PLCA-based algorithm by more than 1dB in SAR and SDR and by 6dB in SIR, on our experimental database. Moreover, the ratios are not so far from the ones provided by the oracle-like PLCA-based algorithm. In table 2, the performance importantly deteriorates for both algorithms, probably due to the reverberation produced by M2. Our algorithm has slightly weaker results than the PLCA-based one.

It is worth noting that:

• the isolated samples of violin and viola from M1 contain important vibrato although there is only very slight vibrato in the sounds from M2. As our model explicitly takes vibrato into account, this phenomenon is more accurately modeled with our algorithm which probably partially explains the results of table 1;

• our approach is less supervised than the PLCA approach since it does not require to synthesize the tracks. The PLCA performance strongly depends on the synthesizer used to generate the training tracks: for instance the signals synthesized with M2 do not have the same decay times than the ones from M1 and have reverberation, probably resulting in incorrect activation priors in the PLCA. Thus, our method performs similarly to the PLCA-based method but has the strength of being less supervised.

SIR SAR SDR

CONCLUSION

In this paper, we presented a new method to perform score informed source separation on a musical mixture. This method is based on a parametric model of non-negative spectrogram which uses harmonic atoms. We designed an evaluation database and our method reached a performance comparable to that of a PLCA-based algorithm while not requiring any priors.

As our model is limited to harmonic instruments, future works will include percussive instruments. Moreover, the spectrogram model can be further refined by, for instance, including other timbral parameters such as those proposed in [START_REF] Hennequin | NMF with time-frequency activations to model nonstationary audio events[END_REF]. The supervised learning of harmonic templates should also be explored.

Fig. 1 .

 1 Fig. 1. Parametric atom defined in (4).

Table 1 .

 1 Results on the signal synthesized with M1.

	Parametric decomposition algorithm	20.2	7.7	7.2
	PLCA-based algo (training set: M2)	14.2	6.3	4.8
	PLCA-based algo (training set: M1)	22.8 11.5	11.1
		SIR SAR SDR
	Parametric decomposition	12.6	3.3	2.1
	PLCA-based algo (training set: M1)	12.4	4.5	3.1
	PLCA-based algo (training set: M2)	20.1 10.6	10.1

Table 2 .

 2 Results on the signal synthesized with M2.

http://timidity.sourceforge.net/

http://www.bismutnetwork.com/10Music/Crisis/Soundfont3.0.php


Holger Kirchhoff

Roland Badeau

Simon Dixon

TOWARDS COMPLEX MATRIX DECOMPOSITION OF SPECTROGRAMS BASED ON THE RELATIVE PHASE OFFSETS OF HARMONIC SOUNDS

Keywords: harmonic signals, relative phase offsets of partials, complex matrix decomposition, nonnegative matrix factorisation

In this paper we study the relative phase offsets between partials in the sustained part of harmonic sounds and investigate their suitability for complex matrix decomposition of spectrograms. We formally introduce this property in a sinusoidal model and visualise the phase relations of a musical instrument. A model of complex matrix decomposition in the time-frequency domain is derived and equations for the estimation of the model parameters are provided in the monophonic case. We illustrate the model with the analysis of a monophonic saxophone signal. The results suggest that the phase offset is able to capture inherent time-invariant phase properties of harmonic sounds and outline its potential use for complex matrix decomposition.

INTRODUCTION

Spectrogram factorisation techniques -particularly nonnegative matrix factorisation (NMF) [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF] -have proven to be useful methods for the analysis of instrument sounds and have been successfully applied for tasks such as music transcription [START_REF] Smaragdis | Non-negative matrix factorization for polyphonic music transcription[END_REF], source separation [START_REF] Virtanen | Monaural sound source separation by nonnegative matrix factorization with temporal continuity and sparseness criteria[END_REF] and instrument recognition [START_REF] Heittola | Musical instrument recognition in polyphonic audio using source-filter model for sound separation[END_REF]. Most spectrogram factorisation techniques rely on the assumption that magnitude spectra of sound mixtures can be approximated by the superposition of the magnitude spectra of the sound sources. Although this assumption provides reasonable analysis results in practice, the linearity only holds for the complex coefficients of the short-time Fourier transform (STFT).

Phase information is often discarded in the analysis stage, not only because the human auditory system is considered insensitive to absolute phase shifts of harmonic partials [START_REF] Risset | Exploration of timbre by analysis and synthesis[END_REF], but also because the magnitude spectrogram is often considered more intuitive and easier to model. For all applications in which sounds have to be synthesised from a time-frequency representation, however, the correct estimation of phase values is crucial in order to avoid artefacts due to phase-incoherent overlap of consecutive time frames. For the task of instrument separation, for example, the phase information for the synthesis of each sound source either has to be estimated from the magnitude spectrogram [START_REF] Griffin | Signal reconstruction from shorttime Fourier transform magnitude[END_REF], or the phases of the original mixture have to be employed for each source [START_REF] Virtanen | Monaural sound source separation by nonnegative matrix factorization with temporal continuity and sparseness criteria[END_REF]. Using the mixture phases can lead to reasonable results when the number of sources is small and when most time-frequency bins are mainly influenced by a single source. For higher numbers of sources and larger time-frequency overlap, however, it can lead to cross-talk artefacts.

The STFT analysis phase of harmonic partials is not constant over time. Therefore the instantaneous phase cannot be simply embedded in a matrix factorisation framework, where the basis functions include prototypical properties of the underlying instrument spectra that vary little over time. Nevertheless, several approaches have been proposed to consider phase information in matrix factorisation frameworks. Parry and Essa [START_REF] Parry | Phase-aware non-negative spectrogram factorization[END_REF] propose a phase-aware nonnegative matrix factorisation. The authors model the STFT bins as complex random variables, thereby assuming the phase to be uniformly distributed, and derive iterative update rules based on this assumption. The update rules, however, still estimate the matrices based on the magnitude spectrogram only. In a similar way, Févotte et al. [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis[END_REF] show that Itakura-Saito NMF is equivalent to a maximumlikelihood parameter estimation of a sum of complex Gaussian components. The Gaussian components have zero mean and a diagonal covariance matrix, which also assumes a uniformly distributed phase. An attempt to explicitly estimate the phase values of the individual sources was made by Kameoka et al. [START_REF] Kameoka | Complex NMF: A new sparse representation for acoustic signals[END_REF]. Their complex NMF algorithm combines the outer product of each NMF basis function and gain vector with a phase spectrogram with the same dimensions as the original spectrogram. In [START_REF] King | Single-channel source separation using complex matrix factorization[END_REF], complex NMF was shown to outperform NMF for speech separation. Complex NMF is not a complex matrix factorisation technique, but a combination of NMF with time-frequency phase estimates. The algorithm is heavily overparameterised and it can be shown that an initialisation with the original phase values leaves the phase parameters unaltered (up to ±π). Lastly, a high resolution NMF framework was introduced in [START_REF] Badeau | Gaussian modeling of mixtures of nonstationary signals in the time-frequency domain (HR-NMF)[END_REF][START_REF] Badeau | High resolution NMF for modeling mixtures of non-stationary signals in the time-frequency domain[END_REF], in order to model both the magnitude and phase of complex or real-valued time-frequency representations. However this framework does not take the phase relations of the partials into account.

In this paper, we exploit the relative phase offsets between partials in the sustained part of the sounds of harmonic instruments as a step towards complex matrix decomposition. The concept will be reviewed and illustrated in Section 2, where we also present a mathematical formulation. In Section 3 we derive the model for complex matrix decomposition and present the parameter estimation equations for the monophonic case. An example analysis of a monophonic signal is provided in Section 4, and we conclude this work in Section 5.

PHASE RELATIONS OF HARMONIC PARTIALS

Concept

Pitched musical instruments generally produce harmonic sounds which can be represented by a superposition of P sinusoids at in-teger multiples of a fundamental frequency. Each harmonic partial can be described by its angular frequency ωp > 0, its amplitude ap ≥ 0 and an absolute phase shift φp ∈ [-π, π): ∀t ∈ Z, s(t) = P p=1 ape j [ωpt+φp] .

(1)

For strictly harmonic sounds, the frequency of each harmonic is given as the p-th multiple of the fundamental frequency: ωp = p ω1.

Complex exponentials are used here to reflect the fact that we only consider the baseband of the DFT in our model.

In this paper we are interested in the relation between the absolute phase shifts of the harmonic partials. To capture this relation, we express the phase shift of each partial in relation to the instantaneous phase of the fundamental frequency ω1:

s(t) = P p=1 ape j[p•θ(t)+∆φp] , (2)
where θ(t) = ω1t + φ1 denotes the instantaneous phase of the fundamental and ∆φp = φpp φ1 represents the phase offset between the p-th partial and the fundamental (with ∆φ1 = 0).

Figure 1 shows a graphical illustration of the parameters in Eq. ( 2). The upper part (Fig. 1a) displays the waveform of the first three partials of a harmonic sound, and the lower part (Fig. 1b) the instantaneous phases. The phase offsets ∆φp correspond to the instantaneous phases of the partials at the time where θ(t) = 0. Modifying ∆φp translates the p-th partial relative to the fundamental along the time axis. Since a translation by ∆φp is equivalent to a translation by ∆φp + c • 2π with c ∈ Z, ∆φp is uniquely defined in the range [-π, π). Given all phase offsets ∆φp of the partials, the instantaneous phase of each partial can be computed at any given time t0 based on the instantaneous phase of the fundamental θ(t0) at that time. Note that even though we can only measure the wrapped phase of θ(t) (i.e. in the interval [-π, π)), the correct wrapped phase of each partial can still be calculated.

Example

To illustrate the phase relations, we display the phase offsets ∆φp of the partials with indices p ∈ {2, 3, 4} in a monophonic saxophone recording of "Summertime" by G. Gershwin in Figure 2. The first four bars of this small excerpt are displayed in Fig. 2a, and Fig. displays the fundamental frequency of the saxophone performance measured by the YIN algorithm [START_REF] De Cheveigné | YIN, a fundamental frequency estimator for speech and music[END_REF]. In Fig. 2c, the partial offsets ∆φp are plotted over time. Phase offsets were obtained from the STFT and computed as the wrapped difference between the measured instantaneous phases of each partial and p times the instantaneous phase of the fundamental. It can be seen that the partial offsets exhibit little variation during the steady state of each note -which is not surprising given the fact that the sound is harmonic. In addition to that, however, the same phase offsets occur at different notes of the same pitch. The area shaded in dark grey highlights all renditions of the note E4 and the light grey area highlights all occurrences of the note D4. These observations make this property suitable for the use in a complex matrix decomposition framework as we will illustrate in the next section. It should be noted that the relative phase offsets can only be defined if the partial frequencies are in a strictly harmonic relation. For instruments with slightly inharmonic frequency relations -such as the piano -a constant phase offset does not exist.

PARAMETER ESTIMATION

Frequency domain model

We aim at estimating the parameters of the model in Eq. ( 2) from the STFT which is given by The STFT of the signal s(t) from Eq. ( 2) is given by

X(n, k) = K-1 t=-K x(t + n•m) • h(t) • e -jΩ k t , (3)
S(n, k) = P p=1 apH(Ω k -pω1)e j[pΘ(n)+∆φp] . (4)
In this equation, H(Ω) = K-1 t=-K h(t) • e -jΩt denotes the Fourier spectrum of the window function h(t) and Θ(n) = θ(n • m).

To simplify the monophonic model in Eq. ( 4), we assume that each partial can be represented by the main lobe of the window function only. This assumption holds fairly well if the side lobe attenuation of the window spectrum H(Ω) is sufficiently high and if the frequency resolution of the STFT is high enough so that the main lobes of adjacent partials do not overlap. We denote the partial index belonging to frequency bin k by p k . We set p k = 0 for all k that lie outside the main lobes of the partials, and set a0 = 0. This allows us to drop the sum in Eq. ( 4):

S ′ (n, k) = ap k H(Ω k -p k ω1)e j[p k Θ(n)+∆φp k ] . (5)
We additionally introduce a real time-varying gain factor g(n) > 0 that enables a uniform scaling of the magnitudes in order to accommodate loudness variations (similar to the gains in NMF):

B(n, k) = g(n) • S ′ (n, k). (6)
Scaling ambiguities between g(n) and ap can be resolved by normalising ap. Finally, the model can be extended to incorporate multiple harmonic sounds. We denote the index of each harmonic sound by r and append it to the quantities in Eq. ( 6):

V (n, k) = R r=1 gr(n) • S ′ r (n, k) (7)
By substituting Eq. ( 5) into Eq. ( 7), we finally obtain:

V (n, k) = R r=1 wr(k) • hr(n, k) (8)
where wr(k) = ap k,r H(Ω kp k,r ω1,r)e j∆φp k,r and hr(n, k) = gr(n)e jp k,r Θr (n) . The term wr(k) is not time-dependent and is therefore referred to as a complex basis function. Accordingly, the term hr(n, k) is referred to as a complex activation. Note that hr(n, k) is a 2-dimensional function. Eq. ( 8) is therefore not a complex matrix factorisation, but a decomposition of a complex spectrogram V (n, k) into a matrix of complex basis functions wr(k), a matrix of real-valued gain factors gr(n) and a matrix of real-valued instantaneous phases of the fundamentals Θr(n).

In this paper we will not investigate the case of multiple concurrent sounds. Our goal is to prove that phase offsets between partials are a viable concept for sound analysis purposes. The model parameters will thus be estimated in the monophonic case of Eq. ( 6) only.

Parameter estimation

The parameters in Eq. ( 6) can be estimated by minimizing the error between the original complex spectrogram B(n, k) and the model approximation B(n, k) for all n ∈ [1 . . . N ] and k ∈ [1 . . . K] with N > 0 and K > 0. We choose to minimise the following cost function: [START_REF] King | Single-channel source separation using complex matrix factorization[END_REF] where ∠B(n, k) denotes the argument of the complex number B(n, k). The term q(n, k) ∈ Z stems from the fact that the logarithm of a complex number has an infinite number of solutions which are obtained by adding integer multiples of 2π to the imaginary part of the solution [START_REF] Sarason | Complex Function Theory[END_REF]. The integer q(n, k) is here treated as an additional parameter that has to be estimated. In Eq. ( 10), H(Ω) is assumed positive, since we only consider the main lobe of the window function. The model parameters are estimated by means of a coordinate descent (J is minimized w.r.t. each parameter):

J = N n=1 K k=1 ln(B(n, k)) -ln( B(n, k)) 2 (9) = N n=1 K k=1 ln |B(n, k)| g(n) • ap k H(Ω k -p k ω1) 2 + [∠B(n, k) -∆φp k -p k Θ(n) + 2πq(n, k)] 2
g(n) = K k=1 |B(n, k)| ap k H(Ω k -p k ω1) 1 K (11) ap =   N n=1 {k|p k =p} |B(n, k)| g(n)H(Ω k -pω1)   1 N •#{k|p k =p} (12) Θ(n) = K k=1 p k [∠B(n, k) -∆φp k + 2πq(n, k)] K k=1 p 2 k (13) ∆φp = N n=1 {k|p k =p} ∠B(n, k) -pΘ(n) + 2πq(n, k) N • #{k|p k = p} (14) q(n, k) = round -1 2π [∠B(n, k) -∆φp k -p k Θ(n)] (15)
In these equations, the expression {k|p k = p} denotes the set of frequency indices k at which p k = p and the operator #{. . .} denotes the cardinality of the set. The function round() denotes the rounding of a real number to the nearest integer.

ANALYSIS OF AN EXAMPLE SIGNAL

In this section we apply the estimation method to an example signal and illustrate how it can be used for a transcription task. The signal is the same monophonic saxophone recording of "Summertime" that we used to illustrate the phase relations in Figure 2. It has a sample rate of 44.1 kHz and we use the first eight bars of the recording. The recording is split into two parts. The first part contains the first four bars (cf. Fig. 2a) and is used as training material in which prototypical partial amplitudes ap and phase offsets ∆φp are learned for different pitches ω1. The spectrogram with K = 2049 frequency bins and N = 5380 time frames is manually segmented in time into the different notes. All spectrogram parts with the same nominal pitch are concatenated, the fundamental frequency is estimated by employing the YIN algorithm and ω1 is computed as the average across all frames for each nominal pitch. g(n) is estimated from the original spectrogram by taking the mean of the magnitudes in each time frame. In order to compute ap we alternately apply Eq. ( 12) and (11) for 10 iterations. For the computation of ∆φp, Θ(n) is initialised by the instantaneous phase value of the frequency bin corresponding to the fundamental in each frame. An initial estimate for ∆φp is obtained by replacing the terms in the summation in the numerator of Eq. ( 14) by wrap(∠B(n, k) -pΘ(n)), where wrap(α) calculates the principal argument of α. q(n, k) is computed according to Eq. ( 15) and we iterate over Eq. ( 13)-( 15) until q(n, k) converges. ∆φp is eventually given by result of Eq. ( 14) in the last iteration.

The second part of the recording contains the remaining four bars (cf. Fig. 3a) and is used as test material. The learned prototype amplitudes ap and ∆φp for all pitches ω1 that occurred in the first part are employed to estimate g(n) and Θ(n) in the following way. First, Θ(n) is initialised with the instantaneous phase values at the frequency bins corresponding to ω1. Then q(n, k) is estimated according to Eq. ( 15). Finally, g(n) and Θ(n) are estimated according to Eqs. ( 11) and [START_REF] De Cheveigné | YIN, a fundamental frequency estimator for speech and music[END_REF].

Active pitches can be estimated from both g(n) and Θ(n). While for g(n) this is obvious -high values indicate activity, low values indicate inactivity -, the instantaneous phase Θ(n) of the fundamental can also be used as an activity detector. We here use a measure inspired by the phase-based onset detection function described in [START_REF] Bello | A tutorial on onset detection in music signals[END_REF]. The measure is based on the unwrapped phase, which can be assumed to be linear when the note is active and non-linear when the note is inactive. We will denote the unwrapped phase of Θ(n) by Θu(n). The second phase difference can be used as a measure of phase-linearity. It is given by

∆Θu(n) = Θu(n) -2Θu(n -1) + Θu(n -2). (16)
If the unwrapped phase is strictly linear, ∆Θu(n) will be close to zero, if it is non-linear ∆Θu(n) is likely to take on values with larger magnitudes. Additionally, ∆Θu(n) is likely to take on low values in several consecutive active frames and more random values in consecutive inactive frames. We therefore compute the mean square of ∆Θu(n) over a sliding window as

σ(n) = 1 M M 2 -1 n ′ =-M 2 ∆Θ 2 u (n + n ′ ), (17)
and define the phase-based activity measure as f (n) = -ln(σ(n)) .

In our simulations a window length of 50 ms (M = 37) was used. The results of the estimation are displayed in Figure 3. In Fig. 3b, the measured fundamental frequencies of the four bar excerpt are shown. Fig. 3c shows the gains g(n) and Fig. 3d the results for the phase-based activity measure f (n). The gains clearly show the activity of the different pitches and are very much reminiscent of activity measurements in NMF analyses. The results of the phase-based activity measure also reveal the active pitches very well, which confirms that the phase relations between harmonic partials can actually be used to characterise pitches of certain instruments and distinguish between them. Note that the pitch E3 does not occur in bars 5-8, and that the note B3, the last note in Fig. 3a, is missing because it did not occur in bars 1-4.

CONCLUSIONS

In this paper we have investigated the relative phase relations within the sustained part of harmonic sounds and their potential use for complex matrix decomposition. The phase relations between harmonic partials have been expressed as relative phase offsets of the partials w.r.t. the fundamental. Equations for the estimation of the model parameters have been presented based on a complex logarithmic cost function between the original spectrogram and the model approximation. With the analysis of an example signal, we demonstrated the potential of the phase coupling property to capture inherent time-independent phase characteristics of harmonic sounds.

In future work the method should be extended to deal with mixtures of harmonic sounds in order to obtain a complex matrix decomposition that can be used to unmix spectral components in the complex domain. A formulation of such a complex matrix decomposition framework has been provided in Section 3.1. For the monophonic case, the complex logarithmic cost function proved to be useful, not only because logarithmic amplitudes better match the human perception than linear amplitudes, but also because it separates the modulus and argument of the model, which allowed us to treat them separately. In the polyphonic case however, a complex matrix decomposition framework would need to deal with magnitudes and phases jointly, since the sum of two complex time-frequency components depends on both their modulus and phase.

  (a) Waveform of the first three harmonics. (b) Instantaneous phases of the first three harmonics.

Fig. 1 :

 1 Fig. 1: Illustration of the model parameters.

  2b

  First four bars of "Summertime" by G. Gershwin. (b) Fundamental frequencies of saxophone performance. (c) Phase differences ∆φp over time.

Fig. 2 :

 2 Fig. 2: Visualisation of the phase relations of a musical instrument.

  where x(t) is the signal under analysis and n and k represent the time frame and frequency index, respectively. h(t) denotes the analysis window of time support [-K . . . K -1]. The distance between consecutive audio frames in samples (hopsize) is denoted by m.Ω k = 2πkN is the normalised angular frequency of the k-th frequency index.

  Bars 5-8 of "Summertime" by G. Gershwin. (b) Fundamental frequencies for bars 5-8. (c) Gains g(n) for bars 5-8 for each pitch from bars 1-4. (d) Phase activations f (n) for bars 5-8 for each pitch from bars 1-4.

Fig. 3 :

 3 Fig. 3: Example analysis of a monophonic saxophone example.

This work was funded by a Queen Mary University of London CDTA studentship.

† This work was undertaken while Roland Badeau was visiting the Centre for Digital Music, partly funded by EPSRC Platform Grant EP/K009559/1.


Le Carrou

François ; Gautier

Roland Badeau
email: roland.badeau@enst.fr

THEORETICAL AND EXPERIMENTAL INVESTIGATIONS OF HARP'S SYMPATHETIC MODES PACS: 43

The harp is composed of a soundboard, a cavity with soundholes and 47 strings. When one string is plucked, other strings are excited via the soundboard. This phenomenon, called sympathetic vibrations, was investigated in a recent paper [1], both theoretically and experimentally using an ersatz of string instrument, composed of 2 strings connected to a beam clamped at both ends. In order to extend the analysis to the case of the harp, an analytical model of a set of 35 strings coupled to a beam, whose characteristics are designed in order to be equivalent to the harp's soundboard, has been developed. The identification of sympathetic modes from the analytical modal basis of the beam-strings assembly is performed using a criterion, the Kinetic Energy Ratio. For some particular tuning conditions, these modes are present in the sound measured on a Camac concert harp. In a same partial, several components, whose frequencies can be really very close one to another are present but cannot be separated using classical Fourier transform. The identification of these components associated to sympathetic modes is achieved using specific signal processing techniques, called High Resolution methods. Identified experimental sympathetic modes are found in good agreement with those theoretically obtained.

INTRODUCTION

The harp is a complex vibroacoustic system, composed of a soundboard, a cavity with five soundholes and 47 strings. When one string is plucked, other strings are excited via the soundboard of the instrument. This phenomenon, called sympathetic vibrations, involves many couplings between the soundboard and the strings. Because of these couplings, the assembly composed by the strings and the soundboard has several modes whose eigenfrequencies are very close one to another. As a consequence, each partial of the musical sound contains multiple components which cannot generally be separated by standard Fourier analysis. The existence of such multiple components in one partial leads to beatings in the musical sound or to non-uniform decrease of the partials [START_REF] Le Carrou | Analyse des modes de cordes couplées d'une harpe par une méthode à haute resolution (Analysis of coupled modes of a harp by a high resolution method)[END_REF]. The sensation caused by this phenomenon of sympathetic vibrations is sometimes described by the harp player as unpleasant and therefore, the harp maker would like to limit them. The aim of this paper is to investigate this phenomenon both experimentally and theoretically. In the first part, an adequate experimental setup and an analysis method are carried out to identify frequency components present in one partial of the musical sound when one string is plucked. In the second part, the sympathetic vibrations are investigated by the use of a theoretical model of the concert harp. In the final part, a comparison between theoretical and experimental results allows us to explain the origin of the sympathetic phenomenon in the instrument.

EXPERIMENTAL INVESTIGATION OF THE HARP'S SYMPATHETIC MODES

Experimental setup

In order to experimentally highlight the presence of sympathetic modes in the concert harp, the following experimental protocol has been set up: the harp is plucked whereas an accelerometer measures the vibratory signal on the soundboard between the D 3 -string and the C 3 -string, respectively labeled 30 and 31. The played string is string 31 whereas the other strings are either damped or stopped during oscillations. The 4 studied experimental configurations are the following: (1) all strings free to vibrate, (2) all strings damped except string 31, (3) strings 31 and 38 stopped during oscillations and (4) string 31 and 42 stopped during oscillations. All measured signals are sampled to 4096 Hz and last 8 seconds. The string is stopped a few seconds after the plucking of string 31 by the harp player.

Fourier analysis

The analysis of the spectral content of the signal's first partial is performed using an appropriate filtering. The magnitude of its Fourier transform is presented in Figure 1, showing that several sinusoidal components are present in the signal. Indeed, since the system is in free oscillations, the response of the instrument to the plucking action, measured by the accelerometer, corresponds to the superposition of several modes whose frequencies are very close one to another. In the spectrum zoomed in the vicinity of the first partial, peaks of every component cannot be clearly separated. This shows the limit of the Fourier analysis. The presence of several vibration modes having close frequencies is a common characteristic for numerous free stringed instruments.

Method for extracting parameters

Various methods, applied to musical instruments, allow the identification of spectral components, like methods based on the Hilbert transform [START_REF] David | Acoustical characterization of vibrating structures in near vacuum. Estimation methods for eigenfrequencies and damping factors[END_REF] or high resolution methods [START_REF] David | Acoustical characterization of vibrating structures in near vacuum. Estimation methods for eigenfrequencies and damping factors[END_REF], [START_REF] Laroche | The use of the matrix pencil method for the spectrum analysis of musical signals[END_REF], [START_REF] Badeau | High resolution methods for estimating and tracking modulated sinusoids[END_REF]. For the identification of modes with close frequencies in the first partial, we choose to use a high resolution method: the ESPRIT method (Estimation of Signal Parameters via Rotational Invariance Techniques) [START_REF] Roy | ESPRIT-A subspace rotation approach to estimation of parameters of cisoids in noise[END_REF]. In this method, the signal model [ ]

x t is supposed to be the sum of complex exponentials:

[ ] [ ] (2 ) 1 , 1,
k k k K t i f t k k x t a e e t N δ π φ + = = ∈ ∑ (Eq. 1)
where each frequency

1 1 , 2 2 k f     ∈ -        
is associated to a real magnitude ( )

0 k a > , a phase [ ] ( ) , k φ π π ∈ -
and a damping or amplification factor ( ) k δ ∈ ℝ . The whole number K is the number of complex exponentials, also called model order, and N is the number of the signal's samples. The ESPRIT method consists in analyzing the signal measured by using a Kdimensional signal subspace and n -length data vectors with n K > . It can be shown that parameter n is ideal when

3 n N = or 2 3 n N = [5]
. The different steps for determining all parameters of the signal model (Eq. 1) are detailed in [START_REF] Badeau | High resolution methods for estimating and tracking modulated sinusoids[END_REF] and [START_REF] Carrou | Vibroacoustics of the concert harp[END_REF].

The main difficulty of the method consists in the evaluation of the number K of components present in the signal. The technique commonly used is the over-estimation of this number and the discrimination of spurious results by mean of an indicator such as the components energy or the error between the measured signal and the model. Other more effective methods exist for estimating the model order K such as the ESTER method (ESTimation Error) [START_REF] Badeau | High resolution methods for estimating and tracking modulated sinusoids[END_REF], [START_REF] Badeau | A new perturbation analysis for signal enumeration in rotational invariance techniques[END_REF] also used in the study, consisting in the computation of an inverse error function.

In order to minimize the computation time and increase the results accuracy [START_REF] Badeau | High resolution methods for estimating and tracking modulated sinusoids[END_REF], the ESPRIT method implementation is carried out according to the following procedure: after the centering of the studied partial around the null frequency, a finite impulse response (FIR) filter selects the frequency range containing the partial to analyze. The filter is chosen with a linear-phase to keep the signal form. The filter is known to have a finite transitory response which corresponds to the length of its impulse response. We can avoid this problem by cancelling the filtered signal points belonging to its transitory phase [START_REF] David | Acoustical characterization of vibrating structures in near vacuum. Estimation methods for eigenfrequencies and damping factors[END_REF]. The filtered and centered signal is highly decimated to limit the computational time of the ESPRIT method. After the estimation of the model order by the ESTER method, the ESPRIT algorithm is then applied. The final validation of the model order is performed using a comparison between measured and synthesized signals.

Results

The implemented ESPRIT method, as previously explained, is applied to the signal measured on the concert harp in its final part, between 4s and 8s. The estimated components found in the first partial are gathered in Table I. The model order stretches from 1, for configurations (3) and (4), to 4, for configuration (1). In Table I, the components are classified according to the energy type, from top to bottom, and are aligned in order to facilitate the reading of the missing components.

Thanks to the developed experimental protocol and the identification method, the components present in the first partial are obtained. The results show that the estimated parameters are stable following the four experiences. These results are compared in the following section to the modal basis of the simplified harp.

THEORETICAL STUDY OF THE HARP'S SYMPATHETIC MODES

Model of the concert harp

The modal behavior of the concert harp in low frequencies has been studied in a previous paper [START_REF] Le Carrou | Foltête: Experimental study of A0 and T1 modes of the concert harp[END_REF]. The vibrations of the soundboard correspond to bending motions which are symmetric according to the strings. The soundboard vibrations can then be described using an equivalent beam clamped at both ends, one corresponding to the pillar and the other one corresponding to the location of string 11. The vibratory model of the concert harp is thus constituted of an equivalent beam connected to 35 strings as shown in Figure 2.

The mechanical properties of the equivalent beam and of each string have to be determined: for the strings, most parameters are directly measured on the harp [START_REF] Carrou | Vibroacoustics of the concert harp[END_REF] except for the Young modulus and the density which are supposed to equal the data given in [START_REF] Bell | The physical properties of gut musical instrument strings[END_REF] and [START_REF] Bell | An acoustical investigation of the concert harp[END_REF]. Values of the tension are computed from the fundamental frequency of the tones, by considering strings fixed at both ends. For the equivalent beam, the geometrical parameters are directly measured on an isolated soundboard, allowing us to evaluate the density (ρ=553 Kg/m 3

). The area of the cross section A and the second moment of area J are measured at different positions along the axis. The average value of these two parameters is thus retained to characterize the equivalent beam (A eq =38.3 cm 2 and J eq =38.9 cm 4 ). Finally, the Young modulus E of the equivalent beam is determined so that the first bending mode of the beam in the beam-35 strings assembly equals 150 Hz (E= 5.9 Gpa).

Modal basis

The modal basis of the simplified harp is computed using the transfer matrix method detailed in [START_REF] Le Carrou | Modelling of sympathetic string vibrations[END_REF]. This computation is based on a wave guide model for each sub-structure of the beamstrings assembly. In the frequency range [0Hz-500Hz], 151 modes are obtained for the coupled system. Examples of modal shapes, associated to modes 26, 27 and 28 are presented in Figure 3. These modes are selected because their eigenfrequencies are close to the fundamental frequency of string 31, studied in the experimental section. In order to classify the modes of the assembly, a criterion called Kinetic Energy Ratio (KER) (see [START_REF] Le Carrou | Modelling of sympathetic string vibrations[END_REF] for an exact definition) has been defined. It corresponds to the ratio of kinetic energy of one sub-structure divided by the total kinetic energy of the structure. Its value is a percentage and allows us to identify the relative importance of each sub-structure displacement field.

ANALYSIS OF THE SPECTRAL COMPONENTS OF THE FIRST PARTIAL

When string 31 of the concert harp is free to oscillate (configuration (1) in Table I), two sinusoidal components are identified in the first partial. When this string is stopped (configurations (2), ( 3) and ( 4)), these two components disappear. In playing configuration, these two sinusoids have close frequencies, separated by 0.2 Hz. This result thus points out that the string vibrates following its two polarizations. These two polarizations are excited by the harp player in a different way depending on the string plucking, involving different initial amplitudes (showed by the difference of energy of each component). Moreover, as for the piano, one of the polarizations seems to be favored for transmitting its energy to the soundboard implying a rapidly decaying component at the plucking moment [START_REF] Weinreich | Coupled piano strings[END_REF]. This result shows that for a good modeling of the strings modal behavior, the two polarizations of string vibrations have to be taken into account. In the model of simplified harp, only one polarization is considered, which constitutes a limitation.

When the string is stopped during the oscillations, modes present in the instrument response are selected.

CONCLUSION

In order to analyze the sympathetic vibrations phenomenon present in the concert harp, experimental and theoretical investigations have been carried out. It has been shown that sympathetic modes, responsible of the sympathetic vibrations, are present in the instrument response. For the identification of these modes, a high resolution method, called ESPRIT, has been used. It allows the analysis of the spectral components present in the first partial of the plucked string. The components are very close one to another and cannot be identified with a classical Fourier analysis. With a simplified model of the vibratory behavior of the concert harp, composed of a beam connected to several strings, these modes are also found. The comparison between experimental and theoretical results is very concordant showing that our approach for modeling the harp is adequate.

The principal restriction of the model is that only one polarization of each string is taken into account. An extension of that study will be to add a supplementary polarization per string. From an experimental point of view, a study of the repeatability of the results will be investigated. The analysis of partials of higher order, where other phenomena such as octave vibrations can be present, will also be performed.

Figure 1 .

 1 Figure 1.-Spectrum of the first partial of the vibratory signal measured on the soundboard when string 31 is plucked. All strings are free to vibrate.

Figure 2 .

 2 Figure 2.-Model of the concert harp: beam-35 strings assembly.

Table I .

 I -Frequencies of identified components in the first partial of the accelerometer signals in the 4 configurations: (1) all free strings, (2) all strings damped except string 31, (3) strings 31 and 38 stopped during oscillations, (4) strings 31 and 42 stopped during oscillations.

			Configurations
		(1)	(2)	(3)	(4)
	Spectral Components [Hz]	123.59 123.35 123.08 123.78	123.53 123.66	123.28	123.08

Table II .

 II -Eigenfrequencies and Kinetic Energy Ratio for each sub-structure (strings 42, 38, 31 and the beam). KER expressed in % and rounded to the nearest whole number.In TableII, the eigenfrequencies and KER of each sub-structure of modes 26, 27 and 28 are reported. For these modes, the KER is significant only for three strings: 31, 38 and 42. String 31 (Cb 3 note of fundamental frequency 123.5 Hz) corresponds to the upper octave of string 38 (Cb 2 note at fundamental frequency 61.7 Hz) which is the upper fifth of string 42 (Fb 1 at fundamental frequency 41.2 Hz). The beam's KER equals zero, showing that mode shapes 26, 27 and 28 are dominated by the string's motion since the KER is distributed according to two or three strings, allowing the definition of these modes as string-string modes or sympathetic modes[START_REF] Le Carrou | Modelling of sympathetic string vibrations[END_REF], because they are responsible for the sympathetic vibrations. Actually, by the modal superposition principle, if string 31 is plucked, modes 26 to 28 are set into vibration and so are strings 38 and 42, generating the phenomenon of sympathetic vibrations. Figure 3.-Modal shape associated to modes 26, 27 and 28. The number indicated below the harp's arm points out the string number.

	42 38		42 38		42 38
		31	31		31
	Mode 26	Mode 27		Mode 28
	Mode	Eigenfrequency	KER for each string [%]	KER [%]
	label	[Hz]	String 42	String 38	String 31	Beam
	26	122,95	0	99	1	0
	27	123,29	97	1	2	0
	28	123,48	1	2	97	0

  By stopping strings 31 and 38, one more vibrating component disappears and by stopping strings 31 and 42, another one is absent from the instrument's response. With configurations (3) and (4), we can deduce that strings 38 and 42 participate in the sound radiated by the instrument. When the harpist plucks string 31, four modes are set into vibrations at the same time: two modes involving string 31 (one mode per polarization), one mode involving strings 31-38 and one mode involving strings 31-42. These last two modes are thus sympathetic modes. The stopping of string 31 does not necessary lead to the weakening of these modes, proving that the kinetic energy present in strings 42 and 38 is definitely more important than in string 31. Moreover, these results show that the sympathetic modes coming from each polarization of string 31 are not visible in the response. This fact can be explained by the weakening energies of these modal components, not allowing them to emerge from the noise.These experimental results can be compared to those obtained from the theory. Experimentally, we found two modes which mostly involve string 31 at 123.59 Hz and at 123.78 Hz. With the vibratory model, the string mode associated to string 31 is found at 123.48 Hz. Note that the model takes only one string polarization into account. The agreement between model and measurement is thus very good. For sympathetic modes 31-42 and 31-38, their eigenfrequencies are measured at 123.28 Hz and at 123.08 Hz. It should be noticed that the theoretical eigenfrequencies of these modes are found at 123.29 Hz and at 122.95 Hz, coinciding almost perfectly with experimental results. Apart from the fact that the vibratory model of the concert harp does not take the two polarizations of the string into account, results obtained from the model are in very good agreement with those obtained during experiments, thus validating our approach.


J L De Bougrenet De La Tocnaye

B Cochener

S Ferragut

D Iorgovan

Y Fattakhova

M Lamard

on 3D TV monitors

Keywords: physiological factors, 3D displays, binocular vision, stereopsis, stereotest Supervised stereo visual acuity tests implemented

In this paper we discuss under which conditions standard stereo visual acuity tests can be implemented on 3D TV monitors. In particular, we emphasize the role of environmental lighting conditions, on the measurement of the stereo visual acuity, when using conventional 3D tests, such as Wirt stereotests. We investigate the impact of parameters such as luminance, backlight and contrast when these tests are implemented on 3D TV monitors. We demonstrate that some deviations are observed when modifying the room luminance and the type of displays used (e.g. plasma (PDP) or liquid crystal (LCD) displays). Our measurements carried out on an human sample are supervised by pupil size measurements, using an eyes-tracker, enabling a better interpretation of the results. Finally, we discuss the benefit of using 3D tools to implement stereo visual acuity measurements.

I. INTRODUCTION

UMEROUS investigations have been performed to measure the effects of luminance and contrast on monocular or binocular visual acuity [START_REF] Johnson | Effects of luminance, contrast and blur on visual acuity[END_REF]- [START_REF] Chiu | Effect of screen contrast ratio and luminance level on visual lobe shape[END_REF]. In parallel, the recent revival of interest for stereoscopic vision and the related visual discomfort which can result [START_REF] Lambooija | Visual Discomfort in Stereoscopic Displays: A Review[END_REF]- [START_REF] Shibata | The zone of comfort: predicting visual discomfort with stereo displays[END_REF], has focused the attention on stereoscopic vision, justifying a fine analysis of 3D perception dysfunctions, definition of appropriate tests and protocols to measure them as well as recommendations such as ITU-R BT1428 [START_REF]Subjective assessment of stereoscopic television pictures[END_REF]. However, these recommendations do not take enough into account the display type and the environmental lighting conditions which could influence the measurement of the depth perception acuity. The main difference between monocular or binocular tests and depth perception acuity ones is the stimulation of our ability to determine the relative position of objects in space. The retina provides a two-dimensional image. The brain integrates these images from each eye to create a sense of three dimensions.

Factors such as excessive binocular parallax (apparent change in the direction of an object due to a change in observational position that provides a new line of sight) and accommodationconvergence conflicts are no longer of minor importance when disparity values surpass one degree limit of visual angle, so that monocular or binocular visual acuity tests are no longer sufficient to understand what occurs when 3D images are displayed. Factors like excessive demand of accommodationconvergence linkage (fast motion in depth), viewed at short distances, and unnatural amounts of blur, impact strongly this perception. In order to adequately characterize such mechanisms, it is needed to assess the stereo visual acuity in a way close to how 3D contents are displayed to observers. Before that, it is necessary, to investigate the impact of environment lighting on these measurements, with respect to the technology used to display the tests. Implementing conventional depth perception tests on TV displays will enable a first scaling of the test to establish a correspondence with conventional tests and then to analyse the impact of different display technology on these measurements due to their different luminance and correlation between contrast, luminance and backlight. This will enable us to establish a link with the above mentioned works and visual acuity dependence to these features, but applied to 3D perception.

II. RADIOMETRY, PHOTOMETRY AND LIGHTING SOURCES

When display systems are used to implement acuity tests, radiometric units are used, whereas photometry concerns the eyes' sensitivity (both are concerned here). Different lighting sources would influence the measurements making this point critical because it impacts the mydriasis, thus modifying the visual acuity. Fig. 1 summarizes the different lighting sources, being measured using a photometer (ILM-1335 from ISO-TECH). The most critical parameter is the display choice, being the main source of lighting in the direct sight of light. Two options are possible: plasma display panel (PDP) and liquid crystal display (LCD). Plasma has a very advanced brightness control technology minimizing the visual fatigue [START_REF] Meesters | A survey of perceptual evaluations and requirements of three-dimensional TV[END_REF], is gentle on the eyes and minimizes the chance of visual dysfunctions even after watching a long movie. Its brightness depends on the white area size. In contrast, LCD has a better brightness dynamic, due to the balance between direct luminance and backlight. Table I presents a comparison between these two options.

The contrast corresponds here to the luminance contrast which is the ratio between the higher luminance, L H and the lower luminance, L L . This ratio (or the Log of this ratio), often called contrast ratio (CR) is used for electronic visual display devices, like here [START_REF] Chen | Handbook of Visual Display Technology[END_REF]. This is a dimensionless number, often indicated by adding ":1" to the value of the quotient (e.g. CR = 100:1 with 1 ≤ CR <

). A CR = 1 means no contrast. When these displays are used outside a completely dark room, e.g. in the living room (around 100 lx) or in an office (around 300 lx minimum), ambient light is reflected from the display surface, adding to the luminance of the dark state and thus reducing the contrast. This point will be confirmed later on.

III. BINOCULAR DEPTH PERCEPTION AND STEREO-ACUITY

MEASUREMENTS USING A 3D TV Eyes are horizontally separated by the inter-pupillar distance (IPD) between 50 to 70 mm. Therefore, both eyes receive slightly different retinal images; stereopsis is a depth perception resulting on binocular retinal disparity between these images (difference in image location on the retina of an object seen by the left and right eyes); the brain fuses left and right images and form retinal disparity, from the distance, the brain extracts depth information. When eyes are focused on a certain point, eyes optical axes converge on that point intersecting at an angle called the parallactic angle. Smallest differences in this angle which can be interpreted as an impression of depth, determine the stereo-acuity. A usual way to measure it, is to use Randot or Wirt stereotests [START_REF] Bennett | Clinical visual optics[END_REF]. The Randot test is made up of 10 rectangles of 3 circles from 400" to 20", the Wirt, of 9 rhombs of 4 circles, from 800" to 40" (Fig. 2). One in each pattern is designed to appear standing forward. At a 40 cm distance, the range covered is from 800" to 40" of arc. Turning the Wirt upside down reverses the depth, the test figure receding from the background. For some patients (e.g. with fixation disparity), impressions of depth are affected by the direction in which the depth is presented. The depth is obtained using vectographs (Polaroid recording the stereoscopic pair), requiring polarized passive glasses to present left and right images to the corresponding eye. The idea is to implement this test on a TV, as it has been done for monocular tests [START_REF] Lenne | Automated visual acuity testing[END_REF], but on a TV 3D ready, generating left and right images. The first step consists in scaling the Wirt stereotest, at the correct size, to obtain the same depth, when using a 3D TV, as when using Wirt vectographs. The watching distance for a TV is between 1 and 2 m. Similar tests could be implemented for near vision at 40 cm, with a 3D laptop. For the considered 42' screen diagonal TV, we have chosen 1 m. For such a distance, stereo-acuities are divided by a factor 2.5 as shown Table II. Bold arc values correspond to values not resolved at 1 m, due to the pixel size of the considered display. This limit is here 25" but can be improved by increasing the distance between the observer and the display. Using active glasses (alternated vision by optical shutters) instead of passive ones provides better separation of the video pair than passive [START_REF] De Bougrenet | Liquid crystal actives glasses for 3-D cinema[END_REF], limiting possible image ghosting. This solution has been chosen here.

Because the parallax angle for infinity equals zero, objet (here the grey rings of Fig. 2) appears in relief as far as this angle reaches the stereo-acuity threshold. When this angle is determined, for a given IDP, one can derive the correct observation distance, which is then used to scale the Wirt stereotest size.

IV. PATIENTS AND METHOD

The testing sample has been obtained from 16 different observers (7 females and 9 males) ranging in the age from 22 to 57, (average age 32), serving as subjects. All subjects were informed about the nature of the study prior to experimentation and agreed to participate. Approval for the publication of subject data was obtained from Brest University hospital's institutional review board, according to the tenets of the Declaration of Helsinki. To perform a good calibration of the 3D Wirt tests it is required to operate with normal subjects, with normal acuity and good fusional range. Therefore, each volunteer had healthy eyes, had a good visual acuity (< 0.05 logMAR for a contrast > 100:1), and did not present visual dysfunctions after refractive corrections. They had good image merging capabilities : between 25 to 30 D for the far convergence, 8 to 10 D for the far divergence, and 35 to 40 for the near convergence, 12 to 14 for the near divergence, and they exhibit no objective obstacles to good stereopsis.

To calibrate our experiments a Standard Wirth stereotest has been used in a lighted room, with light sources directed to the test, to provide a better contrast. The test was located 40 cm far from the observer eyes, who is asked to tell which of four circles is highlighted in each rhomb, with a difficulty level increasing at each step. At this distance (corresponding to near vision) it is established that effects such as heterotropia, ocular accommodation asymmetry, oculomotor imbalance, fusional disparity (function of the age and population [START_REF] Bennett | Clinical visual optics[END_REF]) can bias the calibration. In addition, TV 3D ready is not appropriate at such a distance, because, first its resolution is not good enough and second the display being an important source of light, at the patient can be dazzled. This is why we decided to locate the patient at 1 m far from the screen. Calibrations performed, the test would be then transferred to a laptop, according to Table III, for near vision measurements. Concerning the test, practice trials are displayed before the testing phase so as to familiarize the observer with the procedure. A 2 minutes break (involving accommodation relaxation by fixing far field) was given every 10 minutes and the whole experiment lasted around less than one hour. The subjects were required to tell whether the randomly displayed depth is observed or not. Since the depth kept changing randomly, the chance of making a correct hit without seeing is very low. Therefore, in order to trigger a stimulus for testing, the subject had to keep staring at the centre fixation point. After a valid response, a neutral scene is displayed for a preset duration. A post-exposure masking stimulus with all the background object and target positions with no depth is presented immediately after stimulation. Subjects were required to indicate the estimated target location, by indicating the number of the Wirt stereotest. Because the visual acuity (and the stereo-acuity) is greatly affected by the level of background luminance [START_REF] Lovasik | Effects of aniseikonia, anisometropia, accommodation, retinal illuminance, and pupil size on stereopsis[END_REF], we considered two cases: the first where the tests are carried out in a dark room (lighting only due to the TV screen) and the second where the room is lighted by an indirect lighting of 300 lx Ne tubes (corresponding to home lighting conditions and luminance measured on a table in the middle of a room). To take into account factors affecting the dark adaptation, we measured intensity and duration of the pre-adapting light. Finally, we used two different test sets to de-correlate the impact of contrast and luminance on stereo-acuity measurements. First, for a given luminance value, we varied the contrast, and second for a contrast value, we varied the luminance. These tests have been carried out for two different 3D TV types a PDP and a LCD.

They are carried out for lighted and dark room (dark room does not mean scotopic vision < 10-2 cd/m 2 , but, in our case, a value close to 1 cd/m 2 , due to the display lighting in the darkest case, corresponding to mesopic conditions). We test three screen brightness values, for which we consider 4 different contrasts. For lighted room the adaptation time (for retina complete adaptation) is 5 minutes, for dark room 10 minutes. First we present the Wirt test with nine different depths (800'', 400'', 200'', 100'', 70'', 50'', 40'', 30'', 20'') to determine the threshold value for each observer. Then we present successively 4 Wirt stereotests with each nine rhomb of the same depth. One above the threshold, one at threshold and two below the threshold to check if the value found in the first test is the real observer limit. For instance, if 40" is the threshold, we test : 45", 40", 35", 30", if it is 100" we test : 110", 100", 95", 90". It means 4 contrasts, for 3 background luminances and 5 different stereotests (i.e. 60 patterns in the worst case). Then, the observer rests almost 2 minutes, and we change the screen brightness and increment. The protocol is implemented for both LCD and PDP displays. In order to control the impact of the luminance and contrast variations on the stereo-acuity measurements, we measure the pupil size, using an eye-tracker (Seeing Machines, FaceLab). Small pupil sizes reducing the angular resolution, they impact the stereoacuity threshold. The Eye tracker has been used particularly in transition areas (between depth perception and no perception). Dedicated software has been developed to determine the raw data of the eye tracker. This application works as a plug-in for the ObserverTM, developed in Java language [START_REF] Diverrez | Development of a procedure to examine the behavioural patterns relating to the perception of a virtual 3D image[END_REF]. Its main task is to recover, display, and store the data logged from the eye tracker. The Eye-tracker provides information such as sight tracking, eyes closing and blinking, PERCLOS (fatigue factor), concentration and jerks and head orientations, and pupil size measurement. This information is particularly

A. Experimental results

Tables III present SVA results obtained from a LCD TV, for various contrast and brightness values. Values are percentage of the population perceiving the displayed depth, with respect to whole, for both lighting conditions. The depth perception depending on different contrast, brightness, and lighting conditions is presented here. The individual features of depth perception of each testers were taken into account. Hence, values from tables allow us to select optimal conditions in which most of all testers (or all testers) can perceive the depth. Because large backlight values result in spatial non homogeneities (impacting the contrast dynamic), we operated with a reduced backlight value. Thus, the LCD contrast dynamic was between 1:1 and 1:100 assuming the contrast variation linear, except at boundaries.

To calibrate the results we compare them with ones obtained with a Wirt stereotest under the same room lighting conditions: 471 cd/m 2 (light reflected by Wirt stereotests is 368 cd/m 2 ). All subjects tested were able to distinguish the minimum depth of 40". Wirt stereotests require bright room lighting (because passive). Then, we substituted the PDP to the LCD and carried out the same experiments (Table IV A and B). We observed, in average, better performance with the PDD (probably due to backlight issues). This is noticeable in particular for the smallest depth columns (30" and 25"). In contrast, changes in environmental lighting conditions do not impact significantly the results, both for LCD and PDP.

To confirm our analysis, we measured the pupil size (Table V) in extreme lighting cases (for minimum and maximum brightness). Pupil contour and diameter have been obtained using an Eyes-tracker Facelab 5.0, from Seeing Machines (Fig. 3). As already shown in various papers (e.g. [START_REF] Lee | Experimental investigations of pupil accommodation factors[END_REF]), the pupil size changes with the intensity of the videos displayed on the screen, as well as with the depth fixation due to pupil accommodation. Although the second effect is less significant, it should be considered in our analysis (despite psychological effects are neglected). As reported in [START_REF] Lee | Experimental investigations of pupil accommodation factors[END_REF], it is shown clearly that, if the pupil diameter increases as the depth value increases, the correlation between the pupil size and depth fixation remains weak. Table V provides measurements of the pupil size for a male and a female representative sample, in dark room conditions. We do not notice significant differences between LCD and PDP according to this parameter. Values displayed in Table V give the average pupil sizes measured over the duration of depth fixation sequence. For 0.9 cd/m 2 brightness and a minimum contrast, in the dark case, we observe the mydriasis (around 7 mm) whereas for 62 cd/m 2 brightness and a maximum contrast the myosis is clearly visible (around 4 mm). In contrast, the evolution of the pupil size with the depth value, for a given brightness and contrast, is not in accordance with what was expected from [START_REF] Lee | Experimental investigations of pupil accommodation factors[END_REF]. The pupil size seems to increase as the depth value decreases. An explanation is the relatively small variation in the considered depths, so that the behaviour noticed in [START_REF] Lee | Experimental investigations of pupil accommodation factors[END_REF] is not relevant here. Furthermore, one can notice that these values are at the limit of detection for some patients (see Tables IIIA, IVA). Both considered male and female samples have been able to detect them after a more or less complex oculo-motor exploration. Hence, we suspect the pupil size increase is mainly due to the fixation strength in detecting small depths (at the limit of the SVA) in the picture. This is confirmed by a fine analysis of ocular fixation trajectories with the eyes-tracker (Fig. 4). Fixation points are depicted by large disks (durations > 100 ms). We observed situations where depth perception is not immediate, but obtained after complex eye motion strategies. Our protocol does not take into account this parameter, which probably plays an important role in the SVA threshold determination, as well as in the decision protocol (see [START_REF] Lenne | Automated visual acuity testing[END_REF]).

Fig. 4a shows ocular trajectories, for a 65" depth on each rhomb, but at different locations on the rhomb. This corresponds to a case where all locations have been detected by the observer. Fig. 4b shows ocular trajectories for a 35" depth on each rhomb. Answers show many mistakes (apart for 1 and 4) : left for 1 (left), up for 2 (right), left for 3 (up), down for 4 (down) and down for 5 (up), right for 6 (up), down for 7 (left), left for 8 (right) and left for 9 (up). We have displayed fixation disks, for a 200 ms duration threshold instead of 100ms. Depicted values correspond to the total fixation duration τ spent for each rhomb, involving several fixation points. We notice that rhombs 2, 3, 5 and 6 require intensive research (τ = 6.4, 6.4 , 9.1 , 7.3 s respectively). Fixation is focused, for rhombs 3 and 5, on left-down, whereas the depth (located up) is ignored. Finally, we observed that, over a given fixation duration, perception is no longer possible. The fatigue operates and the observer needs to rest and move the fixation to other parts of the test to recover some depth acuity.

Concerning the display type, we observed a slight difference at the eye-tracker between the PDP and the LCD, in the central part of the screen. Fixation durations are longer in case of LCD, probably due to the backlight which generates contrast non uniformity between the centre and the top and bottom of the screen, which is not the case with PDP. This point emphasizes the critical role of the eyes-tracker in the result interpretation, as well as the importance of taking into account various parameters in the determination of the SVA threshold. This use, of course, should be further investigated to clarify, for instance, small pupil size variations.

V. DISCUSSION

Results confirm the impact of brightness and contrast on the stereo visual acuity. SVA improves, when stereotest brightness and contrast increase (for both displays). This is hampered, for large brightness values, by possible dazzling or Modulation Transfer Function (MTF) reduction due to pupil closing [START_REF] Campbell | Optical quality of human eye[END_REF]. This explains why pupil size measurements should be used to validate the SAV and fix the right operating conditions of the display. Whatever the monitor choice, best results are obtained, for large contrast and brightness values (to prevent mydriasis), with a weak impact of the environmental lighting conditions. It makes a difference with fixed Wirt stereotest measurements, very sensitive to environmental lighting, which can vary substantially and therefore bias the test. With TV 3D ready, lighting is controlled by the display parameters (brightness and contrast). Determintaion of mydriaisis and myoisis is however required to fix the right contrast and brightness operating range, for a given display. About technology, PDP should be preferred to LCD, because providing more homogeneous lighting (confirmed both in the table and in the ocular strategy analysis). This is due to the LCD backlight systems which generate screen brightness and contrast non homogeneities (darker band in the middle of the screen for low brightness). This effect is mitigated with LCD laptops, well suited to implement 3D Wirt tests for near vision. Furthermore, PDP, even in the dark state, maintains a minimum lighting, smoothing the transition between dark and bright screen, during the test. Although these restrictions, we demonstrated the benefit of using a 3D TV to implement SVA tests: first because enabling more reliable and robust analysis, thanks to a better calibration and control of the lighting and contrast parameters, second due to the fact that smaller depths can be displayed randomly, offering more freedom degrees for practitioners to determine the depth sequence displayed to the patient. Another interesting result is the correlation with pupil size and eye motion analysis, providing a consolidation of diagnosis. Such a work could be extended to near vision (40 cm) using a 3D LCD laptop required at such a distance and by adjusting display parameters according to the above analysis. The goal is to demonstrate that 3D tools are very powerful to implement more sophisticated and new tests and protocols to provide efficient, reliable and faster binocular system analysis and diagnosis, but for which SVA implementation appears to be a necessary calibration step.

Fig. 1 .

 1 Fig. 1. Different sources of light which could impact the binocular acuity measurements.

Fig. 2 .

 2 Fig. 2. Wirt stereotest made up of circularly polarized images

Fig. 3 .

 3 Fig. 3. Pupil contour and diameter extracted by the eyes-tracker

Fig. 4b .

 4b Fig. 4b. Ocular trajectories and total fixation durations τ corresponding to bottom table 6, 35" depths are located on left for 1 (τ = 4.5 s) ,right for 2 (τ = 6.4 s), up for 3 (τ = 6.4 s), down for 4 (τ =4.5 s), up for 5 (τ = 9.1 s) and 6 (τ = 7.3 s), left for 7 (τ = 5.5 s), right for 8 (τ = 4.5 s), up for 9 (τ = 5.5 s), (brightness of 62 cd/m 2 )

TABLE I MEASURED

 I

	LUMINANCE AND CONTRAST DYNAMICS COMPARISON BETWEEN
		PDP AND LCD
	LCD LG42LX6500 [8] :
	Max brightness	Max backlight	Contrast = 707:1
	Min brightness	Max backlight	Contrast = 237:1
	Max Brightness	Min backlight	Contrast = 78:1
	Min Brightness	Min backlight	Contrast = 36:1
	PDP Panasonic TXP42GT20 [9] :
	Max brightness		Contrast = 187:1
	Min brightness		Contrast = 76:1

TABLE II STEREO

 II

	rhomb n° 1	800	320
	rhomb n° 2	400	160
	rhomb n° 3	200	80
	rhomb n° 4	140	56
	rhomb n° 5	100	40
	rhomb n° 6	80	32
	rhomb n° 7	60	24
	rhomb n° 8	50	20
	rhomb n° 9	40	16

-ACUITY (IN ARC-SEC) FOR TWO OBSERVATION DISTANCES 40 CM AND 1 METER TEST Patient at 40 cm Patient at 1 meter

TABLE III

 III

					A
				DARK ROOM
	Brightness	Contrast	65" 55" 45" 35" 30"	25"
	cd/m 2
	0.9	33:1	0	0	0	0	0	0
	1.2	55:1	31	31	13	13	0	0
	2	67:1	63	50	25	19	6	0
	3.5	100:1	88	75	38	25	25	0
	2.3	33:1	56	44	25	25	19	19
	9	55:1	88	75	25	19	19	0
	13	67:1	88	75	38	19	19	0
	23.5	100:1	94	88	50	38	19	0
	9	33:1	94	88	31	25	19	19
	26.5	55:1	88	88	50	31	31	0
	31	67:1	94	88	75	69	31	0
	62	100:1	100	94	75	75	25	0
	0.9 cd/m 2 measured in front of eyes for the minimum LCD brightness.

TABLE III

 III

						B
				BRIGHT ROOM
	Brightness	Contrast	65" 55" 45" 35" 30" 25"
	cd/m 2
	471	33:1	0	0	0	0	0	0
	479	55:1	38	19	6	6	0	0
	485	67:1	81	44	6	6	6	0
	494	100:1	88	81	13	13	13	0
	488	33:1	63	50	19	19	13	6
	497	55:1	94	94	38	31	13	0
	500	67:1	88	88	75	50	25	6
	503	100:1	88	88	63	50	19	0
	500	33:1	94	94	38	25	25	13
	518	55:1	94	94	63	44	31	0
	532	67:1	94	94	69	62	38	0
	556	100:1	94	94	69	69	31	0
	441 cd/m 2 measured in front of eyes for the minimum LCD brightness.

TABLE IV

 IV

						A
					DARK ROOM
	Brightness	Contrast	65" 55" 45" 35" 30" 25"
	cd/m 2
	0.9	20:1	0	0	0	0	0	0
	5	33:1	100	71	42	14	0	14
	8	40:1	100	100	42	42	14	0
	10	60:1	100	100	71	42	42	0
	4	20:1	100	100	42	0	0	0
	10	33:1	100	100	42	42	42	0
	17	40:1	100	100	42	42	42	14
	43	60:1	100	100	57	57	57	14
	9	20:1	100	100	57	14	14	0
	21	33:1	100	100	57	57	57	29
	29	40:1	100	100	57	57	42	14
	53	60:1	100	100	57	57	57	14
	0.9 cd/m 2 measured in front of eyes for the minimum PDP brightness.

TABLE IV

 IV

						B
				BRIGHT ROOM
	Brightness	Contrast	65" 55" 45" 35" 30" 25"
	cd/m 2
	441	20:1	0	0	0	0	0	0
	447	33:1	86	86	42	29	42	0
	468	40:1	100 100	57	42	42	42
	491	60:1	100 100	57	57	42	14
	447	20:1	100	71	14	0	0	0
	462	33:1	100 100	57	57	42	0
	488	40:1	100 100	57	57	57	14
	497	60:1	100 100	71	57	42	14
	471	20:1	100	86	57	14	14	14
	500	33:1	100 100	57	57	42	0
	503	40:1	100 100	71	71	42	29
	509	60:1	100 100	71	71	42	14
	441 cd/m 2 measured in front of the eye for minimum PDP brightness.

TABLE V

 V Top: a representative male and bottom a representative female sample. Each value is the pupil diameter in mm (np means no perception of depth, the pupil size remains constant).

	Brightness,	Contrast	65"		55"		45"		35"		30"		25"
	cd/m 2		L	R	L	R	L	R	L	R	L	R	L	R
	0.9	33:1	5.9	6.1	np		np		np		np		np
	1.2	55:1	4.7	5.0	np		np		np		np		np
	31	67:1	3.4	3.4	3.7	3.6	3.6	3.6	3.7	3.8	3.8	3.8	3.9	3.9
	62	100:1	3.5	3.4	3.5	3.5	3.8	3.7	3.9	3.8	4.0	3.9	3.8	3.8
	Brightness,	Contrast	65"		55"		45"		35"		30"		25"
	cd/m 2		L	R	L	R	L	R	L	R	L	R	L	R
	0.9	33:1	6.5	6.2	np		np		np		np		np
	1.2	55:1	6.1	5.9	np		np		np		np		np
	31	67:1	5.5	5.4	5.4	5.3	5.6	5.5	5.7	5.6	5.9	5.7	5.7	5.5
	62	100:1	4.6	4.5	4.9	4.8	5.0	5.0	5.3	5.2	5.3	5.1	5.2	5.1