Table of Contents
International Journal of Spectroscopy
Volume 2017, Article ID 9265084, 29 pages
https://doi.org/10.1155/2017/9265084
Research Article

Fourier Spectroscopy: A Bayesian Way

1CCFE, Culham Science Centre, Abingdon OX14 3DB, UK
2Max-Planck-Institut für Plasmaphysik, Teilinstitut Greifswald, Wendelsteinstraße 1, 17491 Greifswald, Germany

Correspondence should be addressed to Stefan Schmuck; ti.rnc.pfi@kcumhcs

Received 6 February 2017; Revised 11 July 2017; Accepted 27 July 2017; Published 31 October 2017

Academic Editor: Wei Kong

Copyright © 2017 Stefan Schmuck and Jakob Svensson. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

The concepts of standard analysis techniques applied in the field of Fourier spectroscopy treat fundamental aspects insufficiently. For example, the spectra to be inferred are influenced by the noise contribution to the interferometric data, by nonprobed spatial domains which are linked to Fourier coefficients above a certain order, by the spectral limits which are in general not given by the Nyquist assumptions, and by additional parameters of the problem at hand like the zero-path difference. To consider these fundamentals, a probabilistic approach based on Bayes’ theorem is introduced which exploits multivariate normal distributions. For the example application, we model the spectra by the Gaussian process of a Brownian bridge stated by a prior covariance. The spectra themselves are represented by a number of parameters which map linearly to the data domain. The posterior for these linear parameters is analytically obtained, and the marginalisation over these parameters is trivial. This allows the straightforward investigation of the posterior for the involved nonlinear parameters, like the zero-path difference location and the spectral limits, and hyperparameters, like the scaling of the Gaussian process. With respect to the linear problem, this can be interpreted as an implementation of Ockham’s razor principle.

1. Introduction

Fourier spectroscopy is a diagnostic application which reveals information about spectral quantities like refractive index, absorption, and transmission of a medium under test. In addition, the characterisation in absolute terms is possible for broadband spectra, for example, emitted by electrons of a high-temperature plasma, being magnetically confined [1].

Commonly, an interferometer diagnostic, let us say of Michelson [2] or Martin-Puplett [3] design type, probes the Fourier transform of a spectral quantity. The corresponding interferometric data is a discrete set with finite length and includes noise contributions. Standard Fourier data analysis techniques [46] have been developed. These techniques lack describing and capturing properly several fundamental aspects, like the noisy nature of measured data and possible spectral limits and their impact on the spectral quantity to be inferred.

One misconception, arising from the standard formulation, is that certain spectral information must be lost inherently, because only a finite amount of data is acquired. This is proven in standard literature by evaluating the convolution function which has a finite full width at half maximum (FWHM), implying that only a finite amount of Fourier coefficients is accessible via measurements. While this conclusion remains valid when a continuous spectrum is probed, the reasoning does not hold in general for a discrete spectrum. This fact was exploited to develop a (self-)deconvolution procedure, so that some discrete lines which were separated by less than the FWHM width of the convolution function have been inferred [7].

Opposed to the standard data analysis techniques, a probabilistic ansatz was introduced to estimate model parameters, like amplitude spectra and frequencies, in the field of Fourier spectroscopy [8, 9]. For a spectroscopic problem, a model is formulated via Bayes’ theorem which allows to state prior information about model parameters. Furthermore, a Gaussian likelihood connects functionally the parameters with the noisy interferometric data. Then, after having measured a noisy data set and framing the reality by a certain model, the knowledge about model parameters is expressed in probabilistic terms by the posterior. This approach [8, 9] demonstrated that the uncertainty on the posterior mean of a frequency to be estimated for a single-frequency problem can be orders of magnitude below the FWHM width of the equivalent convolution function. On top of that, a criterion has been derived how far frequencies need to be separated, so that the Bayesian approach is still able to make a distinction. This separation can be well below the FWHM width of the convolution function. These findings are a direct consequence of the use of probabilistic theory.

One of the advantages of the Bayesian approach is that different models and, hence, their assumptions can be compared with each other also known as Ockham’s razor. This enables the identification of the best, that is, the most likely, model, complying with the data. Given this context, the (self-)deconvolution procedure mentioned above is interpreted here as an optimisation to find a minimising set of discrete frequencies to describe the data sufficiently. In general, the most fundamental issue is whether the spectrum to be inferred is discrete or continuous. If a discrete spectrum is more likely or follows by a physics model, how many discrete frequencies are involved and what are their estimates including uncertainties? Each frequency is associated with an amplitude and a phase which need to be estimated as well. If a continuous spectrum is at hand, then the spectral limits are of interest. In addition, in case, the underlying physics process is understood, what is the uncertainty on the spectrum and the phase following from experimentally inaccessible regions in the data domain.

