Indoor acoustic simulator for ultrasonic broadband signals with Doppler effect

Applied Acoustics 97 (2015) 140–151

Contents lists available at ScienceDirect

Applied Acoustics journal homepage: www.elsevier.com/locate/apacoust

Technical Note

Indoor acoustic simulator for ultrasonic broadband signals with Doppler effect Daniel F. Albuquerque a,⇑, José M.N. Vieira b, Sérgio I. Lopes b, Carlos A.C. Bastos b, Paulo J.S.G. Ferreira b a b

CI&DETS, ESTGV, Polytechnic Institute of Viseu, Portugal Signal Processing Lab, IEETA/DETI – University of Aveiro, Portugal

a r t i c l e

i n f o

Article history: Received 1 October 2014 Received in revised form 6 February 2015 Accepted 14 April 2015

Keywords: Ultrasonic signals Acoustic simulator Impulse response Ultrasonic reflection Doppler effect

a b s t r a c t In this paper we present an indoor acoustics simulator that is able to deal with high bandwidth ultrasonic signals. The simulator can handle the multiple reflections present in a room, which could be produced by the room walls and by the objects inside it. Moreover, it takes into account the transducer’s beam patterns as well as the ultrasonic wave propagation losses. Additionally the simulator may be used to simulate the Doppler effect for any type of movement. The simulator was implemented using Matlab and is freely available, so that it can be easily modified by the user. Field tests have shown that the simulator is very accurate, turning it into a valuable tool to study and compare different ultrasonic pulses or ultrasonic location systems under controlled conditions. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction An ultrasonic (US) room simulator is a valuable tool to study and compare the performance of different location algorithms and acoustic pulses in a controlled environment. By using a simulator it is very easy to change the parameters that affect the response of the acoustic system and to test extreme situations that are not easy to replicate in the real field. Moreover, it is much faster to use a simulator than to setup and carry out real field tests in order to test different algorithms and ultrasonic signals under the same circumstances. Commercially available acoustic simulators, like Ramsete, ODEON or CATT-Acoustic, are very expensive and complex and since they are not open source, changes are not possible [1]. However, there are examples of freely available simulators, [2–4] but they are implemented to simulate very specific situations, most of them in the audible range. Moreover, the majority of them do not take into account the Doppler effect or implement it for very simple situations such as uniform movement. Note that the Doppler effect has a considerably stronger impact on US signals than on radio-frequency (RF), as it will be seen in Section 5. Indoor ultrasonic acoustic propagation presents problems similar to those experienced in electromagnetic waves propagation, such as multiple reflections, Doppler effect and propagation ⇑ Corresponding author. Tel.: +351 232480500; fax: +351 232424651. E-mail address: [email protected] (D.F. Albuquerque). http://dx.doi.org/10.1016/j.apacoust.2015.04.010 0003-682X/Ó 2015 Elsevier Ltd. All rights reserved.

attenuation [4,5]. Due to these problems, a good model for the ultrasonic channel is needed in order to properly design ultrasonic communication or location systems. The proposed simulator is based on changing and simplifying some techniques that were used in the acoustics and radio-frequency fields. The approach presented here allows to produce a very fast and easy to use simulator to evaluate and compare different ultrasonic location systems and ultrasonic communication pulses. 2. Simulator architecture The problem that we want to address is schematically presented in Fig. 1. Each source s receives an electric signal xs ðnÞ and transforms it into an acoustic wave that travels through the Room Channel. Moreover, this Room Channel represents a room where the acoustic waves can propagate and reflect in the walls before arriving at the receiver position. Therefore, the receiver block combines each received wave into the signal yðnÞ. The simulator architecture can be represented by the block diagram of Fig. 2, where hs ðnÞ represents the impulse response (IR) that simulates the effect of the transducer and the acoustic channel on the input signal xs ðnÞ. This approach brings the first two simplifications to the models to be used in the simulator: 1. all signals are digital and therefore discrete (time and value); 2. the source, the receiver and the room effects are considered linear (represented by an IR).

141

D.F. Albuquerque et al. / Applied Acoustics 97 (2015) 140–151

Fig. 1. The ultrasonic room simulator architecture.

Fig. 3. The division of hs ðnÞ in three main IRs. Fig. 2. Block diagram of the simulator.

