Method for estimating the frequency-wavenumber
resolved power spectrum density using the
maximum entropy method for limited spatial
pointsメタデータ

言語: English
出版者:
公開日: 2021-07-14
キーワード (Ja):
キーワード (En):作成者: Kobayashi, Tatsuya, Nishizawa, Takashi, Sasaki,
Makoto, YOSHINUMA, Mikiro, IDA, Katsumi
メールアドレス:所属:

URL

http://hdl.handle.net/10655/00012586
This work is licensed under a Creative Commons
Attribution-NonCommercial-ShareAlike 3.0
International License.

1

Method for estimating the frequency-wavenumber
resolved power spectrum density using the maximum
entropy method for limited spatial points
T. Kobayashi1,2 , T. Nishizawa3 , M. Sasaki4 , M. Yoshinuma1 , andK. Ida1
1

National Institute for Fusion Science, National Institutes of Natural Sciences, Toki 509-5292,

Japan
2

The Graduate University for Advanced Studies, SOKENDAI, Toki 509-5292, Japan

3

Max Planck Institute for Plasma Physics, Boltzmannstr. 2, 85748 Garching, Germany

4

Research Institute for Applied Mechanics, Kyushu University, Kasuga 816-8580, Japan

E-mail: kobayashi.tatsuya@nifs.ac.jp

Abstract.

Combination of the Fourier transform and the maximum entropy method for

estimating the frequency-wavenumber resolved power spectrum density is proposed.

Afterillustrating the physical insight of the maximum entropy method by using synthetic test data,capability of the proposed method is tested using numerical simulation data. The method is also
applied to experimental data obtained by the beam emission spectroscopy in the Large Helical
Device. All of those examinations show that the proposed method provides more plausible results
than conventional methods when the available spatial points are limited.

1. Introduction

The frequency-wavenumber resolved power spectrum density is a powerful analysis
tool widely used in the magnetically confined plasma community. Such an analysis

2
is often performed for the poloidal or azimuthal direction in low temperature basic
plasma devices by using the Fourier transform (FT) [1, 2, 3, 4, 5, 6, 7]. This analysis
allows the direct comparison between the experimental observation and the dispersion
relation theoretically predicted. A poloidal or azimuthal probe array that covers thefull circumference [8, 9, 10] plays a crucial role for realizing the spatial FT with theperiodic boundary condition. Availability of such a probe array is an advantage of basic
devices. Even in high temperature torus devices, the frequency-wavenumber resolved
power spectrum density is calculated as well with a limited coverage of the diagnostics
by using the spatial FT [11, 12, 13].
An alternative way to compute the frequency-wavenumber resolved spectrum is theso-called two-point statistical method [14]. Due to its applicability in a wide variety
of situations, the two-point statistical method is frequently used for various occasions
both in linear plasma columns [15, 16, 17] and in torus devices [18, 19, 20]. However,this method is in principle inapplicable for two or more wave components at a frequency
range, for instance, the standing wave [21]. In such a case, at least three probe tips
are required to resolve the amplitude and the wavenumber of the constituents. A longdischarge duration is also necessary for performing ensemble averaging to achieve thestatistical convergence.
In a condition where a multiple channel fluctuation diagnostic that partly covers
a particular spatial range is available, which is often realized in fusion devices throughthe beam emission spectroscopy, the Doppler back scattering, and others, a diﬀerent

3path to reach the frequency-wavenumber resolved power spectrum density is possible byusing the maximum entropy method (MEM) [22, 23, 24, 25, 26, 27]. MEM can estimate
the power spectrum density with a very high resolution even when the input wave istruncated in a part of one cycle [24, 25]. MEM is occasionally applied to fusion plasma
data as well, aiming for improving the resolution in the wavenumber estimation witha limited number of measurement points [11, 28, 29, 30, 31, 32, 33]. In this paper, we
systematically document a method for estimating the frequency-wavenumber resolved
power spectrum density by applying FT to the time domain and MEM to the spatial
domain. In section 2, the analysis principle of MEM is described by contrasting thediﬀerence between FT with a window function and MEM on a set of test data. Here,the physical insight into MEM is focused upon, therefore readers are recommended
to refer to specialized literatures (e.g., [24, 25]) for the comprehensive review ofmathematical backgrounds. Section 3 introduces the two-dimensional FT-MEM method
for estimating the frequency-wavenumber resolved power spectrum density. Applicationof the proposed method is presented in sections 4 and 5, to a numerical fluid turbulencesimulation data and to a multichannel beam emission spectroscopy (BES) data in the
Large Helical Device (LHD), respectively. A summary is given in section 6.

2. Comparison of Fourier transform and maximum entropy method

