Wei-Hsiang Liao
email: wliao@ircam.fr

Axel Roebel
email: roebel@ircam.fr

Alvin W Y Su
email: alvinsu@mail.ncku.edu.tw

ON THE MODELING OF SOUND TEXTURES BASED ON THE STFT REPRESENTATION

Sound textures are often noisy and chaotic. The processing of these sounds must be based on the statistics of its corresponding time-frequency representation. In order to transform sound textures with existing mechanisms, a statistical model based on the STFT representation is favored. In this article, the relation between statistics of a sound texture and its time-frequency representation is explored. We proposed an algorithm to extract and modify the statistical properties of a sound texture based on its STFT representation. It allows us to extract the statistical model of a sound texture and resynthesise the sound texture after modifications have been made. It could also be used to generate new samples of the sound texture from a given sample. The results of the experiment show that the algorithm is capable of generating high quality sounds from an extracted model. This result could serve as a basis for transformations like morphing or high-level control of sound textures.

INTRODUCTION

Environmental sounds such as wind, rain or crowd noises are common in our daily life. In contrast to instrumental sounds, these sounds are often inharmonic, chaotic and possess stronger noise components. They are often known as sound textures.

Recently, sound texture processing has become more important in the entertainment industry and various other areas. High quality sound texture transformations are essential for work in these fields. State-of-the-art algorithms are capable of achieving highquality transformations with sinusoidal sounds [START_REF] Dolson | The phase vocoder: A tutorial[END_REF] [2] [3] [START_REF] Bloit | Towards morphological sound description using segmental models[END_REF], while preserving the envelope and naturality of transients [START_REF] Roebel | A new approach to transient processing in the phase vocoder[END_REF] [6] [START_REF] Roebel | Shape-invariant speech transformation with the phase vocoder[END_REF]. However, these algorithms yield less satisfying results when applied to sound textures. The underlying reason is that sound textures do not fit the common assumptions of most state-of-the-art algorithms, such as harmonicity and sinusoidality.

Sound textures are often driven by a stochastic process. They usually exhibit some type of structure over a moderate time period. Julesz's conjecture [START_REF] Julesz | Visual pattern discrimination[END_REF] and [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] suggests that the structure can usually be described by a series of statistical properties. Moreover, in [START_REF] Mcdermott | Sound texture perception via statistics of the auditory periphery: Evidence from sound synthesis[END_REF] and [START_REF] Mcdermott | Summary statistics in auditory perception[END_REF], it is stated that humans distinguish between sound textures based on their statistical properties. Based on these assumptions, it is possible to establish a parametrized model which describes a sound texture using a set of statistical properties. This model will allow the analysis, synthesis and transformation of sound textures.

There are several previous works on sound textures. Saint-Arnaud [START_REF] Saint | Computational auditory scene analysis[END_REF] proposed an analysis/synthesis scheme based on modeling the atomic events in sound textures. In [START_REF] Hanna | Time scale modification of noises using a spectral and statistical model[END_REF], stochastic noises are synthesized with randomized sinusoids. Schwarz [START_REF] Schwarz | Data-Driven Concatenative Sound Synthesis[END_REF] proposed a descriptor-driven, corpus-based approach to synthesize sound textures. In addition, several parametric models have been proposed for image textures. For example, Portilla [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] proposed a texture model based on the statistics of wavelet coefficients. Bruna [START_REF] Bruna | Classification with invariant scattering representations[END_REF] proposed a new wavelet transform which provides a better view of textures while capturing high-order statistics. For sound textures, McDermott [START_REF] Mcdermott | Sound texture perception via statistics of the auditory periphery: Evidence from sound synthesis[END_REF] proposed a parametrized model, adapted from [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF], which characterizes target sound textures with high order statistics and correlation between subband envelopes. In his article, statistical properties are applied to Gaussian noises to generate different sound textures.

A common approach of signal processing is to use the STFT (Short Time Fourier Transform) representation. The idea is to use STFT to obtain the time-frequency representation of the signal and apply the transformation. However, conventional transformation techniques, such as phase vocoder, cannot be applied to sound textures [START_REF] Liao | On stretching gaussian noises with the phase vocoder[END_REF]. To apply transformations to a sound texture over its STFT representation, a statistical model based on the signal's timefrequency representation is required. Also, establishing a model based on the STFT representation can allow us to utilize existing efficient implementations of DFT.

In this article, we would like to propose a model which can be used to analyse and synthesise sound textures. It's based on the statistics of time-frequency representations of sound textures. From the model, one could generate new samples of the input texture. It could serve as a basis for high-level transformations of sound textures. An experiment has been conducted to investigate the quality of sound generated from the model.

The paper is organized as follows: The first section describes the present state of sound texture processing and the importance of establishing a statistical model. The second section discusses perceptually important statistics of sound textures and the set of properties to be used in the model. The third section describes the process of modeling a sample texture and resynthesising it with the resulting model. The fourth section is an experiment of texture generation with the model. The fifth section are the conclusion and future works.

SOUND TEXTURE STATISTICS

Mcdermott's work [START_REF] Mcdermott | Sound texture perception via statistics of the auditory periphery: Evidence from sound synthesis[END_REF] suggests that a proper model of a sound texture is composed of the statistical properties of individual subband signals. These subband signals can be formed from the timefrequency representation of the original signal. Furthermore, from his experiment, we know that the perceptually important statistical properties can be summarized into three categories : Subband Moments, Temporal Correlation and Spectral Correlation.

Subband Moments

Subband moments are the statistical moments of the spectral coefficients of subbands. They describe the shape of the histogram of a subband signal.

Temporal Correlation

Temporal correlation refers to the autocorrelation of a subband signal. It can be characterized by the autocorrelation function.

Spectral Correlation

Spectral correlation means the cross-correlation between different subband signals.

While the autocorrelation function and cross-correlation characterize the horizontal and vertical relationships in the time-frequency representation, slant relationships should also be considered. This can be characterized by the delayed cross-correlation, which is the cross-correlation function. As such, the three categories could be summarized in two, where moments characterize the property of one subband signal, and correlation functions characterize the relationships between subband signals in all directions.

Here, the proposed model uses STFT as the underlying timefrequency representation. Each bin of the STFT representation is treated as a single subband, and the coefficients of each subband form the subband signal. Since humans mainly perceive envelope statistics [START_REF] Mcdermott | Sound texture perception via statistics of the auditory periphery: Evidence from sound synthesis[END_REF], the statistics of the subband signals and their envelopes are both included. In order to establish a model which contains all the statistics from the categories described above, the following statistics are included: the first four subband moments, the auto-correlation function of each subband, and the cross-correlation function of each pair of subbands.

MODELING SOUND TEXTURES

In this section we describe the proposed modeling mechanism. The first subsection describes the process by which we extracts the model from a given sound texture. The second subsection describes how to apply the model to a STFT representation and the steps to follow to resynthesize sounds.

Model Extraction

Assume a given input signal x, along with a window function w, hop size h and the size of Fourier transform N . Applying STFT to the signal will obtain the analysis frames. The analysis frame X l which is centered at time point l is:

X l (k) = ∑ t x(t)w(t -lh)e -j2πtk N (1)
These analysis frames X l form the STFT representation. Thus, the subband signals are formed as:

S k (n) = Xn(k)e j2πlhk N (2)
In ( 2), a demodulation term is applied to remove the carrier frequency in each subband. The carrier frequency is introduced by the STFT and will affect the evaluation of correlation functions. Next, we extract the envelope E of the subband signal, which can be done by any envelope extraction algorithm. One could also use the amplitude of the subband signal for simplicity.

E k (n) = envelope(S k (n)) (3)
At last, statistics are evaluated for every S k and E k . Here, only the moments of the envelopes E k are evaluated. This is because currently we have not yet found a proper way to define complex moments for the proposed model. Finally, the extracted model Φ of the input texture can be written as:

Φ ≡ {AS, CSS, ME, AE, CEE} (4)
AS k (l) = { ∑ t S k (t)S k (t + l)} (5)
CS x Sy (l) = { ∑ t Sx(t)Sy(t + l)} (6)
ME k = {µ1(E k ), µ2(E k ), η(E k ), κ(E k )} (7) AE k (l) = { ∑ t E k (t)E k (t + l)} (8)
CE x Ey (l) = { ∑ t Ex(t)Ey(t + l)} (9)
In (4), A denotes the autocorrelation function, C denotes the cross-correlation function and M denotes the first four moments, namely mean, variance, skewness and kurtosis. Where mean and variance are central moments, skewness and kurtosis are normalized central moments. After the model has been extracted, certain transformations could be applied to it. For example, timestretching is equal to stretching the autocorrelation and cross-correlation functions.

Texture Resynthesis

Sound textures can be generated from the model extracted in the previous subsection. The idea is to apply the statistical properties to a randomized STFT representation. Usually, it requires immense iterative computation to estimate the output signal. However, this can be avoided by utilizing the relationship between the correlation and the spectrum.