At the end the simulator provides to the user the IRs for each individual source s; hs ðnÞ. The user may compute a set of impulse responses at different receiver’s positions and compare the effect of each source and input signal xs ðnÞ in the received signal yðnÞ. Each hs ðnÞ can be decomposed in three main IRs which represent the effect of the source, ts ðnÞ, the receiver, r s ðnÞ, and the sound propagation, ps ðnÞ, see Fig. 3. The source s IR, t s ðnÞ, models the effect of the transmission chain (sound card; cables; amplifier; speaker; etc.), in other words, it models all the path from the digital world up to the acoustic wave generated by the speaker. In contrast, the receiver IR, r s ðnÞ, models the effect of the reception chain (microphone; amplifier; cables; sound card; etc.), in other words, it simulates the path from the received acoustic wave by the microphone up to the digital world. Finally, the sound propagation IR, ps ðnÞ, models the acoustic wave propagation in a specific room for a pair of receiver and source positions (attenuation; distance delay; etc). 2.1. Transmitter/receiver model In order to obtain the source and the receiver IR two approaches can be employed. If the hardware manufacturer gives the IR (or frequency response) of their components, then it will be possible to obtain the IR for all the chain. If not, measuring the IR of the system with calibrated hardware could be the only way. Note that, for a complete hardware description (e.g.: speaker and microphone) IR measures for each azimuth and elevation angle must be obtained. However, this needs to be carried out just once because the measured IR are considered to be independent of the room conditions. Both IR must be provided to the simulator, a source IR for a given source output direction and a receiver IR for a given receiver input direction. The receiver or the source direction is defined by the azimuth angle h, and the elevation angle / as presented in Fig. 4. 2.2. Sound propagation model

Fig. 4. Receiver or source direction.

for the source and the receiver. Therefore, the proposed simulator, uses a physic based model instead of a statistically based model [1,6,7]. In contrast to the statistically based models the physic based models allow to simulate the behavior of a system in a particular room or in a particular situation. In order to reduce the simulation time and the user intervention some approximations and simplifications were made in the sound propagation model, wherever possible, without compromising the closeness to reality. On that account, the simulator incorporates the following simplifications and specifications, that will be detailed later: 1. The room ambient conditions (temperature, pressure and humidity) are constant over time and space. 2. The wave propagation attenuates due only to dispersion and energy absorption and it is considered to be a function of the signal frequency. 3. The wind effect and the gradients of temperature were not taken into account. The sound speed, c, depends on temperature, pressure and density of the air [8]. Nevertheless, at typical room conditions, the sound speed may be assumed to be independent of small pressure and density variations [8]. Therefore, c can be written only as a function of the temperature:

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi T þ1 273:15

c ¼ 331:5641 This is the most crucial element of the simulator. However, it is not possible to measure the IR for all the ambient conditions and for all possible source/receiver distances, being impractical to model the sound propagation by applying the same approach used

ðm=sÞ

ð1Þ

where T is the air temperature in °C. This speed implies that the ultrasonic wave takes some time to travel the distance, d, from the source up to the receiver. Due to the

142

D.F. Albuquerque et al. / Applied Acoustics 97 (2015) 140–151

LD ðnÞ ¼ Source

m¼0 m–n

d

d0

Air Attenuation

0.25 0.2 0.15 0.1

0

10

20

30

40

50

60

70

nm

ð2Þ

;

where D is a real number that represents the filter delay and M the filter order. Lagrange fractional delay filters only operate properly for frequencies up to a quarter or a third of the sampling frequency, depending on the filter order, which can be pointed out as the main Lagrange fractional delay filter’s disadvantage [9]. For that reason, the simulator must use a sampling rate around four times higher than the maximum frequency of the signal in order to reduce the delay and amplitude errors. Considering that the majority of the ultrasonic transducers manufacturers only provide the pressure amplitude, P0 , at a certain known distance d0 and ultrasonic waves propagate as spherical waves [8], as shown in Fig. 5, the pressure amplitude at distance d (in meters) relates with P0 as:

Fig. 5. Ultrasonic wave propagation.

0.05

M Y Dm

80

90

100

Frequency (kHz) Fig. 6. Total attenuation in air for a receiver at 4 m normalized by the attenuation at 1 m with the following ambient conditions: 20 °C of temperature, 30% of relative humidity and 1 atm of pressure.

PðdÞ ¼

d0  P0 : d

ð3Þ

This equation only takes into account the wave attenuation due to dispersion. However, in [10], it is pointed out that the propagation losses in ultrasonic frequencies have a strong impact in the wave attenuation, especially for distances above a few meters. This attenuation can be expressed as [10]: aðdd0 Þ

need to accurately simulate this delay, the input signal xs ðnÞ may need to be delayed by a non integer number of sampling periods. For that reason, Lagrange fractional delay filters [9] were used to simulate the fractional delay. The choice of this particular type of filters relies on the fact that they produce a finite IR, LD ðnÞ, that can be easily integrated with the other IRs of the simulator. Furthermore, they were also chosen because their filter coefficients can be readily updated by a close-form equation [9] given by:

PðdÞ ¼ P0  10

where a is the air absorption coefficient in dB/m and it can be computed by the 9613-1 ISO norm [11] using only parameters easily obtained, such as, the temperature, the pressure, frequency and relative humidity. Combining the attenuation due to dispersion (3) and due to propagation losses (4), the resultant acoustic wave pressure amplitude at distance d (in meters) is given by:

0

−3

10

Error Amplitude

Amplitude

10

−2

10

−4

10

−6

10

ð4Þ

20

−4

10

−5

10

−6

10

−7

−30

−20

−10

0

10

20

30

Samples (a) Impulse response with 64 samples length.

10

0

20

40

60

80

100

Frequency (kHz) (b) Error obtained by the truncation process.

Fig. 7. Impulse response with 64 samples length that represents the total air attenuation and the error obtained by the truncation process.