The investigation of the fundamental issues listed above is quite challenging from numerical point of view. However, the computational effort is largely reduced when even and odd amplitudes are used instead of the phase and amplitude. This ensures a linear dependence between the even and odd amplitude parameters and the data. Formulating the prior information about all even and odd amplitudes as a multivariate normal with a specific prior mean and covariance gives straightforwardly the posterior mean and covariance. Furthermore, the marginalisation can be carried out analytically for these linear parameters. The remaining posterior quantity carries information about the nonlinear model parameters like frequencies or spectral limits and so-called hyperparameters, entering merely in the prior.

In Section 2, the basic equations for Fourier spectroscopy and their implications are generally investigated, and the main concepts of the standard analysis and their drawbacks are pointed out. Section 3 of this paper presents a Bayesian formalism, so that it lends itself to applications in the field of Fourier spectroscopy for Gaussian (white) noise. The fundamental information about the spectral quantity, that it must vanish at a lower and upper spectral limit, can be stated by the covariance function of a Brownian bridge process as described in Section 4. This covariance is used as prior in the example application of the Bayesian approach presented in Section 5 to infer continuous but band-limited spectral quantities given an actually measured interferometric data set. Thereby, some diagnostic imperfections like a drifting signal offset, the zero-path difference, and a nonuniform spatial sampling are also taken into account. Section 6 discusses the strategy for a plausibility study of models, using different priors for the spectral quantities, and attempts to compare results and computational efforts obtained with the Bayesian model and with a standard model. The last section presents the conclusions.

2. Fourier Spectroscopy

2.1. General Definitions

Commonly, the complementary coordinates used for Fourier transformations are the frequency and time , or the wavenumber and a spatial coordinate . For the moment, the latter pair is used to state the basic operations of Fourier transformation. Afterwards, the wavenumber is replaced via for convenience.

The real-valued continuous functions and form a Fourier transformation pair stated byThe inverse operation readsNote that so far holds. To find a relation, one inserts (1) in the above expression. After applying trigonometric identities, the spatial integral becomesbecause the sinusoidal contribution vanishes. Hence,follow, because the two wavenumber coordinates equal each other. In fact, the delta distribution occurs, because is a distribution.

Replacing with gives in the spectral domain the quantity . This givesand rescales the inverse operation liketo match the units /Hz.

2.2. Selection of Representation

The spectral domain over which the integration is performed in (5) includes the whole negative and the positive ranges. While the cosine transform acts on the even part of , the sine transform is linked only to the odd part . Thus, an alternative formulation readsand it becomes clear that the even and the odd parts of and are connected.

Another representation uses the amplitude and the phase by settingwhich givesBecause the model representations (5) and (7) are linear in , , and , both are favoured over the formulation (9), when linear inversion techniques are to be applied. Of both favoured representations, the one, using even and odd parts, has preferred properties. Since the orders of magnitude of the amplitudes can be quite different for the even and odd functions, a separation is logical.

2.3. Finite Bandwidth

If the spectral domain is band-limited, such that and are finite for the range from and with the bandwidth and the centre frequency , (7) becomesHere, the assumption must be mentioned that both functions have the same spectral limits. This is assumed in the following but not mandatory. Furthermore, for any combination of and the relation must be fulfilled to be meaningful.

2.3.1. Relation: Fourier Transform-Fourier Coefficients

When the bandwidth is finite, one can express the spectral functions by Fourier coefficients multiplied each with the associated sinusoidal basis function of order (). These coefficients are defined here by the integralswhich carry the unit as and . The coefficients label the mean values of and in the spectral domain covered. Then, one can replace the even and odd functions withwhich allows performing the spectral integration in (10) analytically. Since the resultwithfollows, it becomes clear that the Fourier transform of a band-limited function can be expressed by Fourier coefficients scaled with the bandwidth and multiplied each with the associated continuous basis function in the spatial domain. These basis functions have two contributions. The first is a sum/difference of two functions which depend on the order , the spatial coordinate, and the bandwidth. The latter quantity determines the spatial width of . Furthermore, the localisation is permitted at , where a coefficient for a given mainly acts. Hence, increasing the order implies the localisation at a larger distance from the spatial origin. This explains the occurrence of factor 2 for for which both sinc functions coincide.

The second contribution causes a modulation of and is given by a sine/cosine with the centre frequency and spatial coordinate in the argument. This dependency makes the basis function for vanish at the spatial origin.

With respect to the spatial origin, the transformed basis functions of the coefficients for and are symmetric and antisymmetric, respectively.

Some basis functions in the spatial domain are shown in Figures 1(a) and 1(b) for = 500 GHz and = 1000 GHz.

Figure 1: Basis functions in spatial domain of Fourier coefficients (of order = 0, 1, and 12) , for even spectral function (a) and , for odd spectral function (b). The spectral functions are finite for the interval from 0 GHz to 1000 GHz for which the centre frequency = 500 GHz and bandwidth = 1000 GHz follow. The basis functions in the spatial domain are given by a sum or difference of two sinc functions, depending each on the optical path difference , , and . Furthermore, the sinc functions are modulated by sine or cosine with and in the argument. The sinc functions for fixed become unity at , obviously increasing in .
2.3.2. Embedding into Larger Spectral Domain

