Generating Sequences of Acoustic Scintillations

Summary Spatial and temporal inhomogeneities in temperature and wind velocity a ﬀ ect sound propagation resulting in amplitude and phase ﬂuctuations called scintillations. A computationally e ﬃ cient method is presented to generate sequences of scintillations. The method, already used in the ﬁeld of wireless communication to predict the performance of wireless communication links, could be used in the ﬁeld of acoustics to create more perceptually valid auralizations. A Gaussian spectrum and a spherical wavefront is considered, but the method can also be used in combination with other spectra like the Von Karman spectrum as well as plane waves. Two examples are given, one is a pure tone a ﬀ ected by the scintillations and the other is an auralization of an aircraft ﬂy-over. The e ﬀ ect of the transverse speed of the source is demonstrated as well. terms


Introduction
In an outdoor situation, spatial and temporal variations in temperature and wind velocity cause small changes in the refractive index. As wavesp ass through the atmosphere, the index-of-refraction variations in effect cause scintillations, i.e. fluctuations in the receivedi ntensity of the wave.Scintillations affect both sound and electromagnetic waves. Theya re am ajor limitation for astronomical observations using Earth-based telescopes and also reduce performance of wireless communication systems. Scintillations can also be noticed when hearing sound emitted by asource at alarge distance, likefor example by an aircraft or distant windturbine [1].
Within the SONORUS project [2], atool wasdeveloped to auralize the sound of aircraft flying overurban areas [3]. Scintillations affect the audible sound, and it is presumed that including scintillations increases the perceivedr ealism of auralizations. Earlier work introduced acoherence factor to include the loss of coherence due to phase fluctuations [4,5,6]. Fluctuations due to turbulence have also been included in auralizations by simulating the amplitude modulations that were observed in measurements [1,7].
In this paper am ethod is presented to generate time series of sound pressure fluctuations caused by line-ofsight propagation through aweakly turbulent atmosphere. Novelinthis field, the method includes both log-amplitude Received12October 2016, accepted 2January 2017. and phase fluctuations. The presented method is based on [8] where it wasused to predict the performance of wireless communication links, and [9] where it wasused to determine the influence of turbulence on the performance of abarrier.The method wasused by the authors to increase the perceptual validity of auralizations [10,11].

Rytova pproximation
Va riations in temperature and wind in both position r and time t cause variations in the refractive indexfield n(r,t). We are interested in howthese variations affect wave propagation and follow [12] and [8], butinstead of electromagnetic waves, we consider sound waves. We consider for nowspatial variations only,and as astarting point we use the Helmholtz equation with pressure p,wavenumber k,and refractive-indexfield with mean value n 0 = E[n(r)] = 1and first-order perturbation n 1 (r) 1. Merging equations (1) and (2) results in ∇ 2 + k 2 (n 0 + n 1 (r)) 2 p(r) = 0.
Forw eak fluctuations, an approximation to equation (3) for small n 1 is used. The Rytovsolution to equation (3) is where ψ 0 is the complexp hase of the unperturbed wave in free space, and ψ 1 and ψ 2 respectively first-order and second-order complexphase perturbations. We are interested in the effect of first-order perturbations n 1 ,o nt he sound pressure, and therefore write ψ = ψ 0 + ψ 1 .T he refractive index n is written in terms of an average n and fluctuation n 1 ,with As derivedin [12], ψ 1 satisfies the following integral equation where G(r − r )isthe free field Green'sfunction. By iteration aseries solution can be obtained. Forthe first iteration we set ψ 1 = 0inside the integral and obtain the first Rytov solution where p 0 (r)i st he unperturbed sound pressure field. The sound pressure after the first iteration is then p(r) = e (ψ 0 +ψ 10 ) = p 0 (r)e (ψ 10 ) .
The first-order complexphase perturbation ψ 1 can be understood as asum of waves, generated at various points r throughout the scattering volume V .The strength of each of these wavesisproportional to the product of the unperturbed field term p 0 ,and the refractive-indexperturbation δn at apoint r [8].

Amplitude and phase fluctuations
We nowwant to findexpressions for the log-amplitude and phase fluctuations, and will use Rytov'sfirst solution. We approximate the refractive-indexfluctuation as δn = 2n 1 + n 2 1 2n 1 (9) and write where A and S are respectively the amplitude and phase of the fluctuating field p(r), and obtain for the first order perturbations In this expression χ and S represent respectively the logamplitude fluctuation and phase fluctuation. By applying the central limit theorem to the first Rytovs olution, it follows that the complexp hase follows a normal probability distribution [8]. This is an important result to keep in mind when generating sequences of fluctuations.