Fig. 8. Spectrum band error reduction algorithm.

143

D.F. Albuquerque et al. / Applied Acoustics 97 (2015) 140–151 0

−3

10

Error Amplitude

Amplitude

10

−2

10

−4

10

10−6

−4

10

−5

10

−6

10

−7

−30

−20

−10

0

10

20

30

Samples (a) Impulse response with 64 samples length.

10

0

20

40

60

80

100

Frequency (kHz) (b) Error obtained by the truncation process.

Fig. 9. Impulse response with 64 samples length that represents the total air attenuation and the error obtained by the truncation process using the proposed algorithm.

PðdÞ ¼ P0 

aðdd0 Þ d0  10 20 : d

ð5Þ

As the result of considering that the room ambient conditions are constant over the simulation, the absorption coefficient in air is only a function of the wave frequency f. Therefore, the atmospheric attenuation faced by an ultrasonic wave at d meters from the source is given by:

AT ðd; f Þ ¼

aðf Þðdd0 Þ d0  10 20 d

ð6Þ

As a result of this, the AT ðd; f Þ will be only a function of the wave frequency f for a fixed distance d and therefore can be considered as a linear filter, which results into an impulse response, aðnÞ. To exemplify the process used in the simulator to obtain the impulse response aðnÞ, a receiver at 4 m from the source and the following ambient conditions: a temperature of 20 °C; a relative humidity of 30%; and a pressure of 1 atm was considered. Under this constraints the total attenuation in air, AT ð4; f Þ, is depicted in Fig. 6. Our approach to obtain the AT ðd; f Þ IR is to use the inverse Fourier transform. Which leads us to an infinite impulse response. Therefore, the resultant IR must be truncated in order to obtain an IR with only a few samples. Fig. 7 presents the IR with 64 samples obtained by this process for the example presented in Fig. 6, it also shows the error produced in the frequency domain. Note that this IR must be delayed due to causality however this delay can be easily subtracted from the total delay resulted from propagation (in this example the IR must be delayed of 32 samples which represents a distance of around 5.5 cm for a sampling frequency of 200 kHz). We can see in Fig. 7b that the resultant error is very low, moreover it can be considerably reduced by increasing the IR length. However, this approach will increase the minimum

distance between the source and the receiver that can be simulated. Nevertheless, the error can be reduced if we discard the frequency response accuracy in some of the frequency components, for example in the audible range (up to 20 kHz). In order to increase the accuracy of the frequency response in one part of the spectrum while giving up elsewhere, we can use the algorithm presented in Fig. 8 based on the Papoulis–Gerchberg algorithm [12,13]. The algorithm starts by computing the total attenuation frequency response AT ðd; f Þ using (6) in K frequency points f k , obtaining the discrete Fourier transform AðkÞ, Note that, the frequency points f k must be computed for the discrete Fourier transform bins, moreover the number of frequency points K must be greater than the desired IR length L to ensure proper frequency response interpolation. After that, the IR ai ðnÞ is computed by the inverse discrete Fourier transform of Ai ðkÞ, resulting into a length K IR. The IR ai ðnÞ is then truncated using a rectangular window wL ðnÞ with length L and the corresponding frequency response is computed using the discrete Fourier transform. This process is repeated until the error in the desired band B is less than a given value . At the beginning of each iteration the frequency points in the desired band B of Ai ðkÞ are restored to the proper values AðBÞ. By using the proposed algorithm it is possible to considerably reduce the error for a predetermined band. Using the same example and applying the algorithm to the band from 20 kHz to 50 kHz and running it one hundred times, it is possible to increase the accuracy of the total attenuation frequency response ten times for the same IR length (see Fig. 9). 2.3. Surfaces model For the surface modeling (e.g.: walls, floor and objects) the simulator uses a polygonal based modeling, widely used in 3D

Fig. 10. Modules and methods interaction in the implemented simulator.

144

D.F. Albuquerque et al. / Applied Acoustics 97 (2015) 140–151