Applying autocorrelation

To apply the autocorrelation functions, we could use the Wiener-Khinchin theorem. It suggests that the amplitude spectrum of a signal is also the spectrum of its auto-correlation function.

F(As) = F (s(t) * s(-t)) = |F(s)| 2 (10)
That is, we can assign any autocorrelation function to a subband signal by constraining its amplitude spectrum. To avoid circular aliasing, a Fourier transform with twice the length of the autocorrelation function is used when calculating the spectrum of the subband signal. This introduces a windowing effect to the signal. We use Griffin Lim's [START_REF] Griffin | Signal estimation from modified short-time fourier transform[END_REF] technique to re-establish the window shape after the inversion.

Applying cross-correlation

The cross-correlation functions can be dealt with in a similar fashion. Here we introduce an extension of [START_REF] Mcdermott | Sound texture perception via statistics of the auditory periphery: Evidence from sound synthesis[END_REF], the cross-correlation theorem:

F(Cx,y) = F (x(t) * y(-t)) = F (x(t))F (y(t)) = |X(k)||Y (k)|e j(θ(X(k))-θ(Y (k))) (11)
Equation ( 11) can be rewritten as:

θ( F (Cx,y) |X(k)||Y (k)| ) = θ(X(k)) -θ(Y (k)) (12)
Since we already constrained the amplitude spectrum by ( 10), we could control the cross-correlation function between x and y by adjusting their phase differences without affecting the autocorrelation functions. To establish full cross-correlation across all subband signals, choose one subband signal as the starting point, initialize phases for the spectrum of the chosen subband signal, assign phases to neighboring subbbands with the phase differences extracted from the original input texture, as in [START_REF] Saint | Computational auditory scene analysis[END_REF]. Here we choose the zeroth bin ,S0, as the starting point. Since S0 is real, the phase of it's spectrum must be Hermitian, that is, θ(FS 0 (k)) = -θ(FS 0 (-k)).

Applying Moments

The moments only apply to envelopes. Since mean and variance are already contained in the autocorrelation function, the remaining moments are the skewness and kurtosis. A gradient descent algorithm is used to apply the skewness and kurtosis, the detail of which can be found in [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF]. The following formula is used to apply the moments to a signal x:

x ′ = x + γ ∂f (x) ∂x (13)
The target signal x ′ can be obtained by solving γ with respect to f (x ′ ) = t. If there is no real root for γ, a real-value approximation which is closest to the target value will be chosen. In [START_REF] Hanna | Time scale modification of noises using a spectral and statistical model[END_REF], f is either skewness or kurtosis, where ∂f (x) ∂x is the derivative of the moment function respected to x, and t is the target value of the desired moment. For skewness η and kurtosis κ, their derivatives with respect to x are:

∂η(x) ∂x ≡ x • x -µ 1/2 2 (x)η(x) -µ2(x) (14) ∂κ(x) ∂x ≡ x • x • x -µ2(x)κ(x)x -µ3(x) (15)
With the gradient descent algorithm, we could apply the moments to the subband signals. Unfortunately, the gradient descent algorithm would change the correlation functions of the subband signals, so, an iterative approach is used to generate the target STFT representation.

The whole resynthesis step is depicted in fig. 1. It begins with a randomly initialized STFT representation. Next, the STFT representation is decomposed into the envelope part and the remaining part. Envelope statistics {AE, CEE, ME} are applied by the methods described in previous sections. Then, reconstruct the STFT representation by combining the modified envelope and the remaining part, and apply the subband signal statistics {AS, CSS} accordingly. If the iteration has converged or reached maximum iteration number, output the STFT representation, otherwise, repeat the iteration process. At last, apply the remodulation, inverse STFT and overlap-add to generate the resynthesized signal x:

x ′ (n) = ∑ k S k (l)e j2πk(n-lh) N x(n) = ∑ n w(n -lh)x ′ (n) ∑ w 2 (n -lh) (16)
Figure 1: The flow diagram of the resynthesis procedure

EXPERIMENT

To demonstrate the texture sound generated from the model, we used some texture samples from [START_REF] Mcdermott | Sound texture perception via statistics of the auditory periphery: Evidence from sound synthesis[END_REF]. Statistical models were extracted from the samples of sound textures, then new samples were regenerated from the model. To generate a texture longer than the original, the model is used twice and connected with overlap-andadd. The discontinuation between generated blocks will be handled in the future. The detailed configuration is listed below:

• number of samples : 6

• original signal : 16bit mono wav, 5-7 secs

• sampling rate : 20khz

• hop size of input frame : 16 samples

• size of Fourier transform : 256, hanning window

• generated sample length : 10-14 secs (2 times long) The generated sound samples can be found in [18]. From the quality of generated sounds, we could see the proposed model is capable of capturing most of the perceptually important features. However, the iterative process is not guaranteed to converge to a good solution, as the case of fig. 2f. It may require more constraints to improve the quality of the click sounds.

CONCLUSION

In this paper, we proposed a parametric statistical model based on the STFT representation. It can be used to analyse and synthesise sound textures and generate high quality sounds. Since it's based on the STFT, it could benefit from many existing related implementations. Possible applications of this model are timestretching, texture morphing and providing high-level control of sound textures. Integrating the model with the phase vocoder could also be a future goal. The future works seek a way to apply all the statistics at once (eliminating the iterative process), connect the model with semantic controls and investigate an efficient way to apply only specific parts of correlation functions.

(a)

  Bubbling Water (b) Insects In A Swamp (c) Heavy Rain Falling And Dripping (d) Bees

  (e) Wind Whistling (f) Rustling Papers

Figure 2 :

 2 Figure 2: The spectrograms of original(left) and generated(right) textures.

ACKNOWLEDGMENTS

This research is funded by the French National Research Agency(ANR) project PHYSIS. Thanks to Alex Baker for proof-reading the manuscript.


Wei-Hsiang Liao
email: wliao@ircam.fr

Axel Roebel

Alvin W Y Su
email: alvinsu@mail.ncku.edu.tw

ON STRETCHING GAUSSIAN NOISES WITH THE PHASE VOCODER

Recently, the processing of non-sinusoidal signals, or sound textures, has become an important topic in various areas. In general, the transformation is done by the phase vocoder techniques. Since the phase vocoder technique is based on a sinusoidal model, it's performance is not satisfying when applied to transform sound textures. The following article investigates into the problem using as example the most basic non-sinusoidal sounds, that are noise signals. We demonstrate the problems that arise when time stretching noise with the phase vocoder, provide a description of some relevant statistical properties of the time frequency representation of noise and introduce an algorithm that allows to preserve these statistical properties when time stretching noise with the phase vocoder. The resulting algorithm significantly improves the perceptual quality of the time stretched noise signals and therefore it is seen as a promising first step towards an algorithm for transformation of sound textures.

INTRODUCTION

In present time, there exist many algorithms that have been developed to achieve high quality transformation of sounds. These state-of-the-art algorithms mainly focused on sinusoidal sounds [START_REF] Dolson | The phase vocoder: A tutorial[END_REF] [2] [3] [START_REF] Bloit | Towards morphological sound description using segmental models[END_REF]. Such as instruments and tonal sounds. While the properties of sinusoidal sounds are well-researched, a large class of sounds cannot be dealed with the way dealing with sinusoidal sounds. These sounds often driven by a stochastic process. They often called sound textures.

Sound textures are common in the environment. They could be either natural of artificial. For example, the sounds of wind blowing and rain dropping are natural textures, and the noise of crowd and the sound of train passing are artificial. Although there's no clear definition of sound texture [START_REF] Saint-Arnaud | Computational auditory scene analysis[END_REF], this class of sounds often exploits certain stable structure over a time period. The structure can be described by a series of statistic properties. There are some different approaches dealing with sound textures. For example, Saint-Arnaud [START_REF] Saint-Arnaud | Computational auditory scene analysis[END_REF] proposed an analysis/synthesis scheme which based on atomic features in sound textures, Hanna [START_REF] Hanna | Time scale modification of noises using a spectral and statistical model[END_REF] used randomized sinusoids to synthesize stochastic noises. Schwarz [START_REF] Schwarz | Data-Driven Concatenative Sound Synthesis[END_REF] proposed a descriptor-driven, corpus-based approach to synthesize sound textures. But to enable the further analysis, synthesis and manipulation of sound textures, a parametrized modeling which describes a texture with a set of statistical properties is desired. Fortunately, some previous works on parametrization has been done for image textures, Portilla [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] proposed a texture model based on the statistics of wavelet coefficients. Bruna [START_REF] Bruna | Classification with invariant scattering representations[END_REF] proposed a new wavelet transform which provides a better view for textures while capturing high-order statistics. For sound textures, McDermott [START_REF] Mcdermott | Sound texture perception via statistics of the auditory periphery: Evidence from sound synthesis[END_REF] pro-posed a parametrized model, adapted from [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] and which characterizes target sound texture with high order statistics and correlation between subband envelopes. An algorithm is proposed in his article to convert Gaussian noise into different target textures.