The functions and may be finite in the spectral domain with limits and or centre frequency and bandwidth . Embedding this domain in a larger one with limits and (), another set of Fourier coefficients and is obtained with associated basis functions for the spectral domain. Without going into more detail here, these coefficients can be evaluated from and and the scalar products of the basis functions labeled with and . For instance, one finds for the ratio of the means . Basically, and maximise the information per coefficient when and are known. For example, a function which is constant inside a spectral domain and zero outside appears as a boxcar function from outside this domain. Thus, the only coefficient is mapped to an infinite number of coefficients and which are mandatory to capture both discontinuities.

In the spatial domain, the basis functions for and behave differently than the ones for and . Important to mention is the effect which the larger bandwidth has; are spatially narrower than . In addition, the number of coefficients per spatial domain increases which is expressed by .

For = 1873.7 GHz and Figure 2 shows some basis functions (, 1, and 2) in the spatial domain. Indeed, the basis function for with = 500 GHz and = 1000 GHz (see Figure 1(a)) is broader.

Figure 2: Basis functions in spatial domain of Fourier coefficients , with order , 1, and 2 for even spectral function, being finite in spectral domain with centre frequency = 1873.7 GHz and bandwidth . For this domain, the sinc functions, contributing to the basis functions, are narrower than for = 500 GHz and = 1000 GHz (see Figure 1(a)). As an example, the basis function for (dashed-black) is given. Since holds, approximately 4 times more Fourier coefficients locate inside a given spatial domain. was chosen to equal the Nyquist frequency , having set the spatial sampling increment μm (black dots).
2.4. Parseval’s Theorem

Parseval’s theorem states abstractly that the length of a function in the spectral domain equals the length of its Fourier transform counterpart in the spatial domain. The length of the band-limited function is stated bybecause only the term remains odd and cancels by the integration. Furthermore, the scaling by the factor appears which originates in . Replacing with the expression (12) and exploiting that the basis functions are perpendicular in the spectral domain leavesThus, the length is given by the sum of the Fourier coefficients squared. According to Parseval’s theorem,the length in the spectral and spatial domain remains unchanged. Inserting (13) in the above expression implies that the spatial basis functions must be orthogonal for , and the spatial integral yields for and for . Analytically, this is hard to prove; however, this was numerically investigated and is considered to be valid.

2.4.1. Square-Integrable Functions

The function is said to be square-integrable, when the condition holds. Furthermore, if is square-integrable, then the Fourier series representations in (12) converges towards and almost everywhere in the spectral domain as the order grows [10]. Hence, the requirement on to be square-integrable seems reasonable.

2.5. Interferometric Data and Basic Model
2.5.1. Ideal and Real-World Interferometer

The Fourier transform can be performed by an interferometer, achieving an optical path difference between two partial beams, and the real-valued function can be sampled. From theoretical point of view, with an ideal interferometer diagnostic, a purely symmetric and noiseless interferogram is acquired. However, a real-world interferometer suffers from diagnostic imperfections like, for example, dispersion of any kind and/or misalignment. As a consequence, any acquired interferogram is to some degree asymmetric, and, hence, an odd feature is inherent due to the measuring principle. Furthermore, a measurement involves noise, always.

2.5.2. Spatial Sampling and Implications

In the spatial domain, is sampled at a finite set of optical path difference locations with , and marks the number of sample points. Usually, the sampling with constant increment between subsequent locations is preferred which puts constraints on the diagnostic design. Furthermore, the spatial origin is most likely missed by the sampling, and, thus, the absolute value of might be unknown. If so, it is mandatory to introduce the zero-path difference which is in the following set that holds.

The finite spatial sampling leaves undetermined between the sampling nodes and outside the limits and . Assuming holds, the Nyquist theorem states that the maximum frequency accessible is given by the Nyquist frequency . Hence, for the spectral quantities and to be inferred, a maximum for the upper limit follows from sampling theory. To prevent aliasing, needs to be chosen small enough so that and vanish below . In case, can be acted on by reducing the diagnostic throughput via optical filters, the transmission line, the detector sensitivity, and postdetection amplifier settings. In fact, solely by these precautions, one can make sure that no other band well above contributes to . If and only if no such band exists, then the interferogram is smooth with respect to the chosen sampling nodes, and missing to sample exactly at has no profound impact.

A diagnostic limitation is that the distance is finite, and, thus, no sampling is achieved below and above these limits. To gain information about and or the phase (see (8)), needs to be sampled on both sides of the spatial origin, so that the asymmetric feature in the interferogram is captured. Hence, in the following and, thus, are set to be negative, and the double-sided region is identified for the locations . This diminishes the maximum optical path difference achievable, being positive, and, thus, the length of the single-sided domain is identified by the relation . Because scales with the order of the Fourier coefficients (see Section 2.3), only a finite number of coefficients can be probed. According to Parseval’s theorem (see Section 2.4), information about the total length is missing. Furthermore, the Gibbs phenomenon, that is, a ringing, is present when and are inferred. Hence, should be maximal, so that as many as possible coefficients can be probed to decrease the loss of information. However, a trade-off between lengths of the single-sided and double-side domains is inevitable, depending on the level of the asymmetric imperfection.