computer graphics [14]. Due to the typical geometric shape of the room surfaces we chose quadrilaterals instead of the conventional triangles as the basic element. Therefore, most of the room surfaces can be represented by making use of only one quadrilateral. However, some of them must be represented has a group of quadrilaterals connected to each other to form a particular geometric shape. Furthermore, many of the room models produced by specific software are based in polygonal modeling which may be easily adapted to be used into the simulator. Nevertheless, this type of geometry considerably simplifies the process of computing the surface reflections. However, to properly model curved surfaces, it is necessary to use a huge number of quadrilaterals. 2.4. Multiple reflections model At this point, it is possible to model the transmission, propagation and reception chain using separate impulse responses and without considering the reflections. As a result of this, it is already possible to properly model the direct path from the sent signal xs ðnÞ to the received signal ys ðnÞ. Moreover, it is also possible to represent all the surfaces in the room that may produce wave reflections. Therefore, all the main requirements to introduce the reflections into the simulator are fulfilled. However, the theory of sound propagation in closed spaces can be very hard to apply and model and it may bring an important drawback associated to it: the computation time. Nevertheless, ultrasonic waves have a huge advantage when compared to sound, which is that its maximum wave length is about 1.7 cm (for 20 kHz). As most of the obstacles usually present in a room are larger that this wavelength, we can consider only specular reflections. Moreover, the ultrasonic scattering energy is very low when compared to scattering energy in the audible range. Due to its complexity, it was also discarded by the simulator. These simplifications considerably reduce the simulation time and the user intervention without compromising the closeness to reality. On that account, the simulator presents the following characteristics and specifications: 1. It is considered that the sound will propagate along sound rays. 2. All possible reflections up to a predefined order are computed and they are considered to be only specular in a plane surface with a specific IR (previously measured) associated to it. Diffraction and scattering are ignored due to its modeling and computational complexity. 3. All the possible received sound rays are subsequently validated by applying the ray-tracing method. When an ultrasonic wave is emitted in a common room, the wave travels radially in all directions. When this wave finds an obstacle (e.g. a wall), the propagation direction changes and the ultrasonic wave reflects on the obstacle [8]. However to simulate the ultrasonic propagation it is better to consider that it travels like rays and the wave front is flat [15]. In other words, the ultrasonic wave behaves, in the presence of an object, as a ray of light behaves in the presence of a mirror. This means that it is possible to apply the well known Snell Law, from the optics theory, to determine the reflection wave direction

hi ¼ hr :

ð7Þ

It is important to notice that, in the simulator, the transmitted wave is discarded due to its low energy, therefore it is considered that the transmitted wave is completely absorbed by the object [15]. For most of the solid materials (e.g. wood, metal, plastic, glass, etc.) the reflected wave has almost the same energy as the incident wave. Nevertheless, some soft materials in a common room can absorb some wave energy and consequently the reflected wave has less energy than the incident wave (e.g. carpets, cork, curtains,

etc.). In addiction to these, if the energy lost in the reflection depends on the frequency of the incident wave, an IR of the wall reflection must be considered. Nevertheless, similar to what happens with the source and the receiver, this has to be carried out only once and it is considered to be independent of the room conditions. In order to model the reflection in the surfaces there are two classical geometric methods: the Ray Tracing Method and the Image Source Method [16]. In both methods, the accuracy increases as the frequency increases and they present very good results for high audible frequencies [16]. Therefore, at the ultrasonic range it is expected that they produce very accurate results. An hybrid method that combines the Ray-Tracing Method and the Virtual Sources Method was used in the simulator. Firstly, all the virtual sources are created and then each incoming ray is validated by using the ray-tracing method. 3. Implementation The implementation of the simulator was carried out in Matlab.1 Since Matlab is proprietary software an effort was made, during the simulator implementation, to only use typical Matlab functions that can be found in several compatible freeware packages (such as GNU Octave2). This way, the implemented simulator can be used almost for free. Moreover, this tool will allow the user to have complete control over the simulation and to understand its behavior. The simulator was implemented as a modular program in order to simplify the implementation of future extensions and the substitution of code by any user. The modular structure allows the user to easily modify any module in order to obtain the desired behavior. The simulator can be divided into six different modules and six different methods for interaction between the modules. A module can be a simple data structure or a complex object depending on its goal. The methods are designed to obtain the information from the modules and/or create new modules. The interaction between the different modules and methods is schematically represented in Fig. 10. 3.1. Walls array module The Walls_Array module is a simple data structure that contains the information of each surface of the room. Each surface in the array is represented by a reflection coefficient and by the four edges of the quadrilateral as discussed previously in Section 2.4. This module is created and modified by the Add_Wall method, which creates the array and fills it with the data. 3.2. Transfer function module The Transfer_Func module is an object that, for a particular source or receiver, outputs to the user the IR for a given azimuth and elevation angle. Depending on the origin of the information this module can have different implementations, however, this must me transparent to the user. In the present simulator implementation, this module needs three inputs: – the direction: ðh; /Þ; – the operation frequency band: B; – the sampling frequency: f s .

1 Matlab is a numerical computing environment developed by MathWorks, it allows the implementation of complex algorithms providing several tools to reduce the implementation time. 2 GNU Octave is a freeware software very similar to Matlab and most of the algorithms implemented in Matlab can be run in GNU Octave without almost any changes.

145

D.F. Albuquerque et al. / Applied Acoustics 97 (2015) 140–151

(a) No reflections.

(b) Up to three reflections.

Table 1 Characteristics of the computer and software used in the simulator performance tests. Description

Value

Matlab version Operating system Processor RAM memory Graphic card

2012b (64 bits) Microsoft Windows 7 SP1 (64 bits) Dual Core 2.5 GHz (Intel T9300) 4 GBytes Nvidea Geforce 8600 M GT

This module is used by the Add_Source and Add_Receiver methods.

The Sources_Array module is a simple data structure that contains the information of each source placed in the room. Each source in the array is represented by: – the position in the room (x; y; z) coordinates; – the direction in Cartesian coordinates (ex ; ey ; ez ); – the Transfer_Func module associated to the source.

