ap-5-11.dvi Acta Polytechnica Vol. 51 No. 5/2011 Synthesis of Room Impulse Responses for Variable Source Characteristics M. Kunkemöller, P. Dietrich, M. Pollow Abstract Every acoustic source, e.g. a speaker, a musical instrument or a loudspeaker, generally has a frequency dependent characteristic radiation pattern, which is preeminent at higher frequencies. Room acoustic measurements nowadays only account for omnidirectional source characteristics. This motivates a measurement method that is capable of obtaining room impulse responses for these specific radiation patterns by using a superposition approach of several measurements with technically well-defined sound sources. We propose a method based on measurements with a 12-channel indepen- dently driven dodecahedron loudspeaker array rotated by an automatically controlled turntable. Radiation patterns can be efficiently described with the use of spherical harmonics representation. We propose a method that uses this representation for the spherical loudspeaker array used for the measurements and the target radiation pattern to be used for the synthesis. We show validating results for a deterministic test sound source inside in a small lecture hall. Keywords: spherical harmonics, adjustable source directivity, room impulse response, linear room acoustics. 1 Introduction In order to determine roomacoustic parameters, e.g. reverberation time, clarity index or evenbinaural pa- rameters (IACC), room impulse responses are mea- sured with an omni-directional sound source, as re- quired by the ISO 3382 standard. These sound sources in general consist of several single loud- speaker chassis placed on a spherical array and ex- cited in a coherent way with the exact same sig- nal. Measured impulse responses in the room under test entirely describe the linear behavior for the ex- act combination of sound source position and micro- phonepositions. This canbeusedafterwardstomake the acoustic situation in this particular roomaudible (concept of auralization) but lacking the characteris- tic effect of the source radiation pattern. Methods were therefore developed to drive the single loudspeaker chassis of these compact spherical loudspeaker arrays with individual signals in order to directly approximate certain radiation patterns of target sound sources [7,8,10]. This directly implies measuring the room impulse responses with approxi- mated radiationpatterns, e.g. ofa speaker, an instru- ment, even though only a technical source is present during themeasurements. Ifweare interested in syn- thesizing the sound of various sources it becomes ob- vious that themeasurement time rises with each tar- get sourceand, of course, all target sourceshave tobe specified and measured in advance. This motivates a novel measurement and synthesis method that al- lows us to measure universal sets of room impulse responses that can be used to synthesize arbitrary radiation patterns after the measurement has been completed. The number of available loudspeaker chassis, and therefore the number of different basis radiation pat- terns, can be artificially increased by using several rotation and tilting angles of the loudspeaker ar- ray. The proposedmethod requires that the acoustic transfer characteristics in a room can be assumed as linear and time-invariant in order to use a superpo- sition approach. Hence, the reasonable limits will be studied and discussed. The proposed synthesis method is based on a de- scription of radiation patterns in the spherical har- monic domain. This enables us to model the radi- ation patterns of the source to be approximated as well as the spherical loudspeaker array used formea- surements on the same basis. 2 Method The proposedmethod can be divided into two parts, measurement and synthesis, which can also be en- tirely separated from each other. Both parts use the same calculus, and the inverse problem is also for- mulated in the same way. The spherical harmonics representation is used throughout. 2.1 Measurement of room impulse responses The core of the measurement consists of well knownimpulse responsemeasurementsof linear time- invariant (LTI) systems. This assumption holds for 69 Acta Polytechnica Vol. 51 No. 5/2011 most acoustical systems within certain limits. A de- tailed overview of these methods can be found in [6]. For each loudspeaker chassis of the array, the impulse response h(t) or its frequency representation H(ω)1 is obtained by using exponentially swept sines (chirps, sweeps) as excitation signals and using proper de- convolution techniques for the signal recorded by the microphones [2]. The chosen signal is advantageous for the given taskbymeansofnon-linearbehaviorde- tection possibilities. Furthermore, we employ a time saving approach that uses interleaved excitation sig- nals allowing several loudspeaker chassis to run at the same time, but also allowing us to perfectly sep- arate the responses is used as proposed by Madjak et. al [4]. For each of the M orientation angles of the loud- speaker array the N impulse responses aremeasured, one for each loudspeaker chassis, leading to a set of L = M · N room impulse responses hl(t) or Hl(ω) with l =1 . . . L. (1) Each response corresponds to a different radiation pattern. Figure 1 illustrates the method schemati- cally for an array of three loudspeakers: the impulse responses of each driver are measured in two orien- tations and are subsequently superposed. Fig. 1: Superposition of several, differently oriented loud- speakers chassis 2.2 Synthesis of target responses In order to approximate the target radiation pat- tern by the spherical loudspeaker array, complex and frequency dependent weighting factors wl are deter- mined to obtain the room impulse response hT(t) or the transfer function HT(ω) of the approximated tar- get radiation pattern by superposition. HT(ω) ≈ L∑ l=1 Hl(ω) · wl. (2) The superposition approach is only applicable if the roomcanbeconsideredasanLTI system. Linear- ity is in general not problematic for air-borne sound paths as in the room for moderate sound pressures, as in the case of standard room acoustic measure- ments. However, time-variances become problematic if the room changes significantly during a measure- ment session. Time variances in rooms are caused by temperature shifts, changes in humidity or light winds and circulations. In order to detect these vari- ances leading to errors in the ongoing synthesis, a concept is used as described in section 4.1. The radiationcharacteristicsof anacoustic source can be described by the directivity factor Γ [5]: Γ(θ, φ) := p(r, θ, φ) p(r, θ0, φ0) . (3) This gives the complex factor between the pressure p in a reference radiation angle (θ0, φ0) and the sound pressure in any direction (θ, φ). (r, θ, φ) are the ra- dius, the vertical and the horizontal angle of the com- mon spherical coordinate system. In general, p andΓ are complex and frequency dependent, but for better readability they are used without subscripts in the following. Since the directivity value can be regarded as a function which only depends on the radiation angle (θ, φ) and which is furthermore continuously inte- grable, it can be represented by a set of spherical harmonic coefficients Γ̂n,m, as shownbyWilliams [9]: Γ(θ, φ)= ∞∑ n=0 m∑ m=−n Γ̂n,m · Y mn (θ, φ) (4) where Γ̂n,m are frequency dependent and complex valued spherical harmonic coefficients, and Y mn are spherical harmonic base functions, which can be de- fined as: Y mn (θ, φ)= √ 2n +1 4π (n − m)! (n + m)! P mn (cosθ) e imφ. (5) Indices n and m denote the spatial periodicity of the function Y mn (θ, φ). They are called order n ∈ N0, and degree m ∈ Z : −n ≤ m ≤ n. P mn (μ) is an asso- ciated Legendre polynomial of the first kind. A de- tailed work explaining the characterization of acous- tic sources and radiation pattern with spherical har- monics is given by Zotter [10]. The radiation pattern of a real source has a fi- nite roughness over the surface. Therefore its char- acterization in the spherical domain can be limited to a maximum order Nmax, and the spherical har- monic coefficients can be summarized in a column vector [10]: Γ̂= [ Γ̂n,m ] (6) where 0 ≤ n ≤ Nmax and −n ≤ m ≤ n. 1The two representations have a fixed mathematical relationship. Hence they are used as synonyms in this work. 70 Acta Polytechnica Vol. 51 No. 5/2011 Eachof the above-mentioned L measured impulse responseswith the spherical loudspeaker array corre- spond to a certain source radiation pattern, which can be also written in such a vector d̂l. Let Γ̂T be the radiation pattern of the target to be synthesized, andwecan formulatebyanalogywith equation (2): Γ̂T ≈ L∑ l=1 d̂l · wl. (7) The vectors d̂l can be summarized in a matrix characterizing the radiationpatterns of the entire ex- tended array: D̂= [ d̂1 . . . d̂L ] . (8) Hence equation (7) can be extended towards a matrix formulation: Γ̂T ≈ D̂ ·w and w= ⎡ ⎢⎢⎣ w1 ... wL ⎤ ⎥⎥⎦ . (9) For the optimum weighting vector w we formulate, || D̂ ·w − Γ̂T || −→ min (10) leading to an inverse problem, that can be solved by using the Moore-Penrose pseudo inverse D̂ + [3]: w= ( D̂H D̂ )−1 D̂H︸ ︷︷ ︸ D̂ + · Γ̂T . (11) All quantities in equations (2) and (11) aremeasured quantities and are therefore subject to noise. In or- der to suppress the influence of these measurement uncertainties in the synthesis result, the rangeof pos- sible solutions is limited by Tikhonov regularization methods [3]: w= ( D̂H D̂+ I · ν )−1 D̂H︸ ︷︷ ︸ D̂ ⊕ · Γ̂T , (12) where I is the unit matrix of dimension L × L and ν is the so called Tikhonov parameter. 3 Implementation and instrumentation The measurement methods were implemented using MATLAB and the ITA-Toolbox. This toolbox is de- veloped at the Institute of Technical Acoustics and provides various tools for acoustics measurements and post-processing. Hence, the calculus for the inverse problem and the synthesis was also imple- mented in MATLAB. A complex calibrated instrumentation setup was used, consisting of the following elements (the num- bers correspond to the numbers in Figure 2). Fig. 2: Instrumentation setup • Measurement PC (1): (Windows XP, 32-Bit) ASIO driver, MATLAB and ITA-Toolbox. • D/A and A/D converter (2, 6): RME-Multiface II connected via RME-HDSP and Behringer ADA8000PRO8 connectedviaADATwithMul- tiface II. • Power amplifier (3): two 8-channel Stageline- IMG STA-1508. • Microphones (4): Radiation Pattern: Brüel & Kjær. Room impulse responses: Sennheiser KE4 elec- tret condensor. • Signal conditioning (5): Behringer ADA8000 PRO8, Preamplification and Phantom Power Supply for condensor microphones. Our spherical loudspeaker arrayconsists of amid- tone dodecahedron loudspeaker developed at the In- stitute of Technical Acoustics. The single loud- speaker chassis can be driven independently and the radiationpattern of each chassiswasmeasuredunder free-field conditions in the anechoic chamber with a controlled measurement scan unit. The results were transformedto the sphericalharmonicsdomain. This loudspeaker array was used along with a computer- ized turntable to allow arbitrary horizontal orienta- tion of the array in the room for measurements as shown in Figure 3. The array was inclined at an an- gle so that the elevation angles of the single chassis were equally distributed. Fig. 3: Dodecahedron array, devices for tilting and rotat- ing 71 Acta Polytechnica Vol. 51 No. 5/2011 The superposition method based on the Moore- Penrose pseudo inverse D̂ ⊕ of the radiationpatterns of the arraywas introduced in section 2.2. The error of this inversion is a good measure of the quality of the method to expect for ongoing calculations with reasonable target responses. In the ideal case, the residual matrix R̂= | Î − D̂ ⊕ · D̂| (13) would be the zero-matrix. The energy of its columns ε̂n,m corresponds to the error that arises when syn- thesizing several spherical harmonics Y mn . Figure 4 shows this error over frequency2. Frequency in Hz O rd e r o f S p h e ri ca l H a rm o n ic s Mean square error of synthesis in dB 400 630 1000 1600 2500 4000 6300 20 18 16 14 12 10 8 6 4 2 0 −60 −50 −40 −30 −20 −10 Fig. 4: Error of synthesizing spherical harmonics As canbe seen, the possible order of the spherical harmonics for synthesized target sources rises with frequency, and the error of synthesis rises as well. The low number of possible reproducible orders for low frequencies can be explained by the fact that the single loudspeakers donot have a dominant radiation pattern for low frequencies themselves. The synthesis error is causedby limited resolutionof theorientation angles in vertical direction of the single chassis, as this angle could not be adjusted automatically with the given measurement setup. 4 Experimental results In order to evaluate the proposedmethod a compar- ativemeasurementwas conducted. The roomchosen for the measurements was an easily accessible lec- ture hall in the Institute of Technical Acoustics with a mean reverberation time of approx. 0.9 seconds at mid frequencies. Two main measurements were conducted in this room: one measurement with the spherical loudspeaker array, and the other measure- mentwith a loudspeaker of a certain target radiation pattern, which was also used as target response for synthesis. Figure5 illustrates themeasurement setup inside the lecture hall. The upper picture shows the spherical loudspeaker array on the left side and the bottom picture shows the target loudspeaker in the same position in the room. Additionally, the refer- ence dodecahedron loudspeaker (right side) was used in afixedposition, as canalsobe seen in thepictures. 4.1 Detection of time variances As mentioned above measurements with a reference loudspeaker are conducted for each orientation angle of the spherical loudspeaker array. The results of the reference loudspeaker are used for a correlation anal- ysis of the impulse responses in the timedomainwith a mean impulse response. Figure 6 shows this cor- relation coefficient. The dotted line marks the time when the room was briefly entered to replace the ar- ray by the target loudspeaker. It is obvious that the acoustic behavior of the room changes over time. At the beginning, after the personnel has left the room, the changes are greater than at the end. This can be explained by the fact that the room still needs some time to completely settle down after objects havemoved. Measurements arechosenatmeasurementtimeswhere the timevari- ances are low. In this case we chose 100 measure- ments as input for the ongoing synthesis. Fig. 5: Measurement setup for comparative meausre- ments 20 40 60 80 100 120 140 160 180 200 0.7 0.75 0.8 0.85 0.9 0.95 1 Correlation coefficient (reference: all measurements) Index of measurement position C o rr e la tio n c o e ff ic ie n t 5e+002 Hz [1] 1e+003 Hz [1] 2e+003 Hz [1] 4e+003 Hz [1] Fig. 6: Correlation analysis to detect time variances 2Only the maximum error of all degrees in one order is shown. 72 Acta Polytechnica Vol. 51 No. 5/2011 4.2 Results The upper image inFigure 7 shows themeasurement with the real source and the synthesis result in the time domain, and the lower image zooms into the range of the first reflections in the room impulse re- sponse. The results lookvery similar in this represen- tation. The frequency domain is plotted in Figure 8. The results show good agreement in the range from 300Hz to 1.5 kHz. The chosen cut off frequency lim- its the result at the low end. The deviation above 1.5 kHz grows over frequency, which is in correspon- dence with the results shown in Figure 4. Fig. 7: Measured and synthesized impulse response (time domain) Fig. 8: Measured and synthesized impulse response (fre- quency domain) 5 Conclusion We have proposed a measurement method for a spe- cial set of room impulse responses and synthesis in a post-processing step for room impulse responses of arbitrary target radiation patterns. In addition an approach was introduced to fully separate measure- mentandsynthesisby transformingthemeasurement results into a universal representation. Since the method assumes linear time-invariant systems, this assumption was studied and a measure to quantify and monitor time variances was used based on mea- surementswith a reference loudspeaker in a fixed po- sition. The method was validated in a small lecture hall using a 12-channel dodecahedron spherical loud- speaker arraywith automatically adjustable orienta- tionangles to virtually increase thenumber of drivers and therefore the number of different radiation pat- terns. The results obtained by the synthesis of the proposed method were compared to measurements with the sourcewhichwas also used as target for the synthesis. The frequency range is limited towards higher frequencies at around 3 kHz with the given measurement setup. As themain idea of this workwas to develop such a measurement method, there are still some limita- tions to overcome in future work. In order to cover the entire audible frequency range from 20 Hz to 20kHztwomodifications seemreasonable: the spher- ical array should be replaced by a high-tone version for the higher frequency range, and the vertical res- olution of the spherical array needs to be increased. Acknowledgement Research described in this paper was supervised by Prof. Dr. Michael Vorländer. The authors would like to thank the electrical and mechanical workshop at the Institute of Technical Acoustics for their support for the special measurement devices. References [1] Dietrich, P., Masiero, B., Mueller-Trapet, M., Pollow, M., Scharrer, R.: MATLAB Toolbox for theComprehension ofAcousticMeasurement and Signal Processing. DAGA, 2010. [2] Farina, A.: Advancements in impulse response measurements by sine sweeps. AES 122nd Con- vention, Vienna, Austria, 2007. [3] Kress, R.: Inverse Problems. Institute for Nu- merical and Applied Mathematics, TU Goettin- gen. [4] Majdak,P.,Balazs,P.,Laback,B.: MultipleEx- ponential Sweep Method for Fast Measurement of Head-Related Transfer Functions, J. Audio Eng. Soc., 2007. [5] Mechel, F. P.: Formulas of Acoustics. Springer, 2002. 73 Acta Polytechnica Vol. 51 No. 5/2011 [6] Müller, S., Massarani, P.: Transfer-Function Measurement with Sweeps, Journal of the Au- dio Engineering Society, 2001. [7] Pollow, M., Behler, G.: Variable Directivity for Platonic SoundSourcesBased onSphericalHar- monics Optimization, Acta Acustica, 2009. [8] Warusfel, O., Derogis, P., Causse, R.: Radiation Synthesis with Digitally Controlled Loudspeak- ers. AES, 1997. [9] Williams, G.: Fourier Acoustics, Sound Radia- tion and Nearfield Acoustical Holography. Naval Research Laboratory, Washington, D.C. [10] Zotter, F.: Analysis and Synthesis of Sound- Radiation with Spherical Arrays Dissertation, Institut für Elektronische Akustik, Graz, 2009. About the authors Martin Kunkemöller received his diploma degree in electrical engineering and information technology at RWTH Aachen, Germany, in 2011. He was em- ployedas researchassistantat the InstituteofTechni- calAcoustics ofRWTHAachenuntil April 2011with the focus on developing and evaluating the proposed method. He is member of the German Acoustical Association (DEGA). Pascal Dietrich received his diploma degree from the faculty of electrical engineering andcomputer sci- ence at RWTH Aachen University in 2006. In 2007 he joined the Institute of Technical Acoustics as a PhD student and researcher in the field of electro- acoustics, transfer-path analysis and structure-borne sound prediction always with the focus on measure- ment and modeling uncertainty. He is a member of the Audio Engineering Society (AES) and the Ger- man Acoustical Association (DEGA). Martin Pollow received his diploma degree in elec- trical engineering from RWTH Aachen University, Germany, in 2007. He is currently employed as a researcher and enrolled as PhD student at the Insti- tute of Technical Acoustics of RWTH Aachen Uni- versity. His focus of research encompasses complex sound fields, airborne sound radiation, array systems and analytical models. He is member of the German Acoustical Association (DEGA). Martin Kunkemöller Pascal Dietrich Martin Pollow E-mail: martin.kunkemoeller@akustik.rwth-aachen.de Institute of Technical Acoustics RWTH Aachen University, Germany 74