2.5.3. Noise Contribution

Since any measurement has a noise contribution, the noisy data value can be written as , and the actual interferometric data is expressed by the vector . As spectral quantities are investigated, photons are involved in the measuring principle, and, hence, a part of has a Poissonian origin. However, the diagnostic under investigation later probes broadband spectra in the microwave and far-infrared range, and, thus, a large number of photons are present. Hence, the central limit theorem suggests that is a sample of a normal distribution with vanishing mean and a certain variance given by the squared noise level . In any case, dedicated diagnostics tests are mandatory to characterise for a given interferometer.

2.5.4. Basic Model

The combination of the relation (10) with the interferometric data, being noisy and sampled in a finite spatial domain, gives the most basic model. Formally, this model is stated here bywith the additional information: the zero-path difference location , the spatial limits are restricted by , and , the sampling increment is constant and known, and have the same spectral limits which obey and , and in case of Gaussian noise . To be precise, the unknowns of the model are , , , and , and .

The basic model is a starting point and must be amended by diagnostic imperfections and specifics to the interferometer design type.

2.6. Inferring Spectra by Standard Analysis Techniques

To infer the spectral quantities and from an interferometric data set , the standard techniques rely on a noiseless model and follow a hierarchical ansatz. After making assumptions on the spectral limits and , the zero-path difference location is estimated. The next step evaluates a phase which is a measure of the ratio , relying on the data located in the double-sided region. Given the model, the spectral limits, the spatial origin, and the phase, and are estimated up to the Nyquist frequency from the whole data set. To reduce the Gibbs phenomenon on the inferred spectral quantities, window functions are multiplied to the interferometric data. In the following, the weak points of the standard analysis techniques are described.

2.6.1. Noiseless Model

The model used and stated by (18) lacks the noise contribution by definition, and, hence, treats the noisy data as being not noisy. Strictly speaking, the model is not applicable for the problem at hand. Only by repeating the measurement sufficiently often so that the total noise contribution given by the vector becomes small, implying , the model applies. Because the integration time remains finite or only a single measurement is possible, the influence of the noise on the inferred quantities cannot be derived from the noiseless model.

2.6.2. Hierarchical Ansatz

Several steps are carried out to deduce the quantities of main interest and . Each step relies on model assumptions, which are not questioned or tested in any way, and results of previous steps, which carry an unstated uncertainty. This hierarchical ansatz lacks the uncertainty propagation onto and entirely.

2.6.3. Spectral Limits: Nyquist Assumptions

Two fundamental assumptions, called Nyquist assumptions in the following, are made by setting the spectral limits to 0 and the Nyquist frequency (see Section 2.5.2). Hence, the chosen spatial sampling would determine the bandwidth of the spectrum which is a misconception. Furthermore, the associated Fourier coefficients are located apart via their basis functions in the spatial domain, and the maximum order probed is artificially blown up to due to the Nyquist assumptions (see Section 2.3). However, if the functions and are finite in the spectral domain with limits and , then the embedding of the smaller domain into the domain leads to the reduction of the information content per Fourier coefficient as discussed in Section 2.3.2. Figure 2 compares the narrower spatial basis functions for the Nyquist assumptions with μm, that is, = 0 GHz and = 3747.4 GHz ( = 1873.7 GHz, ) with the wider basis function for the absolute term for = 0 GHz and = 1000 GHz ( GHz, = 1000 GHz).

The uncertainty of a Fourier coefficient, relying on the Nyquist assumptions, scales like (noise level/maximum bandwidth) which follows from the linear uncertainty propagation for (13). But for the band-limited case, the uncertainty would scale like , where the square root term states that more than one data point is related to one coefficient. Hence, if a band-limitation exists but is not taken into account, then the uncertainty is maximised on the inferred spectral quantities.

2.6.4. Estimation of Spatial Origin

The spatial origin or zero-path difference is most likely missed by the spatial sampling. One of the standard approaches to estimate fits a parabola to the main interferogram peak without any information about the even and odd spectra itself. However, as one can see from (13), the basis functions for the even and odd absolute terms ( and ) are of leading order close to the spatial origin. Hence, information about the zeroth-order coefficients and the spectral limits should be at hand for the estimation of .

With available, the double-sided and single-sided regions are identified. Though, a systematically affected estimate of the origin causes an additional asymmetry in the interferometric data which would result in an increase of and a decrease in which is usually interpreted as a phase ramp feature. Hence, the origin should be determined with the criterion that it minimises the odd spectral function.

2.6.5. Windowing

