A measurement method for the linear and ... - Semantic Scholar
SIGNAL PROCESSING, vol. 64, no. 1, 1998, pp. 49-60
A measurement method for the linear and nonlinear properties of electro-acoustic transmission systems F. Heinle b, R. Rabenstein a and A. Stenger a a Lehrstuhl
fur Nachrichtentechnik, Universitat Erlangen-Nurnberg, Cauerstr. 7, D-91058 Erlangen, Germany b Philips Semiconductors, Technology Centre for Mobile Communications, Stromerstr. 5-7, D-90443 Nurnberg, Germany
Abstract Competitive audio consumer products require not only cheap signal processing hardware but also low-cost analog equipment and sound transducers. The nonlinear distortions produced by these electro-acoustic transmission systems cannot be described and analyzed by standard methods based on linear systems theory alone. In order to take the nonlinear properties into account, we present a measurement method for the linear and nonlinear transmission characteristics of almost arbitrary systems and show its application to the analysis of electro-acoustic systems. Examples demonstrate the measurement of the impulse response of a loudspeaker-enclosuremicrophone-system with cheap analog equipment in the presence of background noise.
1 Introduction Real-time signal processing used to be the most expensive part of acoustic echo and noise control applications. The additional cost of complementing the digital hardware with good quality audio equipment was tolerable and consequently, nonlinear e ects of sound transducers could be neglected. The situation is changing with the availability of cheap computing power for real-time applications. Competitive audio consumer products require not only cheap signal processing hardware but also low-cost sound transducers. Software-only solutions of speech communication features for desktop computers must rely on Preprint submitted to Elsevier Preprint
12 February 1998
built-in microphones and speakers, whatever their quality may be. Therefore, nonlinear distortions have to be taken into account in the design of low-cost electro-acoustic systems. Furthermore, the digital transmission of speech or audio signals are subject to nonlinear e ects. Subband-coding with the least possible number of bits assigned to each band is a standard technique. Also low-cost analog-digital-converters produce distortions, which cannot be described by linear e ects alone. Many acoustic echo and noise control applications require measuring or estimating the properties of the loudspeaker-enclosure-microphone-system (LEMS) including any digital pre- and post-processing. The nonlinear e ects described above may show up in many places of this electro-acoustic transmission chain, so that it is usually impossible to consider them by a proper theoretical analysis. Also measurements of the distortion factor of isolated components do not give a complete picture of the nonlinear behaviour of the overall system. Moreover, common measurement methods for the room impulse response record only the linear transmission characterics and are blind to nonlinearities. This contribution presents a measurement method for the linear and nonlinear transmission characteristics of almost arbitrary systems and shows its application to the analysis of electro-acoustic systems. The measurement principle described here is based on an analysis method for nite wordlength implementations of time invariant digital lters [1,11,12] and periodically time-varying multirate systems [5,6,8]. It has been successfully applied to the simultaneous determination of the frequency response, the alias components, and the quantization noise of implemented lter banks. The following section 2 describes the electro-acoustic transmission system under consideration. The theoretical foundations of the measurement method are covered by section 3 along with two alternate interpretations given in section 4. Finally, some examples for the measurement of room inpulse responses with cheap sound transducers and background noise are presented in section 5.
2 Problem description 2.1 Overview
Fig. 1 shows an electro-acoustic transmission system for the measurement of an LEMS, consisting of ampli ers, digital-to-analog and analog-to-digital converters. A discrete measurement system MS excites the transmission system with a discrete input signal v(k) and records its response y(k). The design of the measurement system is the topic of this contribution. 2
D
v(k)
A
LEMS
MS y(k)
A
D
Fig. 1. Electro-acoustic transmission system for the measurement of a loudspeaker-enclosure-microphone-system (LEMS)
2.2 Linear approximation
The transmission of acoustic signals from the loudspeaker to the microphone is described by the linear wave equation [17]. However, if other sound sources are present in the LEMS, the measurement signal may be superimposed by additive background noise. Ampli ers and sound transducers can be modelled as linear systems if the distortions produced by them are negligible. This assumption is only valid for good quality (expensive) equipment. The converters are nonlinear systems by their very nature. They can be approximated by linear systems only for a high wordlength of the digital signals. If background noise, distortions of the analog equipment, and quantization e ects are neglected, the electro-acoustic system from Fig. 1 can be replaced by a linear model with an overall impulse response hlin(k) between input v(k) and output y(k) (see Fig. 2). v(k)
hlin (k)
y(k)
Fig. 2. Linear model
An estimation h^ lin(k) of this impulse response can be obtained with the Discrete Fourier Transformation (DFT), e.g. in the most simple case as ^ lin(k) = DFT ?1 h M
(
DFTM fy(k)g DFTM fv(k)g
)
(1)
A number of techniques exist, which di er in the choice of the input signal v(k), the DFT length M , and in more elaborate details of obtaining the estimate ^ lin(k) [16]. However, all these techniques rely on the assumption of linearity h and are blind to any of the nonlinear e ects described above. 3
2.3 Weakly nonlinear approximation
A detailed account of all possible nonlinear contributions in Fig. 1 leads to a general nonlinear model according to Fig. 3. Due to the variety and di erent nature of the nonlinear components, an exact nonlinear model S of the transmission system from Fig. 1 would be too complex to be handled eciently.
S
v(k)
y(k)
Fig. 3. Nonlinear model SL
yL (k)
v(k)
y(k) SN
n(k)
Fig. 4. Weakly nonlinear model
In order to circumvent this diculty, we approximate the nonlinear system
S by a weakly nonlinear model [11] as shown in Fig. 4. It consists of a parallel arrangement of a linear system SL and a nonlinear system SN. Again, it is not feasible to derive SL and SN analytically from the components of the transmission system. Instead, the characteristics of both systems have to be obtained experimentally from suitably chosen input signals v(k) and the corresponding responses y(k). A method for the determination of SL and SN from measurements of v(k) and y(k) is described in the following.
3 Description of the measurement method 3.1 Frequency response of the linear subsystem SL
The frequency response of the linear subsystem SL is determined such that the contribution of the nonlinear subsystem SN to the total output signal y (k ) becomes a minimum. This is achieved by expressing the output yL(k ) of SL in Fig. 4 by a convolution with the impulse response h(k) of SL: yL(k) = h(k ) v (k ). The minimization of the quadratic mean value of the output of the nonlinear system "(k) = E fj n(k) j2g is equivalent to determining the impulse response h(k) such that the expected value E fg n
( ) = E j y(k) ? h(k) v(k) j2
" k
4
o
(2)
becomes a minimum. By expanding the squared magnitude in (2) we obtain ( ) = E fy(k)y(k)g ?
" k
? +
+ 1 X
?1
=
+ 1 X
+ 1 X
?1
=
(
)E fy(k)v(k ? )g ?
h
( )E fy(k)v(k ? )g +
h
+ 1 X
( ) ( )E fv(k ? 1)v(k ? 2)g
h 1 h 2
?1 2 =?1 + 1 + 1 X X = y2 ? h ()'yv () ? h()'vy (?) + =?1 =?1 + 1 +X 1 X + h (2 )h(1 )'vv (2 ? 1 ): 1 =?1 2 =?1 1 =
(3)
where 'vv (k); 'vy (k) and 'yv (k) are the auto- and cross-correlation sequences between input and output. The unknown values h(k) of the impulse response follow by di erentiation of "(k ) with respect to h () from the minimum condition @"(k ) = 0; 8: (4) @h() The di erentiation of a function f (x) with respect to its conjugate complex variable x can be carried out in a formally correct way by the so-called Wirtinger calculus [2,13]. The result is ( )=
'yv
+ 1 X
?1
=
()
( ? ):
h 'vv
(5)
Finally, we obtain by Fourier transformation F an estimation H (ej ) = Ffh(k)g of the frequency response of the linear subsystem yv (ej ) : j
H (e ) = (6) vv (ej ) vv (ej ) = F f'vv ()g and yv (ej ) = F f'yv ()g are the power spectral density (PSD) and the cross power density of v(k) and y(k) respectively. For simplicity, we furthermore assume that vv (ej ) 6= 0 for all frequencies of interest. Since these power densities are not known in advance, they have to be determined from the signals v(k) and y(k). In general, a wide range of possi5
ble input signals and a variety of corresponding spectral estimation methods could be utilized. Here, we will focus on periodic signals v(k) and y(k) with period M . Then it is sucient to determine H (ej ) at M discrete frequencies = 2=M; = 0; ; M ? 1. The values H (ej ) of (6) can now be expressed by V () = DFTM fv(k)g and Y () = DFTM fy(k)g [1] H ej
(
) = E fY ()V (2)g ; = 0; ; M ? 1 : E jV ()j o
n
(7)
The expected values E fg are calculated from the measured input and output sequences by applying L sample sequences v(k); = 0; : : : ; L ? 1 of the same stochastic process to the transmission system and by averaging the results. Thus, the estimated value H^ (ej ) for the frequency response of the linear subsystem is obtained ?
LP1
^ (ej ) = H
Y =0 LP1
?
=0
()V()
jV()j2
;
= 0; ; M ? 1 :
(8)
The procedure is simpli ed, if the DFT-spectrum of the sample sequences has a constant magnitude jV()j = jV ()j. This can be achieved by deriving the sample sequences v(k) from an arbitrary periodic, deterministic signal v(k) by adding a stochastic phase term ( ) = v(k + M ) = v(k)ej' :
(9)
v k
The phase ' is a stochastic variable, equally distributed in the interval [?; +). For input signals of this kind, the computation of the estimated value simpli es to ^ (ej ) = 1 DFTM H LV ()
(
?
L X1 =0
y k e?j'
()
)
:
(10)
Fig. 5 shows the structure of the resulting measurement arrangement. In our application, a complex chirp signal ( ) = ej(=M )k2
v k
(11)
has been used as the deterministic component v(k) of the input. This signal achieves the minimum possible crest factor. 6
e
j
v(k)
1L L
S
1
DFT
ˆ j ) H(e
0
DFT
Fig. 5. Measurement arrangement for the frequency response of the linear subsystem SL. ( ) denotes the conjugate value. 3.2 Power spectral density at the output of the nonlinear subsystem SN
Once the frequency response of the linear subsystem SL is estimated, we can use this result to obtain the PSD nn (ej ) of the nonlinear subsystem SN. For the estimated output of the nonlinear system follows from Fig. 4 n ^ (k) = y(k) ? ^h(k) v(k) ; (12) where ^h(k) is the estimated impulse response of SL. The estimation of nn (ej ) is performed in the frequency domain by evaluation of L?1 jY()j2 ? H^ (ej ) 2 jV()j2 ^ nn (ej ) = M (L1? 1) X
=0
(13)
3.3 A new windowing technique
The proposed method for estimating the PSD was originally developed for the investigation of xed-point digital lters under the assumption of periodic quantization noise. However, it is well-known that averaged periodograms L?1 ^ nn (ej ) = M (L1? 1) jN()j2 X
=0
(14)
are biased in case of non-periodic noise sequences n(k), i.e. the estimate does not equal the true PSD even for an in nite number L of single measurements [7]. Due to the implicit windowing of n(k) with a rectangular window ( )=
r k
8 > < > :
1; k = 0 : : : M ? 1 0; else 7
(15)
the expected value of (14) yields
E ^ nn (e n
j
+ 2 1 ) = 2M nn (ej( ?)) R(ej ) d: ? Z
o
(16)
Thus, the PSD estimate is smeared since the spectrum R(ej ) of r(k) has high sidelobes. Many spectral estimation procedures such as the Welch-Bartlett method or the Blackman-Tukey method try to circumvent this problem by using di erent windows, e.g. the Bartlett, Hamming, or Hanning window. Nevertheless, all these procedures only use possibly overlapping signal blocks of length M for the estimation. As the window length M remains unchanged the sidelobes can only be reduced on the expense of the width of the mainlobe. Recently, in [9] a new windowing technique was proposed which overcomes most of these problems. Here, we will only give a brief survey of some major points of this new method. Obviously, the spectral properties of a window sequence can be improved by increasing the window length. Therefore, the method uses N subsequent blocks of n(k) and one common window f (k ) of length N M for estimating the PSD at M discrete frequencies = 2=M; = 0 : : : M ? 1. Thus, the PSD estimate (14) has to be modi ed as follows: L?1 ^ nn () = M (L1? 1)
?1
M X NX
=0
1
= M (L ? 1)
?
k =0
() ()
?
?
L X1 X1 N X1 M =0
L?1 = M (L1? 1)
k =0 =0
?1
=0
k =0
( + M )n (k + M ) e?j
f k
|
X X M
2
2k
f k n k e?j M
{z
2k
M
2
}
~()
n k
~()
n k e
2 ?j 2k M :
(17)
With F (ej ) = Fff (k)g, the PSD estimate then yields
E ^ nn (e n
j
+ 2 1 ) = 2M nn (ej( ?)) F (ej ) d: ? o
Z
(18)
Now, it is possible to deduce the following desirable properties of the window sequence f (k) from the above equations:
The PSD estimate at = 2=M should comprise all spectral components of the continous PSD in the frequency bin 2( ? 1=2)=M < 2( + 8
1=2)=M . The optimum window therefore has the spectrum Fopt ej
( )=
8 >
:
:
(19)
Consequently, the PSD estimate is the integral over the respective frequency bin at . Obviously, Fopt(ej ) cannot be realized with a nite window length. Therefore, we have to design a window f (k) which approximates the optimum window under the given constraints. The proposed method is capable of measuring the frequency response and the PSD simultaneously. For eciency reasons, it is therefore necessary that the windowing does not cause any interference of periodic output signal components which result from di erent periodic input signal components. This can be easily achieved by prescribing the following zeroes of the window spectrum: F ej
(
) = 0;
= 1 : : : M ? 1:
(20)
Finally, the windowing should preserve the total power, i.e. 1 M
+ ! 1 j ^ E nn (e ) = 2 nn (ej ) d : =0 ?
?
M X1
n
Z
o
(21)
From this constraint, the last condition immediately results: ?
M X1 =0
F
2
(ej( ?2=M )) =! M 2:
(22)
Evidently, the window sequence has to be a so-called M th-band lter which is well-known from lter bank theory [15]. In [9] a highly ecient method for designing such windows was introduced. Our PSD estimation procedure which is based on these window sequences shows an excellent performance even in case of a modest N . This will also be illustrated by the examples at the end of this paper.
4 Alternate interpretations The derivation of the frequency response of the linear subsystem in section 3.1 is quite straightforward and self-contained. In order to show interrelations with other results, we give here two alternate interpretations of the obtained estimate from the point of view of optimal ltering theory and estimation theory. 9
4.1 Wiener lter
The minimization of a quadratic error term according to (2) is also used for the derivation of the optimal reconstruction lter (Wiener lter). However, the scope of Wiener ltering is di erent from the measurement method considered here. Nevertheless, the measurement problem and the optimal lter problem can be discussed in the same framework. Both lead to the same minimization problem and consequently the solutions are formally identical. At rst we consider the general optimal reconstruction problem as shown in Fig. 6. For later reference to the measurement problem, the signals in Fig. 6 have been denoted with the same letters as in the weakly nonlinear model from Fig. 4. However, the meaning of the signals in the optimal lter problem is di erent. y(k) is an unknown input signal to the system S. Only the output v (k ) is available by observation. The problem is to nd a linear system H such that its output yL(k) is a reconstruction of the unknown signal y(k). The frequency response H (ej ) of the system H is determined by least squares minimization of the di erence ( ) = y(k) ? yL(k) = y(k) ? h(k) v(k)
n k
(23)
with respect to the coecients h(k) of the impulse response of H. The result H (ej ) is expressed by the power spectral density vv (ej ) of the observed signal v(k) and the cross power density yv (ej ), which re ects the available knowledge of the system S as [3,4] yv (ej ) : (24) H (ej ) = vv (ej ) y(k)
S
v(k)
H
yL (k)
n(k)
Fig. 6. Optimal lter problem
In order to establish a relation to the measurement problem considered previously, we redraw Fig. 6 and arrive at the modi ed arrangement of Fig. 7. Note that the problem description is identical to Fig. 6. The shaded area indicates that the available knowledge of the reconstruction problem consists of the correlation properties of the signals v(k) and y(k) represented by vv (ej ) and yv (ej ). Now we return to the measurement problem as described in section 2. Fig. 8 shows the arrangement for the detection of the linear subsystem of the weakly 10
vv
v(k)
(e j ),
(e j )
y(k)
S H
yv
n(k)
yL (k)
Fig. 7. Modi ed arrangement of the optimal lter problem
nonlinear model according to Fig. 4. The signal designations in Figs. 8 and 4 correspond to each other. The derivation of the frequency response H (ej ) of the linear subsystem by minimization of ( ) = y(k) ? yL(k) = y(k) ? h(k) v(k)
n k
(25)
has been shown in section 3.1. Here, we focus on the similiarities and differences of the measurement problem according to Fig. 8 and the modi ed arrangement of the optimal lter problem according to Fig. 7. An obvious di erence is that the role of v(k) and y(k) as input and output signals of S are interchanged. However, both problems share the property that the system S is described by the cross power spectrum yv (ej ). A second di erence is the purpose of the linear system H. It serves to reconstruct an unknown signal in the optimal lter problem and it describes the linear part of an unknown system in the measurement problem. However, the minimization problems (23) and (25) are identical for both cases, irrespective of the di erent meaning of the signals v(k) and y(k) in Figs. 7 and 8. vv
v(k)
(e j ),
S H
yv
(e j )
y(k)
n(k)
yL (k)
Fig. 8. Measurement problem
Since both the way of describing the system S and the minimization problems (23) and (25) formally agree, we can expect, that also the solutions for the frequency responses H (ej ) of the unknown systems H coincide. This is indeed the case, as is seen by comparison of the frequency response of the optimal lter according to (24) and the frequency response of the linear subsystem of (6). This shows that the foundation of the presented measurement method can also be interpreted by formal agreement with optimal lter theory. 11
4.2 Maximum likelihood method
The frequency response estimate can even be derived in a quite di erent way which ensures optimality in terms of estimation theory. This approach is called maximum likelihood (ML) estimate. We will derive this estimate for a single measurement and subsequently generalize the results. We start from a frequency domain description of the weakly nonlinear system according to Fig. 9 (compare Figs. 4 and 8). It can be formulated in terms of the DFT values of blocks of measured data with V () = DFTfv(k)g and Y (), N () alike. H () is the DFT of the impulse response of the unknown linear subsystem. ( ) = H ()V () + N ()
(26)
Y
v(k)
H
y(k) n(k)
Fig. 9. Additive noise description of the weakly nonlinear system
Under the assumption that the probability density function (PDF) of the disturbing signal component N () is known, we try to determine a frequency response estimate H^ ML () which most probably caused the observed output Y () of the system under test. For convenience we introduce vectors n, N, V, Y, H, and H^ ML which consist of the corresponding time or frequency domain samples. The ML estimate of the frequency response then can be found by maximizing the conditional PDF pYjH(YjH) with respect to the frequency response H, i.e. we have to determine ^ ML ) = max pYjH(YjH) : (27) pYjH (YjH n
o
H
Due to equation (26) the output Y () is distributed around the deterministic value H ()V () with the PDF of the disturbing signal N () [14], i.e. pYjH
(YjH) = pN(Y ? diag fH ()g V) = pN(N)
(28)
Thus, we have to maximize the noise PDF pN(N) with respect to the frequency response estimate. The required frequency domain PDF pN(N) can be determined from its time domain equivalent in a straightforward manner. For this purpose, the joint PDF of M subsequent samples of the underlying noise-like process n(k) = nr (k) + jni(k) is needed. Under quite general 12
conditions, this process is normally distributed and has an autocorrelation ( ) = 'nr nr () + 'ni ni () with 'nr nr () = 'nini (). With the autocorrelation matrix nn = E nnH the required joint PDF then can be written as 1 H ?1 n: pn (n) = M exp ? n (29) nn jdet nn j 'nn
o
n
Now, the PDF in the frequency domain can be determined easily using the DFT matrix WM and n = M1 WMH N. We obtain 1 NH W ?1 WH N: 1 exp ? (30) pN (N) = M M nn M jdet nn j M2 =: C |
{z
}
Instead of maximizing pN (N) it is convenient to introduce the so-called loglikelihood function L(H) = ln M jdet nn j pN (N) (31)
and maximize it with respect to H (). With N () = Y () ? H ()V () the unknown noise sequence can be eliminated and the log-likelihood function 1 M ?1 M ?1 c (Y ( ) ? H ( )V ( )) (Y () ? H ()V ())(32) L(H) = ? X
M2
X
=0 =0
follows. From the maximum condition 1 M ?1 c V () @L(H) = @H () H=H ^ ML M 2 =0
X
( ) ? H^ ML()V () =! 0 (33)
Y
the ML estimate immediately results: ^ ML ()V (): Y () = H
(34)
If we furthermore restrict ourselves to signals with a deterministic magnitude, i.e. (35) jV ()j2 = E jV ()j2 ; n
o
this estimate can be generalized for multiple measurements and we obtain ^ ML () = 1 2 E fY ()V ()g : HML () = E H (36) jV ()j n
o
13
Obviously, this is equivalent to the aforementioned estimate (7) if (35) is taken into account. Furthermore, it can be shown that the estimation error () = E fN ()=V ()g is given by () =
?
M X1 =0
a
@
ln pN(N) @H ( )
(37)
with constant factors a . This is a generalization of the Cramer-Rao bound [10,14] for the joint estimation of multiple parameters. Due to the equality in (37), the ML estimate (34) meets this bound. Thus, it is an ecient estimate in terms of estimation theory. Nevertheless, it has to be emphasized that this result only holds for multiple measurements if input signals with a deterministic magnitude according to (35) are used.
5 Results This section presents measurements of an LEMS both with high quality and low quality equipment. The LEMS was an anechoic chamber containing the measurement gear which caused some re ections. The high quality equipment was a studio microphone, a preampli er (Bryston Mod. PB4) and an active loudspeaker (GENELEC Mod. 1013A). For the low quality measurements, preampli er and loudspeaker were replaced by a one-chip ampli er (TBA 820M) with a cardboard mounted noname speaker (6 cm diameter). 20 log |H(f)|
0 −10 −20 −30 0
1
2 f [kHz]
3
4
Fig. 10. Frequency response using high quality audio equipment and standard measurement method
Fig. 10 shows the frequency response of the LEMS, measured with high quality audio equipment. Here, we used a standard measurement procedure according to (1) which assumes a strictly linear system behaviour. Fig. 11 shows the same measurement using low-cost equipment and again a standard measurement method. Due to this equipment a di erent frequency response has to be expected even though the room impulse response was not a ected by the change. Nevertheless, the measured frequency response shows 14
20 log |H(f)|
0 −10 −20 −30 0
1
2 f [kHz]
3
4
Fig. 11. Frequency response using low-cost audio equipment (standard measurement method)
an unexpectedly strong noise-like behaviour. This is due to the fact that conventional measurement procedures are not able to distinguish between linear and nonlinear system components, the latter being caused by distortions of the ampli er and speaker. Thus, the result is a superposition of both components. 20 log |H(f)|
0 −10 −20 −30 0
1
2 f [kHz]
3
4
Fig. 12. Frequency response using low-cost audio equipment with background noise (standard measurement method)
The frequency response measurement shown in Fig. 12 was performed with an additional broadband background noise 4dB below the signal level. It was produced by a noise generator fed into a guitar ampli er. The detrimental distortion of the measurement result by the superimposed noise source is obvious. In contrast, Fig. 13 shows the measurement results with our proposed method for the low-cost equipment with background noise as in Fig. 12. This method is capable of recording linear and nonlinear e ects separately. The method yields the frequency response of the linear subsystem as well as the PSD of nonlinear distortions including noise. Obviously, the true frequency response deviates from the original one (Fig. 10) due to additional linear distortions of the low-cost equipment. The remaining set of results demonstrates the suitability of the windowing technique described in section 3.3 for the discrimination between LEMS response and narrowband background noise. To this end, we repeat the measurement with high quality equipment similar to Fig. 10. During the measurement, an additional noise source was present in the LEMS which emitted a narrowband random signal with a signal-to15
20 log |H(f)|
0 −10 −20 1
2 f [kHz]
3
4
0 −10 −20 −30 −40 −50 −60 0
1
2 f [kHz]
3
4
PSD [dB]
−30 0
Fig. 13. Frequency response and PSD using low-cost audio equipment with background noise (proposed measurement method)
noise-ratio (SNR) of 3.5 dB. The noise signal was obtained from a broadband noise source by ltering with a digital FIR lter with a passband between 1 and 1.5 kHz. Fig. 14 shows the frequency response and the PSD obtained with the proposed measurement method and a rectangular data window. As a reference, the squared magnitude of the FIR form lter which was used for the generation of the narrowband noise is shown. Although averaging with L = 100 sample sequences was used (see (8) and (10) ), the measured PSD still di ers signi cantly from the true one. This deviation can be attributed to the non-periodicity of the noise signal. 20 log |H(f)|
0 −10 −20 1
2 f [kHz]
3
4
0 −10 −20 −30 −40 −50 −60 0
1
2 f [kHz]
3
4
PSD [dB]
−30 0
Fig. 14. Frequency response and PSD using high quality audio equipment with bandpass noise (proposed measurement method, rectangluar window, L = 100 sample sequences)
The same system was measured with the M th band window described in section 3.3 and with L = 10 and L = 100 sample sequences, respectively. Two 16
e ects can be observed from the results in Figs. 15 and 16. First, the corruption of the frequency response is evidently restricted to the actual frequency range of the bandpass noise. Second, the PSD of the nonlinear disturbances closely resembles the true PSD of the noise signal. The latter is given again by the squared magnitude of the FIR form lter. The obvious improvement in performance is due to the use of an appropriately designed data window. 20 log |H(f)|
0 −10 −20 −30 0
2 f [kHz]
3
4
3
4
form filter
PSD [dB]
0 −10 −20 −30 −40 −50 −60 0
1
1
2 f [kHz]
Fig. 15. Frequency response and PSD using high quality audio equipment with bandpass noise (proposed measurement method, M th band window, L = 10 sample sequences), squared magnitude of form lter 20 log |H(f)|
0 −10 −20 −30 0
2 f [kHz]
3
4
3
4
form filter
PSD [dB]
0 −10 −20 −30 −40 −50 −60 0
1
1
2 f [kHz]
Fig. 16. Frequency response and PSD using high quality audio equipment with bandpass noise (proposed measurement method, M th band window, L = 100 sample sequences), squared magnitude of form lter
Finally, we repeat the same experiment with the bandpass noise replaced by a sine wave of 1 kHz and an SNR of 4 dB. Fig. 17 shows the results for L = 100 sample sequences. The distortion of the measured frequency response of the linear subsystem by the sine wave is hardly visible. On the other hand, 17
we see that the PSD gives a perfect picture of the sinusoidal disturbance of the measurement. The absence of any spectral smearing or leakage clearly demonstrates the suitability of the M th-band window for spectral estimation. 20 log |H(f)|
0 −10 −20 1
2 f [kHz]
3
4
0 −10 −20 −30 −40 −50 −60 0
1
2 f [kHz]
3
4
PSD [dB]
−30 0
Fig. 17. Frequency response and PSD using high quality audio equipment with sinusoidal noise (proposed measurement method, M th band window, L = 100 sample sequences)
6 Conclusion We have presented a method for measuring the properties of the loudspeakerenclosure-microphone-system in the presence of nonlinear distortions and background noise. These disturbances have to be considered in the evaluation of the performance of low-cost audio products in oce or home environments. The measurement method is based on the analysis of implemented digital lters and multirate systems. Here, it has been extended to the investigation of the linear and nonlinear properties of electro-acoustic transmission systems. The fundamentals of the measurement method can be derived from a weakly nonlinear system model by minimization of a quadratic error criterion. Nonlinear e ects are characterized by the PSD of an additive noise signal. A specially adapted windowing technique provides extremely sharp estimates also for narrowband noise sources. The interrelations with optimal ltering theory and estimation theory have been shown by alternate interpretations of the measurement principle using the Wiener lter and the maximum likelihood approach. The results presented show the ability to clearly distinguish between the linear and nonlinear properties of an electro-acoustic transmission system. The 18
frequency response of the linear subsystem and the PSD which describes nonlinear distortions and background noise can be obtained separately.
Acknowledgement The authors wish to thank Eberhard Hansler for valuable hints, Thomas Klinger for preparing the results and Ivan Selesnick for his help with the manuscript.
References [1] Y. Dong. Ein neues Verfahren zur Messung der Eigenschaften schwach nichtlinearer Systeme. PhD thesis, Universitat Erlangen-Nurnberg, Lehrstuhl fur Nachrichtentechnik, Erlangen, 1990. [2] W. Fischer and I. Lieb. Ausgewahlte Kapitel aus der Funktionentheorie. Vieweg Verlag, Braunschweig, 1988. [3] B. Girod, R. Rabenstein, and A. Stenger. Einfuhrung in die Systemtheorie. B. G. Teubner, Stuttgart, 1997. [4] E. Hansler. Statistische Signale. Springer Verlag, Berlin, 1997. [5] F. Heinle and H. W. Schuler. An enhanced method for measuring the performance of multirate systems. In Proc. Int. Conf. on DSP, pages 182{187. Limassol, Cyprus, 1995. [6] F. Heinle and H. W. Schuler. Measuring the performance of implemented multirate systems. In Proc. Int. Conf. Acoustics, Speech and Signal Proc. (ICASSP 96), pages 2754{2757, Atlanta, USA, 1996. IEEE. [7] S.M. Kay. Modern Spectral Estimation { Theory & Application. Prentice-Hall, Englewood Cli s, USA, 1988. [8] R. Reng and H. W. Schuler. Measurement of aliasing distortions and quantization noise in multirate systems. In Proc. Int. Symp. on Circuits and Systems, pages 2328{2331, 1992. [9] H. Repp. Fensterung beim Rauschklirrmeverfahren. PhD thesis, Universitat Erlangen-Nurnberg, Lehrstuhl fur Nachrichtentechnik, Erlangen, 1998. [10] L.L. Scharf. Statistical Signal Processing: Detection, Estimation, and Time Series Analysis. Addison-Wesley, Reading, USA, 1991. [11] H. W. Schuler and Y. Dong. A new method for measuring the performance of weakly nonlinear and noisy systems. FREQUENZ, 44:82{87, 1990.
19
[12] H. W. Schuler and F. Heinle. Measuring the properties of implemented digital systems. FREQUENZ, 48:3{7, 1994. [13] H. Tietz. Funktionentheorie. In: Mathematische Hilfsmittel des Ingenieurs, vol. I, R. Sauer, I. Szabo (eds.). Springer Verlag, Berlin, 1970. [14] H.L. Van Trees. Detection, Estimation, and Modulation Theory, Part I. John Wiley & Sons, New York, 1968. [15] P.P. Vaidyanathan. Multirate Systems and Filter Banks. Prentice Hall, Englewood Cli s, USA, 1993. [16] N. Xiang. Using M-sequences for determining the impulse responses of LTIsystems. Signal Processing, 28:139{159, 1992. [17] L. J. Ziomek. Fundamentals of Acoustic Field Theory and Space-Time Signal Processing. CRC Press, Boca Raton, 1995.
20
A measurement method for the linear and nonlinear properties of electro-acoustic transmission systems F. Heinle b, R. Rabenstein a and A. Stenger a a Lehrstuhl
fur Nachrichtentechnik, Universitat Erlangen-Nurnberg, Cauerstr. 7, D-91058 Erlangen, Germany b Philips Semiconductors, Technology Centre for Mobile Communications, Stromerstr. 5-7, D-90443 Nurnberg, Germany
Abstract Competitive audio consumer products require not only cheap signal processing hardware but also low-cost analog equipment and sound transducers. The nonlinear distortions produced by these electro-acoustic transmission systems cannot be described and analyzed by standard methods based on linear systems theory alone. In order to take the nonlinear properties into account, we present a measurement method for the linear and nonlinear transmission characteristics of almost arbitrary systems and show its application to the analysis of electro-acoustic systems. Examples demonstrate the measurement of the impulse response of a loudspeaker-enclosuremicrophone-system with cheap analog equipment in the presence of background noise.
1 Introduction Real-time signal processing used to be the most expensive part of acoustic echo and noise control applications. The additional cost of complementing the digital hardware with good quality audio equipment was tolerable and consequently, nonlinear e ects of sound transducers could be neglected. The situation is changing with the availability of cheap computing power for real-time applications. Competitive audio consumer products require not only cheap signal processing hardware but also low-cost sound transducers. Software-only solutions of speech communication features for desktop computers must rely on Preprint submitted to Elsevier Preprint
12 February 1998
built-in microphones and speakers, whatever their quality may be. Therefore, nonlinear distortions have to be taken into account in the design of low-cost electro-acoustic systems. Furthermore, the digital transmission of speech or audio signals are subject to nonlinear e ects. Subband-coding with the least possible number of bits assigned to each band is a standard technique. Also low-cost analog-digital-converters produce distortions, which cannot be described by linear e ects alone. Many acoustic echo and noise control applications require measuring or estimating the properties of the loudspeaker-enclosure-microphone-system (LEMS) including any digital pre- and post-processing. The nonlinear e ects described above may show up in many places of this electro-acoustic transmission chain, so that it is usually impossible to consider them by a proper theoretical analysis. Also measurements of the distortion factor of isolated components do not give a complete picture of the nonlinear behaviour of the overall system. Moreover, common measurement methods for the room impulse response record only the linear transmission characterics and are blind to nonlinearities. This contribution presents a measurement method for the linear and nonlinear transmission characteristics of almost arbitrary systems and shows its application to the analysis of electro-acoustic systems. The measurement principle described here is based on an analysis method for nite wordlength implementations of time invariant digital lters [1,11,12] and periodically time-varying multirate systems [5,6,8]. It has been successfully applied to the simultaneous determination of the frequency response, the alias components, and the quantization noise of implemented lter banks. The following section 2 describes the electro-acoustic transmission system under consideration. The theoretical foundations of the measurement method are covered by section 3 along with two alternate interpretations given in section 4. Finally, some examples for the measurement of room inpulse responses with cheap sound transducers and background noise are presented in section 5.
2 Problem description 2.1 Overview
Fig. 1 shows an electro-acoustic transmission system for the measurement of an LEMS, consisting of ampli ers, digital-to-analog and analog-to-digital converters. A discrete measurement system MS excites the transmission system with a discrete input signal v(k) and records its response y(k). The design of the measurement system is the topic of this contribution. 2
D
v(k)
A
LEMS
MS y(k)
A
D
Fig. 1. Electro-acoustic transmission system for the measurement of a loudspeaker-enclosure-microphone-system (LEMS)
2.2 Linear approximation
The transmission of acoustic signals from the loudspeaker to the microphone is described by the linear wave equation [17]. However, if other sound sources are present in the LEMS, the measurement signal may be superimposed by additive background noise. Ampli ers and sound transducers can be modelled as linear systems if the distortions produced by them are negligible. This assumption is only valid for good quality (expensive) equipment. The converters are nonlinear systems by their very nature. They can be approximated by linear systems only for a high wordlength of the digital signals. If background noise, distortions of the analog equipment, and quantization e ects are neglected, the electro-acoustic system from Fig. 1 can be replaced by a linear model with an overall impulse response hlin(k) between input v(k) and output y(k) (see Fig. 2). v(k)
hlin (k)
y(k)
Fig. 2. Linear model
An estimation h^ lin(k) of this impulse response can be obtained with the Discrete Fourier Transformation (DFT), e.g. in the most simple case as ^ lin(k) = DFT ?1 h M
(
DFTM fy(k)g DFTM fv(k)g
)
(1)
A number of techniques exist, which di er in the choice of the input signal v(k), the DFT length M , and in more elaborate details of obtaining the estimate ^ lin(k) [16]. However, all these techniques rely on the assumption of linearity h and are blind to any of the nonlinear e ects described above. 3
2.3 Weakly nonlinear approximation
A detailed account of all possible nonlinear contributions in Fig. 1 leads to a general nonlinear model according to Fig. 3. Due to the variety and di erent nature of the nonlinear components, an exact nonlinear model S of the transmission system from Fig. 1 would be too complex to be handled eciently.
S
v(k)
y(k)
Fig. 3. Nonlinear model SL
yL (k)
v(k)
y(k) SN
n(k)
Fig. 4. Weakly nonlinear model
In order to circumvent this diculty, we approximate the nonlinear system
S by a weakly nonlinear model [11] as shown in Fig. 4. It consists of a parallel arrangement of a linear system SL and a nonlinear system SN. Again, it is not feasible to derive SL and SN analytically from the components of the transmission system. Instead, the characteristics of both systems have to be obtained experimentally from suitably chosen input signals v(k) and the corresponding responses y(k). A method for the determination of SL and SN from measurements of v(k) and y(k) is described in the following.
3 Description of the measurement method 3.1 Frequency response of the linear subsystem SL
The frequency response of the linear subsystem SL is determined such that the contribution of the nonlinear subsystem SN to the total output signal y (k ) becomes a minimum. This is achieved by expressing the output yL(k ) of SL in Fig. 4 by a convolution with the impulse response h(k) of SL: yL(k) = h(k ) v (k ). The minimization of the quadratic mean value of the output of the nonlinear system "(k) = E fj n(k) j2g is equivalent to determining the impulse response h(k) such that the expected value E fg n
( ) = E j y(k) ? h(k) v(k) j2
" k
4
o
(2)
becomes a minimum. By expanding the squared magnitude in (2) we obtain ( ) = E fy(k)y(k)g ?
" k
? +
+ 1 X
?1
=
+ 1 X
+ 1 X
?1
=
(
)E fy(k)v(k ? )g ?
h
( )E fy(k)v(k ? )g +
h
+ 1 X
( ) ( )E fv(k ? 1)v(k ? 2)g
h 1 h 2
?1 2 =?1 + 1 + 1 X X = y2 ? h ()'yv () ? h()'vy (?) + =?1 =?1 + 1 +X 1 X + h (2 )h(1 )'vv (2 ? 1 ): 1 =?1 2 =?1 1 =
(3)
where 'vv (k); 'vy (k) and 'yv (k) are the auto- and cross-correlation sequences between input and output. The unknown values h(k) of the impulse response follow by di erentiation of "(k ) with respect to h () from the minimum condition @"(k ) = 0; 8: (4) @h() The di erentiation of a function f (x) with respect to its conjugate complex variable x can be carried out in a formally correct way by the so-called Wirtinger calculus [2,13]. The result is ( )=
'yv
+ 1 X
?1
=
()
( ? ):
h 'vv
(5)
Finally, we obtain by Fourier transformation F an estimation H (ej ) = Ffh(k)g of the frequency response of the linear subsystem yv (ej ) : j
H (e ) = (6) vv (ej ) vv (ej ) = F f'vv ()g and yv (ej ) = F f'yv ()g are the power spectral density (PSD) and the cross power density of v(k) and y(k) respectively. For simplicity, we furthermore assume that vv (ej ) 6= 0 for all frequencies of interest. Since these power densities are not known in advance, they have to be determined from the signals v(k) and y(k). In general, a wide range of possi5
ble input signals and a variety of corresponding spectral estimation methods could be utilized. Here, we will focus on periodic signals v(k) and y(k) with period M . Then it is sucient to determine H (ej ) at M discrete frequencies = 2=M; = 0; ; M ? 1. The values H (ej ) of (6) can now be expressed by V () = DFTM fv(k)g and Y () = DFTM fy(k)g [1] H ej
(
) = E fY ()V (2)g ; = 0; ; M ? 1 : E jV ()j o
n
(7)
The expected values E fg are calculated from the measured input and output sequences by applying L sample sequences v(k); = 0; : : : ; L ? 1 of the same stochastic process to the transmission system and by averaging the results. Thus, the estimated value H^ (ej ) for the frequency response of the linear subsystem is obtained ?
LP1
^ (ej ) = H
Y =0 LP1
?
=0
()V()
jV()j2
;
= 0; ; M ? 1 :
(8)
The procedure is simpli ed, if the DFT-spectrum of the sample sequences has a constant magnitude jV()j = jV ()j. This can be achieved by deriving the sample sequences v(k) from an arbitrary periodic, deterministic signal v(k) by adding a stochastic phase term ( ) = v(k + M ) = v(k)ej' :
(9)
v k
The phase ' is a stochastic variable, equally distributed in the interval [?; +). For input signals of this kind, the computation of the estimated value simpli es to ^ (ej ) = 1 DFTM H LV ()
(
?
L X1 =0
y k e?j'
()
)
:
(10)
Fig. 5 shows the structure of the resulting measurement arrangement. In our application, a complex chirp signal ( ) = ej(=M )k2
v k
(11)
has been used as the deterministic component v(k) of the input. This signal achieves the minimum possible crest factor. 6
e
j
v(k)
1L L
S
1
DFT
ˆ j ) H(e
0
DFT
Fig. 5. Measurement arrangement for the frequency response of the linear subsystem SL. ( ) denotes the conjugate value. 3.2 Power spectral density at the output of the nonlinear subsystem SN
Once the frequency response of the linear subsystem SL is estimated, we can use this result to obtain the PSD nn (ej ) of the nonlinear subsystem SN. For the estimated output of the nonlinear system follows from Fig. 4 n ^ (k) = y(k) ? ^h(k) v(k) ; (12) where ^h(k) is the estimated impulse response of SL. The estimation of nn (ej ) is performed in the frequency domain by evaluation of L?1 jY()j2 ? H^ (ej ) 2 jV()j2 ^ nn (ej ) = M (L1? 1) X
=0
(13)
3.3 A new windowing technique
The proposed method for estimating the PSD was originally developed for the investigation of xed-point digital lters under the assumption of periodic quantization noise. However, it is well-known that averaged periodograms L?1 ^ nn (ej ) = M (L1? 1) jN()j2 X
=0
(14)
are biased in case of non-periodic noise sequences n(k), i.e. the estimate does not equal the true PSD even for an in nite number L of single measurements [7]. Due to the implicit windowing of n(k) with a rectangular window ( )=
r k
8 > < > :
1; k = 0 : : : M ? 1 0; else 7
(15)
the expected value of (14) yields
E ^ nn (e n
j
+ 2 1 ) = 2M nn (ej( ?)) R(ej ) d: ? Z
o
(16)
Thus, the PSD estimate is smeared since the spectrum R(ej ) of r(k) has high sidelobes. Many spectral estimation procedures such as the Welch-Bartlett method or the Blackman-Tukey method try to circumvent this problem by using di erent windows, e.g. the Bartlett, Hamming, or Hanning window. Nevertheless, all these procedures only use possibly overlapping signal blocks of length M for the estimation. As the window length M remains unchanged the sidelobes can only be reduced on the expense of the width of the mainlobe. Recently, in [9] a new windowing technique was proposed which overcomes most of these problems. Here, we will only give a brief survey of some major points of this new method. Obviously, the spectral properties of a window sequence can be improved by increasing the window length. Therefore, the method uses N subsequent blocks of n(k) and one common window f (k ) of length N M for estimating the PSD at M discrete frequencies = 2=M; = 0 : : : M ? 1. Thus, the PSD estimate (14) has to be modi ed as follows: L?1 ^ nn () = M (L1? 1)
?1
M X NX
=0
1
= M (L ? 1)
?
k =0
() ()
?
?
L X1 X1 N X1 M =0
L?1 = M (L1? 1)
k =0 =0
?1
=0
k =0
( + M )n (k + M ) e?j
f k
|
X X M
2
2k
f k n k e?j M
{z
2k
M
2
}
~()
n k
~()
n k e
2 ?j 2k M :
(17)
With F (ej ) = Fff (k)g, the PSD estimate then yields
E ^ nn (e n
j
+ 2 1 ) = 2M nn (ej( ?)) F (ej ) d: ? o
Z
(18)
Now, it is possible to deduce the following desirable properties of the window sequence f (k) from the above equations:
The PSD estimate at = 2=M should comprise all spectral components of the continous PSD in the frequency bin 2( ? 1=2)=M < 2( + 8
1=2)=M . The optimum window therefore has the spectrum Fopt ej
( )=
8 >
:
:
(19)
Consequently, the PSD estimate is the integral over the respective frequency bin at . Obviously, Fopt(ej ) cannot be realized with a nite window length. Therefore, we have to design a window f (k) which approximates the optimum window under the given constraints. The proposed method is capable of measuring the frequency response and the PSD simultaneously. For eciency reasons, it is therefore necessary that the windowing does not cause any interference of periodic output signal components which result from di erent periodic input signal components. This can be easily achieved by prescribing the following zeroes of the window spectrum: F ej
(
) = 0;
= 1 : : : M ? 1:
(20)
Finally, the windowing should preserve the total power, i.e. 1 M
+ ! 1 j ^ E nn (e ) = 2 nn (ej ) d : =0 ?
?
M X1
n
Z
o
(21)
From this constraint, the last condition immediately results: ?
M X1 =0
F
2
(ej( ?2=M )) =! M 2:
(22)
Evidently, the window sequence has to be a so-called M th-band lter which is well-known from lter bank theory [15]. In [9] a highly ecient method for designing such windows was introduced. Our PSD estimation procedure which is based on these window sequences shows an excellent performance even in case of a modest N . This will also be illustrated by the examples at the end of this paper.
4 Alternate interpretations The derivation of the frequency response of the linear subsystem in section 3.1 is quite straightforward and self-contained. In order to show interrelations with other results, we give here two alternate interpretations of the obtained estimate from the point of view of optimal ltering theory and estimation theory. 9
4.1 Wiener lter
The minimization of a quadratic error term according to (2) is also used for the derivation of the optimal reconstruction lter (Wiener lter). However, the scope of Wiener ltering is di erent from the measurement method considered here. Nevertheless, the measurement problem and the optimal lter problem can be discussed in the same framework. Both lead to the same minimization problem and consequently the solutions are formally identical. At rst we consider the general optimal reconstruction problem as shown in Fig. 6. For later reference to the measurement problem, the signals in Fig. 6 have been denoted with the same letters as in the weakly nonlinear model from Fig. 4. However, the meaning of the signals in the optimal lter problem is di erent. y(k) is an unknown input signal to the system S. Only the output v (k ) is available by observation. The problem is to nd a linear system H such that its output yL(k) is a reconstruction of the unknown signal y(k). The frequency response H (ej ) of the system H is determined by least squares minimization of the di erence ( ) = y(k) ? yL(k) = y(k) ? h(k) v(k)
n k
(23)
with respect to the coecients h(k) of the impulse response of H. The result H (ej ) is expressed by the power spectral density vv (ej ) of the observed signal v(k) and the cross power density yv (ej ), which re ects the available knowledge of the system S as [3,4] yv (ej ) : (24) H (ej ) = vv (ej ) y(k)
S
v(k)
H
yL (k)
n(k)
Fig. 6. Optimal lter problem
In order to establish a relation to the measurement problem considered previously, we redraw Fig. 6 and arrive at the modi ed arrangement of Fig. 7. Note that the problem description is identical to Fig. 6. The shaded area indicates that the available knowledge of the reconstruction problem consists of the correlation properties of the signals v(k) and y(k) represented by vv (ej ) and yv (ej ). Now we return to the measurement problem as described in section 2. Fig. 8 shows the arrangement for the detection of the linear subsystem of the weakly 10
vv
v(k)
(e j ),
(e j )
y(k)
S H
yv
n(k)
yL (k)
Fig. 7. Modi ed arrangement of the optimal lter problem
nonlinear model according to Fig. 4. The signal designations in Figs. 8 and 4 correspond to each other. The derivation of the frequency response H (ej ) of the linear subsystem by minimization of ( ) = y(k) ? yL(k) = y(k) ? h(k) v(k)
n k
(25)
has been shown in section 3.1. Here, we focus on the similiarities and differences of the measurement problem according to Fig. 8 and the modi ed arrangement of the optimal lter problem according to Fig. 7. An obvious di erence is that the role of v(k) and y(k) as input and output signals of S are interchanged. However, both problems share the property that the system S is described by the cross power spectrum yv (ej ). A second di erence is the purpose of the linear system H. It serves to reconstruct an unknown signal in the optimal lter problem and it describes the linear part of an unknown system in the measurement problem. However, the minimization problems (23) and (25) are identical for both cases, irrespective of the di erent meaning of the signals v(k) and y(k) in Figs. 7 and 8. vv
v(k)
(e j ),
S H
yv
(e j )
y(k)
n(k)
yL (k)
Fig. 8. Measurement problem
Since both the way of describing the system S and the minimization problems (23) and (25) formally agree, we can expect, that also the solutions for the frequency responses H (ej ) of the unknown systems H coincide. This is indeed the case, as is seen by comparison of the frequency response of the optimal lter according to (24) and the frequency response of the linear subsystem of (6). This shows that the foundation of the presented measurement method can also be interpreted by formal agreement with optimal lter theory. 11
4.2 Maximum likelihood method
The frequency response estimate can even be derived in a quite di erent way which ensures optimality in terms of estimation theory. This approach is called maximum likelihood (ML) estimate. We will derive this estimate for a single measurement and subsequently generalize the results. We start from a frequency domain description of the weakly nonlinear system according to Fig. 9 (compare Figs. 4 and 8). It can be formulated in terms of the DFT values of blocks of measured data with V () = DFTfv(k)g and Y (), N () alike. H () is the DFT of the impulse response of the unknown linear subsystem. ( ) = H ()V () + N ()
(26)
Y
v(k)
H
y(k) n(k)
Fig. 9. Additive noise description of the weakly nonlinear system
Under the assumption that the probability density function (PDF) of the disturbing signal component N () is known, we try to determine a frequency response estimate H^ ML () which most probably caused the observed output Y () of the system under test. For convenience we introduce vectors n, N, V, Y, H, and H^ ML which consist of the corresponding time or frequency domain samples. The ML estimate of the frequency response then can be found by maximizing the conditional PDF pYjH(YjH) with respect to the frequency response H, i.e. we have to determine ^ ML ) = max pYjH(YjH) : (27) pYjH (YjH n
o
H
Due to equation (26) the output Y () is distributed around the deterministic value H ()V () with the PDF of the disturbing signal N () [14], i.e. pYjH
(YjH) = pN(Y ? diag fH ()g V) = pN(N)
(28)
Thus, we have to maximize the noise PDF pN(N) with respect to the frequency response estimate. The required frequency domain PDF pN(N) can be determined from its time domain equivalent in a straightforward manner. For this purpose, the joint PDF of M subsequent samples of the underlying noise-like process n(k) = nr (k) + jni(k) is needed. Under quite general 12
conditions, this process is normally distributed and has an autocorrelation ( ) = 'nr nr () + 'ni ni () with 'nr nr () = 'nini (). With the autocorrelation matrix nn = E nnH the required joint PDF then can be written as 1 H ?1 n: pn (n) = M exp ? n (29) nn jdet nn j 'nn
o
n
Now, the PDF in the frequency domain can be determined easily using the DFT matrix WM and n = M1 WMH N. We obtain 1 NH W ?1 WH N: 1 exp ? (30) pN (N) = M M nn M jdet nn j M2 =: C |
{z
}
Instead of maximizing pN (N) it is convenient to introduce the so-called loglikelihood function L(H) = ln M jdet nn j pN (N) (31)
and maximize it with respect to H (). With N () = Y () ? H ()V () the unknown noise sequence can be eliminated and the log-likelihood function 1 M ?1 M ?1 c (Y ( ) ? H ( )V ( )) (Y () ? H ()V ())(32) L(H) = ? X
M2
X
=0 =0
follows. From the maximum condition 1 M ?1 c V () @L(H) = @H () H=H ^ ML M 2 =0
X
( ) ? H^ ML()V () =! 0 (33)
Y
the ML estimate immediately results: ^ ML ()V (): Y () = H
(34)
If we furthermore restrict ourselves to signals with a deterministic magnitude, i.e. (35) jV ()j2 = E jV ()j2 ; n
o
this estimate can be generalized for multiple measurements and we obtain ^ ML () = 1 2 E fY ()V ()g : HML () = E H (36) jV ()j n
o
13
Obviously, this is equivalent to the aforementioned estimate (7) if (35) is taken into account. Furthermore, it can be shown that the estimation error () = E fN ()=V ()g is given by () =
?
M X1 =0
a
@
ln pN(N) @H ( )
(37)
with constant factors a . This is a generalization of the Cramer-Rao bound [10,14] for the joint estimation of multiple parameters. Due to the equality in (37), the ML estimate (34) meets this bound. Thus, it is an ecient estimate in terms of estimation theory. Nevertheless, it has to be emphasized that this result only holds for multiple measurements if input signals with a deterministic magnitude according to (35) are used.
5 Results This section presents measurements of an LEMS both with high quality and low quality equipment. The LEMS was an anechoic chamber containing the measurement gear which caused some re ections. The high quality equipment was a studio microphone, a preampli er (Bryston Mod. PB4) and an active loudspeaker (GENELEC Mod. 1013A). For the low quality measurements, preampli er and loudspeaker were replaced by a one-chip ampli er (TBA 820M) with a cardboard mounted noname speaker (6 cm diameter). 20 log |H(f)|
0 −10 −20 −30 0
1
2 f [kHz]
3
4
Fig. 10. Frequency response using high quality audio equipment and standard measurement method
Fig. 10 shows the frequency response of the LEMS, measured with high quality audio equipment. Here, we used a standard measurement procedure according to (1) which assumes a strictly linear system behaviour. Fig. 11 shows the same measurement using low-cost equipment and again a standard measurement method. Due to this equipment a di erent frequency response has to be expected even though the room impulse response was not a ected by the change. Nevertheless, the measured frequency response shows 14
20 log |H(f)|
0 −10 −20 −30 0
1
2 f [kHz]
3
4
Fig. 11. Frequency response using low-cost audio equipment (standard measurement method)
an unexpectedly strong noise-like behaviour. This is due to the fact that conventional measurement procedures are not able to distinguish between linear and nonlinear system components, the latter being caused by distortions of the ampli er and speaker. Thus, the result is a superposition of both components. 20 log |H(f)|
0 −10 −20 −30 0
1
2 f [kHz]
3
4
Fig. 12. Frequency response using low-cost audio equipment with background noise (standard measurement method)
The frequency response measurement shown in Fig. 12 was performed with an additional broadband background noise 4dB below the signal level. It was produced by a noise generator fed into a guitar ampli er. The detrimental distortion of the measurement result by the superimposed noise source is obvious. In contrast, Fig. 13 shows the measurement results with our proposed method for the low-cost equipment with background noise as in Fig. 12. This method is capable of recording linear and nonlinear e ects separately. The method yields the frequency response of the linear subsystem as well as the PSD of nonlinear distortions including noise. Obviously, the true frequency response deviates from the original one (Fig. 10) due to additional linear distortions of the low-cost equipment. The remaining set of results demonstrates the suitability of the windowing technique described in section 3.3 for the discrimination between LEMS response and narrowband background noise. To this end, we repeat the measurement with high quality equipment similar to Fig. 10. During the measurement, an additional noise source was present in the LEMS which emitted a narrowband random signal with a signal-to15
20 log |H(f)|
0 −10 −20 1
2 f [kHz]
3
4
0 −10 −20 −30 −40 −50 −60 0
1
2 f [kHz]
3
4
PSD [dB]
−30 0
Fig. 13. Frequency response and PSD using low-cost audio equipment with background noise (proposed measurement method)
noise-ratio (SNR) of 3.5 dB. The noise signal was obtained from a broadband noise source by ltering with a digital FIR lter with a passband between 1 and 1.5 kHz. Fig. 14 shows the frequency response and the PSD obtained with the proposed measurement method and a rectangular data window. As a reference, the squared magnitude of the FIR form lter which was used for the generation of the narrowband noise is shown. Although averaging with L = 100 sample sequences was used (see (8) and (10) ), the measured PSD still di ers signi cantly from the true one. This deviation can be attributed to the non-periodicity of the noise signal. 20 log |H(f)|
0 −10 −20 1
2 f [kHz]
3
4
0 −10 −20 −30 −40 −50 −60 0
1
2 f [kHz]
3
4
PSD [dB]
−30 0
Fig. 14. Frequency response and PSD using high quality audio equipment with bandpass noise (proposed measurement method, rectangluar window, L = 100 sample sequences)
The same system was measured with the M th band window described in section 3.3 and with L = 10 and L = 100 sample sequences, respectively. Two 16
e ects can be observed from the results in Figs. 15 and 16. First, the corruption of the frequency response is evidently restricted to the actual frequency range of the bandpass noise. Second, the PSD of the nonlinear disturbances closely resembles the true PSD of the noise signal. The latter is given again by the squared magnitude of the FIR form lter. The obvious improvement in performance is due to the use of an appropriately designed data window. 20 log |H(f)|
0 −10 −20 −30 0
2 f [kHz]
3
4
3
4
form filter
PSD [dB]
0 −10 −20 −30 −40 −50 −60 0
1
1
2 f [kHz]
Fig. 15. Frequency response and PSD using high quality audio equipment with bandpass noise (proposed measurement method, M th band window, L = 10 sample sequences), squared magnitude of form lter 20 log |H(f)|
0 −10 −20 −30 0
2 f [kHz]
3
4
3
4
form filter
PSD [dB]
0 −10 −20 −30 −40 −50 −60 0
1
1
2 f [kHz]
Fig. 16. Frequency response and PSD using high quality audio equipment with bandpass noise (proposed measurement method, M th band window, L = 100 sample sequences), squared magnitude of form lter
Finally, we repeat the same experiment with the bandpass noise replaced by a sine wave of 1 kHz and an SNR of 4 dB. Fig. 17 shows the results for L = 100 sample sequences. The distortion of the measured frequency response of the linear subsystem by the sine wave is hardly visible. On the other hand, 17
we see that the PSD gives a perfect picture of the sinusoidal disturbance of the measurement. The absence of any spectral smearing or leakage clearly demonstrates the suitability of the M th-band window for spectral estimation. 20 log |H(f)|
0 −10 −20 1
2 f [kHz]
3
4
0 −10 −20 −30 −40 −50 −60 0
1
2 f [kHz]
3
4
PSD [dB]
−30 0
Fig. 17. Frequency response and PSD using high quality audio equipment with sinusoidal noise (proposed measurement method, M th band window, L = 100 sample sequences)
6 Conclusion We have presented a method for measuring the properties of the loudspeakerenclosure-microphone-system in the presence of nonlinear distortions and background noise. These disturbances have to be considered in the evaluation of the performance of low-cost audio products in oce or home environments. The measurement method is based on the analysis of implemented digital lters and multirate systems. Here, it has been extended to the investigation of the linear and nonlinear properties of electro-acoustic transmission systems. The fundamentals of the measurement method can be derived from a weakly nonlinear system model by minimization of a quadratic error criterion. Nonlinear e ects are characterized by the PSD of an additive noise signal. A specially adapted windowing technique provides extremely sharp estimates also for narrowband noise sources. The interrelations with optimal ltering theory and estimation theory have been shown by alternate interpretations of the measurement principle using the Wiener lter and the maximum likelihood approach. The results presented show the ability to clearly distinguish between the linear and nonlinear properties of an electro-acoustic transmission system. The 18
frequency response of the linear subsystem and the PSD which describes nonlinear distortions and background noise can be obtained separately.
Acknowledgement The authors wish to thank Eberhard Hansler for valuable hints, Thomas Klinger for preparing the results and Ivan Selesnick for his help with the manuscript.
References [1] Y. Dong. Ein neues Verfahren zur Messung der Eigenschaften schwach nichtlinearer Systeme. PhD thesis, Universitat Erlangen-Nurnberg, Lehrstuhl fur Nachrichtentechnik, Erlangen, 1990. [2] W. Fischer and I. Lieb. Ausgewahlte Kapitel aus der Funktionentheorie. Vieweg Verlag, Braunschweig, 1988. [3] B. Girod, R. Rabenstein, and A. Stenger. Einfuhrung in die Systemtheorie. B. G. Teubner, Stuttgart, 1997. [4] E. Hansler. Statistische Signale. Springer Verlag, Berlin, 1997. [5] F. Heinle and H. W. Schuler. An enhanced method for measuring the performance of multirate systems. In Proc. Int. Conf. on DSP, pages 182{187. Limassol, Cyprus, 1995. [6] F. Heinle and H. W. Schuler. Measuring the performance of implemented multirate systems. In Proc. Int. Conf. Acoustics, Speech and Signal Proc. (ICASSP 96), pages 2754{2757, Atlanta, USA, 1996. IEEE. [7] S.M. Kay. Modern Spectral Estimation { Theory & Application. Prentice-Hall, Englewood Cli s, USA, 1988. [8] R. Reng and H. W. Schuler. Measurement of aliasing distortions and quantization noise in multirate systems. In Proc. Int. Symp. on Circuits and Systems, pages 2328{2331, 1992. [9] H. Repp. Fensterung beim Rauschklirrmeverfahren. PhD thesis, Universitat Erlangen-Nurnberg, Lehrstuhl fur Nachrichtentechnik, Erlangen, 1998. [10] L.L. Scharf. Statistical Signal Processing: Detection, Estimation, and Time Series Analysis. Addison-Wesley, Reading, USA, 1991. [11] H. W. Schuler and Y. Dong. A new method for measuring the performance of weakly nonlinear and noisy systems. FREQUENZ, 44:82{87, 1990.
19
[12] H. W. Schuler and F. Heinle. Measuring the properties of implemented digital systems. FREQUENZ, 48:3{7, 1994. [13] H. Tietz. Funktionentheorie. In: Mathematische Hilfsmittel des Ingenieurs, vol. I, R. Sauer, I. Szabo (eds.). Springer Verlag, Berlin, 1970. [14] H.L. Van Trees. Detection, Estimation, and Modulation Theory, Part I. John Wiley & Sons, New York, 1968. [15] P.P. Vaidyanathan. Multirate Systems and Filter Banks. Prentice Hall, Englewood Cli s, USA, 1993. [16] N. Xiang. Using M-sequences for determining the impulse responses of LTIsystems. Signal Processing, 28:139{159, 1992. [17] L. J. Ziomek. Fundamentals of Acoustic Field Theory and Space-Time Signal Processing. CRC Press, Boca Raton, 1995.
20
