2.1.1. Continuous wavelet transform

The wavelet transform (Grossmann & Morlet Reference Grossmann and Morlet1984) unfolds any signal (e.g. in time) or any field (e.g. in three-dimensional space) into both space (or time) and scale (or time scale), and possibly directions (for dimensions higher than one). The building block of the wavelet transform is the ‘mother wavelet’, ${\it\psi}(x)\in L^{2}(\mathbb{R})$ with $x\in \mathbb{R}$ , that is a well-localized function with fast decay at infinity and at least one vanishing moment (i.e. zero mean). It is also smooth enough in order that its Fourier transform, $\hat{{\it\psi}}(k)$ , exhibits fast decay for $|k|$ tending to infinity. From the mother wavelet one then generates a family of wavelets, translated by $b\in \mathbb{R}$ , the position parameter, dilated (or contracted) by $a\in \mathbb{R}^{+}$ , the scale parameter, and normalized in $L^{2}$ -norm (i.e. $\Vert {\it\psi}_{a,b}\Vert _{2}=1$ ) to obtain the set

(2.1) $$\begin{eqnarray}{\it\psi}_{a,b}(x)=\frac{1}{\sqrt{a}}{\it\psi}\left(\frac{x-b}{a}\right).\end{eqnarray}$$

The wavelet transform of $f\in L^{2}(\mathbb{R})$ is the inner product of $f$ with the analysing wavelets ${\it\psi}_{a,b}$ , and the wavelet coefficients that measure the fluctuations of $f$ at scale $a$ and position $b$ given by

(2.2) $$\begin{eqnarray}\widetilde{f}(a,b)=\langle \,f,{\it\psi}_{a,b}\rangle =\int _{\mathbb{R}}f(x){\it\psi}_{a,b}^{\star }(x)\,\text{d}x,\end{eqnarray}$$

with $^{\star }$ denoting the complex conjugate. The function $f$ is reconstructed from its wavelet coefficients, as the inner product of $\widetilde{f}$ with the set of analysing wavelets ${\it\psi}_{a,b}$

(2.3) $$\begin{eqnarray}f(x)=\frac{1}{C_{{\it\psi}}}\int _{\mathbb{R}^{+}}\int _{\mathbb{R}}\widetilde{f}(a,b){\it\psi}_{a,b}(x)\,\frac{\text{d}a\,\text{d}b}{a^{2}},\end{eqnarray}$$

where $C_{{\it\psi}}=\int _{\mathbb{R}^{+}}|\hat{{\it\psi}}(k)|^{2}k^{-1}\,\text{d}k$ is a constant that depends on the wavelet ${\it\psi}$ . Similarly to the Fourier transform, the wavelet transform corresponds to a change of basis (from physical space to wavelet space) and, since it is an isometry, it preserves the inner product $(\langle \,f,g\rangle =\langle \widetilde{f},\widetilde{g}\rangle )$ (Plancherel’s theorem) and conserves energy (Parseval’s identity), therefore

(2.4) $$\begin{eqnarray}\int _{\mathbb{R}}|\,f(x)|^{2}\,\text{d}x=\frac{1}{C_{{\it\psi}}}\int _{\mathbb{R}^{+}}\int _{\mathbb{R}}|\widetilde{f}(a,b)|^{2}\,\frac{\text{d}a\,\text{d}b}{a^{2}}.\end{eqnarray}$$

Note that the wavelet coefficients of the continuous wavelet transform are redundant and therefore correlated. This could be illustrated by the patterns one observes within the continuous wavelet coefficients of white noise, which correspond to the correlation between the dilated and translated wavelets (white noise being decorrelated by construction), which visualise the ‘reproducing kernel’ of the continuous wavelet transform. Due to the fact that wavelets are well localized in physical space, the behaviour of the signal at infinity does not play any role. Therefore both wavelet analysis and wavelet synthesis can be performed locally, in contrast to the Fourier transform which is intrinsically non-local (Fourier modes are spread all over space). One can also construct peculiar wavelets on a dyadic grid ${\it\lambda}=(j,i)$ (i.e. scale is sampled by octaves $j$ and space by positions $2^{-j}i$ ) that are orthogonal to each other and are used to construct wavelet orthonormal bases. In contrast to the continuous wavelet coefficients equation (2.2) that are redundant and correlated, the orthogonal wavelet coefficients are decorrelated and non-redundant (i.e. a signal sampled on $N$ points is perfectly represented by $N$ orthogonal wavelet coefficients only). As for the Fourier transform, there exists a fast wavelet transform (FWT) that is even faster than the fast Fourier transform (FFT) whose operation count for a one-dimensional signal sampled at $N$ points is proportional to $N$ , instead of $N\log _{2}N$ for the FFT.

Figure 1.

Wavelet representation. Physical space (a) and spectral space (b). Note that ${\rm\Delta}x{\rm\Delta}k>C$ is due to the Fourier uncertainty principle.

2.1.2. Orthogonal wavelet transform

A discrete wavelet representation is obtained by sampling dyadically the scale $a$ and the position $b$ introducing $a_{j}=2^{-j}$ and $b_{ji}=\text{i}a_{j}$ with $i,j\in \mathbb{Z}$ . The resulting discrete wavelets

(2.5) $$\begin{eqnarray}{\it\psi}_{ji}(x)=a_{j}^{-1/2}{\it\psi}\left(\frac{x-b_{ji}}{a_{j}}\right)=2^{j/2}{\it\psi}(2^{j}x-i),\end{eqnarray}$$

generate orthogonal bases for peculiar wavelets. Figure 1 shows five discrete wavelets ${\it\psi}_{ji}$ for $j=3,\ldots ,7$ and their corresponding Fourier transforms, the modulus $|\widehat{{\it\psi}}_{ji}|$ . Note that the scale $2^{-j}$ is related to the wavenumber $k_{j}$ as

(2.6) $$\begin{eqnarray}k_{j}=k_{{\it\psi}}2^{j},\end{eqnarray}$$

where $k_{{\it\psi}}=\int _{0}^{\infty }k|\hat{{\it\psi}}(k)|\,\text{d}k/\!\int _{0}^{\infty }|\widehat{{\it\psi}}(k)|\,\text{d}k$ is the centroid wavenumber of the chosen wavelet. In figure 1 we observe the duality between physical and spectral space, namely small-scale wavelets are well localized in physical space and badly localized in spectral space, and vice versa. Denoting the support of a wavelet in physical space by ${\rm\Delta}x$ and the one in spectral space by ${\rm\Delta}k$ the Fourier uncertainty principle requires that the product ${\rm\Delta}x{\rm\Delta}k$ is bounded from below. In this case, the orthogonal wavelet coefficients of a function $f\in L^{2}(\mathbb{R})$ are given by

