Wavenumber spectrum determination for aeroacoustic applications using FISTA

Hans-Georg Raumer 1

German Aerospace Center (DLR) Bunsenstr. 10, 37073 Göttingen

Carsten Spehr German Aerospace Center (DLR) Bunsenstr. 10, 37073 Göttingen

ABSTRACT This conference paper deals with computational methods to determine the wavenumber spectrum of acoustic data measured by a phased microphone array. Such problems occur e.g. within the analysis of pressure fluctuations due to a turbulent boundary layer on a surface such as a wind-tunnel wall or the skin of an aircraft. The problem is closely related to the deconvolution of dirty beamforming maps in wavenumber domain, which seeks to determine the wavenumber spectrum by removing the influence of the shift-invariant point spread function from the beamforming result. Firstly, we recall how this task can be formulated as a minimization problem and then discuss a specific solver for this problem, provided by the framework of the generalized FISTA algorithm. The resulting method includes regularization with L 1 and L 2 penalties as well as a nonnegativity constraint. By exploiting convolutional structures, the computation can be further accelerated. Finally, the presented algorithmic framework is demonstrated with numerical examples.

1. INTRODUCTION

The analysis of turbulent boundary layers is an important task for various applications in experimental aeroacoustic, such as the prediction of cabin noise due to structural vibrations on the fuselage or denoising of aeroacoustic windtunnel measurement data. If a sensor array is used as measurement device, the frequency wavenumber spectrum of the raw pressure fluctuations data is an essential quantity that needs to be computed. A rough estimator of the frequency wavenumber spectrum can be obtained by beamforming in the wavenumber domain. Subsequent postprocessing of the "dirty" beamforming map (also known as deconvolution) further improves the accuracy of the result. This can be done, for instance, by the DAMAS 2.1 scheme (cf. [1]). However, without further modifications, this approach often converges very slowly and it does not guarantee a unique solution. Therefore we discuss a framework that has a well-defined unique solution and provable optimal convergence properties. The optimization scheme relies on the generalized FISTA algorithm (cf. [2]).

1 hans-georg.raumer@dlr.de

a slaty. inter.noise 21-24 AUGUST SCOTTISH EVENT CAMPUS O ¥, ? GLASGOW

2. DISCRETE PROBLEM FORMULATION

We consider a planar region beneath a turbulent boundary layer of a flow field. For the entire analysis we assume that the flow is aligned with the x -direction of the spatial coordinate system. Moreover, assume that the correlation of pressure fluctuations in frequency domain between two points ( x , y ) , ( x ′ , y ′ ) in the xy -plane can be described by a function that depends only on their spatial separation. More precisely, for a given frequency f and a complex pressure signal p we have

E  p ( x , y , f ) p ( x ′ , y ′ , f )  = Φ x − x ′ , y − y ′ , f
.

To simplify the notation, we introduce the separation coordinates ξ (separation in flow direction) and η (separation in cross flow direction). The frequency wavenumber spectrum of the pressure data is then given by the two-dimensional spatial Fourier transform of Φ (cf. [3]) i.e.

∞ Z

∞ Z

b Φ ( k x , k y , f ) = 1 (2 π )

−∞ Φ ( ξ, η, f ) exp  − i ( k x ξ + k y η )  d ξ d η . (1)

−∞

Conversely, Φ can be recovered by the inverse Fourier transform

∞ Z

∞ Z

Φ ( ξ, η, f ) = 1 (2 π )

b Φ ( k x , k y , f ) exp  i ( k x ξ + k y η )  dk x dk y . (2)

−∞

−∞

From now on we will omit the explicit dependency on the frequency f in the notation.

For the application to an experimental setup, one needs to discretize the equations above. We consider a planar array of M microphones

M = { ( x m , y m ) | m = 1 , . . . , M } . (3)

The noise-free correlation data is given by the cross spectral matrix C ∈ C M × M with entries