Mean simulation time (min.)

This module is created and modified by the Add_Source method, which creates the array and fills it with the data.

5 reflections 4 reflections 3 reflections 2 reflections 1 reflection

0

3 2 1 0

0

1

2

3

4

5

6

3.4. Room representation file module The Room_File module is a simple data structure stored in a file that contains:

3.3. Sources array module

10

3 sources 2 sources 1 source

4

Fig. 13. Simulation time in function of the maximum number of successive reflections to consider for different number of sources in the room.

– the impulse response: ts ðnÞ or rs ðnÞ; – the index of the IR’s first sample.

1

5

Maximum number of successive reflections

From these inputs the module provides two outputs:

10

Mean simulation time (min.)

Fig. 11. Output example of the Display_Room method.

−1

10

– the information of all the room surfaces; – all the sources in the room (real and virtual). This module is created by the Get_Room method, which by processing the Walls_Array and the Sources_Array modules can create the file and fill it with the data for a maximum reflection order given by the user. This module has a specific method to show graphically the information stored in the file: the Display_Room method. Fig. 11 presents two example outputs of this method. In Fig. 11a the room without reflections is presented and therefore only the real source marked with a circle is shown. Fig. 11b presents the same example but considering a maximum of three successive reflections. The virtual sources, that produce the reflections are marked with a light gray circles and the real source with a dark gray circle. Note that, this file is independent of the room ambient conditions and the receiver position. Therefore, the resultant file can be used for any room ambient conditions and for any receiver position without the need to create it again.

−2

10

1

2

3

4

5

6

7

8

9

10

11

12

Number of Surfaces Fig. 12. Simulation time as a function of the number of room surfaces and the maximum number of successive reflections to consider.

3.5. Receiver module The Receiver module is similar to the Sources_Array module. The biggest difference is that it only has one receiver instead of having an array of them. Therefore, the data structure contains:

146

x (n)

D.F. Albuquerque et al. / Applied Acoustics 97 (2015) 140–151

...

n-M

n-M+1

...

n-2

n-1

n

n+1

hs(n) (m) hs(n)(M-1)

...

hs(n)(2)

hs(n)(1)

hs(n)(0)

account the ultrasonic wave attenuation presented in Section 2.2. In the present simulator, the implementation of this module needs the following inputs:

...

– – – –

+

y (n)

...

n-M

n-M+1

...

n-2

n-1

n

n+1

...

the the the the

distance between the source and the receiver; room temperature, humidity and pressure; operation frequency band, B; sampling frequency, f s .

From these inputs the module provides the following outputs:

Fig. 14. Block diagram for Doppler effect processing.

– the impulse response, aðnÞ; – the index of the IR’s first sample (causality problem).

– the position in the room (x; y; z) coordinates; – the direction in Cartesian coordinates (ex ; ey ; ey ); – the receiver Transfer_Func module.

3.7. Get impulse response method The Get_Impulse method is the most complex method in the simulator due to the fact that it collects and combines all the information from other modules (Room_File, Receiver and Propagation) to compute the IR from each source, hs ðnÞ. In the first stage, by using the Receiver module, it validates all the sources (real and virtual) in the Room_File module. For each valid source, it obtains the source output direction and the receiver input direction. It computes the total attenuation due the reflection coefficient of each surface and the total traveling distance. In the

This module has a creator method, Add_Receiver, which creates it and fills it with the data. 3.6. Propagation module The Propagation module is an object that for particular ambient conditions provides the user with the IR for a given distance between the source and the receiver. This IR only takes into

Microphone Amplifier

d

Data Acquisition

PC

Microphone

Speaker

Buffer

Fig. 15. Measure system chain for practical testes.

SPL/V @1m re 20µPa

100

80

60

40

20

0

20

40

60

100

80

Frequency (kHz) (a) Speaker.

(b) Frequency response of the speaker. Fig. 16. Speaker Kemo L10 and its frequency response.

Atmospheric measure point

5.5 m

Microphone

Speaker 1.45 m

Floor Fig. 17. Scheme of the setup used to test the simulator for one single reflection.

147

Normalized Amplitude

D.F. Albuquerque et al. / Applied Acoustics 97 (2015) 140–151

1 0.8 0.6 0.4 0.2 0 15

15.5

16

16.5

17

17.5

18

18.5

19

19.5

Normalized Amplitude

Time (ms) (a) Impulse Response obtained in the simulator. 1 0.8 0.6 0.4 0.2 0 15

15.5

16

16.5

17

17.5

18

18.5

19

19.5

Time (ms) (b) Real Impulse Response obtained by correlation. Fig. 18. Impulse response comparison for a single reflection.

second stage, it computes the IR from the Transfer_Func module of the source, ts ðnÞ, the IR from the Transfer_Func module of the receiver, rs ðnÞ, and the IR from the Propagation module, aðnÞ. At the end it groups together all the sources that represent the same real source and computes the resultant impulse response for each real source, hs ðnÞ.

4. Performance

Fig. 19. Simulation environment used to test the simulator in a complex room.