(2.7) $$\begin{eqnarray}\widetilde{f}_{ji}=\langle \,f,{\it\psi}_{ji}\rangle ,\end{eqnarray}$$

and the corresponding orthogonal wavelet series reads

(2.8) $$\begin{eqnarray}f(x)=\mathop{\sum }_{j,i\in \mathbb{Z}}\widetilde{f}_{ji}{\it\psi}_{ji}(x).\end{eqnarray}$$

The integral in the continuous reconstruction forumla, (2.3), can thus be replaced by a discrete sum. In practical applications, the infinite sums of the wavelet series have to be truncated in both scale and position. Limiting the analysis to the largest accessible scale of the domain, $2^{0}=L$ , the scaling function associated to the wavelet has to be introduced and the wavelet series becomes

(2.9) $$\begin{eqnarray}f(x)=\mathop{\sum }_{i\in \mathbb{Z}}\overline{f}{\it\phi}_{0i}(x)+\mathop{\sum }_{j\geqslant 0,i\in \mathbb{Z}}\widetilde{f}_{ji}{\it\psi}_{ji}(x),\end{eqnarray}$$

where ${\it\phi}$ is the scaling function and $\overline{f}=\langle \,f,{\it\phi}_{0i}\rangle$ the corresponding scaling coefficients. The smallest scale $2^{-J}$ is given by the sampling rate of the function $f$ , which determines the number of grid points $N=2^{J}$ . The finite domain size implies that the number of positions also becomes finite and, choosing $L=1$ , we obtain the range $i=0,\ldots ,2^{j}-1$ for $j=0,\ldots ,J-1$ . Figure 2 illustrates for an orthogonal spline wavelet the discrete scale-space representation for three different scales ( $j=6,7,8$ ) and positions. There exists a fast wavelet transform algorithm which computes the orthogonal wavelet coefficients in $O(N)$ operations, therefore it is even faster than the fast Fourier transform whose operation count is $O(N\log _{2}N)$ (Mallat Reference Mallat1998).

Figure 2.

Space-scale representation of an orthogonal spline wavelet at three different scales and positions, i.e. ${\it\psi}_{6,6},{\it\psi}_{7,32},{\it\psi}_{8,108}$ . The modulus of the Fourier transform of three corresponding wavelets is shown in the inset (top, left).

As an example we show in figure 3 the orthogonal wavelet coefficients of an academic function presenting discontinuities. We observe that wavelet coefficients at small scales only have significant values in the vicinity of the discontinuities. Hence only a few coefficients are needed to represent the function after discarding the small wavelet coefficients.

Figure 3.

Academic example: function with two discontinuities and one in its derivative (a), corresponding modulus of orthogonal wavelet coefficients (b) in logarithmic scale using periodic spline wavelets of degree five.

Extension to higher dimensions. The orthogonal wavelet representation can be extended to represent functions in higher space dimensions using tensor product constructions, see e.g. Daubechies (Reference Daubechies1992), Mallat (Reference Mallat1998), Schneider & Farge (Reference Schneider, Farge, Françoise, Naber and Tsun2006). Figure 4 shows two-dimensional orthogonal wavelets constructed by tensor products.

Figure 4.

Two-dimensional orthogonal wavelets. Scaling function (a) and the three associated directional wavelets in the horizontal (b), vertical (c) and diagonal (d) direction.

The wavelet transform can also be generalized to treat vector-valued functions (e.g. velocity or magnetic fields) in $d$ space dimensions by decomposing each vector component into an orthogonal wavelet series. In the following we consider a vector field $\boldsymbol{v}=(v^{(1)},v^{(2)},v^{(3)})$ for $d=3$ sampled at resolution $N=2^{3J}$ with periodic boundary conditions. Its orthogonal wavelet series reads

(2.10) $$\begin{eqnarray}\displaystyle \boldsymbol{v}(\boldsymbol{x})=\mathop{\sum }_{j=0}^{J-1}\mathop{\sum }_{{\it\mu}=1}^{7}\mathop{\sum }_{i_{1},i_{2},i_{3}=0}^{2^{j}-1}\widetilde{\boldsymbol{v}}_{j,{\it\mu},\boldsymbol{i}}{\it\psi}_{j,{\it\mu},\boldsymbol{i}}(\boldsymbol{x}), & & \displaystyle\end{eqnarray}$$

using three-dimensional orthogonal wavelets ${\it\psi}_{j,{\it\mu},\boldsymbol{i}}(\boldsymbol{x})$ . The basis functions are constructed by tensor products of a set of one-dimensional wavelets and scaling functions (Daubechies Reference Daubechies1992; Mallat Reference Mallat1998) which have been periodized since the boundary conditions considered here are periodic. The scale index $j$ varies from $0$ to $J-1$ , the spatial index $\boldsymbol{i}=(i_{1},i_{2},i_{3})$ has $2^{3j}$ values for each scale $2^{-j}$ and each angle indexed by ${\it\mu}=1,\ldots ,7$ . The three Cartesian directions $\boldsymbol{x}=x^{(1)},x^{(2)},x^{(3)}$ correspond to ${\it\mu}=1,2,3$ , while ${\it\mu}=4,5,6,7$ denote the remaining diagonal directions. The wavelet coefficients measure the fluctuations of $\boldsymbol{v}$ at scale $2^{-j}$ and around position $2^{-j}\boldsymbol{i}$ for each of the seven possible directions ${\it\mu}$ . The contribution of the vector field $\boldsymbol{v}$ at scale $2^{-j}$ and direction ${\it\mu}$ can be reconstructed by summation of $\widetilde{\boldsymbol{v}}_{j,{\it\mu},\boldsymbol{i}}{\it\psi}_{j,{\it\mu},\boldsymbol{i}}(\boldsymbol{x})$ over all positions $\boldsymbol{i}$ :

(2.11) $$\begin{eqnarray}\boldsymbol{v}_{j,{\it\mu}}(\boldsymbol{x})=\mathop{\sum }_{i_{1},i_{2},i_{3}=0}^{2^{j}-1}\widetilde{\boldsymbol{v}}_{j,{\it\mu},\boldsymbol{i}}{\it\psi}_{j,{\it\mu},\boldsymbol{i}}(\boldsymbol{x}).\end{eqnarray}$$

The contribution of $\boldsymbol{v}$ at scale $2^{-j}$ is obtained by