Having only a finite amount of Fourier coefficients probed causes the Gibbs phenomenon to appear for the spectral quantities inferred. To reduce this ringing feature, window functions are applied in the spatial domain to bring the interferometric data smoothly to zero towards the sampling limits. More precise, probed Fourier coefficients of higher orders are damped out, and a window function corresponds to a certain convolution function in the spectral domain. Hence, a weighted averaging of the spectral quantities is carried out which reduces the ringing. This approach can give a good global approximation of and for regions with no significant gradients. However, the damping of Fourier coefficients worsens the convergence of the inferred quantities in regions with considerable gradients.

Implicitly, the application of window functions excludes the investigation of the uncertainty on and introduced by nonprobed Fourier coefficients. Hence, the requirement of square-integrability of the spectral functions is not taken into account.

3. Bayesian Formalism

3.1. Bayes’ Theorem

The joint probability density function (pdf) captures the chance that the outcome , let us say a data value or set, and the outcome , a single model parameter or a set, are realised simultaneously.

The product ruleintroduces the conditional probabilities for finding the outcome , if the outcome were true and vice versa. By the theorem of Bayesone conditional probability can be expressed by the other, when the marginal distributions and are known. Hence, Bayes’ theorem captures the information/knowledge gained for when a certain outcome for has manifested. For the pdfs occurring in Bayes’ theorem, common names are used, that is, the posterior , the likelihood , the evidence , and the prior . The link or functional dependence enters in the likelihood which takes into account known uncertainties like, for example, measurement noise. Any knowledge about before new data is available can be found in the prior .

Bayes’ rule can be extended tointroducing a set of hyperparameter which enters per definition solely in the prior . The additional pdf is called hyperprior which allocates trust in . Apart from having the posterior for the parameters , the marginalisation with respect to reveals the posteriorfor which measures the plausibility of an outcome of the hyperparameter given the data. Since does not depend on , the most likely hyperparameter set is identified by the maximum of , assuming is uniform.

3.2. Formalism for Linear, Nonlinear, and Hyperparameter Problem for Gaussian Noise
3.2.1. Multivariate Normal

Let the joint pdf for the random vector () be a multivariate normal with mean () and covariance matrix (); then the pdf becomeswith the determinant .

3.2.2. Model for Linear Problem

If the dependency between the data and the parameters of interest is linear, and the likelihood and the prior can be expressed by multivariate normals, then the evaluation of the posterior is analytically straightforward. Such a model is the starting point for investigating a more complex model which includes parameters with a nonlinear mapping to the data domain and/or hyperparameters.

(a) Gaussian Likelihood. The data may be represented by the vector (). The parameters of interest () map linearly to the data domain likewhere the dimensional matrix M encodes the linear mapping, and captures the random noise contribution. When the data is acquired independently, and the noise is independent for each datum and follows a Gaussian with vanishing mean and standard deviation (noise level), then the Gaussian likelihoodcan be found with the covariance matrix .

(b) Gaussian Prior. The prior information about may be expressed by the multivariate normalwith the prior mean and the prior covariance .

(c) Gaussian Posterior and Evidence. Formally, Bayes’ theorem states the posterior byAfter some algebra, one can show that the posterioris a multivariate normal with posterior meanand covariancewhich are both analytically obtained. Furthermore, the evidence readswhere the first part depends explicitly on the measured data, and the second part, being dimensionless, incorporates the ratio dependent on the means and covariances of the prior and posterior.

3.2.3. Model for Linear, Nonlinear, and Hyperparameter Problem

The linear model is amended by hyperparameters, entering in some way in the prior, and parameters with a nonlinear connection to the data domain. Such a model is then applicable in the field of Fourier spectroscopy.

(a) Gaussian Likelihood. The linear mapping of the parameters to the data domain, as stated by (24), should remain valid. However, the mapping itself may depend on the parameters in a nonlinear way, so that . This leaves the Gaussian likelihood in (25) formally unchanged but is symbolically stated as .

(b) Priors. The Gaussian prior for should be given by , where the prior mean and covariance depend on some of the hyperparameters . Similarly, a prior follows for the nonlinear parameters. Finally, the hyperparameters have an assigned prior .

(c) Posteriors and Evidence. According to Bayes’ theorem, one can write the joint posterior likeand the conditional amplitude posterior for becomes a multivariate normalThus, both, the conditional posterior mean and covariance evaluated by (29) and (30), depend on the nonlinear parameters and hyperparameters. After the trivial marginalisation with respect to , the joint posterior for and remains. By expressing the posterior, named settings posterior in the following, likethe evidence is identified withNote, that the dimensionless constant , and, thus, the evidence depend on the chosen model, including likelihood and priors. Hence, is of importance, when the model is even further abstracted or compared with alternative models.

(d) Role of Settings Posterior. The optimisation, that is, the finding of the maximum of the settings posterior , can be interpreted as an implementation of the Ockham’s razor principle and/or as a regularisation procedure. This is essential when the number of parameters exceeds the number of data points.