To evaluate the simulator performance several simulations were carried out, in the same computer, with different parameters. For these tests a computer with the characteristics presented in Table 1 was used. In order to minimize external interference from other programs running in the computer, the Matlab process priority was set to High. In addition to this, at least 10 trials were carried out for each test.

Fig. 20. Measurement points in the test of a complex room.

148

D.F. Albuquerque et al. / Applied Acoustics 97 (2015) 140–151

Simulator IR

Real IR Normalized Amplitude

P1

Normalized Amplitude

1 0.8 0.6 0.4 0.2 0

0

5

10

15

20

25

30

35

40

1 0.8 0.6 0.4 0.2 0

0

5

10

Time (ms) Normalized Amplitude

Normalized Amplitude

0.8 0.6 0.4 0.2 0

0

5

10

15

20

25

30

35

40

Normalized Amplitude

Normalized Amplitude

1

0.6 0.4 0.2 0

5

10

15

20

25

30

35

40

Normalized Amplitude

Normalized Amplitude

0.8 0.6 0.4 0.2 0

5

10

15

20

25

30

35

40

Normalized Amplitude

Normalized Amplitude

0.8 0.6 0.4 0.2 0

5

10

15

20

25

35

40

30

35

40

30

35

40

30

35

40

0.4 0.2 0

0

5

10

15

20

25

1 0.8 0.6 0.4 0.2 0

0

5

10

15

20

25

1 0.8 0.6 0.4 0.2 0

0

5

10

15

20

25

Time (ms)

1

0

30

0.6

Time (ms)

P5

40

Time (ms)

1

0

35

0.8

Time (ms)

P4

30

Time (ms)

0.8

0

25

1

Time (ms)

P3

20

Time (ms)

1

P2

15

30

35

40

Time (ms)

1 0.8 0.6 0.4 0.2 0

0

5

10

15

20

25

Time (ms)

Fig. 21. Impulse response comparison for a complex room shape. Where Px are the points marked in Fig. 20.

4.1. Room complexity impact In this test, the simulation time was measured using one to twelve surfaces and considering a maximum of one to five successive reflections. The simulation times are presented in Fig. 12. From this figure it can be seen that the simulation time is approximately exponential with the number of surfaces (>3) in the room for any number of maximum successive reflections. Considering three

maximum successive reflections the simulation time starts at about 2 ms for three surfaces and goes up to 50 ms. Another important point is that the simulation time increases considerably with the maximum number of successive reflections. It starts around 2 ms for one reflection and ends at 20 s for five reflections for a room with twelve surfaces. Nevertheless, for the majority of the ultrasonic applications, a maximum of three or four successive reflections is more than enough to present good results.

149

D.F. Albuquerque et al. / Applied Acoustics 97 (2015) 140–151

(7.5; 3; z)

START/END POINT

(0; 0; 0)

(2; 0; 0)

Source (0; -1; 0)

y

z

x

(-7.5; -3; z) Fig. 22. Simulation environment to test the Doppler effect for a receiver with uniform and circular movement. All the coordinates are in meters.

In this test, the simulation time was measured using one to three sources and considering a maximum successive reflection number of one to six. The results from this test are presented in Fig. 13. From this figure it can be seen that the simulation time is approximately exponential with the number of reflections and is approximately linear with the number of sources. These results are consistent with the previous test and the simulator implementation, where each source is treated separately and independently. 5. Simulating the Doppler effect The Doppler effect is produced by the receiver and/or source movement. Moreover, the Doppler effect can be characterized by a compression or expansion of the acoustic signal, in the time domain, proportional to the source/receiver relative speed [7]. The Doppler effect has a stronger impact on ultrasonic propagation than in electromagnetic wave propagation due to the lower speed of sound. For example, a person traveling at 1 m/s produces a maximum Doppler shift [7] of 118 Hz on a 40 kHz carrier. For FM radio reception to have an equivalent Doppler shift, the person needs to travel at more than 330 m/s (1188 km/h). The Doppler effect can be simulated by resampling the channel impulse response over time. Therefore, the impulse response will be different for each instant n: ðnÞ

hs ðn; mÞ ¼ hs ðmÞ

ð8Þ

ðnÞ

where the hs ðmÞ is the impulse response hs ðmÞ for the instant n for source s. Fig. 14 presents the block diagram of the Doppler effect implementation in the simulator. For each instant n the source and the receiver position must be determined and the IR must be computed using the simulator. After that the IR can be applied to the input signal xðnÞ to obtain the output sample yðnÞ. For instant n þ 1 the source and the receiver position must be updated in the simulator and a new IR must computed and applied to xðn þ 1Þ to obtain the output sample yðn þ 1Þ. As a result of this, for each output sample an IR must be computed with the simulator which leads to a extremely high computation time. However, it produces a very accurate result with multiple reflections and it can be applied to any kind of movement and room scenarios. 6. Simulator results The system measurement chain presented in Fig. 15 will be used to perform the practical tests. The distance between the microphone and the speaker could be fixed or variable depending on the type of measurement. An PC with Matlab was used to receive and process the ultrasonic waves and to communicate with