(2.12) $$\begin{eqnarray}\boldsymbol{v}_{j}(\boldsymbol{x})=\mathop{\sum }_{{\it\mu}=1}^{7}\boldsymbol{v}_{j,{\it\mu}}(\boldsymbol{x}).\end{eqnarray}$$

For more details on wavelets, we refer the reader to several review articles, e.g. Farge (Reference Farge1992), Farge & Schneider (Reference Farge, Schneider, Françoise, Naber and Tsun2006), Schneider & Farge (Reference Schneider, Farge, Françoise, Naber and Tsun2006), Schneider & Vasilyev (Reference Schneider and Vasilyev2010) and textbooks, e.g. Daubechies (Reference Daubechies1992), Mallat (Reference Mallat1998).

2.2.1. Scale-dependent moments

To study the scale-dependent directional statistics we consider the component $v^{\ell }$ with $\ell =1,2,3$ of a generic vector field $\boldsymbol{v}$ . First we define the $q$ th-order moment of the scale-dependent vector $\boldsymbol{v}_{j}(\boldsymbol{x})=(v_{j}^{(1)},v_{j}^{(2)},v_{j}^{(3)})$ , which is here either the vector field at scale $2^{-j}$ and direction ${\it\mu}$ , $v_{j,{\it\mu}}^{(\ell )}$ , or the vector field at scale $2^{-j}$ , $v_{j}^{(\ell )}$ ,

(2.13) $$\begin{eqnarray}M_{q}[v_{j}^{(\ell )}]=\langle (v_{j}^{(\ell )})^{q}\rangle .\end{eqnarray}$$

By construction the mean value satisfies $\langle v_{j}^{(\ell )}\rangle =0$ and hence the moments are automatically centred. These scale-dependent moments are related to the $q$ th-order structure functions, as shown, e.g. in Schneider et al. (Reference Schneider, Farge, Kevlahan, Tongring and Penner2004). In the following, we consider the second-order moment $M_{2}[v_{j}^{(\ell )}]$ , which is a scale-dependent quadratic mean intensity of $v_{j}^{(\ell )}$ , and the fourth-order moment $M_{4}[v_{j}^{(\ell )}]$ , which contains the scale-dependent spatial fluctuations. Both moments are related via the flatness factor.

In anisotropic turbulence, typically a preferred direction can be defined, e.g. for low magnetic Reynolds number turbulence, or rotating turbulence. These flows have statistical symmetries, given here with respect to the $x_{3}$ -axis. For the remaining perpendicular components, $\ell =1,2$ , the average of these two components is taken as, $M_{q}[v_{j}^{\bot }]=\{M_{q}[v_{j}^{(1)}]+M_{q}[v_{j}^{(2)}]\}/2$ , and the superscript $\bot$ represents the perpendicular contribution. The parallel contribution $v_{j}^{(3)}$ is denoted by $v_{j}^{\Vert }$ .

The wavelet energy spectrum for $v_{j}^{(\ell )}$ is obtained using $M_{2}[v_{j}^{(\ell )}]$ and (2.6),

(2.14) $$\begin{eqnarray}E[v_{j}^{(\ell )}]=\frac{1}{2{\rm\Delta}k_{j}}M_{2}[v_{j}^{(\ell )}],\end{eqnarray}$$

where ${\rm\Delta}k_{j}=(k_{j+1}-k_{j})\ln 2$ (Meneveau Reference Meneveau1991; Addison Reference Addison2002). It is thus directly related to the Fourier energy spectrum and yields a smoothed version (Meneveau Reference Meneveau1991; Farge Reference Farge1992). The orthogonality of the wavelets with respect to scale and direction guarantees that the total energy is obtained by direct summation, $E=\sum _{\ell ,j}E[v_{j}^{(\ell )}]=\sum _{\ell ,j,{\it\mu}}E[v_{j,{\it\mu}}^{(\ell )}]$ .

The standard deviation of the energy spectrum at a given wavenumber $k_{j}$ quantifies the spatial variability

(2.15) $$\begin{eqnarray}{\it\sigma}[v_{j}^{(\ell )}]=\frac{1}{2{\rm\Delta}k_{j}}\sqrt{M_{4}[v_{j}^{(\ell )}]-(M_{2}[v_{j}^{(\ell )}])^{2}}.\end{eqnarray}$$

The ratio of the fourth- and second-order moments defines the scale-dependent flatness factor,

(2.16) $$\begin{eqnarray}F[v_{j}^{(\ell )}]=\frac{M_{4}[v_{j}^{(\ell )}]}{(M_{2}[v_{j}^{(\ell )}])^{2}},\end{eqnarray}$$

which quantifies the flow intermittency at scale $2^{-j}$ .

The scale-dependent flatness is related to the energy spectrum (2.14) and the standard deviation (2.15),

(2.17) $$\begin{eqnarray}F[v_{j}^{(\ell )}]=\left(\frac{{\it\sigma}[v_{j}^{(\ell )}]}{E[v_{j}^{(\ell )}]}\right)^{2}+1,\end{eqnarray}$$

as shown in Bos, Liechtenstein & Schneider (Reference Bos, Liechtenstein and Schneider2007). This relation illustrates that the spatial variability of the energy spectrum is directly reflected by the scale-dependent flatness.

2.2.2. Scale-dependent directional statistics

To quantify scale-dependent spatial flow anisotropy and anisotropic flow intermittency we introduce wavelet-based measures. Both component-wise anisotropy and directional anisotropy of the flow are considered in the following. For the scale-dependent mean energy, $E[v_{j}^{(\ell )}]$ , the anisotropy measure can be defined similarly to the classical Fourier representation. Analogously, this can be extended for its spatial fluctuations, ${\it\sigma}[v_{j}^{(\ell )}]$ . Using the relation between the scale-dependent flatness with the energy spectrum and its spatial fluctuations, (2.17), various measures of anisotropic flow intermittency can be defined.

Component-wise anisotropy. The scale-dependent component-wise anisotropy is defined by the ratio of perpendicular to parallel energy, and its fluctuation, at a given scale $2^{-j}$ , respectively,

(2.18a,b ) $$\begin{eqnarray}c_{E}(k_{j})=\frac{E[v_{j}^{\bot }]}{E[v_{j}^{\Vert }]},\quad c_{{\it\sigma}}(k_{j})=\frac{{\it\sigma}[v_{j}^{\bot }]}{{\it\sigma}[v_{j}^{\Vert }]}.\end{eqnarray}$$