Amplitude and phase covariance
The log-amplitude and phase fluctuations are considered to be the result of arandom temperature fluctuation field. Ac haracteristic of ar andom function or field is its correlation function [13]. The spatial correlation function of ar andom field f (r), as function of distance r = r 2 − r 1 between observation points r 1 and r 2 ,isdefined as In ah omogeneous and isotropic random field the correlation function C(r)d epends only on the distance r = r 2 − r 1 between the observation points and not the path r = r 2 − r 1 [14]. Note that at this point, the atmosphere is assumed frozen in time, i.e., variations are only spatially, not temporal.
We would liketoobtain expressions for the covariance functions of the log-amplitude and phase fluctuations. A specificp art of the turbulence spectrum can be approximated with aGaussian correlation function where σ 2 µ is the variance of the dynamic refractive index, x = r 1 − r 2 the distance between twopoints in space and L the correlation distance or length [12].
We shall nowconsider aline-of-sight situation where d is the distance between the source and areceiverpair along the wave propagation direction, and ρ the spatial separation of the receivers transverse to the wave propagation direction.
If the correlation length L is much smaller than the Fresnel zone size of the sound √ λd,then the log-amplitude and phase variance scale with σ 2 χ = σ 2 S ∼ k 2 d [12] and the variances of the fluctuations are givenby [15] Fors pherical wavest he covariances of the fluctuations, B χ (ρ)and B S (ρ), normalized to their variances, are given by where and erf is the error function. The covariance functions of the fluctuations B χ (ρ)and B S (ρ)are thus Still assuming Taylor'shypothesis regarding frozen turbulence, we can perform as pace-to-time conversion of the correlation and covariance functions to obtain C(τ)a nd B(τ)respectively.The turbulence correlation time is given by and the transverse time lag by where v ⊥ is the transverse speed, corresponding to e.g. the mean speed at which the field is carried by the wind transverse to the wave propagation direction. In this paper we will continue to use the covariance for spherical wavesand aGaussian spectrum, because the Gaussian spectrum is the simplest model to work with and computationally least demanding, butthe method can also be used with covariance functions that describe other turbulence spectra, likee.g. the Vo nKarman spectrum [16].

Propagation in the turbulent atmospherea sa multichannel
The fluctuations in the atmosphere are temporal and/or spatial. Therefore, to model sound propagation of asignal x(t)t hrough such an atmosphere at ime-varying channel is necessary.The receivedsignal y(t)consists of aline-ofsight contribution and additional contributions from scattering, together forming amultichannel [8]. Ignoring beam spreading, we can write this as where α sc n (t)i st he time-varying scintillation sequence representing the effect of the pressure fluctuations on the nth-multipath component, and τ n the propagation delay of the nth component relative to the propagation delay in an undisturbed atmosphere. Assuming the spread in propagation delay overthe channels is small compared to the inverse of the signal bandwidth, so that Because the attenuation is very similar for the different multipaths, we can write that Therefore, the receiveds ignal is obtained by shifting the emitted signal x(t)with atime-varying propagation delay τ(t), and multiplying the result with at ime-varying gain α sc (t).

Design of scintillation sequence
We will nowd esign as equence of scintillations and for that we need to knowt he statistical distribution and power spectral density |H B (f )| 2 of the desired sequence.
As mentioned before the fluctuations are Gaussian-distributed. We can therefore generate asequence with the correct distribution and power spectral density by convolving ar andom unit variance Gaussian signal z(t)w ith the impulse response h B (t)ofthe designed filter.
The power spectral density |H B (f )| 2 of ar andom sequence forms aFourier pair with its autocorrelation function R(τ)through the Wiener-Khinchin theorem 1 This filter has zero phase and is non-causal. To create a causal filter with constant group delay α we can shift the peak in the impulse response by adding alinear-phase factor corresponding to 90 degrees The impulse response of the filter is finally obtained through the Inverse Fourier Transform

Discrete time
We nowconvert from continuous to discrete time The linear-phase factor exp(−j2πfα)isthen givenby where M 1 is the length of the desired impulse response. The frequencyresponse of the filter is where F {} is the Discrete Fourier Transform (DFT). The desired impulse response is real-valued. Therefore we have because of Hermitian symmetry Scintillations are finally obtained through the convolution of the impulse response h B [n]w ith Gaussian white noise z[n] The first M 1 /2samples would have to be dropped because of the filter delay.Ab lock diagram of the steps is shown in Figure 1.