In this section, FT and MEM are compared in estimating the power spectrum densityof a time-dependent test signal. The power spectrum density S(f ) as a function of the

4frequency f is defined by the Wiener-Khinchin (W-K) theorem,
1 ∫∞S(f ) =C(τ )e−i2πf τ dτ
2π −∞

(1)

where C(τ ) is the autocorrelation function, which depends on the time delay τ , and i is
the imaginary unit. In practice, any recorded experimental data have their finite data
length, therefore the integral in Eq. (1) cannot be rigorously performed. Instead, thepower spectrum density is estimated by giving certain assumptions for data outside therecord length.
The most commonly used method for estimating the power spectrum density is FT
via the fast FT (FFT) algorithm. In using FT for the power spectrum density estimation,
the data with a finite length is assumed to repeat itself outside the record length. This
assumption is reasonable for the data which are actually periodic, for example for the
spatial distribution of a physical quantity in the poloidal or the toroidal direction in
torus plasmas. However, in the time domain, the data is usually not periodic, therefore
discontinuities appear at both ends of the data. In order to avoid this problem, various
window functions, such as the Hanning window or the Tukey window, are often used.
Applying a window function on the data reduced the relative importance of the data
points near both ends. To compensate this eﬀect, an overlapping moving average ofthe power spectrum density is performed for data with a long record length. However,FT with a window function does not necessarily perform well when the number of data
points is limited, which is often the case for spatially sampled data.
A completely diﬀerent methodology called MEM was introduced by Burg [22] in51967 for estimating the power spectrum density particularly from data with a limited
record length.

MEM estimates the power spectrum density by extrapolating the

autocorrelation function to an infinite lag in such a way that the entropy density is
maximized. The entropy density of a fluctuation data in a stationary Gaussian process
is defined as
[

]

1 ∫∞
h∝
log S(f )df =log
C(τ )e−i2πf τ dτ df,
2π −∞
−fN
−fN
∫ fN

∫ fN

(2)

where fN is the Nyquist frequency. On the one hand, maximizing the entropy density
holding consistency of S(f ) in the recorded portion is regarded as a variational problem
that can be directly solved by the method of Lagrange multiplier [24, 25]. On the other
hand, maximizing the entropy density in extrapolating the autocorrelation function
corresponds to maintaining the maximal ambiguity for the unavailable information
but being consistent with the known information [24, 25, 27]. This extrapolation isequivalent to the linear prediction (LP) by autoregressive model fitting for the raw
data, i.e.,M∑(M )aj xk−j ,
x̂k =
j=1
(M )

where x̂k

(3)

is the linearly predicted value of a time series xk , aj is the LP coeﬃcients,and M is the number of data points used for LP. The MEM spectrum estimate is givenby the LP coeﬃcients aj as(M )SMEM (f ) =

|1 −

∆tPM
,
2
j=1 aj exp(−i2πf j∆t)|

∑M −1

(4)

6
where ∆t is the data sampling time and PM is the variance of the estimation error, i.e.,
PM = N −1

N
∑

(M )

|xk − x̂k |2 .

(5)

k=1

The LP coeﬃcients aj are obtained so as to minimize PM , which is equivalent tomaximizing the entropy density.
The windowed FT and MEM are compared using test data. The test data is defined
as a simple cosine wave, cos(2πf t), where f = 10.5 and t continues infinitely. For thedata analysis, 100 points is taken with 0.01 sampling interval and labeled as t = 0
to 0.99. At the sampling, the white Gaussian noise with the standard deviation of
0.2 is superposed. In order to perform ensemble averaging, 100 datasets are created.
One of the test data series is shown in Fig. 1 (a). Results of the windowed FT withthree diﬀerent window functions (rectangular, Hanning, and Tukey) and MEM with thedata points for LP of M = 9 are displayed in Fig. 1 (b). The frequency resolutionof MEM, which can be arbitrarily set, is 0.1 Hz. MEM provides the power spectrum
density that is sharpest and closest to the true spectrum, i.e., the Dirac delta function ofS(f ) = δ(f − 10.5). Meanwhile the power spectrum densities obtained by the windowed
FTs show peaks broadened to some extent. In particular, FT with the rectangular
window gives a very wide tail much above the noise floor at O(∼ 10−5 ) likely due to thediscontinuity of the repeated data as discussed below. The power spectrum density with
the Tukey window appears with two sidebands at both sides of the peak. The resultwith the Hanning window is better than the other two window function cases, but stillshows a substantial peak width.

1

ls

y

E

1

ft

M
tos z
2

7
r

(a)

T
0

5

2
w

H

l
H
1=
SFH

R

(b)

M
H

g

F