Unfortunately, a general analytical expression is not available for this posterior, and, thus, it needs to be investigated numerically for the problem at hand. In order to do so, the quantityis of interest, because it is numerically accessible. In case, has a well distinguishable global maximum, can be approximated by a multivariate normal which is estimated by evaluating the Hessian matrix. Thus, one finds the posterior means and with the associated posterior covariance. This allows an approximate marginalisation with respect to and/or . This can be understood as a propagation of the posterior uncertainties in and to the marginalised posterior for the parameters of interest.

(e) Simplifications. For the remainder of this paper, some simplifications are made which modify (34), (35), and (36) accordingly. The prior mean is set to 0, and the priors and are chosen to be uniform. Then, one can set which modifies (37) to

4. Brownian Bridge Covariance

The continuous even and odd spectral functions to be inferred can be modelled each by a Gaussian process [11]. Thereby, the Brownian bridge process is a good starting point, because it exploits a fundamental condition to prevent aliasing for Fourier spectroscopy applications. This condition states that the spectrum and, thus, and must vanish at the spectral origin and at an upper limit which is smaller than the Nyquist frequency (see Section 2.6.3). However, this information is usually not taken into account any further in the analysis. On the contrary, a Brownian bridge and its associated covariance function fulfil the boundary conditions for any lower and upper limit. Hence, the covariance can be used in the Gaussian prior for and . In addition, this process has only one scaling hyperparameter which makes it attractive from data analysis point of view. This scaling can be estimated as well from the Fourier coefficients probed. In fact, this reveals information about the nonprobed coefficients and gives an additional uncertainty on and . After presenting some properties of the Brownian bridge covariance, it is used as prior covariance in the example application (see Section 5).

4.1. Standard Definition

The Brownian bridge is a continuous stochastic process for an interval, say from 0 to . This bridge is constructed by tying-down a Brownian motion process to 0 at the end of the interval in question. Furthermore, the tie-down at the beginning of the interval is inherited from the Brownian motion process. The covariance function for the bridge is defined in standard literature byfor and .

4.2. Adapted Definition

The even and odd spectral functions, being finite for the interval , demand some adaption of the standard covariance expression (39) when modelled by a Brownian bridge process. To keep the same properties on spectral scale, a shift by and the interval length need to be set. Thus, one gets

with the unit = Hz. With the normalisationthe modified covariance becomeswhere = Hz−2 has now the proper unit with respect to the spectral scale.

The parameters and of unit are introduced which are defined each as a scaling factor for the associated process. With these scalings the covariances are obtained, so that the units = V2/Hz2 match.

4.3. Covariance for Fourier Coefficients

The Brownian bridge covariance function for the spectral domain can be studied in the domain of the Fourier coefficients via the coordinate transform stated in (11) [11]. Compactly written, one finds the infinite-dimensional covariance matrix for the Fourier coefficients analytically by for all , and similarly follows. The only finite off-diagonal elements occur for the absolute term in connection with the higher order terms for the even coefficients captured by the infinite-dimensional row and column vectors and , respectively. This is caused by the condition that vanishes at the spectral boundaries where the sine vanishes intrinsically for any but the cosine takes values either 1 or −1 for even and odd orders, respectively. Hence, the covariance imposes the boundary condition.

For orders greater than 0, has no off-diagonal elements due to the Kronecker delta , meaning that these coefficients are independent on each other. Furthermore, for the infinite-dimensional matrices holds, and the amplitudes for even and odd coefficients drop equally with the square of the order.

4.4. Square-Integrable Property

According to Parseval’s theorem (see (16)), is evaluated by summing the squares of the Fourier coefficients. Because the entries of the main diagonal in the covariances and drop with the order squared, the Brownian bridge process ensures square-integrability of and as long as the scalings and remain finite.

4.5. Signal Envelope

For the even process, the signal level can be estimated by the envelope in the data domain. Starting point is the square root of the main diagonal of . Since the argument of (see (14)) localises the even and odd Fourier coefficient at a fixed in the same data domain, the even and odd contributions of must be added for . As can be seen by (13), the mapping of the absolute term to the data domain includes already the factor 2. In addition, the mapping comprises the bandwidth . In total, one finds the envelope asIn the above equation, factor 2 in front of was chosen, so that captures most of the signal. An approximation might be convenient, because 1.

For the envelope of the odd process, the same reasoning can be applied with one modification. The mapping demands that the contribution of the absolute term at the spatial origin vanishes (see (13)). Hence, one findsBoth envelops drop with , and, thus, most of the signal associated with each process would in the data domain.

5. Example Application

5.1. Formulation of Model
5.1.1. Martin-Puplett Interferometer at JET