Apply scintillations to signal
Nowt hat we can generate sequences of fluctuations, we need to apply these to ac arrier signal x(t)o rx [ n ]. According to equation (22) the log-amplitude fluctuations can be applied in amultiplicative manner,which is indeed the case of as ignal with as u ffi ciently small bandwidth, likef or example am onochromatic signal. However, the variance of the fluctuations is frequency-dependent, and since in practice broadband signals are commonly used, a method is sought for applying the fluctuations to abroadband signal.
One possible method would be to decompose the input signal in complexe xponentials with the Inverse Discrete Fourier Transform, and apply per complexexponential the log-amplitude and phase fluctuations which, when combined, can be written as ac omplexe xponential as shown in equations (10) and (11).
Acomputationally less demanding method is to instead filter the carrier signal with multiple band-pass filters, and proceed with applying fluctuations to each of the bandpass filtered signals and combining the resulting signals. Scintillations would be computed for the center frequencies of the bands. The amplitude fluctuations could be applied through multiplication and the phase fluctuations could be converted to time-delay fluctuations (see equation (35))and applied with aV ariable Delay Line (VDL).

Scintillations as time-variant filter
Am ore efficient method to takei nto account the frequency-dependence of the scintillations, would be to instead convolvethe carrier signal with atime-variant filter that is updated fast enough to represent the fluctuations.
The fluctuations χ and S,a te ach time step computed for M 2 frequencies, can be merged into ac omplexp hase (see equation (11)) Adding al inear-phase factor as wasd one with h B [n]a nd taking the Inverse Discrete Fourier Transform will result in an impulse response for every instance in time n.C onvolution of the carrier signal x[n]with this time-variant filter will result in alog-amplitude and phase modulated signal. The first M 2 /2samples will have to be dropped to correct for the filter delay. Because the fluctuations are relatively low-frequent they can be sampled at amuch lower sample frequencythan the carrier signal. Foracorrelation length L of 10 meters and a transverse speed v ⊥ of 2meters per second the correlation time would be 5seconds. If we would sample the sequence of fluctuations at 5t imes the maximum bandwidth of the filter,which is proportional to the inverse of the correlation time τ 0 , f s = 5/τ 0 (34) the required sample frequencywould be 1Hertz. In practice one might want to interpolate the impulse response as it changes overt ime, butb ecause the phase also changes with time this can be problematic. Aside from converting to am inimum-phase representation we can also create al inear-phase filter for which the magnitude changes with time butthe phase does not. If the phase fluctuations are linear-phase, then theyc an be applied to the carrier signal x[n]w ith aV ariable Delay Line. In the Gaussian model the phase fluctuations scale as σ 2 S ∼ f 2 . As catterer affects all frequencies, and therefore the fluctuations of different frequencies move together,r esulting in alinear-phase. We therefore write the phase fluctuation dS(t)asapropagation delay fluctuation With the Gaussian turbulence spectrum it is straightforward to separate the covariance B(τ)i nto the correlation Vol. 103 (2017) C(τ)a nd variance σ as defined in respectively equations (16) and (14).T hat wayw ec an compute afi lter h C [n] that depends solely on the correlation C(τ), and scale this afterwards with the correct variances σ χ and σ S .T he advantage is having to compute only one sequence of fluctuations, which is determined entirely by the covariance and z [n], and can then be scaled with the (frequencydependent)variances σ χ and S.

Moving source
Thus fari tw as assumed neither the source nor the receiverwere moving, and that the fluctuations arise due to sampling aturbulent field that is moving by with constant speed at afixedreceiverposition.
If the field is transported faster,t he spatial fluctuations are sampled faster,and therefore the scintillations χ[n] = S[n]are compressed in time corresponding to fluctuations with relatively higher frequencies.
In case the turbulent field is transported by wind and neither source nor receiverare moving, the source-receiver distance remains constant. If we instead consider amoving source of which the emitted sound samples the turbulent field, then the source-receiverdistance d generally changes and thus also the variances σ χ and σ S .
We will nowconsider amoving source. The transverse velocity is in this case the velocity component of the moving source, perpendicular to the wave propagation direction. This is the vector rejection of the source velocity vector v s from the unit vectorô that represents the orientation from source to receiver. The transverse speed is the norm of this vector The vector projection of the source velocity on the wave propagation direction is the component that is relevant for the Doppler shift. Therefore, as the Doppler shift is maximum the modulations are relatively low-frequent, whereas if the Doppler shift is zero, the modulations are relatively high-frequent. When computing twos equences of scintillations, but for different transverse speeds, one does not obtain acompressed version of the other using the method as described thus farinthis section, even when using the same seed for the Pseudo-Random Number Generator.T his is because, by adjusting v ⊥ ,weeffectively applied the compression on the filter h C [n]b efore convolution with z[n]. Therefore, the twor esulting sequences will not look compressed in time, butinstead just filtered differently.Evenso, the generated sequence will have the desired statistical properties.
We can obtain the desired compression, i.e., scaling in time, if we instead perform the convolution between z