the data acquisition hardware, namely the DT9836S-6-2-BNC from Data Translation. The Brüel & Kjær 4954-A microphone was used to receive the ultrasonic signals, it is calibrated, operates up to 100 kHz and its manufacture provides a well-known frequency response. The speaker that was used is a Kemo L10 (see Fig. 16a), which has an excellent frequency response between 15 kHz and 45 kHz, nevertheless it can produce waves up to 80 kHz, as can be seen in Fig. 16b. 6.1. Single reflection To test the simulator behavior a very simple example with a single reflection will be described. The simulator impulse response will be compared with the real impulse response obtained under similar circumstances. Fig. 17 shows the scheme of the measurement setup system used. The microphone and the source were placed at fixed positions at 1.45 m from the floor and 5.5 m apart. An environment as similar as possible to the real one was created in the simulator. The real impulse response was obtained by correlation. The resultant IR is compared in Fig. 18 with the impulse response obtained by the simulator. One can see that the impulse responses are very similar, proving that the simulator presents a very high degree of similarity to the real environment described in Fig. 17. 6.2. Simulation of a complex room In order to test the simulator in a complex room, the room environment presented in Fig. 19 was used. The room has the following external dimensions 6  10  3 m and in the middle of the room there is a 3 m long wall. The walls and the ceiling are made of medium-density fiberboard (MDF) and the floor of concrete. The speaker was placed in the inside wall corner, near the ceiling pointing in the direction of the room corner (3, 0, 0) as shown in Fig. 20. The impulse response was measured in five different places in the room, as presented in Fig. 20, at 1.5 m from the floor. A measuring

Relative Speed (m/s)

4.2. Number of sources impact

<

>|<

first turn

>

second turn

5 2.5 0 -2.5 -5 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

2.2

2.4

Time (s) Fig. 23. Relative speed between the receiver and source for a receiver with uniform and circular movement (moving at 10 m/s).

150

D.F. Albuquerque et al. / Applied Acoustics 97 (2015) 140–151

(a) Single sinusoidal signal with 30 kHz of frequency.

(b) 5 sinusoidal signals with frequencies equally spaced from 29 to 31 kHz.

(c) 1 Kbps BPSK modulated signal with a central frequency of 30 kHz. Fig. 24. Simulation of the Doppler effect without reflections (in the first half of the time) and with reflections (in the second half of the time) in three different signals and for a receiver with uniform and circular movement.

tape was used to set the positions of the microphone. Therefore the microphone position presents errors of the order of a centimeter. The first four points are inside the same room division of the speaker and the last point is out of the speaker’s line-of-sight. Nevertheless, this point can receive reflected waves from walls. The results for this five points are presented in Fig. 21, where the left column presents the impulse response obtained from the simulator and the right column the real impulse response obtained by correlation. In spite of some differences, especially in the impulses amplitude, it can be said that the impulse response obtained in the simulator matches very well the impulse response obtained in the field. A careful analysis of the results reveals that the main echos of the simulator IR (the ones with a considerable amplitude in

comparison to the direct wave) have a similar amplitude and arrive at the same time that the ones in the real IR. The amplitude difference could be related to some reflection coefficient difference over the wall material and the roughness of the concrete floor. Another important issue that the simulator does not take it into account can be observed in the results for P5, which is the ultrasonic propagation through wall. In P5 all the impulse response spikes before 25 ms are related to ultrasonic waves that cross the inside wall. The distance between the speaker and the microphone is around 2.55 m resulting into a 7.5 ms time delay therefore the first spike in the real IR arises from the direct path. Nevertheless, these small differences do not compromise the realism of the simulator, on the other hand, they reinforce the idea

D.F. Albuquerque et al. / Applied Acoustics 97 (2015) 140–151

that the simulator is a very important tool to estimate the ultrasonic signal propagation in a particular room. 6.3. Doppler effect To test the Doppler effect a receiver moving with an uniform and circular movement was used. For that purpose the simulation setup shown in Fig. 22 was used. In this setup there are six surfaces with a reflection coefficient of 0.5 (four walls, floor and ceiling), the floor was placed at z ¼ 1 m and the ceiling at z ¼ 2 m. The receiver performs a circular movement with constant speed of 10 m/s and 2 m of radius. Moreover, to evaluate and understand the reflections impact, the receiver makes two complete turns, which takes around 2.51 s, in the first turn the simulator considers no reflection and in the second turn the simulator considers a maximum of three successive reflections. The source was placed at 1 m from the coordinate system origin in (0; 1; 0) and its direction was set to the y coordinate. The same speaker of the previous tests (Fig. 16a) was used. As a consequence of the source position, the relative speed will change with time according to Fig. 23 which presents a speed variation of ±5 m/s. Fig. 24a presents the Doppler effect on a single sinusoidal signal, Fig. 24b on a signal composed by 5 sinusoidal signals and Fig. 24c on a binary phase shift keying (BPSK) modulated signal. When multiple reflections are introduced, the resultant signal will be a mixture of different sinusoids with different frequencies due to reflected waves arriving to the receiver with different relative speeds. At 10 m/s this mixture can have sinusoids with frequencies shifts in the range of ±840 Hz for 29 kHz, ±870 Hz for 30 kHz and ±900 Hz for 31 kHz due to the Doppler effect. This effect can be seen in Fig. 24a where the signal frequency changes between 29.13 kHz and 30.87 kHz and Fig. 24b where the signal frequency changes between 28.16 kHz and 31.90 kHz. Doppler effect was not compared with field tests due to the complexity of the setup, the measuring process and the lack of specific hardware.