C ml = E  p ( x m , y m ) p ( x l , y l )  = Φ ξ ml , η ml
,

where ξ ml = x m − x l and η ml = y m − y l . The wavenumber domain of interest is discretized by a grid of wavenumber focus points K given by the cartesian product of two 1D grids K x , K y i.e.

K = { ( k n , x , k n , y ) | n = 1 , . . . N } = { ( k x , k y ) | k x ∈K x , k y ∈K y } . (4)

We wish to reconstruct the frequency wavenumber spectrum at ( k n , x , k n , y ) for n = 1 , . . . N , represented by a vector ˆ φ ∈ C N with entries ˆ φ n = b Φ ( k n , x , k n , y ) . (5)

The relation between the cross spectral matrix and the discrete wavenumber spectrum ˆ φ can be approximated as

∞ Z

∞ Z

C ml = 1 (2 π )

b Φ ( k x , k y ) exp  i ( k x ξ + k y η )  dk x dk y

−∞

−∞

N X

≈ area( K )

(6)

n = 1 ˆ φ n exp  i ( k n , x ξ ml + k n , y η ml ) 

(2 π ) N

n = 1 exp  − i ( k n , x x l + k n , y y l )  ˆ φ n exp  i ( k n , x x m + k n , y y m )  ,

N X

= area( K )

(2 π ) N

where area( K ) =  max n k n , x − min n k n , x  ·  max n k n , y − min n k n , y  . With the wavenumber focus matrix

E mn = exp  i ( x m k n , x + y m k n , y )  (7)

and the scaling factor

s = area( K )

N 2 π (8)

Equation 6 can be compactly reformulated by the discrete forward model

C = s · ED ( ˆ φ ) E H = : T ( ˆ φ ) , (9)

where D ( φ ) = diag  ˆ φ 1 , . . . , ˆ φ N  and E H denotes the Hermitian transpose. The discrete adjoint forward operator T ∗ is given by

T ∗ ( K ) = s · diag  E H KE  ,

for an arbitrary matrix K ∈ C M × M . Hence, setting K = T ( ˆ φ ) we get

T ∗  T ( ˆ φ )  = s 2 · diag  E H ED ( φ ) E H E  . (10)

In practice, the exact (noise-free) data C is not accessible but only a noisy approximation C obs ≈ C . An estimator of the discrete wavenumber spectrum can be obtained by the following regularized non- negative least squares problem

min ˆ φ ≥ 0 1 2

F + α 2

ˆ φ 2

T (ˆ φ ) − C obs 2

2 + α 1 ˆ φ 1 , (11)

2

with regularization parameters α 1 , α 2 . For α 2 > 0 the problem is strictly convex and there exists a unique minimizer.

3. FISTA FOR WAVENUMBER SPECTRUM DETERMINATION

The minimization problem from Equation 11 consists of a smooth part and a convex non-smooth part. Moreover, the later has a computable proximal mapping. Problems of such type can be e ffi ciently solved by the generalized FISTA algorithm [2, p. 291 ff .] (also known as fast proximal gradient method).

3.1. Vectorized FISTA formulation A straightforward implementation of the FISTA algorithm for problem 11 is given below. For guaranteed convergence, the stepsize τ must be chosen such that

− 1 .

∥T ∗ T ( x ) ∥ 2



τ <  sup x ∈ R N , x , 0

∥ x ∥ 2

The upper bound can be su ffi ciently estimated by a few steps of the power method (cf. [4, p. 239]). The expensive step in Algorithm (1) is performed in line 6. Using the specific structure of those matrix products, the computation can be accelerated compared to a straightforward matrix multiplication. For a detailed discussion of this aspect within the context of acoustic source power reconstruction, we refer to [5].

However, for the wavenumber domain problem considered here, the vectorized formulation can be even further accelerated. This will be discussed in the next subsection.

Algorithm 1: FISTA in wavenumber domain (standard formulation)