The scale-dependent mean energy, $c_{E}(k_{j})$ is a smoothed version of the Fourier counterpart $c(k)$ . The component-wise anisotropy of the spatial fluctuations is quantified by $c_{{\it\sigma}}(k_{j})$ . These measures are directly related to the component-wise flatness factors of $v_{j}^{(\ell )}$ , i.e. $F[v_{j}^{\bot }]$ and $F[v_{j}^{\Vert }]$ , as shown in Okamoto et al. (Reference Okamoto, Yoshimatsu, Schneider and Farge2014). Combining (2.17) and (2.18) results in

(2.19) $$\begin{eqnarray}{\it\Lambda}_{j}^{C}\equiv \left\{\frac{c_{{\it\sigma}}(k_{j})}{c_{E}(k_{j})}\right\}^{2}=\frac{F[v_{j}^{\bot }]-1}{F[v_{j}^{\Vert }]-1},\end{eqnarray}$$

which yields a scale-dependent measure of component-wise anisotropic intermittency.

Directional anisotropy. Scale-dependent measures for directional anisotropy can be defined using ratios of perpendicular to parallel energy and fluctuations in longitudinal or transverse directions,

(2.20a,b ) $$\begin{eqnarray}\displaystyle d_{E}^{L}(k_{j})=\frac{E[v_{j,L}^{\bot }]}{E[v_{j,L}^{\Vert }]},\quad d_{{\it\sigma}}^{L}(k_{j})=\frac{{\it\sigma}[v_{j,L}^{\bot }]}{{\it\sigma}[v_{j,L}^{\Vert }]}, & & \displaystyle\end{eqnarray}$$

(2.21a,b ) $$\begin{eqnarray}\displaystyle d_{E}^{T}(k_{j})=\frac{E[v_{j,3}^{\bot }]}{E[v_{j,T}^{\bot }]},\quad d_{{\it\sigma}}^{T}(k_{j})=\frac{{\it\sigma}[v_{j,3}^{\bot }]}{{\it\sigma}[v_{j,T}^{\bot }]}. & & \displaystyle\end{eqnarray}$$

The longitudinal direction is denoted by the index $L$ , i.e. $L={\it\mu}=\ell$ . The subscript ${\it\mu}=3$ corresponds to a transverse direction of the perpendicular components, while $T$ represents the other transverse direction of the perpendicular components, i.e. $T={\it\mu}=1$ for $v_{j,{\it\mu}}^{(2)}$ or $T={\it\mu}=2$ for $v_{j,{\it\mu}}^{(1)}$ .

Three principal directions, i.e. ${\it\mu}=1,2$ and $3$ , out of the seven possible directions have been selected for the directional statistics.

The measures $d_{E}^{L}(k_{j})$ and $d_{E}^{T}(k_{j})$ are smoothed versions of the Fourier representation $2e^{(3)}(k_{3})/\{e^{(1)}(k_{1})+e^{(2)}(k_{2})\}$ and $\{e^{(1)}(k_{3})+e^{(2)}(k_{3})\}/\{e^{(1)}(k_{2})+e^{(2)}(k_{1})\}$ , respectively, following the interpretation of the directional statistics proposed in Bos et al. (Reference Bos, Liechtenstein and Schneider2007). Furthermore, these quantities can be related to second-order structure functions defined in physical space, and we have:

(2.22a,b ) $$\begin{eqnarray}\frac{2D^{(3)}(r\hat{\boldsymbol{l}}_{3})}{\{D^{(1)}(r\hat{\boldsymbol{l}}_{1})+D^{(2)}(r\hat{\boldsymbol{l}}_{2})\}}\quad \text{and}\quad \frac{\{D^{(1)}(r\hat{\boldsymbol{l}}_{3})+D^{(2)}(r\hat{\boldsymbol{l}}_{3})\}}{\{D^{(1)}(r\hat{\boldsymbol{l}}_{2})+D^{(2)}(r\hat{\boldsymbol{l}}_{1})\}}.\end{eqnarray}$$

Structure functions are defined as the spatial average of velocity increments, $D^{(\ell )}(\boldsymbol{r})=\langle \{v^{(\ell )}(\boldsymbol{x}+\boldsymbol{r})-v^{(\ell )}(\boldsymbol{x})\}^{2}\rangle$ . Here $v^{(\ell )}$ consists of contributions of $v^{(\ell )}$ to scales larger than $2^{-j}$ . which are obtained by low pass filtering using the three-dimensional scaling function at scale $2^{-j}$ . The unit vector of the Cartesian direction $x_{\ell }$ is denoted by $\hat{\boldsymbol{l}}_{\ell }$ .

Combining (2.17) and (2.20)–(2.21), yields directional anisotropy measures (Okamoto et al. Reference Okamoto, Yoshimatsu, Schneider and Farge2014):

(2.23) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\Lambda}_{j}^{L}\equiv \left\{\frac{d_{{\it\sigma}}^{L}(k_{j})}{d_{E}^{L}(k_{j})}\right\}^{2}=\frac{F[v_{j,L}^{\bot }]-1}{F[v_{j,L}^{\Vert }]-1}, & \displaystyle\end{eqnarray}$$

(2.24) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\Lambda}_{j}^{T}\equiv \left\{\frac{d_{{\it\sigma}}^{T}(k_{j})}{d_{E}^{T}(k_{j})}\right\}^{2}=\frac{F[v_{j,3}^{\bot }]-1}{F[v_{j,T}^{\bot }]-1}, & \displaystyle\end{eqnarray}$$

$\boldsymbol{B}_{0}$

2.2.3. Scale-dependent topological statistics

Considering the velocity field $\boldsymbol{u}$ and the corresponding vorticity ${\bf\omega}=\boldsymbol{{\rm\nabla}}\times \boldsymbol{u}$ , the kinetic helicity, $H(\boldsymbol{x})=\boldsymbol{u}\boldsymbol{\cdot }{\bf\omega}$ , can be defined. The helicity yields a measure of the geometrical statistics of turbulence. Integrating the helicity over space one obtains the mean helicity $\overline{H}=\langle \boldsymbol{u}\boldsymbol{\cdot }{\bf\omega}\rangle$ . The scale-dependent helicity $H_{j}$ was introduced in Yoshimatsu et al. (Reference Yoshimatsu, Okamoto, Schneider, Kaneda and Farge2009a ) and is defined by

(2.25) $$\begin{eqnarray}H_{j}(\boldsymbol{x})=\boldsymbol{u}_{j}\boldsymbol{\cdot }{\bf\omega}_{j}.\end{eqnarray}$$

It preserves Galilean invariance, though the kinetic helicity itself does not. The corresponding mean helicity is obtained by summing $H_{j}$ over scale, $\overline{H}=\sum _{j}\langle H_{j}\rangle$ due to the orthogonality of the wavelet decomposition.