To manipulate sounds, such as time stretching or transposition, using phase vocoder [START_REF] Dolson | The phase vocoder: A tutorial[END_REF] [2] [3] is a common approach. Phase vocoder could stretch a signal while preserving it's envelope and naturality of transient [START_REF] Roebel | A new approach to transient processing in the phase vocoder[END_REF] [12] [START_REF] Roebel | Shape-invariant speech transformation with the phase vocoder[END_REF]. With phase vocoder, one could perform a large factor stretching which still yields a natural result. While the theory basis of phase vocoder works well for sinusoidal signal, it's not suitable for sound textures. Because of the constraint of phase continuity could no longer be used to estimate output phases. To stretch a sound texture with phase vocoder, one must maintain perceptive important properties of the texture during the stretching.

In this article, we try to resolve the properties being changed during the transformation. In order to aim the phenomenon appeared during the process, the time stretching of white Gaussian noise is used. Since Gaussian noise exploits independency along both time and frequency axes, it introduces no extra dependency to the analysis frames. With this simplified case, we found that it's crucial to compensate the change of temporal correlation during the time stretching transformation. A primitive algorithm is given to correct the distortion of temporal correlation.

The paper is organized as follows: section 1 describes previous works on the sound texture parametrization and the issues encountered when applying phase vocoder to sound texture. Section 2 lists the perceptually important statistical properties of sound textures. Section 3 formulates the time stretching procedure with phase vocoder and correct the result with iterative approaches. Section 4 conducts an experiment of the correction. Section 5 gives the conclusion and future prospective.

PROPERTIES OF SOUND TEXTURES

According to Julesz's conjecture [START_REF] Julesz | Visual pattern discrimination[END_REF] and portilla's work [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF], the perceptually important properties of a texture can be described in statistical features. In Mcdermott's work [START_REF] Mcdermott | Sound texture perception via statistics of the auditory periphery: Evidence from sound synthesis[END_REF], it suggests that the perceptually important properties of sound texture could be summarized in three categories : marginal statistics, frequency correlation and temporal correlation. These properties are further described in the following subsections:

Marginal Statistics

Marginal statistics are the statistical moments of the spectral coefficients in different time points or different subbands. Though the usage of higher orders are possible, the first four moments are used to form the marginal statistics in common cases [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF].

Frequency Correlation

In general, frequency correlation means the cross-correlation between different frequency components. These components can be simply coefficients, subbands, or features of subbands, for example, the envelope of cochlear bands [START_REF] Glasberg | Derivation of auditory filter shapes from notched-noise data[END_REF] [START_REF] Mcdermott | Sound texture perception via statistics of the auditory periphery: Evidence from sound synthesis[END_REF]. Moreover, in [START_REF] Mcdermott | Sound texture perception via statistics of the auditory periphery: Evidence from sound synthesis[END_REF], each subband is further filtered by a modulation filterbank, and the cross-correlation of the modulation bands are evaluated.

Temporal Correlation

Temporal correlation refers to the autocorrelation of a frequency component along time axis. It can be the autocorrelation of STFT coefficients [START_REF] Mcdermott | Sound texture synthesis via filter statistics[END_REF] or the autocorrelation of a temporal feature. It is also the most important property when performing time stretching by phase vocoder. The detail is described in section 3.

TIME STRETCHING GAUSSIAN NOISES

To investigate the problem encountered during the stretching of Gaussian by the phase vocoder, the formulation of phase vocoder must be described. The transformation procedure of phase vocoder consists of analysis, modification and resynthesis [START_REF] Dolson | The phase vocoder: A tutorial[END_REF] [2] [START_REF] Laroche | New phase-vocoder techniques for real-time pitch shifting, chorusing, harmonizing and other exotic audio modifications[END_REF]. The analysis first applies STFT to the input signal, then, proceed modification in the spectral domain, and last, the resynthesis is done by inverse STFT and overlap-and-add. Assuming a given input signal s, along with a window function w and the size of Fourier transform N . Then the analysis frame S l which centered at time point l can be written as:

S l (k) = n s(n)w(n -l)e -j2πnk N (1)
During the transformation, the coefficients of S l may be modified as well as the centering position l might be moved [START_REF] Laroche | Improved phase vocoder timescale modification of audio[END_REF] [3] [START_REF] Roebel | A new approach to transient processing in the phase vocoder[END_REF]. If f (l) is the function of the modified center of S l , then the resynthesis procedure can be written as:

s ′ (n) = k S l (k)e j2πnk N (2)
s(n) = n w(n -f (l))s ′ (n) w 2 (n -f (l)) (3)
If inconsistency happens between overlapping frames, (3) outputs the signal with the least square error respect to the corresponding analysis frames [START_REF] Griffin | Signal estimation from modified short-time fourier transform[END_REF]. If we stretch a Gaussian noise (1)-( 3) and assign the phase of each frame as the way processing sinusoidal signal, one would find that the noise component disappeared in the resynthesized signal. In fact, short distorted sinusoids will be perceived in the resynthesized signal. In this case, the phase and the time position of analysis frames are modified during the time stretching transformation. While the amplitude of the coefficients is unmodified, the marginal statistics of analysis frames remain the same. Therefore, the possible variation of correlation between STFT coefficients is investigated.

The correlation between two STFT coefficients in the analysis frames is the inner product between one coefficient and the complex conjugate of the other. By using (1), the correlation is:

C Sα(k 1 ),S β (k 2 ) = n s(n)w(n -α)e -j2πnk 1 N * n s(n)w(n -β)e -j2πnk 2 N (4)
The equation in ( 4) can be used to calculate the correlation in both time and frequency directions. Let the terms in (4) be column vectors, replace the window function w(n-l) by w l , and e

-j2πnk

N by e k N , (4) can be rearranged into the summation of element-wise product(Hadamard product) of outer product matrices:

C Sα(k 1 ),S β (k 2 ) = (s • wα • e k 1 N )(s • w β • e k 2 N ) T (5) = ss T • wαw T β • e k 1 N e k 2 N T (6)
Equation ( 6) consists of three parts, the outer product of s, shifted outer product of w, and the outer product between e k 1 N and e k 2 N . In fact, the latter two terms form the shape of a Gabor filter [START_REF] Feichtinger | Gabor Analysis and Algorithms[END_REF]. The correlation respect to a direction is equal to apply the corresponding Gabor filter to the outer product matrix of s.

Since the stretching is done by moving the center of analysis frames to new positions, their relative distances in time are changed. Therefore the temporal correlation between the frames would no longer correct. It can be seen from ( 6), moving the position of frames affects only the second term, wαw T β . Assuming s is time-invariant and independent to w, LS is the length of s. The estimated temporal correlation of STFT coefficients can be written as the form of complex autocorrelation function in time:

AS(k, l) = LS • E[ss T ] • w0w T l • e k N e k N T (7)
The autocorrelation of analysis frames in [START_REF] Schwarz | Data-Driven Concatenative Sound Synthesis[END_REF] should be AS(k, dl) after time stretching. As one could see, even if s has no innate autocorrelation (E[ss T ] = I), there still exists autocorrelation between the coefficients of STFT frames at the same bin index, which is induced by wαw T β . Originally, the autocorrelation introduced by the analysis window would be canceled after resynthesis. But since the position of analysis frames are moved, the cancellation fails. Fig. 1 shows the difference of autocorrelation function of a Gaussian noise before and after a time stretching of factor 2 with an ordinary phase vocoder. Also, in Fig. 1, the autocorrelation is not exactly stretched by twice in time. This is because the phase compensation in phase vocoder cannot fix the autocorrelation of a noise correctly.

To compensate the change of correlation, the coefficients of S l (k) must be altered to make their autocorrelation approaches AS(k, dl). The optimization which proposed in [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] is applied. Let T (l ′ ) = AS(k, dl) and the maximum lag is L, we could have:

Φ = L l T (l) -AS(k, l) 2 (8) ∂Φ ∂si(k) ∝ L l (T (l) -AS(k, l))(s i-l (k) + s i+l (k)) (9)
Apply the gradient projection mentioned in [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF], we could have:

s ′ i (k) = si(k) + L l λ(l, k)(s)(s i-l (k) + s i+l (k)) (10)
s ′ i (k) = si(k) ⊗ h(i, k) (11)
h(i, k) in ( 11) is a zero-phase filter function with length L. There are N/2 different filters for bin indices between 1 N/2, while the other half are the complex conjugates. At last, we use the algorithm in [START_REF] Dennis | Trustregion interior-point algorithms for minimization problems with simple bounds[END_REF] to solve the coefficients of the filters.

EXPERIMENT