input : E ∈ C M × N wavenumber focus matrix , C obs ∈ C M × M observed CSM, ˆ φ (0) ∈ R N starting value, α 1 , α 2 > 0 regularization parameters, τ > 0 stepsize, n iter number of maximum iterations output: ˆ φ ( n iter ) ∈ R N approximate solution of Problem (11)

1 t 0 : = 0; ˆ φ ( − 1) : = ˆ φ (0) ; z : = Re diag ( E ∗ CE )

2 for n = 0 , . . . , n iter − 1 do

2  1 + p

1 + 4 t 2 n 

3 t n + 1 : = 1

4 β n : = t n − 1

t n + 1 5 v ( n ) : = ˆ φ ( n ) + β n ( ˆ φ ( n ) − ˆ φ ( n − 1) )

6 u ( n ) : = s 2 · diag  E ∗ ED ( ˆ φ ) E ∗ E 

7 w ( n ) : = v ( n ) − τ  u ( n ) − z 

α 2 τ + 1  + // ( ... ) + takes the positive part

8 ˆ φ ( n + 1) : =  Re  w ( n ) − τα 1

9 end

3.2. Convolutional FISTA formulation As mentioned before, the computational step that dominates the overall computational cost is the evaluation of T ∗ T ( ˆ φ ) (cf. Equation 10). For an arbitrary but fixed component index j ∈{ 1 , . . . , N } we get

j = s 2 N X

( E ∗ E ) jn  2 ˆ φ n

 T ∗ T ( ˆ φ ) 

n = 1

M X

= s 2 N X

n = 1 ˆ φ n

m , l = 1 exp  − i ( x m − x l )( k j , x − k n , x )  exp  − i ( y m − y l )( k j , y − k n , y ) 

= : s 2 N X

n = 1 ˆ φ n P h k j , x − k n , x , k j , y − k n , y i . (12)

Note that P must be evaluated at all possible di ff erences between two wavenumber grid points. Therefore, P is defined for the extended wavenumber grid

K e x = { k x − k ′ x | k x , k ′ x ∈K x } , K e y = { k y − k ′ y | k y , k ′ y ∈K y } , K e = K e x × K e y

i.e. we have P ∈ R | K e x | × | K e y | . Moreover, we consider the frequency wavenumber spectrum in matrix form denoted by X ∈ R |K x |× | K y | such that

X h k x , k y i = b Φ ( k x , k y ) for k x ∈K x , k y ∈K y .

Using this data representation we can reshape the result in Equation 12 to matrix form, where the evaluation index j is replaced by the corresponding 2D index [ k x , k y ]. This yields

s 2 X

X h k ′ x , k ′ y i P h k x − k ′ x , k y − k ′ y i , (13)

k ′ x ∈K x , k ′ y ∈K y

which is essentially a 2D convolution of X and P . The computation of the convolution in 13 can be e ffi ciently carried out by fast convolution schemes, based on the fast Fourier transform such as SciPy’s fftconvolve . For the FISTA algorithm one has to ensure that the output has the same shape as X . With fftconvolve this is achieved by the Python command

fftconvolve(X, P, mode=’same’) .

A full description of the convolutional FISTA scheme is given below.

Algorithm 2: FISTA in wavenumber domain (convolutional formulation)

input : E ∈ C M × N wavenumber focus matrix, P ∈ R | K e x | × | K e y | convolution kernel matrix, C obs ∈ C M × M observed CSM, X (0) ∈ R |K x |× | K y | starting value, α 1 , α 2 > 0 regularization parameters, τ > 0 stepsize, n iter number of maximum iterations output: X ( n iter ) ∈ R |K x |× | K y | approximate solution of Problem (11) in 2D form

1 t 0 : = 0; X ( − 1) : = X (0) ; Z : = reshape  Re diag ( E ∗ CE )
, |K x | ×  K y  