The Martin-Puplett interferometer diagnostic [12] at the fusion device JET (Culham, UK) probes the spectrum emitted by a broadband source and performs the Fourier transform. The interferogram data is acquired in terms of Volts dependent on the optical path difference . However, two different sources are probed for 20 minutes subsequently to remove a class of diagnostic imperfections not treated here any further. By subtraction of the corresponding two interferograms, the data becomes available in form of the difference interferogram acquired at the spatial grid node . Then, the abstract model for the Martin-Puplett interferometer is stated by using the total amplification of the detection system. Furthermore, the offset marks a diagnostic imperfection which varies with . The Gaussian noise contributes to each data sample described by . The unknown quantities in the diagnostic model are the spatial grid , the lower and upper spectral boundaries and of the Fourier transform integral, the even and odd functions and dependent on frequency , and the offset.

5.1.2. Interferometric Data

The data set , that is, the difference interferogram consists of = 788 values (see Figure 3(a)). Merely for graphical presentation a certain is chosen derived from the standard approach (see Section 5.1.3). Globally, the data shows an upward trend with respect to the zero baseline.

Figure 3: (a) Difference interferogram (black) versus optical path difference . The data set includes = 788 values, and the noise level for each reads = 132.29 μV. The values of the spatial grid are determined by the standard model with , the constant increment = 40 μm, and the origin . This origin is estimated as −5.118 mm by a fit of a quadratic polynomial to the maximum and its neighbouring values to either side. Clearly, deviates from the baseline (red). (b) Difference quotient for different subsets of . Since the zig-zag-pattern in the point-to-point variation appears only for the first case (black), the standard model = const. seems to be inadequate. Instead, a model seems appropriate for which = holds.

The components are measured independently on each other, and the noise level for each is captured by = 132.29 μV, and, thus, the variance of the whole data vector is stated by the matrix .

5.1.3. Optical Path Difference

The diagnostic is set up, so that the sampling of the interferogram is triggered ideally when the optical path difference has changed by the increment = 40 μm. Hence, the standard model is obvious with , and the zero-path difference is a free parameter. However, this model is not accurate. Applying the standard approach which fits a second-order polynomial to the maximum and its two nearby values,  mm, is inferred for the data set shown in Figure 3(a). Furthermore, the difference quotient evaluated by is presented versus the optical path difference in Figure 3(b). The point-to-point variation of the quotient has a zig-zag pattern which implies that the assumption = const. is incorrect. Indeed, each of two the difference quotients and is smoother. Hence, the modelseems more appropriate, making use of two free parameters: the zero-path difference and a shift for every other grid value. The priors and are set to be uniform.

5.1.4. Offset

The upward trend of (see Figure 3(a)) is modelled by the offset . Here, a second-order polynomial is chosen to capture the offset by with the free parameters , , and being summarised by the vector . The corresponding mapping to the data domain can be expressed aswith the matrix . The joint prior is expressed by the factorisable multivariate normal distribution:where are considered as hyperparameters to which a uniform prior is assigned.

5.1.5. Spectral Quantities

The kernel of the spectral integration in (46) is finite only for frequencies     . Thereby, and are free parameters. The spectral domain is dicretised using the constant increment which is considered as a free parameter as well. Then, the spectral domain is represented by the set with , covering the band centred at . To be clear, since , , and are free parameters and will be inferred, the number is variable. This number determines the dimensionality of the vectors and which represent the discretised functions and . The mapping of and to the data domain in (46) is written aswith the two matricesThe joint prior for the two vectorial quantities and is factorised, and each prior is chosen as a multivariate normal distribution with vanishing mean. Since the Brownian bridge covariance (see Section 4) describes functions which vanish at the boundaries and and are square-integrable, and its signal envelope decays with the optical path difference like the interferometric data at hand, the priors are chosen bywith the two hyperparameters and for scaling. For each of these hyperparameters a uniform prior is applied.

Finally, the joint priorshould be constant, so that any combination of , , and has the same probability. Furthermore, the conditions , , and (global upper limit) must be fulfilled. For example, the upper limit of is set here to the Nyquist frequency = 3747.4 GHz.

5.2. Bayes’ Theorem

In the following, the linear parameters are summerised by the set , the nonlinear parameters by , and the hyperparameters by .

The matrix maps the parameters to the data domain. Hence, the Gaussian likelihood is written asThe prior for the full problem takes the formwith the multivariate normal priorfor , using the dimensional covariance matrixOne gets the joint posteriorwith the conditional amplitude posterior given by the multivariate normalusing the posterior covariance matrix and the mean . Furthermore, one obtains the settings posteriorwhere the constant is unknown so far.

5.3. Investigation of Posterior
5.3.1. Conditional Amplitude Posterior for Chosen Settings

To give some insight, the conditional posterior for the amplitudes is evaluated given the specific set of values for , , , , , and , , , and .

The values for  mm and  mm are chosen to form an optical path difference grid for which all subsequent nodes are separated by . Then,  mm identifies the double-sided domain (256 data points), and the single-sided domain is marked by 5.118 mm 26.362 mm (532 data points).

The spectral domain is covered from = 0 GHz to the Nyquist frequency = 3747.4 GHz (=). With the bandwidth and the optical path difference set, the conversion (see Section 2.3.1) to the order of Fourier coefficient follows. Thus, the double-sided and full spatial domain coefficients up to the orders 64 and 330 are located (see Figure 4(b)), respectively. Some examples of the spatial basis functions associated with this set of Fourier coefficients are shown in Figure 2.