To verify the result of section 3, an experiment is conducted. In this experiment, we perform the time stretching on a random generated Gaussian noise. Since the marginal statistics and spectral correlation are unaffected by time stretching, only the temporal correlation has to be corrected after the transformation. The first step is the ordinary phase vocoder analysis as [START_REF] Dolson | The phase vocoder: A tutorial[END_REF]. Then the following step corrects temporal correlation for each frequency bins respecting to the stretch factor. The resynthesis is done as [START_REF] Laroche | New phase-vocoder techniques for real-time pitch shifting, chorusing, harmonizing and other exotic audio modifications[END_REF]. The estimation of target autocorrelation follows [START_REF] Schwarz | Data-Driven Concatenative Sound Synthesis[END_REF], and E[ss T ] is replaced by σ 2 I due to the nature of Gaussian noise. The detailed configuration is listed below:

• original signal : Gaussian noise (1 sec).

• sampling rate : 44.1khz

• stretching factor : 3.0

• hop size of input frame : 16

• size of Fourier transform : 256

• maximum lag of correlation : 64

The error of temporal correlation for each frequency bin before and after correction is shown in Fig. 2. In most of the frequency bins, the error is lower than 48dB, through a better correction may still exist. The spectrogram of the Gaussian noise is shown in Fig. 3. In the spectrogram of 'stretched w/o correction', the energy of the spectrogram is not uniformly distributed. Both frequency components and the gaps between them tend to persist longer in time, thus making the sound more sinusoidal and resulting perceivable artifacts such as short pitches. After the correction is applied, the intensity of short sinusoids are mitigated, thus reducing the arfefacts. The sound files could be found on [21]. The difference between the sounds is perceivable, but still not perfect. This may due to the temporal correlation was distorted after the resynthesis.

A possible reason is that the phase coherence is not considered in the correction of temporal correlation. If so, certain constraint has to be applied to the correction mechanism. From this simplified case, we found that the variation will always in place regardless the innate correlation of the signal. Also, the variation in the temporal correlation is perceptually significant. Therefore, to achieve a successful time stretching, it is necessary to correct temporal correlation respect to the stretching factor and analysis window. The temporal correlation function can be obtained by estimation under different stretch ratio. Currently, the correction procedure is done iteratively, but it's likely that there exists a non-iterative solution. However, this correction does not include the correction of correlation induced by the signal itself and the correction of correlation is not perfect. The discovery from the streching of Gaussian could be served as the basis of general sound texture transformation. This could lead to various applications like non-uniform stretching, texture morphing and parametric texture synthesizing. The future work includes considering the general sound texture and seek an efficient mechanism to correct statistical properties during the transformation by phase vocoder.

Figure 1 :

 1 Figure 1: The change of temporal correlation before(thick solid line) and after(thin dotted line) phase vocoder stretching.(k = 64, -11.6 ≤ timelag ≤ 11.6)

Figure 2 :Figure 3 :

 23 Figure 2: The error of temporal correlation before(solid line) and after(dotted line) the correction.


Gerardo Del Guercio

Carlo Cossu

Gregory Pujals

Optimal perturbations of non-parallel wakes and their stabilizing effect on the global instability

We compute the spatial optimal energy amplification of steady inflow perturbations in a non-parallel wake and analyse their stabilizing action on the global mode instability. The optimal inflow perturbations, which are assumed spanwise periodic and varicose, consist in streamwise vortices that induce the downstream spatial transient growth of streamwise streaks. The maximum energy amplification of the streaks increases with the spanwise wavelength of the perturbations, in accordance with previous results obtained for the temporal energy growth supported by parallel wakes. A family of increasingly streaky wakes is obtained by forcing optimal inflow perturbations of increasing amplitude and then solving the nonlinear Navier-Stokes equations. We show that the linear global instability of the wake can be completely suppressed by forcing optimal perturbations of sufficiently large amplitude. The attenuation and suppression of self-sustained oscillations in the wake by optimal 3D perturbations is confirmed by fully nonlinear numerical simulations. We also show that the amplitude of optimal spanwise periodic (3D) perturbations of the basic flow required to stabilize the global instability is much smaller than the one required by spanwise uniform (2D) perturbations despite the fact that the first order sensitivity of the global eigenvalue to basic flow modifications is zero for 3D spanwise periodic modifications and non-zero for 2D modifications. We therefore conclude that first-order sensitivity analyses can be misleading if used far from the instability threshold, where higher order terms are the most relevant.

I. INTRODUCTION

Two-dimensional wakes behind bluff bodies support robust self-sustained vortex shedding for sufficiently large Reynolds numbers. The onset of self-sustained oscillations is associated to a global instability supported by a finite region of local absolute instability in the near wake. [START_REF] Chomaz | Bifurcations to local and global modes in spatially developing flows[END_REF][START_REF] Monkewitz | The absolute and convective nature of instability in two-dimensional wakes at low Reynolds numbers[END_REF][START_REF] Huerre | Local and global instabilities in spatially developing flows[END_REF] There is a continued interest in controlling vortex shedding because, in addition to inducing unsteady loads on the body, it also leads to an increase of the mean drag.

Spanwise periodic (3D) perturbations of spanwise uniform (2D) wakes, e.g., obtained with periodic modulations of the trailing and/or leading edge of the bluff body [START_REF] Tanner | A method of reducing the base drag of wings with blunt trailing edges[END_REF][START_REF] Tombazis | A study of three-dimensional aspects of vortex shedding from a bluff body with a mild geometric disturbance[END_REF][START_REF] Bearman | Reduction of bluff-body drag and suppression of vortex shedding by the introduction of wavy separation lines[END_REF][START_REF] Darekar | Flow past a square-section cylinder with a wavy stagnation face[END_REF] or spanwise periodic blowing and suction [START_REF] Kim | Distributed forcing of flow over a circular cylinder[END_REF] can attenuate and even suppress vortex shedding and reduce the associated undesired drag and unsteady loads (see, e.g., Ref. 9 for a review). Recently, important progress has been made in the understanding of this stabilizing action from a linear stability perspective: Hwang et al. [START_REF] Hwang | Stabilization of absolute instability in spanwise wavy two-dimensional wakes[END_REF] show that appropriate 3D spanwise periodic perturbations of 2D absolutely unstable wake profiles lead to a reduction of the absolute growth rate. This reduction is observed for a range of spanwise wavelengths that is in accordance with experimental results, just as the fact that varicose perturbations are more stabilizing than sinuous ones. However, the question of the higher efficiency of 3D perturbations when compared to 2D ones in reducing the absolute growth rate was left partially open by this study. It is indeed known that lower rates of 3D blowing and suction, compared to 2D one, are required to suppress shedding in a cylinder wake. [START_REF] Kim | Distributed forcing of flow over a circular cylinder[END_REF] This seems to contrast the fact that the first-order sensitivity of the absolute instability growth rate with respect to 3D spanwise periodic modifications of the basic flow is zero, [START_REF] Hwang | Stabilization of absolute instability in spanwise wavy two-dimensional wakes[END_REF][START_REF] Hwang | Control of absolute instability by basic-flow modification in a parallel wake at low Reynolds number[END_REF] therefore predicting that, at first order, 2D perturbations are more effective than 3D ones in reducing the absolute growth rate. A partial explanation of the higher efficiency of 3D perturbations when compared to 2D ones has been given in another recent study [START_REF] Del Guercio | Stabilizing effect of optimally amplified streaks in parallel wakes[END_REF] where we show that parallel "frozen" 2D wakes can support the large temporal amplification of streamwise streaks from stable spanwise periodic and streamwise uniform streamwise vortices via the lift-up effect. [START_REF] Ellingsen | Stability of linear flow[END_REF][START_REF] Landahl | A note on an algebraic instability of inviscid parallel shear flows[END_REF] The optimal perturbations leading to the optimal amplification of the streaks were computed and it was shown that varicose streaks of relatively small amplitude are able to completely quench the absolute instability. [START_REF] Del Guercio | Stabilizing effect of optimally amplified streaks in parallel wakes[END_REF] It was also shown that the initial amplitude of optimal 3D perturbations necessary to quench the absolute instability is much smaller than the initial amplitude required by 2D perturbations.

Many questions were however left unanswered by the local temporal analysis developed in our previous investigation. [START_REF] Del Guercio | Stabilizing effect of optimally amplified streaks in parallel wakes[END_REF] For instance, can optimal spatial amplifications be large in spatially diffusing wakes? The answer is not a priori clear because the wake diffusion not only reduces the basic flow shear fuelling the transient growth but also increases the local spanwise wavenumber of the perturbation, which is known from local analysis to reduce the growth. Another question is: are 3D optimal perturbations more efficient than 2D ones in stabilizing a global instability? The answer to this question is not obvious because finite downstream distances are needed to attain the maximum energy growth, while the pocket of absolute instability that needs to be controlled is located upstream, and therefore it is not clear how efficient optimal perturbations can be in quenching the absolute instability. The scope of the present study is to answer these questions by considering the spatial optimal perturbations and their influence on the global stability of non-parallel wakes.