2 for n = 0 , . . . , n iter − 1 do

2  1 + p

1 + 4 t 2 n 

3 t n + 1 : = 1

4 β n : = t n − 1

t n + 1 5 V ( n ) : = X ( n ) + β n ( X ( n ) − X ( n − 1) ) // fast convolution V ( n ) ∗ P , ensure that U ( n ) has shape |K x | ×  K y

6 U ( n ) : = s 2 · fastConv  V ( n ) , P 

7 W ( n ) : = V ( n ) − τ  U ( n ) − Z 

α 2 τ + 1  + // ( ... ) + takes the positive part

8 X ( n + 1) : =  Re  W ( n ) − τα 1

9 end

4. NUMERICAL EXAMPLES

To illustrate the discussed algorithm we consider an exemplary problem. The convective velocity, i.e. the speed of propagation of turbulent structures on the array surface, is denoted by u c and the convective wavenumber is given by k c = 2 π f

u c . Moreover, l x , l y denote the correlation lengths of turbulent structures in x and y direction and a denotes an amplitude factor. With these parameters we model the exact correlation data by

l x − | η |

Φ ( ξ, η ) = a · exp ( i ξ k c ) · exp − | ξ |

! . (14)

l y

Given this expression we can explicitly evaluate the spatial Fourier transform, which yields the exact wavenumber frequency spectrum

∞ Z

∞ Z

−∞ Φ ( ξ, η, f ) exp  − i ( k x ξ + k y η )  d ξ d η 2 al x l y

b Φ ( k x , k y , f ) = 1 (2 π )

π  1 + l 2 x ( k x + k c ) 2  1 + l 2 y k 2 y  . (15)

−∞

The physical and algorithmic parameters are summarized in table 1. The microphone array consists of an array with M = 144 microphones aranged in spiral arms. This array was used in the benchmark measurement ’DLR1’ (cf. [6, 7]). Figure 1 shows the microphone positions and their spatial separations (the co-array). The true solution (cf. Equation 15) is shown in Figure 2. The dominant structure is the so-called convective ridge.

For the FISTA computations we consider noisy data

C obs = C + ϵ · a · rr H , (16)

with a multivariate complex standard normal random variable r ∼ [ N C (0 , 1)] M and a noise power factor ϵ . Figure 3 shows the FISTA results for several choices of ϵ .

Parameter Value

convective velocity u c 82 m

s

speed of sound c s 343 m

s

frequency f f = 1347 Hz

coherence length l x l x = u c 0 . 1 · 2 π f m

coherence length l y l y = u c 2 π f m

amplitude a 1000 Pa 2

noise power factor ϵ 0 . 0 , 0 . 01 , 0 . 02 , 0 . 05

acoustic wavenumber k 0 = 2 π f

c s 1000 m − 1

wavenumber grid K x uniform grid on [ − 8 , 8] with 64 grid points

wavenumber grid K y uniform grid on [ − 8 , 8] with 64 grid points

number of FISTA iterations n iter 100

− 1

stepsize τ of gradient step 0 . 99 ·  sup x ∈ R N , x , 0



∥T ∗ T ( x ) ∥ 2

∥ x ∥ 2

L 1 regularization parameter α 1 10 − 3 · C obs 2

F L 2 regularization paramet er α 2 10 − 4 · C obs 2

F

Table 1: Parameter settings for the numerical example.

0 . 6

1

0 . 4

0 . 5

0 . 2

η [m]

y [m]

0

0

− 0 . 2

− 0 . 5

− 0 . 4

− 1 − 0 . 5 0 0 . 5 1 − 1

− 0 . 6 − 0 . 4 − 0 . 2 0 0 . 2 0 . 4 0 . 6 − 0 . 6

ξ [m]

x [m]

Figure 1: Array sensor positions (left) and spatial sensor separations (right).

0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6

5

k 0 [m]

0

k y

− 5

− 5 0 5

