Introduction to the Diagnosis of Magnetically Confined Thermonuclear - - PowerPoint PPT Presentation

Fourier Analysis Tomographic reconstruction

Introduction to the Diagnosis of Magnetically Confined Thermonuclear Plasma

Selected Data Analysis Techniques

• J. Arturo Alonso

Laboratorio Nacional de Fusión EURATOM-CIEMAT E6 P2.10 arturo.alonso@ciemat.es

version 0.1 (February 13, 2012)

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 1 / 20

Fourier Analysis Tomographic reconstruction

Outline

1

Fourier Analysis Fourier transform Some aplications

2

Tomographic reconstruction Pixel methods Singular Value Decomposition

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 2 / 20

Fourier Analysis Tomographic reconstruction

Outline

1

Fourier Analysis Fourier transform Some aplications

2

Tomographic reconstruction Pixel methods Singular Value Decomposition

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 3 / 20

Fourier Analysis Tomographic reconstruction

The Fourier Transform

Given a function f(x) we define its Fourier Transform as

Fourier Transform

F(s) = Ff ≡ ∞

−∞

f(x)e−2πixsdx . We assume this integral exist and can be inverted as

Inverse Fourier Transform

f(x) = F−1F ≡ ∞

−∞

F(s)e2πixsds . Which follows from a definition of the Dirac’s delta (generalised) function δ(x) = F−1{1} = ∞

−∞ exp(2πixs)ds

f(x) = F−1Ff = ∞

−∞

ds ∞

−∞

dx′f(x′)e2πis(x−x′) = ∞

−∞

dx′f(x′)δ(x−x′)

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 4 / 20

Fourier Analysis Tomographic reconstruction

The convolution theorem

The convolution of two functions f and g is a third function f ∗ g defined as (f ∗ g)(x) ≡ ∞

−∞

f(x′)g(x − x′)dx′ , then the convolution theorem states that F{f ∗ g} = (Ff)(Fg) Proof: F{f ∗ g} = ∞

−∞

∞

−∞

f(x′)g(x − x′)dx′

• e−2πisxdx

= ∞

−∞

∞

−∞

g(x − x′)e−2πisxdx

• f(x′)dx′

= ∞

−∞

G(s)e−2πisx′f(x′)dx′ = G(s)F(s) It also follows that F−1{F ∗ G} = fg.

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 5 / 20

Fourier Analysis Tomographic reconstruction

Real functions

For a general complex function f(x) we have F{f ∗(x)} =

• f ∗(x)e−2πisxdx =
• f(x)e2πisxdx

∗ = F∗(−s) . For a real function f = f ∗ F(s) = F(−s)∗ i.e., only half of the (complex) transform of a real function contains non-redundant information.

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 6 / 20

Fourier Analysis Tomographic reconstruction

Energy conservation

From the identity (f ∗ h)(x) = F−1{FH} evaluated at x = 0 we have ∞

−∞

f(x′)h(−x′)dx′ = ∞

−∞

F(s)H(s)ds

• r, writting g∗(x′) = h(−x′) ⇒ G∗(−s) = H(−s)

∞

−∞

f(x)g∗(x)dx = ∞

−∞

F(s)G∗(s)ds know as the Power Theorem. If g = f then one obtains Rayleigh’s Theorem ∞

−∞

|f(x)|2dx = ∞

−∞

|F(s)|2ds

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 7 / 20

Fourier Analysis Tomographic reconstruction

Periodic Functions: Fourier Series

A periodic function f(x) : f(x + nT) = f(x) such as a physical quantity (say B) defined over a periodic domanin (say the poloidal angle θ) has a discrete transform1. f(x + nT) = ∞

−∞

F(s)e2πis(x+nT)ds = f(x) ⇒ F(s) = 0 , except if e2πisnT = 1 → sk = k T In this case the function f(x) can be represented as a Fourier series f(x) =

∞

• k=−∞

ckeiωkx ; ωk = 2π T k ; ck = T/2

−T/2

f(x)e−iωkxdx

1In reality the intergral is non-convergent, but it still illustrates that only

some frequencies are present

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 8 / 20

Fourier Analysis Tomographic reconstruction

The Discrete-time Fourier Transform (DFT)

A physical measurement always consist of a vector with finite number N of discrete samples fτ = f(τ∆t) of an ideal continuos function f(t). We define the

Discrete Fourier Transform

Fν = N−1

N−1

• τ=0

fτe−2πi(ν/N)τ ; fτ =

N−1

• ν=0

Fνe2πi(ν/N)τ The consistency of these transformation follows from the fact that for any two intergers τ, τ ′ N−1

N−1

• ν=0

e−2πi(ν/N)(τ−τ ′) = 1 if τ = τ ′

• therwise .

There are DFT versions of the FT theorems (convolution, power, etc.)

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 9 / 20

Fourier Analysis Tomographic reconstruction

Frequencies in the DFT

The following relations hold for real fτ Fν = Fν+nN Fν = F∗

−ν = F∗ N−ν

F N

2 +ν = F∗ N 2 −ν

The discrete Fourier component Fν is associated with the physical frequency νp =

ν N∆t.

Frequency restolution is then ∆νp = 1/N∆t = 1/T Maximum frequency is that of ν = N/2 or νNy

p

= 1/2∆t

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 10 / 20

Fourier Analysis Tomographic reconstruction

Average over realizations

Estimating real data spectral properties.

Experimental data often consist of 220 ≈ 106 samples. Rather than using all the points to get a single spectrum can be split into realizations and average the spectrum