An "artificial" wake, left free to spatially develop downstream the enforced inflow wake profile, is introduced as reference 2D basic flow in Sec. III. The use of such a basic flow allows us to find results which are independent of the specific body shape generating the wake and of the particular devices used to generate the optimal perturbations. The optimal spatial perturbations of this non-parallel wake are computed in Sec. IV following the procedure described in Sec. II. These optimal perturbations are defined as the perturbation profiles enforced at the inflow station that lead to the optimal energy amplification G(x) at the downstream station x. This definition is quite different from that of optimal initial or inflow conditions leading to the optimal temporal energy amplification 15, 16 G(t). Optimal spatial energy amplifications have already been computed in non-parallel boundary layers by using direct-adjoint methods exploiting the parabolic nature of the boundary layer equations. [START_REF] Andersson | Optimal disturbances and bypass transition in boundary layers[END_REF][START_REF] Luchini | Reynolds-number independent instability of the boundary layer over a flat surface. Part 2: Optimal perturbations[END_REF] Here we choose to specifically design an alternative scalable optimization method (see Sec. II) that does not rely on the parabolic nature of the equations and that does not require the explicit computation of adjoint operators. The influence of forcing optimal perturbations on the global linear stability is investigated in Sec. V. The results of fully nonlinear simulations that validate these results in the nonlinear regime are reported in Sec. V D, while the used numerical methods are summarized in the Appendix.

II. PROBLEM FORMULATION

In this section the mathematical formulation of the analysis performed in the paper is briefly introduced. The formulation is general and can be applied to other non-parallel shear flows. Specific details about the particular wake profile and perturbations used in this study are mentioned in Secs. III-V.

A reference two-dimensional (2D) non-parallel plane basic flow U 2D (x, y) is obtained as a steady solution of the Navier-Stokes equations with inflow boundary condition U = U 0 (y)e x given at x = 0 and free-stream conditions U → U ∞ e x given as y → ±∞. We denote by x, y, and z the streamwise, cross-stream, and spanwise coordinates and by e x , e y , e z the associated unit vectors. The Reynolds number Re = U * re f δ * /ν is based on the characteristic velocity and length associated to U 0 (y) and on the kinematic viscosity ν of the fluid. The non-parallel basic flow U 2D = U (x, y)e x + V (x, y)e y is invariant to translations and reflections in the spanwise coordinate z (it is therefore two-dimensional or 2D). Perturbations u to the reference 2D basic flow are ruled by the Navier-Stokes equations in perturbation form:

∇ • u = 0, (1)
∂u ∂t + (∇U) u + ∇u U + ∇u u = -∇ p + 1 Re ∇ 2 u , (2)
using U = U 2D as basic flow.

In the first part of the study, dealing with optimal spatial perturbations of U 2D , we consider steady perturbations u of U 2D obtained by perturbing the inlet profile U 0 (y) with steady inflow perturbations u 0 (y, z). We are interested in steady perturbations both because they are spatially stable and because they are of interest in passive control applications. In particular, spanwise periodic perturbations of wavelength λ z will be considered in the following. Considering small perturbations, the nonlinear term (∇u )u can be neglected, which makes the perturbation equations linear. Defining the local perturbation kinetic energy density as

e (x) = 1 2δλ z ∞ -∞ λ z 0 u • u dy dz, (3)
the optimal spatial energy amplification of inflow perturbations is defined as

G(x) = max u 0 e (x) e 0 . ( 4
)
Different approaches can be used to compute G(x) and the associated optimal inflow perturbation.

We choose here to decompose the inlet perturbation on a set of linearly independent functions b (m) 0 , in practice limited to M terms, as

u 0 (y, z) = M m=1 q m b (m) 0 (y, z). ( 5
)
Denoting by b (m) (x, y, z) the perturbation velocity field obtained using b (m) 0 (y, z) as inlet perturbation, from the linearity of the operator follows that

u (x, y, z) = M m=1 q m b (m) (x, y, z), (6)
where the coefficients q m are the same used in Eq. ( 5). The optimization problem in Eq. ( 4) can therefore be recast in terms of the M-dimensional control vector q as

G(x) = max q q T H(x)q q T H 0 q , ( 7
)
where the components of the symmetric matrices H(x) and H 0 are

H mn (x) = 1 2δλ z ∞ -∞ λ z 0 b (m) (x, y, z) • b (n) (x, y, z) dy dz, ( 8
)
H 0,mn = 1 2δλ z ∞ -∞ λ z 0 b (m) 0 (y, z) • b (n) 0 (y, z) dy dz. ( 9
)
Within this formulation G(x) is easily found as the largest eigenvalue μ max of the generalized M × M eigenvalue problem μH 0 w = Hw. The corresponding eigenvector is the optimal set of coefficients q (opt) maximizing the kinetic energy amplification at the selected streamwise station x.

The corresponding inlet perturbation is u

(opt) 0 (y, z) = M m=1 q
(opt) m b (m) 0 (y, z). In the limit M → ∞ the approximated solution converges to the exact solution.

In the second part of the study, the influence of forcing finite amplitude optimal perturbations on linear global stability is investigated. A family of 3D streaky nonlinear non-parallel basic flows U 3D (x, y, z; A 0 ) is obtained by looking for steady solutions of the (nonlinear) Navier-Stokes equations with inflow boundary condition U 0 (y, z) = U 0 (y)e x + A 0 u (opt) 0 (y, z) given at x = 0. u (opt) 0 is normalized to unit x-local energy so that e (0) = A 2 0 for u = U 3D -U 2D (see the definition of e in Eq. ( 3)). In general, it is not guaranteed that the linear optimal perturbations are also (nonlinearly) optimal at finite amplitude. However, for the present purpose of open-loop control, this is not a problem as long as they are still largely amplified. Using strictly optimal perturbations is also not critical because it is not likely that strictly optimal perturbations can be forced in a real flow and there is no guarantee that they would be also the optimal ones in reducing the global mode growth rate. The global linear stability of the U 3D basic flows is then analysed by integrating in time the linearized form of the Navier-Stokes equations ( 1) and ( 2) in perturbation form with U = U 3D . After the extinction of transients, the leading global mode emerges inducing an exponential growth or decay of the solution. The global growth rate is then deduced from the slope of the global energy amplification curve.

III. NON PARALLEL 2D REFERENCE WAKE

The 2D reference wake is computed by enforcing as inflow boundary condition the following well studied 2 wake profile:

U 0 (y) = 1 + 2 1 + sinh 2N (y sinh -1 1) -1 , ( 10
) with = (U * c -U * ∞ )/(U * ∞ + U * c )
, where U * c is the centreline and U ∞ * the freestream velocity (dimensional variables are starred). The velocity U 0 is made dimensionless with respect to the reference velocity

U * re f = (U * c + U * ∞ )/2.
The spatial coordinates are made dimensionless with respect to the reference length δ * 0 that is the distance from the centreline to the point where the 2D wake velocity is equal to U * re f , computed at the inflow. We set = -1.35 to ensure a small recirculation in the upstream region of the wake. For = -1.35, the wake is globally unstable when Re 39 (not shown). In the following we will consider the value Re = 50 for which numerical simulation (see the Appendix for the numerical details) shows strong self-sustained oscillations in the wake (see Fig. 1(b)). As at Re = 50 the only unstable global mode is sinuous (antisymmetric with respect to the y = 0 axis), the unstable basic flow U 2D is computed by direct temporal integration by enforcing the y-symmetry of the solutions (otherwise a Newton-based continuation method would have been required). The reference basic flow is shown in Fig. 1(a). It can be seen how the basic flow vorticity, which is maximum at x = 0 with peaks at y ≈ ±1, slowly diffuses downstream.

IV. OPTIMAL SPATIAL ENERGY GROWTH

Optimal steady inlet perturbations of U 2D maximizing the spatial energy amplification G(x) are computed following the procedure described in Sec. II. Our previous investigation of the optimal temporal energy growth in parallel wakes [START_REF] Del Guercio | Stabilizing effect of optimally amplified streaks in parallel wakes[END_REF] has shown that the most amplified spanwise periodic and streamwise uniform (corresponding to steady in our spatial framework) perturbations consist in streamwise vortices inducing the growth of streamwise streaks. We therefore consider inlet conditions of the type: u 0 = (u 0 , v 0 , w 0 ) = (0, ∂ψ/∂z, -∂ψ/∂y). Single-harmonic spanwise periodic perturbations can be considered without loss of generality: ψ = f(y)sin (βz). As varicose perturbations (mirror-symmetric with respect to the y = 0 plane) are the most efficient for control, [START_REF] Choi | Control of flow over a bluff body[END_REF][START_REF] Hwang | Stabilization of absolute instability in spanwise wavy two-dimensional wakes[END_REF][START_REF] Del Guercio | Stabilizing effect of optimally amplified streaks in parallel wakes[END_REF] even if they are slightly less amplified than sinuous ones, [START_REF] Del Guercio | Stabilizing effect of optimally amplified streaks in parallel wakes[END_REF] we enforce f (-y) = -f(y) which leads to varicose streaks. The set of linearly independent inflow conditions used in Eq. ( 5) is chosen as b (m) 0 (y, z) = (0, ∂ψ (m) /∂z, -∂ψ (m) /∂ y) with ψ (m) = f m (y)sin (βz) and f m (y) = -f m (y) for m = 1, . . . , M. We have found well suited the set f m (y) = sin (2mπ y/L y ), where the numerical box extends from -L y /2 to L y /2 in the y direction.