k x k 0 [m]

Figure 2: True solution

The computations were carried out on a standard notebook. Figure 4 compares the computational e ff ort of the vectorized FISTA formulation (Algorithm 1) and the convolutional formulation (Algorithm 2). We observe that in the chosen parameter range, the convolutional implementation is about a factor 4 faster than the vectorized method.

5. CONCLUSIONS

We presented a computational framework that computes an estimator of the frequency wavenumber spectrum of pressure fluctuation measurements on a microphone array. The algorithm employed the setup of the well-known FISTA optimization scheme. For this particular problem, the computationally most expensive step has a convolutional structure. Therefore, it is strongly beneficial in terms of e ffi ciency, to use fast convolutions for the implementation.

FUNDING

This research has received funding from the "Flight-LAB / OVAL" research project within the aerospace research program (LuFo V2; Support code 20K1511C) supported by the Federal Ministry for Economic A ff airs and Energy.

0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5

0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5

5

5

k 0 [m]

k 0 [m]

0

0

k y

k y

− 5

− 5

− 5 0 5

− 5 0 5

k x k 0 [m]

k x k 0 [m]

(a) ϵ = 0 . 0, L2-error = 0 . 758

(b) ϵ = 0 . 01, L2-error = 1 . 10

0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5

0 0 . 2 0 . 4 0 . 6

5

5

k 0 [m]

k 0 [m]

0

0

k y

k y

− 5

− 5

− 5 0 5

− 5 0 5

k x k 0 [m]

k x k 0 [m]

(c) ϵ = 0 . 02, L2-error = 1 . 52

(d) ϵ = 0 . 05, L2-error = 2 . 40

Figure 3: FISTA results for the parameter setup from Table 1.

standard convolutional

8 · 10 − 2

Time per iteration [s]

6 · 10 − 2

4 · 10 − 2

2 · 10 − 2

20 40 60 80 100 120 0

|K| x number of grid points per dimension

Figure 4: Computational cost for Algorithm 1 and 2. For the time measurements, a quadratic wavenumber grid was used (i.e. |K x | = |K x | ). The number of grid points for each direction x , y was increased gradually from 8 (i.e. a 8 × 8 grid) to 120 (i.e. a 120 × 120 grid).

REFERENCES

[1] S. Haxter. Extended version: Improving the damas 2 results for wavenumber-space beamforming. In Paper BeBeC-2016-D8; Proceedings of the 7th Berlin Beamforming Conference , February 2016. [2] A. Beck. First-Order Methods in Optimization . Society for Industrial and Applied Mathematics, oct 2017. [3] K. Ehrenfried and L. Koop. Pressure fluctuations beneath a compressible turbulent boundary layer. In 14th AIAA / CEAS Aeroacoustics Conference (29th AIAA Aeroacoustics Conference) . American Institute of Aeronautics and Astronautics, may 2008. [4] H. W. Engl, M. Hanke, and A. Neubauer. Regularization of Inverse Problems . Springer Netherlands, 1996. [5] G. Chardon, J. Picheral, and F. Ollivier. Theoretical analysis of the DAMAS algorithm and e ffi cient implementation of the covariance matrix fitting method for large-scale problems. Journal of Sound and Vibration , 508:116208, sep 2021. [6] T. Ahlefeldt. Aeroacoustic measurements of a scaled half-model at high reynolds numbers. AIAA Journal , 51(12):2783–2791, December 2013. [7] C. J. Bahr, W.M. Humphreys, D. Ernst, T. Ahlefeldt, C. Spehr, A. Pereira, Q. Leclère, C. Picard, R. Porteous, D. Moreau, J. R. Fischer, and C. J. Doolan. A comparison of microphone phased array methods applied to the study of airframe noise in wind tunnel testing. In 23 rd AIAA / CEAS Aeroacoustics Conference , Reston, VA, June 2017. American Institute of Aeronautics and Astronautics.