Figure 4: Examination of Gaussian prior for even and odd spectra modelled by Brownian bridge process. The spectral discretisation = 5.68 GHz is used, and the lower and upper spectral limits are set to = 0 GHz and = 3747.4 GHz (Nyquist assumptions), respectively. The scalings of the processes are chosen by  V2. (a) Prior samples and for even and odd spectra. The spectra vanish towards the spectral limits. (b) Prior samples mapped to data domain to give even (odd) data-like quantity for optical path difference obtained for  mm and = 0. With increasing distance to the origin, the amplitude drops which is well bounded by the signal envelopes for the chosen bandwidth . Fourier coefficients up to order = 64 and 330 settle in the double-sided domain and single-sided domain , respectively.

To match the maximum order (330) of probed Fourier coefficient, the discretisation is set to .68 GHz which reflects the classical increment. Then, the spectral grid has = 660 elements.

The hyperparameters are set to large values like  V2 to prevent from a prior determined posterior. This unwanted feature would be present, if one would choose too small, so that the signal envelops do not include the measured data (see next paragraph).

With the chosen spectral priors, sample functions/vectors and are drawn and shown in Figure 4(a). The function values are of the order of  V/Hz decaying towards the spectral boundaries. The mapping of the spectral prior samples to the spatial domain via shows that can be of the order of some Volts close to the spatial origin (see Figure 4(b)) which exceeds the measured data by about two orders of magnitude. Furthermore, the amplitude of drops with the distance from the origin. All these data samples are well bounded by the signal envelopes obtained by which is the scaled version of the expressions (44) and (45).

For the prior covariance is determined by setting  V,  V/m and  V/m2.

(a) Posterior Mean. For the above values, the mean of the amplitude posterior is given by .41 nV, 5.75 μV/m, 2.42 μV/m2, and as presented in Figure 5(a). Clearly, peaks below 1000 GHz with the maximum of approximately 6 × 10−18 V/Hz around 400 GHz. Furthermore, is much smaller than below 1000 GHz.

Figure 5: Investigation of conditional amplitude posterior for linear parameters . The posterior is obtained for the settings  mm and = 0, = 0 GHz, = 3747.4 GHz, = 5.68 GHz,  V2,  V,  V/m, and  V/m2. (a) Spectral posterior means and one posterior sample. Below 1000 GHz, the means show a finite amplitude. Each sample deviates strongly from its corresponding mean. (b) Posterior mean and 100 samples mapped to give even, odd, and offset contributions in data domain. While in the double-sided region () most of the data described by , and give same contributions in the single-sided region (). In the single-sided region, the spread of the even (odd) contributions for the posterior samples is large but bounded by the signal envelops . (c) Residuals for (black) and one sample (red) both mapped to data domain. An almost perfect description of the noisy data is achieved by , indicating an overfitting of the data. (d) Measured data and 100 samples of amplitude posterior mapped to data domain. The individual contributions have a much wider spread (see (b)), but the total sum agrees well with and its uncertainties. This is caused by the posterior covariance of with suited off-diagonal entries.

The offset contribution to the data has an upwards trend but is small (see Figure 5(b)). The even and odd contributions behave differently in the double- and single-sided domains (see Figure 5(b)). While in the double-sided domain is much larger than , and equal each other in the single-sided domain. This indicates an underestimation (overestimation) of the Fourier coefficients above the order 64 for (). The cause of this finding relies in the choice () which becomes more clear in Section 5.3.3.

Since the quantity equals almost the data set (see residuals in Figure 5(c)), a nearly perfect match is achieved. Hence, all noise contributions are captured by the mean , and an overfitting of the data is obvious.

(b) Posterior Covariance. The posterior covariance is not characterised in detail. From the elements of the main diagonal of , one finds the values = 73.30 nV, = 12.58 μV/m, and = 446.67 μV/m2. These posterior uncertainties are quite large when compared to the values of the corresponding mean , especially for the quadratic coefficient.

For the spectral quantities, the square root of the main diagonal elements of is of the order of some  V/Hz (for one sample function drawn from the conditional amplitude posterior see Figure 5(a)). Hence, a considerable deviation from the posterior mean is possible.

(c) Posterior Samples. 100 samples are drawn from the amplitude posterior and mapped to the data domain (see Figure 5(b)). Thereby, the contributions , , and are split to investigate their interplay. has a noticeable width but remains in the vicinity of for the whole spatial domain.

For the double-sided region and form a narrow band around the corresponding posterior mean quantities and . On the contrary, for the single-sided region and increase drastically but are restricted by and . The deviation is large between the individual and the posterior mean quantities; though, the sum remains well inside the 2- band surrounding (see Figures 5(c) and 5(d)). Hence, the posterior covariance captures properly the correlations.

With the specific values of the settings one gets the number