Optimal energy growths have been computed for a set of spanwise wavenumbers β increasing M until a precision of 1% or higher on G max was achieved (see also the Appendix for the numerical details of the computations). The computed optimal energy growth curves G(x, β) are reported in Fig. 2(a). It is seen how, consistently with results form the local analysis, [START_REF] Del Guercio | Stabilizing effect of optimally amplified streaks in parallel wakes[END_REF] both the maximum growth G max = max x G(x) and the position x max where it is attained increase with increasing spanwise wavelength λ z = 2π/β, i.e. with decreasing β (see also panel (b) of the same figure). The convergence of the optimal growth curves with increasing M is quite fast, and this for all the considered values of β, as can be seen in Fig. 3. Well converged results, with relative variations below 1% are obtained with only M = 16 terms.

z y -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 FIG. 4.
Cross-stream (y-z) view of the cross-stream v 0 -w 0 components of optimal vortices (arrows) forced at the inflow (x = 0) and of the streamwise u component of the corresponding optimal streaks (contour-lines) at x = x max for Re = 50, β = 1. The 2D basic flow wake streamwise velocity at the inflow U 0 (y) is also reported in grey-scale with white corresponding to the freestream velocity and dark grey the minimum velocity (wake centreline).

The optimal inflow perturbations (x = 0) and the maximum response (x = x max ) associated to the maximum growth G max obtained for β = 1 are reported in Fig. 4. The corresponding velocity profiles are reported in Fig. 5, where additional values of β are also considered. The optimal inflow perturbations consist in two rows of counter-rotating vortices on each side of the y = 0 plane, with opposite rotation on each side. These vortices induce the growth of y-symmetric (varicose) streaks. From Fig. 5 it can be seen how, for increasing spanwise wavelengths λ z (decreasing β), the size of optimal perturbations increases in the normal (y) direction (and of course also in the spanwise z direction).

The observed trends are in agreement with those found in our previous local analysis. [START_REF] Del Guercio | Stabilizing effect of optimally amplified streaks in parallel wakes[END_REF] However, the maximum spatial growth rates obtained in the non-parallel case are smaller than the temporal ones obtained at the same nominal β under the frozen and parallel flow approximation. This is not surprising because the nominal values of β and Re of the non-parallel results are based on the properties of the wake profile at the inflow (x = 0). As the dimensional reference length δ * (x) (the y * value where U * 2D (y * ) = U * re f ) increases with x, a dimensionless wavenumber β = β * δ * based on the local scale would increase going downstream. As the maximum growth rate is a decreasing function of β, it is not surprising that the maximum growth rates are smaller than the ones that would be obtained if the wake was parallel. Therefore, the results of the present analysis should be compared to the ones of the local analysis obtained at larger values of β.

V. STABILIZING EFFECT OF OPTIMALLY FORCED STREAKS ON GLOBAL MODES

In this second part of the study, we investigate the influence of the forcing of optimal perturbations on the linear global stability of the wake. The input parameter of this analysis is the amplitude A 0 of the forcing at the inflow boundary and the output is the linear growth rate s r of the global mode supported by the streaky wake. All the results are obtained for Re = 50 and β = 1. The choice of β = 1 is not completely arbitrary. On the one hand, in order to obtain large energy amplifications, one should choose low values of β. For low β, however, not only would the cross-stream size of the inlet optimal vortices be probably too large to be implemented in practical applications but, even more importantly, large amplitudes of the streaks would be obtained only far downstream (e.g., x max ≈ 120 for β = 0.5). This is a problem because the main scope of the control is to reduce the absolute growth rate in the absolute region which extends up to x ≈ 5. On the other hand, selecting large values of β, in order to have x max in the absolute region, would lead to poor energy amplifications and would exclude any damping in the convective region. The value β = 1 is a good compromise between these two extrema. In particular, as can be seen from Fig. 2(a), the obtained G(x) for β = 1 in the region x 10 are sensibly the same of those obtained for higher β.

A. Streaky wakes basic flows

Non-parallel streaky (3D) wake basic flows U 3D (x, y, z; A 0 ) are computed by enforcing at x = 0 the inflow condition U = U 0 (y)e x + A 0 u (opt) 0 and by then computing the corresponding steady solution of the (nonlinear) Navier-Stokes equations, as explained in Sec. II. The solution, which may be unstable, is obtained by enforcing symmetry with respect to the y = 0 plane, exactly as done to compute U 2D . The local amplitude of the streaks is measured extending the standard definition used in previous studies: [START_REF] Del Guercio | Stabilizing effect of optimally amplified streaks in parallel wakes[END_REF][START_REF] Andersson | On the breakdown of boundary layers streaks[END_REF] A s

(x) = 1 2 max y,z U 3D -U 2D -min y,z U 3D -U 2D max y U 0 -min y U 0 . ( 11
)
In this definition the streak amplitude at the station x is defined as half the maximum deviation of the streaky 3D profile from the reference 2D profile, at the same x station, normalized by the maximum velocity variation of the inflow 2D reference profile. The considered values of A 0 and the obtained values of A s at x max (maximum value of A s ) and in the middle of the absolute instability region (x = 2.7) are reported in Table I, where each considered case is given a literal label. Case A corresponds to the reference two-dimensional wake profile U 2D (no streaks) while cases B, C, D, and E are obtained by increasing the inlet amplitude A 0 of the forced optimal perturbations. The nonlinear streaks amplitude evolution A s (x) associated to the velocity fields U 3D are reported in Fig. 6(a) for the considered cases. From Fig. 7, where streaky basic flows are shown in the symmetry plane y = 0, it is seen how, indeed for increasing A 0 , the wake is increasingly 3D. The effect of nonlinearity is to slightly reduce the maximum energy growth (from ≈20 in the linear small amplitude limit to ≈17-15 for streaks D and E, not shown) and to induce a mean flow distortion that slightly counteracts the effect of the streaks.

B. Linear global stability analysis

The linear global stability analysis is performed via a direct numerical simulation of the Navier-Stokes equations ( 1) and ( 2) linearized upon the 3D streaky basic flows U 3D defined above (Sec. V A). In previous local stability analyses [START_REF] Hwang | Stabilization of absolute instability in spanwise wavy two-dimensional wakes[END_REF][START_REF] Del Guercio | Stabilizing effect of optimally amplified streaks in parallel wakes[END_REF] the dominant absolute mode is subharmonic, i.e., its spanwise wavelength is twice that of the basic flow streaks. This is taken in due consideration by integrating the linearized equations in a domain including two basic flow streaks wavelengths (L z = 2λ z ). Noisy initial conditions on u are given for the reference 2D wake (case A). The unstable global mode emerging for large times in the reference 2D wake is then used, upon normalization of its amplitude, as initial condition on u in simulations with the increasingly streaky basic flows B,...,E.

The temporal evolution of the global perturbation kinetic energy

E = 1 2δL x L z L x 0 L y /2 -L y /2 L z 0 u • u dx dy dz
is reported in Fig. 6(b). After an initial transient extending to t ≈ 70, the dependence of E on time is exponential (a straight line in the lin-log scales used in the figure ), where the rate of growth or of decay is twice the growth rate of the global mode.

As anticipated (see, e.g., Fig. 1(b)) the reference 2D wake (case A) is strongly linearly unstable at Re = 50. The forcing of 3D linearly optimal perturbations of increasing amplitude has a stabilizing effect on the global instability. The growth rate first reduced for low amplitude streaks (cases B and C) is then rendered quasi-neutral (case D) and finally completely stable for sufficiently large streak amplitudes (case E).

In the neutral and stable case the streaks amplitudes A s (x = 2.7), measured in the middle of the absolute instability region of the reference 2D wake, are respectively of ≈5% and ≈7%. These values are not far from the ≈8% value at which the absolute instability was completely quenched in our previous local stability analysis. [START_REF] Del Guercio | Stabilizing effect of optimally amplified streaks in parallel wakes[END_REF] Also remark that, in the present non-parallel case A s is given in terms of the entrance reference maximum U 2D (x = 0), but if it was based on the local value of U 2D (x) which is decreased with x, this would result in even larger downstream values of A s . The stabilization of the global mode therefore appears to be associated to a strong reduction of the pocket of absolute instability that drives the global mode oscillations in the 2D reference case. A local stability analysis of the basic flow profiles extracted at x = 2.7 (not shown) indeed confirms that the local absolute growth rate is reduced with increasing streak amplitudes and that it is completely quenched by streak E.