• ver them to gain some statistical significance.

100 200 300 400 500 600 700 800 900 1000 10

−10

10

−8

10

−6

10

−4

10

−2

10 f (kHz) Power Spectrum (a.u.) FFT of all 217 points Periodogram 210 points x 27 windows Welch 210 points 50% ovelap windows 10 20 30 40 50 60 −50 50 Floating potential (2 MHz sampling) Selected Data Analysis Techniques, A. Alonso, copyleft 2010 11 / 20

Fourier Analysis Tomographic reconstruction

Spectrum evolution

Time (ms) Frequency (kHz) 50 100 150 200 250 200 400 600 800 1000 TJII#18998 Vfloat

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 12 / 20

Fourier Analysis Tomographic reconstruction

Cross power spectrum and coherence

Just as the power spectrum is of a signal x(t) estimated as Px(ω) = XωX∗

ω where · is the average over ‘realisations’,

the cross-power spectrum of x(t) and y(t) is defined as Pxy = XωY∗

ω

The cross-coherence is defined as Cxy(ω) = |Pxy(ω)|2 Px(ω)Py(ω) for stationary spectra Cxy(ω) ∼ |ei(φx(ω)−φy(ω))|2 so it is high for frequencies with approx. constant phase difference.

10

3

10

4

10

5

10

6

0.2 0.4 0.6 0.8

Cxy(ν) ν (Hz)

TJII#18998 Vfloat 1 and 2 coherence 1040 1050 1060 1070 1080 −50 50

x(t), y(t)

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 13 / 20

Fourier Analysis Tomographic reconstruction

The Continuous Wavelet Transform

Wavelets were originally intended to analyse highly non-stationary and inhomogeneous signals (as human speech). FT provides all the frequency information about the signal but no temporal localization of the frequencies. WT give certain degree of localization in both frequency and time (local frequency analyser). It is defined as WTψ{x}(a, b) =

• R

x(t)ψa,b(t) dt with ψa,b(t) = 1 √aψ t − b a

• Selected Data Analysis Techniques, A. Alonso, copyleft 2010

14 / 20

Fourier Analysis Tomographic reconstruction

Outline

1

Fourier Analysis Fourier transform Some aplications

2

Tomographic reconstruction Pixel methods Singular Value Decomposition

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 15 / 20

Fourier Analysis Tomographic reconstruction

Tomography: posing the problem

To recover the spatial distribution of, say, the emissivity g(r) from line-integrated measurements fl =

• Sl

dsg(r) , l = 1 . . . nl Ill-posed problem: never sufficient information in a finite number of chords to determine g(r) Needs to reduce the function search space

1 Expansion in a finite set of

• rtogonal basis functions and fit

coefficients to data (Cormak)

2 Pixel methods (or expasion in

‘pixel’ basis functions)

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 16 / 20

Fourier Analysis Tomographic reconstruction

Tomography: posing the problem

To recover the spatial distribution of, say, the emissivity g(r) from line-integrated measurements fl =

• Sl

dsg(r) , l = 1 . . . nl

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 16 / 20

Fourier Analysis Tomographic reconstruction

Pixel methods (I)

Space is divided into (possibly non-square) pixels of approximatelly constant g Arrange the nl f’s and npixels = nx × ny gij’s in two column vectors f and g Then compute the ‘contribution’ nl × npixel matrix T f = Tg Direct inversion is usually inpractical. Instead one tries to minimise the deviation from the data min

g

1 2ξ⊺ξ , ξ(g) = Tg − f .

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 17 / 20

Fourier Analysis Tomographic reconstruction

Pixel methods (and II)

The set of normal equation are obtained imposing

∂ ∂g(ξ⊺ξ) = 0

T⊺Tg = T⊺f A commons situation is that npixel > nl so g is not uniquelly determined To obtain a unique and ‘sensible’ solution we introduce a regularising functional R(g) which imposes regularity on the solution of min

g

1 2ξ⊺ξ + αR I.e. we can require g to be smooth choosing a first-order regulariser R = ∇xg2 + ∇xg2 where ∇x in the finite differences x-derivative matrix. Normal equations in this case read (T⊺T + αH)g = T⊺f , with H = ∇⊺

x∇x + ∇⊺ y∇y

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 18 / 20

Fourier Analysis Tomographic reconstruction

Singular Value Decomposition (I)

A tomographic reconstruction of the emissivity g can be

• btained for different times. This is composed of an

average, roughly flux-constant emissivity and possibly some short-lived modes superimpossed on it. A common technique used to separate the contributions of g(t) is the Singular Value Decomposition (SVD) which decoposes a matrix G as the product of the matrices G = U × S × V† where U and V are unitary (U⊺ = U−1) and S is a diagonal matrix of singular values that are positive and sorted in descending order. The signal energy EG = M

m

N

n G2 nm = K k S2 k so S2 k/EG

measures the relative importance of the uk, vk basis vectors.

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 19 / 20

Fourier Analysis Tomographic reconstruction

Singular Value Decomposition (and II)

Think of Gij = gi(tj) so columns in G are instantaneous tomographic reconstructions. Columns in U, uk are spatial ‘eigenvectors’ or topos whereas vk are temporal eigenvectors or chronos. The example shows the average soft x-ray emissivity matching the flux surfaces (first topos) and a m = 2 mode (second topos). The third topos is esentially noise.

Selected Data Analysis Techniques, A. Alonso, copyleft 2010 20 / 20