The relative helicity

(2.26) $$\begin{eqnarray}h(\boldsymbol{x})=\frac{H}{|\boldsymbol{u}|\,|{\bf\omega}|}\end{eqnarray}$$

defines the cosine of the angle between the velocity and the vorticity at each spatial position. The range of $h$ lies between $-1$ and $+1$ . The scale-dependent relative helicity can be defined correspondingly

(2.27) $$\begin{eqnarray}h_{j}(\boldsymbol{x})=\frac{H_{j}}{|\boldsymbol{u}_{j}|\,|{\bf\omega}_{j}|}.\end{eqnarray}$$

The Euler equations of hydrodynamics conserve the mean kinetic helicity, while in ideal MHD turbulence, the mean cross-helicity $\overline{H}^{C}=\langle \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{b}\rangle$ and the mean magnetic helicity $\overline{H}^{M}=\langle \boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{b}\rangle$ are conserved quantities. Here $\boldsymbol{a}$ is the vector potential of the magnetic field $\boldsymbol{b}$ . The scale-dependent versions of the relative cross and magnetic helicities have been introduced in Yoshimatsu et al. (Reference Yoshimatsu, Schneider, Okamoto, Kawahara and Farge2011) and are defined respectively by

(2.28) $$\begin{eqnarray}h_{j}^{C}(\boldsymbol{x})=\frac{H_{j}^{C}}{|\boldsymbol{u}_{j}||\boldsymbol{b}_{j}|},\end{eqnarray}$$

with $H^{C}(\boldsymbol{x})=\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{b}$ and

(2.29) $$\begin{eqnarray}h_{j}^{M}(\boldsymbol{x})=\frac{H_{j}^{M}}{|\boldsymbol{a}_{j}||\boldsymbol{b}_{j}|},\end{eqnarray}$$

with $H^{M}(\boldsymbol{x})=\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{b}$ . These quantities define the cosine of the angle between the two vector fields.

3.1.1. Principle

We propose a wavelet-based method to extract coherent structures that emerge out of turbulent flows, both in fluids (e.g. vortices, shock waves in compressible fluids,  $\ldots$ ) and in plasmas (e.g. bursts, blobs,  $\ldots$ ). The goal is to study their role regarding the transport and mixing properties of flows in the turbulent regime.

For this, we use the wavelet representation that keeps track of both time and scale, instead of the Fourier representation that keeps track of frequency only. Since there is not yet an universal definition of the coherent structures encountered in turbulent flows, we use an apophatic method (introduced in Hinduist theology several thousands years ago) where one does not try to define what an entity (e.g. a phenomenon, a noumenon,  $\ldots$ ) is but rather what it is not. We thus agree on the minimal and hopefully consensual statement: ‘coherent structures are not noise’, and propose to define them as: ‘coherent structures are what remains after denoising’.

The mathematical definition of noise states that a signal is noise if it cannot be compressed in any functional basis. As a result, the shortest description of a noise is itself. Note that in most of the cases, the experimental noise generated by a measurement device does not fit the definition of mathematical noise since it can be compressed in at least one functional basis (e.g. parasite frequencies can be removed in the Fourier basis).

This new way of defining coherent structures allows us to process signals and fields and also their cuts or projections (e.g. a probe located at one point provides a one-dimensional cut of a four-dimensional space-time field). Indeed, the algorithms commonly used to extract coherent structures cannot work for cuts or projections because they require a template of the structures to extract (one would need to take into account of how the probe sees all possible translations and distortions of the coherent structures). The strength of our algorithm is that it treats fields and projections in the same way.

Since we assume that coherent structures are what remains after denoising, we need a model, not for the structures themselves, but for the noise. Applying ‘Ockham’s Razor principle’ (or the ‘law of parsimony’), we choose as a first guess the simplest possible model: we suppose the noise to be additive, Gaussian and white (i.e. uncorrelated). We then project the turbulent signal (in one-dimension), or turbulent field (in higher dimensions), into wavelet space and retain only the coefficients having a modulus larger than a given threshold. As a threshold value we follow Donoho and Johnstone’s proposition of a threshold value that depends on the variance of the Gaussian noise we want to remove and on the chosen sampling rate (Donoho & Johnstone Reference Donoho and Johnstone1994). Since the noise variance is not known a priori for turbulent signals (the noise being produced by their intrinsic nonlinear dynamics), we designed a recursive method Azzalini et al. (Reference Azzalini, Farge and Schneider2005) to estimate it from the variance of the weakest wavelet coefficients, i.e. those whose modulus is below the threshold value. After applying our algorithm, we obtain two orthogonal fields: the coherent field retaining all coherent structures and the incoherent field corresponding to the noise. We then check a posteriori that the latter is indeed noise-like (i.e. spread all over physical space), Gaussian and uncorrelated (i.e. also spread all over Fourier space), and thus confirm the hypotheses we have a priori chosen for the noise.

3.1.2. Wavelet denoising

We consider a signal $s(t)$ sampled on $N=2^{J}$ points that we want to denoise, assuming the noise to be additive, Gaussian and white. We first project $s(t)$ onto an orthogonal wavelet basis and then filter out some of the wavelet coefficients thus obtained, $\widetilde{s}_{ij}$ . We retain only the wavelet coefficients whose modulus is larger than a threshold value. The main difficulty is to estimate it a priori and we encounter two possible cases:

1. (i) If we a priori know the noise’s variance ${\it\sigma}^{2}$ , the optimal threshold value is given by Donoho and Johnstone’s formula (Donoho & Johnstone Reference Donoho and Johnstone1994)

   (3.1) $$\begin{eqnarray}{\it\epsilon}=(2{\it\sigma}^{2}\ln N)^{1/2}.\end{eqnarray}$$
   In 1994 they proved (Donoho & Johnstone Reference Donoho and Johnstone1994) that such a wavelet thresholding method is optimal to denoise signals in the presence of additive Gaussian white noise because it minimizes the maximal $L^{2}$ -error (between the denoised signal and the noise-free signal) for functions whose regularity is inhomogeneous, such as bursty or intermittent turbulent signals.