Saturation of the log-amplitude fluctuations
According to equation (14) the fluctuations increase with distance and frequencyi ndefinitely.F or longer path lengths or stronger turbulence, the amplitude fluctuations gradually levelo ff .S aturation of the amplitude fluctuations can be observed when measuring aircraft noise at distances of overafew kilometers. The standard deviation of the fluctuating sound pressure levels is in such cases limited to approximately 6dB [15]. Saturation of the log-amplitude fluctuations can be included based on an analysis by Wenzel [17]. In Wenzel's approach the soundwavei ss plit up in ac oherent and incoherent wave.T he amplitude of the coherent wave decreases overd istance while the incoherent wave obtains the energy from the coherent wave.T he coherent wave p is written as  Wenzel defines the distance to saturation r s as the distance at which the coherent wave is reduced to e −1 of its original strength Saturation of the log-amplitude fluctuations can nowb e included by multiplying the log-amplitude with ac orrection factor.The variance of log-amplitude fluctuations that includes saturation σ χ ,isgiven by

As ingle tone
To demonstrate the model we shall nowconsider atone of 1000 Hz and apply modulations to the tone. The following parameters are used: c = 343.2m/s, σ 2 µ = 3 · 10 −6 , L = 1.1m, v ⊥ = 2.0m/s and d = 500 m. Saturation according to equation (42) is also included. The filter h c is designed to have 8192 taps. The correlation time is then τ = L/v ⊥ = 0.55 sand the sample frequencyatwhich the modulations are generated 5/τ = 9.1Hz. Figure 3s hows the time series of amplitude and phase fluctuations that were generated. Figure 4s hows the tone along with amodulated version of the tone. The tone was sampled at 8000 Hz.
The log-amplitude and phase move in-phase. The value for the log-amplitude is slightly lower than that of the phase because of saturation. As expected the phase fluctuations cause spectral broadening of the tone as can be seen in Figure [

Scintillations as function of transverse speed
By changing the transverse speed the refractive-indexfield is sampled at ad i ff erent speed. This effectively changes the correlation length and shifts the frequencyr ange that is covered by our applied filter [19]. Figure 6s hows the correlation function for different transverse speeds. With ah igh transverse speed the correlation drops faster,and this will result in relatively more high-frequencyfluctuations, as is shown in Figure 7.

Application in auralization of aircraft
The method wasdeveloped in order to create more realistic sounding auralizations of aircraft. An aircraft movesat high speed through the atmosphere. Atransverse speed is computed for each propagation path separately.B ecause the correlation time, distance and propagation varies between the twop aths, decorrelation occurs. Figures 8a nd 9showspectrograms respectively with and without turbulence. The vertical lines that can be seen are the amplitude modulations. Earlier it wasa ssumed that the correlation length was much smaller than the Fresnel zone size. In this auralization the aircraft is moving close to the receiver. When the aircraft is closest, the distance is almost entirely givenby the height which wasapproximately 100 meters. The correlation length wasset at 20 meters. In that case the Fresnel zone is larger than the correlation length for the lower frequencies, and equation (14) is not valid. Instead, the log-amplitude and phase variances should scale with respectively d 3 and 2k 2 d instead of k 2 d [12].

Conclusion
Fluctuations in the refractive-indexfield due to variations in temperature and wind affects sound propagation and causes audible modulations. Am ethod wasp resented for generating sequences of modulations and applying these to monochromatic as well as broadband signals.   AR ytova pproximation to first-order refractive-index flucutuations results in ac omplexp hase which we can write as al og-amplitude χ and phase S fluctuation. The propagating sound is modelled as at ime-varying channel where we consider twos equences, one for the logamplitude fluctuations, and another for the phase fluctuations.
The fluctuations are frequency-dependent and therefore afi lter wasd esigned to taket hat into account. AG aussian turbulence spectrum wasc onsidered, butt he general method can be used with other turbulence spectra as well. The Vo nKarman spectrum describes real turbulence spectra typically better than the Gaussian spectrum, however, the Vo nKarman spectrum is computationally much more demanding.