C. Sensitivity of the global growth rate to the amplitude of 3D optimal structures 1. The first order sensitivity of the global growth rate to streaks is zero

In previous studies based on local stability analyses, [START_REF] Hwang | Stabilization of absolute instability in spanwise wavy two-dimensional wakes[END_REF][START_REF] Hwang | Control of absolute instability by basic-flow modification in a parallel wake at low Reynolds number[END_REF] it was shown that the sensitivity of the absolute growth rate to 3D spanwise periodic modifications of the basic flow is zero. The argument developed in the local absolute instability analysis is easily extended to the global stability analysis of the nonparallel wake and proceeds as follows. Denote by L 2D the Navier-Stokes operator linearized near the 2D basic flow U 2D . For small values of the inflow amplitude of optimal perturbations, the basic flow is modified by a small amount δU = U 3D -U 2D that induces a small change δL in the linear operator. At first order, the change of the leading eigenvalue induced by this small variation is [START_REF] Bottaro | The effect of base flow variation on flow stability[END_REF][START_REF] Chomaz | Global instabilities in spatially developing flows: Nonnormality and nonlinearity[END_REF]

δs = w † 2D , δLw 2D / w † 2D , w 2D
, where w 2D , s 2D , w † 2D are, respectively, the leading global mode, eigenvalue, and adjoint global mode associated to U 2D , and the standard inner product is defined as a, b = 2) it is seen that the variation δL induced by δU consists only in spanwise periodic terms as δU is itself spanwise periodic. As w 2D and w † 2D do not depend on z, it follows 22 that w † , δLw = 0 and therefore that δs = 0. This is not the case for 2D (spanwise uniform) perturbations of U for which the variation of the leading eigenvalue is, in general, non-zero.