2. (ii) If we do not a priori know the variance of the noise, which is the most usual case, one should use the wavelet-based recursive algorithm we proposed in Farge et al. (Reference Farge, Schneider and Kevlahan1999), Azzalini et al. (Reference Azzalini, Farge and Schneider2005). This algorithm first estimates the variance of the noise by considering the variance of the noisy signal ${\it\sigma}_{0}^{2}$ and computes the corresponding threshold

   (3.2) $$\begin{eqnarray}{\it\epsilon}_{0}=(2{\it\sigma}_{0}^{2}\ln N)^{1/2}.\end{eqnarray}$$
   The algorithm splits the wavelet coefficients into two classes: the weak coefficients whose moduli are smaller than the threshold, and the remaining strong coefficients. It then computes the variance of the weak coefficients ${\it\sigma}_{n}$ to obtain a better estimation of the variance of the noise (estimated from the wavelet coefficients using Parseval’s theorem)
   (3.3) $$\begin{eqnarray}{\it\sigma}_{n}^{2}=\frac{1}{N}\mathop{\sum }_{(j,i)\in {\mathcal{I}}^{J},|\widetilde{s}_{ji}|<{\it\epsilon}_{n}}|\widetilde{s}_{ji}|^{2},\end{eqnarray}$$
   where ${\mathcal{I}}^{J}=\{0\leqslant j<J,i=0,\ldots ,2^{j}-1\}$ is the index set of the wavelet coefficients. The algorithm then replaces ${\it\epsilon}_{0}$ by ${\it\epsilon}_{n}=(2{\it\sigma}_{n}^{2}\ln N)^{1/2}$ , which yields a better estimate of the threshold. This procedure is iterated until it reaches the optimal threshold value, when ${\it\epsilon}_{n+1}\approx {\it\epsilon}_{n}$ .

   In Azzalini et al. (Reference Azzalini, Farge and Schneider2005) we proved that this algorithm converges for signals having a sufficiently sparse representation in wavelet space, such as the intermittent signals encountered in turbulence. We also showed that the larger the signal to noise ratio (SNR) is, the faster the convergence. Hence, if the signal $s(t)$ is only noise, it converges in one iteration and retains ${\it\epsilon}_{0}$ as the optimal threshold.

Using the optimal threshold, we then separate the wavelet coefficients $\widetilde{s}_{ij}$ into two contributions: the coherent coefficients $\widetilde{s}_{ij}^{C}$ whose moduli are larger than ${\it\epsilon}$ and the remaining incoherent coefficients $\widetilde{s}_{ij}^{I}$ . Finally, the coherent component $s^{C}(t)$ is reconstructed in physical space using the inverse wavelet transform, while the incoherent component is obtained as $s^{I}(t)=s(t)-s^{C}(t)$ .

3.1.3. Extraction algorithm for one-dimensional signals

We detail the iterative extraction algorithm for the one-dimensional case and quote it from Azzalini et al. (Reference Azzalini, Farge and Schneider2005):

Initialization

1. (i) given the signal $s(t)$ of duration $T$ , sampled on an equidistant grid $t_{i}=\text{i}T/N$ for $i=0,N-1$ , with $N=2^{J}$ ;

2. (ii) set $n=0$ and perform a wavelet decomposition, i.e. apply the FWT (Mallat Reference Mallat1998) to $s$ to obtain the wavelet coefficients $\widetilde{s}_{ji}$ for $(j,i)\in {\mathcal{I}}^{J}$ ;

3. (iii) compute the variance ${\it\sigma}_{0}^{2}$ of $s$ as a rough estimate of the variance of the incoherent signal $s^{I}$ and compute the corresponding threshold ${\it\epsilon}_{0}=(2\ln N{\it\sigma}_{0}^{2})^{1/2}$ , where ${\it\sigma}_{0}^{2}=(1/N)\sum _{(j,i)\in I^{J}}|\widetilde{s}_{ji}|^{2}$ ;

4. (iv) set the number of coefficients considered as noise to $N_{I}=N$ , i.e. to the total number of wavelet coefficients.

Main loop

Repeat

1. (i) set $N_{I}^{old}=N_{I}$ and count the number of wavelet coefficients smaller than ${\it\epsilon}_{n}$ , which yields a new value for $N_{I}$ ;