7. Conclusion and future work This paper describes the design and implementation of an ultrasonic room simulator with the following characteristics: – designed for wideband ultrasonic signals; – the signal bandwidth can be specified in order to improve its accuracy within that band; – takes into account the air attenuation and losses by providing the desired room ambient conditions, such as, temperature and humidity; – any type of flat surface can be implemented and its reflections can be modeled through an IR; – makes use of the speaker and/or microphone beam patterns; – simulates multiple signal sources at the same time; – simulates the Doppler effect for any type of movement. On the other hand it presents the following limitations: – only computes the reflections up to a certain order, initially specified by the user; – does not take into account gradients of temperature nor wind effects; – discards the wave propagation through the walls (transmitted wave); – only simulates a single receiver at a time; – rounded surfaces are difficult to model; – the simulation of Doppler effect is a time consuming process.

151

The results presented in Section 6 show that the models used in the simulator are valid, very accurate and realist. The implemented simulator is already in a mature state with acceptable performance and good accuracy. However, there are several points that can be improved both in terms of performance and accuracy. In the present version of the simulator the virtual source validation by ray tracing takes more than half of the total simulation time for some configurations and therefore this is a key point to address in order to improve the performance. Another important point is the user interface, which forces the user to know how to use Matlab. Therefore some work could be done in this field making the simulator more user friendly. In the simulator, the room is implemented by setting all the surfaces coordinates, which can make the room definition very time consuming. Consequently, the simulator integration with other simple and intuitive room design software (like SketchUp3) will simplify the room implementation process. From the accuracy point of view, it will be important to take into account the time and space variation of some of the ambient variables. Another important point in the accuracy field, especially for low US frequencies is the propagation waves through the room surfaces. This simulator is available at the URL: https://code.google.com/ p/ultrasonic-simulator/. Acknowledgments The Portuguese Foundation for Science and Technology (FCT) through the project PEst-OE/CED/UI4016/2014, and the Center for Studies in Education, Technologies and Health (CI&DETS). References [1] Rindel JH. The use of computer modeling in room acoustics. J Vibroeng 2000;3(4):41–72. [2] Holm S. Simulation of acoustic fields from medical ultrasound transducers of arbitrary shape. In: Nordic symposium in physical acoustics. Ustaoset (Norway); 1995. [3] Holm S. Ultrasim – a toolbox for ultrasound field simulation. In: Nordic matlab conference; 2001. [4] Wabnitz A, Epain N, Jin C, Schaik Av. Room acoustics simulation for multichannel microphone arrays. In: International symposium on room acoustics. Melbourne (Australia); 2010. [5] Sonnek P, Rice SV. Synthesizing the acoustic Doppler effect in software. In: 48th annual southeast regional conference. Oxford (Mississippi): ACM; 2010. p. 1–5. [6] Stojanovic M, Preisig J. Underwater acoustic communication channels: propagation models and statistical characterization. Commun Mag IEEE 2009;47(1):84–9. [7] Fontán FP, Espiñeira PM. Modeling the wireless propagation channel: a simulation approach with MATLAB, vol. 5. Wiley; 2008. [8] Kinsler LE, Frey AR, Coppens AB, Sandersq JV. Fundamentals of acoustics. John Wiley & Sons; 1982. [9] Valimaki V, Haghparast A. Fractional delay filter design based on truncated lagrange interpolation. IEEE Signal Process Lett 2007;14:816–9. [10] Phillips S, Dain Y, Lueptow RM. Theory for a gas composition sensor based on acoustic properties. Meas Sci Technol 2003;14(1):70–5. [11] ISO 9613-1. Acoustics – attenuation of sound during propagation outdoors – Part 1: Calculation of the absorption of sound by the atmosphere; 1993. [12] Papoulis A. A new algorithm in spectral analysis and band-limited extrapolation. Circ Syst 1975;22(9):735–42. [13] Gerchberg RW. Super-resolution through error energy reduction. Opt Acta: Int J Opt 1974;21(9):709–20. [14] Agoston MK. Computer graphics and geometric modeling: implementation and algorithms. London: Springer; 2005. [15] Everest FA. Master handbook of acoustics, 4th ed. Tab Electronics; 2000. [16] MCKerrow PJ, min Zhu S. Modelling multiple reflection paths in ultrasonic sensing. In: International conference on intelligent robots and systems, vol. 1; 1996. p. 284–91.

3 SketchUp from Trimble is one of the most popular and used 3D modeling tools, it was created to be used not only by professionals but also by amateurs with no knowledge in 3D modeling.