L x 0 L y /2 -L y /2 L z 0 a • b dx dy dz. From Eq. (

Effective sensitivity of global growth rate to spanwise periodic basic flow modifications and comparison with 2D modifications

We now consider the observed dependence of the most unstable global mode on the control amplitude. Such a control amplitude is unequivocally defined in terms of inflow optimal perturbation FIG. 8. Dependence of the growth rate of the global eigenvalue s r on the inflow optimal perturbation amplitude A 0 (panel (a)) and on the streak amplitude A s (x = 2.7) measured in the centre of the absolute region of the reference 2D wake (panels (b) and (c) for a zoomed plot). A spanwise uniform perturbation (2D) has been also considered for comparison. Symbols denote data points, while lines are best fits to the data points.

amplitudes A 0 . The dependence of the global growth rate s r on the inflow optimal perturbation amplitude A 0 is displayed in Fig. 8(a). If this dependence is also to be reported in term of streaks amplitudes, for the considered streaks with β = 1, it makes no sense to report it in terms of A s,max because this value is attained far downstream, in the convectively unstable region. We instead take as an indicator of the "useful" streak amplitude the amplitude of the streaks in the middle of the absolute region of the unperturbed flow A s (x = 2.7). The dependence of s r on this amplitude is reported in Fig. 8 (panels (b) and (c) for a zoom).

In the same figures the variation of the growth rate s r induced by a 2D perturbation of the basic flow is also reported for comparison. The 2D perturbation has the same y shape as the optimal streak shape in the middle of the absolute region (x = 2.7) but is uniform instead of periodic in the spanwise direction. For this 2D perturbation, A 0 is unambiguously defined and A s is defined as the maximum associated U taken at x = 2.7.

From the figures it is clearly seen how the first order sensitivities ds r /dA 0 and ds r /dA s computed for A 0 = A s = 0 are zero for the 3D perturbations and non-zero for the 2D perturbations as predicted by the first order sensitivity analysis. According to a first-order sensitivity analysis one would expect the 2D perturbations to be more effective than 3D ones in quenching the global instability, but exactly the opposite is observed. Indeed, 2D perturbations are more effective than 3D ones in reducing s r only for very small perturbation amplitudes, while the opposite is observed for larger amplitudes where the higher order dependence of s r on A 0 and A s induces more important reductions of s r . We indeed find that 3D perturbations stabilize the global mode at a value of A s (x = 2.7) more than five times smaller, and more than ten times smaller in terms of A 0 . A higher efficiency of 3D perturbations was expected for results expressed in terms of A 0 , due to the gain associated with the lift-up of the 3D optimal perturbations. However, such a result was somehow unexpected when expressing the growth rate reductions in terms of A s .

D. Nonlinear simulations

Non-linear simulations of the full Navier-Stokes equations have finally been performed to assess the effect of the inflow forcing of 3D optimal perturbations in the nonlinear regime. The same grid used in linear simulations has been used in the nonlinear ones. In a first simulation, the permanent harmonic self-sustained state supported by the reference 2D wake is allowed to develop. This 2D (spanwise uniform) self-sustained state is then given as an initial condition to simulations in the presence of the optimal perturbations (streaky wakes) of increasing amplitude. As expected from the linear analysis, the global perturbation kinetic energy E associated to the self-sustained oscillations in the wake is reduced when the amplitude of the enforced optimal perturbations is increased (see Fig. 9). A stable steady streaky wake is found for case E, where the oscillations are completely suppressed.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: Temporal evolution E (t) of the total perturbation kinetic energy, integrated over the whole computational box, supported by the streaky wakes, normalized to the E A of the reference 2D wake. The results are issued from nonlinear simulations where the permanent periodic state supported by the reference 2D wake is given as initial condition. In the presence of optimal perturbations of increasing amplitude, the amplitude of self-sustained oscillations is initially reduced (cases B, C, D), up to their complete suppression (case E).

Snapshots of the perturbation streamwise velocity in the y = 0 plane are reported in Fig. 10 for all the considered cases. For the reference 2D wake (case A), the self-sustained state is spanwise uniform (2D) with structures corresponding to standard von Kármán vortices. These vortical structures become increasingly modulated in the spanwise direction for increasing amplitudes A 0 of the forcing. Unsteady structures are completely suppressed for case E, where the basic flow streamwise streaks remain the only visible structures in the wake.

VI. SUMMARY AND DISCUSSION

In this study the optimal amplifications supported by an "artificial" non-parallel unstable 2D wake have been computed and their influence on the stability of the wake have been investigated by a global stability analysis. The main results can be summarized as follows:

r The energy of steady, symmetric spanwise periodic streamwise vortices forced at the inflow boundary can be significantly amplified downstream leading to large amplitude varicose streamwise streaks. r An increase of the spanwise wavelength λ z of the perturbations leads to larger energy amplification and to taller (in y) optimal structures. r The used optimization technique, based on the simulation of the responses to a set of linearly independent inflow forcings, has proved very flexible. Only 16 simulations of independent forcings were needed to obtain accuracies higher than 1% on the optimal energy growths. r The unstable global mode of the reference 2D wake at Re = 50 is completely stabilized when optimal inflow perturbations (vortices) are forced with sufficiently large amplitude. r The results of first order local sensitivity analyses [START_REF] Hwang | Stabilization of absolute instability in spanwise wavy two-dimensional wakes[END_REF][START_REF] Hwang | Control of absolute instability by basic-flow modification in a parallel wake at low Reynolds number[END_REF] are easily extended to the non-parallel case to show that the sensitivity of the 2D global mode eigenvalue to 3D spanwise periodic modifications of the basic flow is zero, while it is in general non-zero for 2D modifications. r 3D optimal perturbations require smaller amplitudes than a reference 2D forcing to quench the global instability, and this both in terms of rms-amplitude of the inflow forcing and in terms of the basic flow distortion amplitude measured in the centre of the absolute instability region of the reference 2D wake. This is in contrast with the prediction of the sensitivity analysis.

The trends observed for the optimal energy amplification and the associated optimal perturbations in the non-parallel case are in qualitative agreement with those found in the local stability analysis. [START_REF] Del Guercio | Stabilizing effect of optimally amplified streaks in parallel wakes[END_REF] In particular, the shapes of the optimal inputs (streamwise vortices) and those of the optimal outputs (streamwise streaks) are very similar for both analyses. When comparing the results, care must however be exerted in, e.g., selecting the appropriate spanwise wavenumbers to compare, as the local dimensionless wavenumber keeps increasing while going downstream in the non-parallel wake. Using the wavenumber made dimensionless with respect to the inflow reference length, the growths in the non-parallel wake appear smaller than the ones that would be predicted by keeping the wake frozen and parallel. The fact that the linear global instability can be suppressed by optimal spanwise perturbations is in agreement with the idea that these 3D perturbations efficiently reduce the local absolute growth in the wave-maker region of the flow. [START_REF] Hwang | Stabilization of absolute instability in spanwise wavy two-dimensional wakes[END_REF][START_REF] Del Guercio | Stabilizing effect of optimally amplified streaks in parallel wakes[END_REF] It also extends to flows with an "oscillator" dynamics [START_REF] Huerre | Local and global instabilities in spatially developing flows[END_REF][START_REF] Chomaz | Global instabilities in spatially developing flows: Nonnormality and nonlinearity[END_REF] the control strategy based on the forcing of optimally amplified streaks that has been successfully used to stabilize convectively unstable waves in non-parallel boundary layers. [START_REF] Cossu | Stabilization of Tollmien-Schlichting waves by finite amplitude optimal streaks in the Blasius boundary layer[END_REF][START_REF] Cossu | On Tollmien-Schlichting waves in streaky boundary layers[END_REF][START_REF] Fransson | Experimental study of the stabilisation of Tollmien-Schlichting waves by finite amplitude streaks[END_REF][START_REF] Fransson | Delaying transition to turbulence by a passive mechanism[END_REF] In this type of control strategy, optimal vortices are forced which then efficiently generate the streaks leading in fine to the stabilization. This control technique is much more efficient than directly forcing the stabilizing streaks because the lift-up effect is used as an amplifier of the control action. Otherwise, much larger forcing amplitudes would be required to directly force the streaks.

The conclusions concerning the sensitivity analysis are, we believe, probably the most relevant of this study. According to a first order sensitivity analysis, the 3D spanwise periodic forcing or modification of the basic flow with amplitude A is less effective than a 2D one with the same amplitude because the sensitivity of 3D perturbations is zero (growth rate reductions ∼A 2 , with zero derivative in A = 0), while the 2D sensitivity is not zero (growth rate reductions ∼A). While these conclusions are correct for very small amplitudes A of the basic flow modifications, they are not correct for larger amplitudes where the parabola-shaped growth rate reductions (3D control) have grown larger than the straight line ones (2D control). In our specific case this cross-over happens at very small amplitudes of the forcing, of the order of 1/10 of the amplitude required for stabilization by 3D modifications and of the order of 1/100 of the one required by 2D modifications. Considering that these results have been obtained at a Reynolds number only 20% above the critical value for global instability, this means that except in very weakly unstable situations, where small control amplitudes suffice to stabilize the perturbations, a first order sensitivity leads to misleading conclusions when the stabilizing efficiency of 3D and 2D perturbations is compared.

An important question, left for future study, is how to force optimal perturbations in the presence of the bluff body. As already mentioned, many ways to modulate wakes in the spanwise direction have been already implemented, among which, the sinusoidal indentation of the leading and/or the trailing edge of the body, [START_REF] Tanner | A method of reducing the base drag of wings with blunt trailing edges[END_REF][START_REF] Tombazis | A study of three-dimensional aspects of vortex shedding from a bluff body with a mild geometric disturbance[END_REF][START_REF] Bearman | Reduction of bluff-body drag and suppression of vortex shedding by the introduction of wavy separation lines[END_REF][START_REF] Darekar | Flow past a square-section cylinder with a wavy stagnation face[END_REF] and the spanwise periodic blowing and suction at the wall surface. [START_REF] Kim | Distributed forcing of flow over a circular cylinder[END_REF] The achieved wake modulations are strikingly similar to the streaky wakes investigated in the present study, suggesting that, just like in boundary layers, the optimal streaks represent a sort of "attractor" of the spanwise modulated solutions. However, it is not clear which of these strategies, if any, predominantly uses vortices to force the streaks instead of directly forcing the streaks. It would also be interesting to investigate the optimal amplification and the control efficiency of periodic inflow perturbations. Additional work is granted on these issues.

APPENDIX: NUMERICAL SIMULATIONS

The Navier-Stokes equations (nonlinear and linearized) have been numerically integrated using OpenFoam, an open-source finite volumes code (see http://www.openfoam.org). The flow is solved inside the domain [0, L x ] × [-L y /2, L y /2] × [0, L z ] that is discretized using a grid with N x and N z equally spaced points in the streamwise and spanwise directions, respectively. N y points are used in the y direction using stretching to densify points in the region where the basic flow shear is not negligible. The fractional step, pressure correction PISO (Pressure Implicit with Splitting of Operators) scheme is used to advance the solutions in time. Different grids have been used to compute optimal linear perturbations of different spanwise wavenumbers, as reported in Table II.

Nonlinear streaky basic flows for β = 1 have been computed using the same grid used for the computation of the linear optimals at the same β. For the linear and nonlinear simulations of the perturbations to the 3D streaky basic flows, however, the domain is doubled in the y direction ([ -L y /2, L y /2] instead of [0, L y /2]) as the y symmetry is no more enforced. Te box is also doubled in the spanwise direction (L z = 2λ z ) in order to include subharmonic perturbations. The corresponding N y and N z are also doubled leading to a grid with L x = 124, L y /2 = 10, L z = 4π , N x = 300, N y = 160, N z = 48 with x = 0.4, z = 0.26 and a minimum y = 0.01 on the symmetry axis and a maximum y = 0.1 near the freestream boundary.

10 FIG. 1 .FIG. 2 .

 1012 FIG. 1. Spanwise vorticity fields ω z (x, y) associated to the reference 2D non-parallel wake at Re = 50. (a) (unstable) Basic 2D flow profile obtained by enforcing the y-symmetry of the solution. (b) Snapshot of the periodic self-sustained state obtained without enforcing the y-symmetry of the solution.

FIG. 3 .

 3 FIG.3. Convergence of the optimal energy growth G(x) for β = 1 (panel (a)) and of the maximum energy growth G max (β) (panel (b)) when the number M of linearly independent inflow conditions is increased at Re = 50. Well converged results are obtained for M = 16.

FIG. 5 .

 5 FIG.5. Normalized amplitude of the v(x = 0, y, z = 0) component (panel (a)) and w(x = 0, y, z = λ z /4) components (panel (b)) of the optimal inflow boundary vortices. The normalized amplitude of the u(x = x max , y, z = 0) streamwise component of the corresponding optimally amplified streaks is plotted in panel (c). Three selected spanwise wavenumbers are considered: β = 0.5 (dashed line, green), β = 1 (solid line, red), and β = 1.5 (dotted line, blue).

FIG. 6 . 5 FIG. 7 .

 657 FIG. 6. Spatial evolution of the streaks amplitudes A s (x) for increasing amplitudes A 0 of the inflow optimal perturbations (panel (a)) and temporal evolution of the global kinetic energy of secondary perturbations E (t) to the considered reference and streaky basic flows (panel (b)). All the results have been obtained for Re = 50 and β = 1.

  FIG.9. Temporal evolution E (t) of the total perturbation kinetic energy, integrated over the whole computational box, supported by the streaky wakes, normalized to the E A of the reference 2D wake. The results are issued from nonlinear simulations where the permanent periodic state supported by the reference 2D wake is given as initial condition. In the presence of optimal perturbations of increasing amplitude, the amplitude of self-sustained oscillations is initially reduced (cases B, C, D), up to their complete suppression (case E).

5 FIG. 10 .

 510 FIG.10. Snapshots from fully nonlinear simulations. Streamwise perturbation velocity u (x, y = 0, z) = u(x, y = 0, z) -U 2D (x, y) in the y = 0 symmetry plane in the permanent regime (t = 250). The reference 2D case A is reported in panel (a), while cases B, C, D, and E, obtained by increasing the amplitude A 0 of the inflow optimal perturbations, are reported in panels (b) to (e) (top to bottom). Case A displays self-sustained periodic oscillations of 2D structures in the wake. These structures become increasingly 3D and of smaller rms value for increasing values of the enforced A 0 (cases B to D). The oscillations are completely suppressed in case E where the stable streaky basic flow is observed after transients are extinguished.

  ACKNOWLEDGMENTS G.D.G. acknowledges the support of ANRT via the convention CIFRE 742/2011.

  This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 195.83.231.74 On: Mon, 24 Feb 2014 17:07:04

  it was shown that for large streaks amplitudes,

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 195.83.231.74 On: Mon, 24 Feb 2014 17:07:04

TABLE I .

 I Considered nonlinear streaky wake basic flows. A 0 is the finite amplitude given, at the inflow, to the linear optimal boundary perturbations (vortices). A s,max is the maximum streak amplitude reached in the nonlinear numerical simulation. Case A corresponds to the reference two-dimensional wake. Cases B, C, D, and E are obtained by increasing A 0 .

	Case	A 0	A s,max (%)	A s (x = 2.7) (%)
	A	0.000	0	0
	B	0.057	10.3	2.5
	C	0.085	15.1	3.7
	D	0.120	20.4	5.2
	E	0.171	27.3	7.4

TABLE II .

 II Numerical grids used for the computation of optimal perturbations. As a y-symmetry is enforced the equations are solved only in the [0, L y /2] half-domain.

	β	L z	Nz	Ly/2	Ny	Lx	Nx
	0.5	12.57	48	20	120	124	300
	0.75	8.38	48	20	120	124	300
	1	6.28	24	10	80	124	300
	1.25	5.03	24	10	80	124	300
	1.5	4.20	12	10	80	124	300
	1.75	3.60	12	10	80	124	300

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 195.83.231.74 On: Mon, 24 Feb 2014 17:07:04