2. (ii) compute the new variance ${\it\sigma}_{n+1}^{2}$ from the wavelet coefficiens smaller than ${\it\epsilon}_{n}$ , i.e. ${\it\sigma}_{n+1}^{2}=(1/N)\sum _{(j,i)\in {\mathcal{I}}^{J}}|\widetilde{s}_{ji}^{I}|^{2}$ , where

   (3.4) $$\begin{eqnarray}\displaystyle \widetilde{s}_{ji}^{I}=\left\{\begin{array}{@{}ll@{}}\widetilde{s}_{ji} & \text{for }|\widetilde{s}_{ji}|\leqslant {\it\epsilon}_{n}\\ 0 & \text{else}\end{array}\right. & & \displaystyle\end{eqnarray}$$
   and the new threshold ${\it\epsilon}_{n+1}=(2\ln N{\it\sigma}_{n+1}^{2})^{1/2}$ ;

3. (iii) set $n=n+1$

until ( $N_{I}=N_{I}^{old}$ ).

Final step

1. (i) reconstruct the coherent signal $s^{C}$ from the coefficients $\widetilde{s}_{ji}^{C}$ using the inverse FWT, where

   (3.5) $$\begin{eqnarray}\displaystyle \widetilde{s}_{ji}^{C}=\left\{\begin{array}{@{}ll@{}}\widetilde{s}_{ji} & \text{for }|\widetilde{s}_{ji}|>{\it\epsilon}_{n}\\ 0 & \text{else}\end{array}\right. & & \displaystyle\end{eqnarray}$$

2. (ii) finally, compute pointwise the incoherent signal $s^{I}(t_{i})=s(t_{i})-s^{C}(t_{i})$ for $i=0,\ldots ,N-1$ .

End

Note that the signal is split into $s(t)=s^{C}(t)+s^{I}(t)$ and its energy into ${\it\sigma}^{2}={\it\sigma}_{C}^{2}+{\it\sigma}_{I}^{2}$ , since the coherent and incoherent components are orthogonal, i.e. $\langle s^{C},s^{I}\rangle =0$ .

We use the FWT (Mallat Reference Mallat1998) that is computed with $(2MN)$ multiplications, $M$ being the length of the discrete filter defining the orthogonal wavelet used. Remark: for all applications presented in this paper, we use Coiflet 12 wavelets (Daubechies Reference Daubechies1992), unless otherwise stated. As long as the filter length $M<(\log _{2}N)/2$ , the FWT is faster than the FFT computed with $N\log _{2}N$ operations. Consequently, the extraction algorithm requires $(2nMN)$ operations, $n$ being the number of iterations, which is small, typically less than $\log _{2}N$ .

This algorithm defines a sequence of estimated thresholds $({\it\epsilon}_{n})_{n\in \mathbb{N}}$ and the corresponding sequence of estimated variances $({\it\sigma}_{n}^{2})_{n\in \mathbb{N}}$ . In Azzalini et al. (Reference Azzalini, Farge and Schneider2005) we proved that this sequence converges after a finite number of iterations by applying a fixed point type argument to the iteration function

(3.6) $$\begin{eqnarray}{\mathcal{F}}_{s,N}({\it\epsilon}_{n+1})=\left(\frac{2\ln N}{N}\mathop{\sum }_{(j,i)\in {\mathcal{I}}^{J}}|\widetilde{s}_{ji}^{I}({\it\epsilon}_{n})|^{2}\right)^{1/2}.\end{eqnarray}$$

The algorithm stops after $n$ iterations, when ${\mathcal{F}}_{s,N}({\it\epsilon}_{n})={\it\epsilon}_{n+1}$ , since the number of samples $N$ is finite. In Azzalini et al. (Reference Azzalini, Farge and Schneider2005) we also proved that the convergence rate depends on the SNR ( $\text{SNR}=10\log _{10}({\it\sigma}^{2}/{\it\sigma}_{I}^{2})$ ), since the smaller the SNR, i.e. the stronger the noise, the faster the convergence, moreover, if the algorithm is applied to a Gaussian white noise, it converges in one iteration only. If it is applied to a signal without noise, the signal is fully preserved. In Azzalini et al. (Reference Azzalini, Farge and Schneider2005) we have also proven the algorithm’s idempotence, i.e. if it is applied several times the noise is eliminated the first time and the coherent signal will remain the same if the algorithm is reapplied several times. This would be the case for a Gaussian filter which, in contrast, is not idempotent.

Figure 10.

Denoising of a piecewise regular signal using iterative wavelet thresholding. (a) Original academic signal $a(t)$ . (b) Noisy signal $s(t)$ with a $\text{SNR}=11\;\text{dB}$ . (c) Denoised signal $s^{C}(t)$ with a $\text{SNR}=28\;\text{dB}$ .

3.1.4. Application to an academic test signal

To illustrate the performance of the iterative algorithm we consider a one-dimensional noisy test signal $s(t)$ sampled on $N=2^{13}=8192$ points (figure 10 b). It is made by adding a Gaussian white noise $w(t)$ , of mean zero and variance ${\it\sigma}_{w}^{2}=25$ , to a piecewise regular academic signal $a(t)$ presenting several discontinuities, in the function or in its derivatives (figure 10 a). The SNR is $\text{SNR}=10\log _{10}({\it\sigma}_{a}^{2}/{\it\sigma}_{w}^{2})=11\;\text{dB}$ . After applying the extraction algorithm we estimate the noise variance to be $25.6$ and we obtain a coherent signal $s^{C}(t)$ very close to the original academic signal $a(t)$ (figure 10 c). The incoherent part $s^{I}(t)$ is homogeneous and noise-like with flatness $3.03$ , which corresponds to quasi-Gaussianity. In figure 10(c) we observe that the coherent signal retains all discontinuities and peaks present in the academic signal $a(t)$ , which is an advantage with respect to standard denoising techniques, e.g. low pass Fourier filtering, which smooth them. In the vicinity of the discontinuities we observe slight overshoots, which are more local than the classical Gibbs phenomena and could for example be removed using the translation invariant wavelet transform (Mallat Reference Mallat1998).

3.1.5. Extension of the algorithm to higher dimensional scalar- and vector-valued fields

The extraction algorithm was described in § 3.1.3 for one-dimensional scalar-valued signals $s(t)$ varying in time. First, it can be extended to higher dimensional scalar fields $s(\boldsymbol{x})$ varying in space $\boldsymbol{x}\in \mathbb{R}^{d}$ where $d$ is the space dimension. To this end the extraction algorithm only requires that the one-dimensional wavelets are replaced by their equivalent $d$ -dimensional wavelets using tensor product constructions, see, e.g. (Daubechies Reference Daubechies1992; Mallat Reference Mallat1998; Schneider & Farge Reference Schneider, Farge, Françoise, Naber and Tsun2006).

Second, the extraction algorithm can also be extended to vector-valued fields $\boldsymbol{v}=(v^{(1)},\ldots ,v^{(d)})$ where each component $v^{\ell },\ell =1,\ldots ,d$ is a scalar valued field. The extraction algorithm is then applied to each component of the vector field. For thresholding the wavelet coefficients we consider the vector $\widetilde{\boldsymbol{v}}_{j,{\it\mu},\boldsymbol{i}}$ in (2.10). Assuming statistical isotropy of the noise, the modulus of the wavelet coefficient vector is computed. The coherent contribution is then reconstructed from those coefficients whose modulus is larger than the threshold defined as ${\it\epsilon}=(2/d\;{\it\sigma}^{2}\ln N)^{1/2}$ where $d$ is the dimensionality of the vector field, ${\it\sigma}$ the variance of the noise and $N$ the total number of grid points. The iterative algorithm in § 3.1.3 can then be applied in a straightforward way.

To extract coherent structures out of turbulent flows we consider the vorticity field, which is decomposed in wavelet space. Applying the extraction algorithm then yields two orthogonal components, the coherent and incoherent vorticity fields. Subsequently the corresponding induced velocity fields can be reconstructed by applying the Biot–Savart kernel, which is the inverse curl operator. For MHD turbulence, we consider in addition the current density and we likewise split it into two components, the coherent and incoherent current density fields. Using Biot–Savart’s kernel we reconstruct the coherent and incoherent magnetic fields.

Note that the employed wavelet bases do not a priori constitute divergence-free bases. Thus the resulting coherent and incoherent vector fields are not necessarily divergence free. However, we checked that the departure from incompressibility only occurs in the dissipative range and remains negligible (Yoshimatsu et al. Reference Yoshimatsu, Kondo, Schneider, Okamoto, Hagiwara and Farge2009b ). Another solution would be to use directly div-free wavelets, however these are much more cumbersome to implement (Deriaz, Farge & Schneider Reference Deriaz, Farge and Schneider2010).

Figure 11.

(a) Position of the reciprocating Langmuir probe in the scrape-off layer of the tokamak Tore Supra in Cadarache. (b) Schematic top view of the probe.

3.3.1. Tomographic reconstruction using wavelet-vaguelette decomposition

Cameras installed in tokamaks acquire images which are difficult to interpret, since the three-dimensional structure of the plasma is mapped onto two spatial dimensions and thus flattened in a non-trivial way. This implies that the received flux cannot be directly related to the volumic emissivity of the plasma, which is a major limitation of such optical diagnostics. The reason is that the photons collected by each pixel on the camera sensor have been emitted along a corresponding ray, rather than out of a single point in space. Nevertheless the three-dimensional radiation can be related to the two-dimensional image using tomographic reconstruction, because the dominant structures in tokamak edge turbulence happen to be field-aligned filaments, commonly known as blobs. They have a higher density than their surroundings, and their structure varies more slowly along magnetic field lines than in their orthogonal directions.

Mathematically, the tomographic reconstruction corresponds to an inverse problem which has a formal solution under the assumed symmetry, but is ill-posed in the presence of noise. Taking advantage of the slow variation of the fluctuations along magnetic field lines in tokamaks, this inverse problem can be modelled by a helical Abel transform, which is a Volterra integral operator of the first kind. In Nguyen van yen et al. (Reference Nguyen van yen, Fedorczak, Brochard, Bonhomme, Schneider, Farge and Monier-Garbet2012) we proposed a tomographic inversion technique, based on a wavelet-vaguelette decomposition coupled with wavelet denoising, to extract coherent structures which allows us to detect individual blobs on the projected movie and to analyse their behaviour. The wavelet-vaguelette decomposition (WVD) was introduced by Tchamitchian (Reference Tchamitchian1987) and used by Donoho (Reference Donoho1995) to solve inverse problems in the presence of localized structures. Tomographic inversion using the wavelet-vaguelette decomposition is as an alternative to singular value decomposition (SVD). Both decompositions regularize the problem by damping the modes of the inverse transform to prevent amplification of the noise, i.e. modes below a given threshold are eliminated. For WVD, the nonlinear iterative thresholding procedure (see § 3.1.3) is applied to the vaguelette coefficients. Here Coiflets with two vanishing moments are used (Daubechies Reference Daubechies1992). However, in contrast to SVD, WVD takes in addition advantage of the spatial localization of coherent structures present in the plasma.

The technicalities of WVD are described in detail in Nguyen van yen et al. (Reference Nguyen van yen, Fedorczak, Brochard, Bonhomme, Schneider, Farge and Monier-Garbet2012), in the following we explain only the principle. The helical Abel transform related the plasma light emissivity $S$ (a scalar-valued field) to the integral of the volume emissivity received by the camera $I=KS$ , where $K$ is a compact continuous operator. The reconstruction of the plasma light emissivity $S$ from $I$ is an inverse problem which becomes very difficult when $S$ is corrupted by noise, since computing $K^{-1}$ is an ill-posed problem which amplifies the noise. The vaguelettes are operator adapted wavelets and a biorthogonal set of basis functions is obtained from the wavelet bases ${\it\psi}_{{\it\lambda}}$ by computing $K{\it\psi}_{{\it\lambda}}$ and ${K^{\star }}^{-1}{\it\psi}_{{\it\lambda}}$ , where ${K^{\star }}^{-1}$ denotes the adjoint inverse operator (Tchamitchian Reference Tchamitchian1987). Note that vaguelettes inherit the localization features of wavelets but may lose the translation and scale invariance, and thus the FWT cannot be applied anymore.

3.3.2. Application to an academic example

To illustrate the method, we first consider an academic test case with an given emissivity map $S$ , having a uniform radiating shell at constant value one and zero elsewhere. A two-dimensional cut in the poloidal plane is shown in figure 16(a). Applying the helical Abel transform we generate the corresponding synthetic image $I=KS$ (figure 16 b). Then we add a Gaussian white noise with standard deviation $0.5$ , which yields the synthetic noisy image (figure 16 c).

Figure 16.

Denoising WVD academic test case with a uniform radiating shell. (a) Source emission intensity $S$ in the poloidal plane. (b) Corresponding noiseless image $I=KS$ in the image plane. (c) Noisy image obtained by adding Gaussian white noise with variance 0.5 (from Nguyen van yen et al. Reference Nguyen van yen, Fedorczak, Brochard, Bonhomme, Schneider, Farge and Monier-Garbet2012).

Applying the WVD reconstruction to the synthetic noisy image (figure 16 c) gives a denoised emissivity map, a poloidal cut is shown in figure 17(a). We observe that the main features are preserved, i.e. the constant emissivity shell is well recovered, besides some spurious oscillations close to discontinuities. The corresponding denoised image $I_{d}=KS_{d}$ (figure 17 b) illustrates that the noise has been successfully removed. A comparison with the standard SVD technique in Nguyen van yen et al. (Reference Nguyen van yen, Fedorczak, Brochard, Bonhomme, Schneider, Farge and Monier-Garbet2012) (not shown here) illustrates the superiority of the wavelet-vaguelette technique.

Figure 17.

Denoising WVD academic test case. WVD inversion results. (a) Reconstructed poloidal emissivity map $S_{d}$ . (b) Denoised image $I_{d}=KS_{d}$ (from Nguyen van yen et al. Reference Nguyen van yen, Fedorczak, Brochard, Bonhomme, Schneider, Farge and Monier-Garbet2012).

3.3.3. Application to fast camera data from tokamaks

Now we present an application to an experimental movie acquired during the Tore Supra discharge TS42967, where the plasma was fully detached and stabilized over several seconds using a feedback control. The movie has been obtained using a fast camera recording at 40 kHz. Moreover, the time average of the whole movie was subtracted from each frame, which helps us to decrease the effect of reflection on the chamber wall. The algorithm is then applied directly to the fluctuations in the signal instead of the full signal. The experimental conditions can be found in Nguyen van yen et al. (Reference Nguyen van yen, Fedorczak, Brochard, Bonhomme, Schneider, Farge and Monier-Garbet2012). One frame of the movie is shown in figure 18(a) and used as input for the WVD reconstruction algorithm. The resulting emissivity map in the poloidal plane, in figure 18(b), shows the presence of localized blobs, which propagate counter-clockwise, as observed in the movies, not shown here. Thus their propagation velocity can be determined. The corresponding denoised movie frame $I_{d}$ (figure 18 c) is obtained by applying the operator $K$ to the inverted emissivity map $S_{d}$ . We observe that the noise has been removed and the local features such as blobs and fronts have been extracted.

Figure 18.

WVD-inversion of a snapshot from a movie obtained from Tore Supra, discharge TS42967. (a) Noisy frame used as input for the WVD algorithm. (b) Reconstructed emissivity map obtained as a result of WVD. (c) Denoised frame obtained by applying the operator $K$ to the reconstructed emissivity map (from Nguyen van yen et al. Reference Nguyen van yen, Fedorczak, Brochard, Bonhomme, Schneider, Farge and Monier-Garbet2012).