Figure 1. (a) Test data and (b) its power spectrum densities obtained by the Fourier transformwith various window functions and the maximum entropy method.

Calculating the autocorrelation function from the obtained power spectrum
densities provides a much clearer insight into the diﬀerence between the windowed FT
and MEM. The autocorrelation function is given byC(τ ) ∼

∫ fN
−fN

S(f )ei2πf τ df,(6)

where the approximately equal sign is due to the integration truncated at the Nyquist
frequency fN . Figures 2 (a-d) correspond to the autocorrelation functions obtained using
Eq. (6) with various values of τ . An envelope modulation with the period of τ = 1 isseen in all the autocorrelation functions for FT. In contrast, MEM provides an almost
complete cosine wave with a constant amplitude of unity. The fluctuation frequency is
consistent with the test data frequency of f = 10.5, as expected of the W-K theorem.
In order to discuss how the autocorrelation function is modulated for the FT cases
and why not for the MEM case, the autocorrelation function is directly calculated from

•

fiT

hTc - 1 1 )
h T N
w

/

1
1

•

# $ ~2 ' +
(

8 q o

(

E

0

ACF from PSD with W-H theorem (solid) &ACF from windowed raw data (dashed)

Raw signal with windows or LP applied(a) FT (Rectangular)

(e)

•

(f)

•
•

(c) FT (Tukey)(g)(d) MEM

(h)

a
Eh

•
(b) FT (Hanning))*

EW ntlrv f
1 / !"h•••hTc pnEWq oTci ntlrv- 0/Ei
hcfiT q o e
h
q o 5
NT
h
i
. 1
eF
h
1 0. .
N ( dN
10 /
NT
WEh
e
0
/
ntlrv
fih
MiT

Figure 2. (a-d) Autocorrelation functions obtained by power spectrum densities in Fig. 1 (b)through the Wiener-Khinchin (W-K) theorem, (e-g) repeated raw data multiplied by windowfunctions (rectangular, Hanning, and Tukey, respectively), and (h) raw data extrapolated by thelinear prediction (LP). Autocorrelation functions directly obtained by data in (e-h) are overlaidin (a-d) by dashed curves.

the test data in the time domain. Figure 2 (e) is the repeated test data with therectangular window function, i.e., no window function. Here, no statistical noise isapplied to the data for simplicity. There are discontinuous points at t = 0 and 1, since
the test cosine wave and the sampled data length are not synchronizing. The period of
discontinuous jumps in the wave phase corresponds to the lifetime of this wave, whichcauses the envelope modulation in the autocorrelation function. Applying the windowfunctions removes discontinuous points as shown in Figs. 2 (f) and (g), but in return
modulates the data amplitude. The autocorrelation functions directly calculated from9the data in Figs. 2 (e-g) are overlaid on Figs. 2 (a-c), which almost completely agree withthose obtained via Eq. (6). The finite lifetime of the wave is reflected to the spectral
peak width shown in Fig. 1 (b), where the half widths at the half maximum of thepeaks are almost identical to the windowed FT results, although the tail shapes arevery diﬀerent. Note that the phase diﬀerence between two curves in Fig. 2 (d) gradually
increases as |τ | increases. This is because the peak frequency of the MEM spectrum has
1 % error.
Input for the direct calculation of the autocorrelation function in the MEM caseis the raw data extrapolated by the LP scheme, which is plotted in Fig. 2 (h). Here,M = 9 is again used as the number of points for LP. By LP, the almost complete cosine
wave is reproduced for t < 0 and t > 1, which therefore maximizes the entropy density
as a result. LP extends the lifetime of the data extremely long, which is reflected in
the sharp spectral peak of MEM as shown in Fig. 1 (b). The autocorrelation function
directly calculated from the extended raw data is overlaid on that obtained via Eq. (6),showing a good agreement here as well.
In the actual application of MEM, the result is sometimes sensitive to the choice ofthe parameter M . Too large value of M can cause an overfitting in the estimated power
spectrum density. There is a scheme to systematically specify the most appropriate
value of M according to the Akaike final prediction error (FPE) criterion [34, 24, 25].
The value of M for the current example is chosen based on this scheme, which hasprovided plausible results. However, the M value provided by the Akaike FPE criterion

10
is not necessarily always most reasonable in practice, and a heuristic determination of M
is occasionally inevitable. In contrast, FT is stable in most of the situations. Therefore,
if there is no specific need in using MEM, e.g., the limitation of the data points, it issafe to use FT.

3. Two-dimensional FT-MEM for the frequency-wavenumber spectral
decomposition

The frequency-wavenumber resolved power spectrum density is calculated by applying
FT and MEM to data in the time-space domain. The analysis procedure is shown in
Fig. 3. First, FT is applied in the time domain to signals at each channel to computethe complex Fourier component as a function of the frequency. The complex MEM[26] is then operated in the space domain for each frequency, providing the complex
frequency-wavenumber resolved FT-MEM component X(f, k). The power spectrum
density is calculated as ⟨X(f, k)X(f, k)∗ ⟩, where ∗ and ⟨⟩ denote the complex conjugate
and the ensemble averaging, respectively. As X(f, k) is a complex quantity, the crossspectrum or the cross coherence can also be calculated, although they are out of thescope in this paper. See Appendix for an example of the python code for the complex
MEM calculation.

) ( (

3

T

11
Real
ImaginaryFTin timeMEM
in spaceof the2two-dimensional
FT-MEM for the frequency-wavenumber• ( ( Figure. 3.- Analysis2procedure
3
3 T
•
EF
3
3
M
D
spectral decomposition.

4. Application I: Numerical simulation data

The proposed analysis method is applied to the results of a numerical simulation for thecylindrical plasma turbulence performed in the Numerical Linear Device (NLD) [35].
The simulation code provides the three dimensional distribution of plasma quantities,
including in the azimuthal direction (the poloidal direction for the torus geometry)
distribution. By using the entire azimuthal circumference for the spatial spectrum
estimation, the true power spectrum density is obtained since the cyclic boundary
condition is naturally satisfied. Meanwhile, by limiting the available azimuthal angle,
the experimental situation is simulated and the capability of the FT-MEM (for thetime-space domain) method is evaluated by contrasting that of the conventional FT-FT
method and the two-point estimation [14].
The model for the numerical simulation is the global three-field Hasegawa-Wakatani
reduced fluid model for the electron density, the electrostatic potential, and the ion
flow velocity parallel to the magnetic field direction (the axial direction) [36]. Theni

N=4 =.
g
N D

12Figure 4. Spatiotemporal evolution of the electron density fluctuation in the Numerical Linear
Device (NLD) fluid simulation.

plasma profiles and the turbulence dynamics are self-consistently evolved by the applied
sources. The spatially inhomogeneous equilibrium profiles in the electron density, theelectrostatic potential, and the parallel flow may drive the resistive drift wave, theKelvin-Helmholtz instability, and the parallel flow shear driven instability (D’Angelomode), respectively. The excited waves saturate nonlinearly and also back-interact withthe equilibrium profiles, forming a variety of the plasma flow structures [36, 37]. The
plasma parameters are given according to the nominal values in the actual cylindricalplasma discharge experiment [38]. In the present case, an azimuthally inhomogeneous
parallel momentum source is applied to observe the nonlinear interaction between
turbulence and the background flow structure.
Figure 4 shows the spatiotemporal evolution of the electron density fluctuation at
the mid-radius of the cylindrical plasma. Turbulence is excited with a large relative
amplitude that occasionally shows spatially isolated intermittent bursts. The phasepropagation direction seems to be in the positive direction of the azimuthal angle,
corresponding to the electron diamagnetic direction.

doah=4 =.

5

. .

13
(a)

(b)

FT-FT
(Full data)

FT-FT(Limited data)

tar
,,1

E (
T
EeL

(c)

h NtarcFM taE

D

FT-MEM(Limited data)(d)Two-point

Figure 5. Power spectrum densities of the electron density fluctuation in NLD obtained withvarious method: (a) FT-FT with the entire azimuthal circumference, (b) FT-FT with 25 % of
the azimuthal angle, (c) FT-MEM with 25 % of the azimuthal angle, and (d) two-point statisticalmethod based on FT.

Here, the mode number m is used for describing the spatial scale of the fluctuation
instead of the wavenumber k. It is defined as m = kr, where r is the radius of the cylinder
of interest and m > 0 corresponds to modes with the electron diamagnetic propagation.
The frequency-mode number resolved power spectrum densities obtained with various
methods are shown in Fig. 5. The original data has a simulation mesh of 128 points inthe azimuthal direction. Figure 5 (a) shows the FT-FT power spectrum density using
all the azimuthal data points, which corresponds to the true power spectrum density.
In order to simulate the experimental situation, where measurement is only performed
in a limited spatial domain, eight data points with an interval of four data points aretaken. By this data sampling, one fourth of the entire azimuthal circumference is used
for estimating the power spectrum density. Figures 5 (b) and (c) show the results of

14# 156774

(a)

ch 49
ch 49

(b)

ch 50
ch 51

ch 56

ch 52
ch 53
ch 54
ch 55
ch 56

Figure 6. (a) Measurement location and (b) the raw data of the BES in LHD.

the FT-FT method with the Hanning spatial window and the FT-MEM method. Here,M = 1 is used as the data points for LP for MEM. Although the Akaike FPE criterion
suggests M = 2 as the most appropriate value but with a marginal diﬀerence against
other values nearby, M = 1 is chosen to avoid having possible mimic spectral peaks. The
FT-MEM method provides a better mode number resolution than the FT-FT methodthat has less pronounced mode number peaks.
In addition to the FT-FT method and the FT-MEM method, the two-point
statistical method based on FT [14] is also applied for comparison. Although the roughfeatures of the spectrum are captured, statistical convergence seems to be insuﬃcientbecause of the lack of simulation data points particularly for high frequency fluctuation
components.

5. Application II: Experimental data

The FT-MEM method is applied to an experimental data in the Large Helical Device(LHD). The target data is the poloidal distribution of the electron density fluctuation

c15# 156774, t = 4.70-4.97 s(a) FT-FT(b) FT-MEMi(c) Two-point (ch49-ch50)a(d) Two-point (ch49-ch54)

M
LE=
E
(
Sv
B

H

Figure 7. Power spectrum densities of the BES data obtained with various methods: (a) FTFT, (b) FT-MEM, and (c,d) two-point statistical method between ch49-ch50 and ch49-ch54,
respectively.

measured by a beam emission spectroscopy (BES) recently installed in LHD [39].
Figure 6 (a) shows the measurement location of the BES in LHD for the shot numberof # 156774. The electron density data are obtained by 8 × 8 radial-vertical sight line
array with a sampling rate of 200 kHz and an analog bandwidth of 100 kHz. One
column of the array slightly inside the plasma edge (the second rightmost column) is
extracted to examine the poloidal distribution of the fluctuation. Figure 6 (b) is thetime evolution of the electron density fluctuation at diﬀerent vertical positions. The
phase of the dominant fluctuating mode propagates from top to bottom, the electron
diamagnetic direction.

16
The frequency-wavenumber resolved power spectrum densities estimated with
diﬀerent methods are shown in Fig. 7. Positive wavenumber corresponds to modes
with the electron diamagnetic propagation. For MEM, M = 1 is again used. Larger M
values were also examined and the results were found to be insensitive to the choice ofM except for exaggerated noise peaks in higher frequency range. As is demonstrated
for the numerical simulation data, the FT-MEM power spectrum density gives the most
sharpest peaks at each frequency. In contrast to the simulation example shown in section
4, the poloidal coverage of the measurement is much smaller fraction, approximately 4
to 5 % of the entire poloidal circumference. Because of this, the wavenumber resolution
in the FT-FT method is extremely poor. The two-point statistical method using two
adjacent channels (Fig. 7 (c)) also exhibits a less focused wavenumber spectrum, eventhough the number of ensembles is large enough (∼ 200 periods of the dominantfrequency components). When the channel interval for the two-point statistical method
is increased, the wavenumber resolution becomes finer and the wavenumber spectrum
peaks are sharpened as shown by Fig. 7 (d). As a result, the spectral shape resembles tothat obtained by the FT-MEM method. However, as the channel interval is increased,
the Nyqvist wavenumber is reduced and the observable range of the wavenumber shrinks
accordingly. The FT-MEM method reconciles the wide dynamic range and the fineresolution for the wavenumber decomposition. The averaged wavenumbers of the peaks
at f = 1 − 2 kHz and f = 3 − 5 kHz are 0.6 m−1 and 7.2 m−1 , which correspond to
the poloidal mode number of m = 0 − 1 and 4 − 5, respectively. Note that m = 0

17
may be excluded because a poloidally symmetric MHD mode is not realistic. Unlike theexamination of the methods using the numerical simulation, the true spectrum cannot
be referred to, therefore, it is impossible to judge the performance of each method in thiscase. However, considering that the fluctuation is driven by MHD instabilities [39], thepower spectrum density with coherent peaks provided by the FT-MEM method seems
to the the most plausible result.

6. Summary

In this paper, we documented the method for estimating the frequency-wavenumber
resolved two-dimensional power spectrum density through combining the Fouriertransform and the maximum entropy method. Aiming to improve the wavenumber
resolution in a situation, in which available spatial points are severely limited, themaximum entropy method was employed for the spatial structure decomposition. The
proposed method was applied to numerical fluid simulation data and experimental dataobtained by the beam emission spectroscopy in LHD. Through those examinations,
it was found that the proposed method can provide more plausible power spectrum
densities than the conventional two-dimensional Fourier transform method or the twopoint statistical method.

18
Acknowledgments

Authors acknowledge T. Ido and T. Morisaki for useful discussion. One of the authors(T. K.) wishes to thank Y. Suzuki and S. Sakakibara for strong support. This work is
partly supported by the National Institute for Fusion Science grants (ULHH033) andby the Grant-in-Aid for Scientific Research of JSPS (17K14898).

Appendix

The python subroutines for obtaining the MEM estimate [26] is shown below. The main
routine mcme_k calculates the complex MEM cross spectrum from a multichannel timeseries of data. The subroutine fx2fk_k performs MEM to the multichannel Fouriercomponents in spacial dimension for each frequency, i.e., the procedure from the second
panel to the third panel in Fig. 3. For this spacial MEM, only the auto spectrum
calculation is used.
#!/ usr / bin / env python# coding : utf -8## Python subroutines for calculating complex MEM .
# Written by Tatsuya Kobayashi ( National Institute for Fusion Science ).
# Email : kobayashi . tatsuya@nifs . ac . jp

import numpy as np

def fx2fk_k ( dx , dk , Ns , sigRe , sigIm ):19# dx : Spatial resolution# dk : Desirable wavenumber resolution# Ns : Number of data points for the linear prediction# sigRe : sigRe ( frequency , space ) , Real part of input spectrum# sigIm : sigIm ( frequency , space ) , Imaginary part of input spectrum
nf = sigRe . shape [0]nx = sigRe . shape [1]mk = dx /2.
kwave = np . arange ( - mk , mk + dk , dk )nk = len ( kwave )y = np . zeros (( nx , 1) , dtype = np . complex )spct = np . zeros (( nf , nk , 1))
for i in range ( nf ):y [: ,0] = sigRe [i ,:] + sigIm [i ,:]*1 jspct [i ,: ,0] = np . abs ( mcme_k ( y [: ,:] , dx , kwave , Ns )[: ,0 ,0])
return kwave , spct

def mcme_k (y , dt , freq , Ns ):# y : y ( time , space ) , Input complex two - dimensional data
# dt : Time resolution# freq : Frequency vector on which MEM is calculated# Ns : Number of data points for the linear prediction
df = len ( freq )z = np . exp ( -2* np . pi *1 j * freq * dt )Nd = y . shape [0]p= y . shape [1]20
Nk = min ( Ns , Nd -1)
spct = np . zeros (( df ,p , p ) , dtype = np . complex )F , B , Pm , Pd = mcmecof_k (y , Ns )for i in range ( df ):iFz = np . linalg . inv ( np . sum ( F [: Nk +1 ,: ,:]* zext ( z [ i ] , Nk , p ) , axis =0))
spct [i ,: ,:] = dt * np . dot ( np . dot ( iFz . conjugate (). T , Pm ) , iFz )return spct

def zext (z , n , p ):return np . repeat ( z ** - np . arange ( n +1) , p **2). reshape ( n +1 ,p , p )def mcmecof_k (y , Ns ):Nd = y . shape [0]p= y . shape [1]Nk = min ( Ns , Nd -1)
F = np . zeros (( Nd ,p , p ) , dtype = np . complex )B = np . zeros (( Nd ,p , p ) , dtype = np . complex )Pm = np . zeros (( p , p ) , dtype = np . complex )F [0 ,: ,:] = np . eye ( p )B [0 ,: ,:] = np . eye ( p )for i in range ( p ):for j in range ( p ):Pm [i , j ] = np . sum ( y [: , i ]* y [: , j ]. conjugate ())/ Nd
Pd = Pm . copy ()21
er_f = y [1: Nd ,:]er_b = y [0: Nd -1 ,:]for N in range (1 , Nk +1):
NN = N
er_f , er_b , Pm , Pd , F , B = mcmecof_calc_k ( er_f , er_b , Pm , Pd , F , B , NN )return (F , B , Pm , Pd )def mcmecof_calc_k ( er_f , er_b , Pm , Pd , F , B , N ):Nd = F . shape [0]p= F . shape [1]M = Nd - N
Wm = 1./ float ( M )Fp = np . zeros (( Nd ,p , p ) , dtype = np . complex )Bp = np . zeros (( Nd ,p , p ) , dtype = np . complex )er_fp = np . zeros (( Nd -1 , p ) , dtype = np . complex )er_bp = np . zeros (( Nd -1 , p ) , dtype = np . complex )Em = np . zeros (( p , p ) , dtype = np . complex )Gm = np . zeros (( p , p ) , dtype = np . complex )Bm = np . zeros (( p , p ) , dtype = np . complex )exCOEF = np . zeros (( p **2 , p **2) , dtype = np . complex )C = np . zeros (p , dtype = np . complex )for i in range ( p ):for j in range ( p ):22
Em [i , j ] = np . sum ( Wm * er_f [: M , i ]* er_f [: M , j ]. conjugate ())
Gm [i , j ] = np . sum ( Wm * er_b [: M , i ]* er_f [: M , j ]. conjugate ())
Bm [i , j ] = np . sum ( Wm * er_b [: M , i ]* er_b [: M , j ]. conjugate ())
iPm = np . linalg . inv ( Pm )for i in range ( p ):for j in range ( p ):for k in range ( p ):for l in range ( p ):exCOEF [ i * p +j , k * p + l ] = Bm [i , k ]* int ( l == j ) + np . sum ( Pd [i , k ]* iPm [l ,:]*
C = -2.* Gm . reshape ( p **2)
Cm = np . dot ( np . linalg . inv ( exCOEF ) , C ). reshape (p , p )Cd = np . dot ( np . dot ( iPm , Cm . conjugate (). T ) , Pd )Fp [0 ,: ,:] = F [0 ,: ,:]Bp [0 ,: ,:] = np . dot ( F [0 ,: ,:] , Cd )for i in range (1 , N ):Fp [i ,: ,:] = F [i ,: ,:] + np . dot ( B [i -1 ,: ,:] , Cm )Bp [i ,: ,:] = B [i -1 ,: ,:] + np . dot ( F [i ,: ,:] , Cd )Fp [N ,: ,:] = np . dot ( B [N -1 ,: ,:] , Cm )Bp [N ,: ,:] = B [N -1 ,: ,:]F = Fp . copy ()B = Bp . copy ()Pmp = Pm - np . dot ( np . dot ( Cm . conjugate (). T , Pd ) , Cm )Pdp = Pd - np . dot ( np . dot ( Cd . conjugate (). T , Pm ) , Cd )

23for i in range ( Nd -N -1):
er_fp [i ,:] = er_f [ i +1 ,:] + np . dot ( Cm . conjugate (). T , er_b [ i +1 ,:])
er_bp [i ,:] = er_b [i ,:]+ np . dot ( Cd . conjugate (). T , er_f [i ,:])

er_f = er_fp . copy ()er_b = er_bp . copy ()Pm= Pmp . copy ()Pd= Pdp . copy ()return ( er_f , er_b , Pm , Pd , F , B )def mcme_fpe_k (y , Ns_max ):Nd = y . shape [0]p= y . shape [1]Ns = min ( Ns_max , Nd -1)
fpe = np . zeros ( Ns )for i in range (0 , Ns ):N = i + 1
F , B , Pm , Pd = mcmecof_k (y , N )fpe [ i ] = np . real ( np . linalg . det ( Pm ))*( float ( Nd +1.+ p * N )/ float ( Nd -1. - p * N ))** p
return fpe

24References[1] C Schröder, T Klinger, D Block, A Piel, G Bonhomme, and V Naulin 2001 Phys. Rev. Lett. 86
5711[2] O Grulke and T Klinger 2002 New J. Phys. 4 67
[3] T Windisch, O Grulke, Volker Naulin, and T Klinger 2011 Plasma Phys. Control. Fusion 53 085001[4] T Yamada, S-I Itoh, T Maruta, N Kasuya, Y Nagashima, S Shinohara, K Terasaka, M Yagi,
S Inagaki, Y Kawai, et al 2008 Nat. Phys. 4 721–725[5] T Kobayashi, S Inagaki, M Sasaki, Y Kosuga, H Arakawa, F Kin, T Yamada, Y Nagashima,
N Kasuya, A Fujisawa, et al 2017 Plasma Fusion Res. 12 1401019[6] U Stroth, F Greiner, C Lechte, N Mahdizadeh, K Rahbarnia, and M Ramisch 2004 Phys. Plasmas
11 2558–2564[7] B Schmid, P Manz, M Ramisch, and U Stroth 2017 Phys. Rev. Lett. 118 055001[8] A Latten, T Klinger, A Piel, and Th Pierre 1995 Rev. Sci. Instrum. 66 3254–3262[9] N Krause, C Lechte, J Stöber, U Stroth, E Ascasibar, J Alonso, and S Niedner 2002 Rev. Sci.
Instrum. 73 3474–3481[10] T Yamada, Y Nagashima, S Inagaki, Y Kawai, M Yagi, S-I Itoh, T Maruta, S Shinohara,
K Terasaka, M Kawaguchi, et al 2007 Rev. Sci. Instrum. 78 123501[11] S Kado, H Nakatake, K Muraoka, K Kondo, F Sano, T Mizuuchi, S Besshou, H Okada, K Nagasaki,
H Funaba, et al 1996 J. Phys. Soc. Jpn. 65 3434–3437[12] T Ido, Y Miura, K Hoshino, K Kamiya, Y Hamada, A Nishizawa, Y Kawasumi, H Ogawa,Y Nagashima, K Shinohara, et al 2006 Nucl. Fusion 46 512[13] Z Yan, P Gohil, G R McKee, D Eldon, B Grierson, T Rhodes, and C C Petty 2017 Nucl. Fusion
57 126015[14] J M Beall, Y C Kim, and E J Powers 1982 J. Appl. Phys. 53 3933–3940[15] M VanZeeland, W Gekelman, S Vincena, and J Maggs 2003 Phys. Plasmas 10 1243–1252

25[16] M J Burin, G R Tynan, G Y Antar, N A Crocker, and C Holland 2005 Phys. Plasmas 12 052320[17] P Manz, M Xu, N Fedorczak, S C Thakur, and G R Tynan 2012 Phys. Plasmas 19 012309[18] A Möller 1998 Phys. Plasmas 5 4321–4327[19] Y Xu, M Van Schoor, R R Weynants, S Jachmich, M Vergote, M W Jakubowski, P Beyer, M Mitri,
B Schweer, D Reiser, et al 2007 Nucl. Fusion 47 1696
[20] A D Liu, T Lan, C X Yu, H L Zhao, L W Yan, W Y Hong, J Q Dong, K J Zhao, J Qian, J Cheng,
et al 2009 Phys. Rev. Lett. 103 095002[21] T Kobayashi, Y Nagashima, S Inagaki, H Arakawa, M Sasaki, T Yamada, M Yagi, N Kasuya,
A Fujisawa, S-I Itoh, et al 2011 Plasma Fusion Res. 6 1401050[22] J P Burg. Maximum entropy spectral analysis. In Proceedings of the 37th Meeting of the Societyof Exploration Geophysicists, Oklahoma City, Okla, 1967.
[23] T J Ulrych 1972 J. Geophys. Res. 77 1396–1400[24] T J Ulrych and T N Bishop 1975 Rev. Geophys. Space Phys. 13 183–200[25] T J Ulrych and Clayton R W 1976 Phys. Earth Planet. Inter. 12 188–200[26] O Strand 1977 IEEE Trans. Autom. Control 22 634–640[27] W H Press, S A Teukolsky, B P Flannery, and W T Vetterling. Numerical recipes in C, chapter 13.
Cambridge University Press, 1992.
[28] Y Amagishi, K Saeki, and I J Donnelly 1989 Plasma Phys. Control. fusion 31 675[29] K Tanaka, C A Michael, L N Vyacheslavov, A L Sanin, K Kawahata, T Akiyama, T Tokuzawa,and S Okajima 2008 Rev. Sci. Instrum. 79 10E702[30] Z Shi, Y Nagayama, S Yamaguchi, D Kuwahara, T Yoshinaga, S Sugito, Y Hirano, H Koguchi,
S Kiyama, H Sakakita, et al 2010 Plasma Fusion Res. 5 S1019[31] K Matsuo, H Iguchi, S Okamura, and K Matsuoka 2012 Rev. Sci. Instrum. 83 013501[32] C A Michael, K Tanaka, L Vyacheslavov, A Sanin, and K Kawahata 2015 Rev. Sci. Instrum. 86
093503[33] M Ono, K Ida, T Kobayashi, M Yoshinuma, G R McKee, Z Yan, K H Burrell, and X Chen.

26
Frequency-Poloidal Wave Number Spectral Analysis of Turbulence in QH-mode plasmas
Measured with BES on DIII-D. In 58th Annual Meeting of the APS Division of Plasma Physics,
page PP10.074, San Jose, California, 2016.
[34] H Akaike 1970 Ann. Inst. Statist. Math. 22 203–217[35] N Kasuya, M Yagi, M Azumi, K Itoh, and S-I Itoh 2007 J. Phys. Soc. Jpn. 76 044501[36] M Sasaki, N Kasuya, K Itoh, S Toda, T Yamada, Y Kosuga, Y Nagashima, T Kobayashi,
H Arakawa, K Yamasaki, et al 2017 Phys. Plasmas 24 112103[37] M Sasaki, Yann Camenen, A Escarguel, S Inagaki, N Kasuya, K Itoh, and T Kobayashi 2019 Phys.
Plasmas 26 042305[38] T Kobayashi, S Inagaki, Y Kosuga, M Sasaki, Y Nagashima, T Yamada, H Arakawa, N Kasuya,
A Fujisawa, S-I Itoh, et al 2016 Phys. Plasmas 23 102311[39] T Kobayashi, M Yoshinuma, and K Ida 2020 Plasma Phys. Control. Fusion 62 125011