# Fourier Spectroscopy: A Bayesian Way

Authors: Stefan Schmuck; Jakob Svensson
Journal: International Journal of Spectroscopy
 (2017)

---

## 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.


---

## Body

## 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 [4–6] 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 Section2, 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 frequencyf and time t, or the wavenumber σ and a spatial coordinate x. For the moment, the latter pair is used to state the basic operations of Fourier transformation. Afterwards, the wavenumber is replaced via σ=f/c for convenience.The real-valued continuous functionsV(x) and Sσ(σ) form a Fourier transformation pair stated by(1)Vx=∫-∞∞Sσcos⁡2πσx+sin⁡2πσxdσ.The inverse operation reads(2)Sσ′=∫-∞∞Vcos⁡2πσ′x+sin⁡2πσ′xdx.Note that so far σ′≠σ holds. To find a relation, one inserts (1) in the above expression. After applying trigonometric identities, the spatial integral becomes(3)∫-∞∞cos⁡2πσ-σ′x+sin⁡2πσ+σ′xdx=δσ-σ′,because the sinusoidal contribution vanishes. Hence,(4)Sσ′=∫-∞∞Sσδσ-σ′dσ,Sσ′=Sσfollow, because the two wavenumber coordinates equal each other. In fact, the delta distribution occurs, because Sσ is a distribution.Replacingσ with f/c gives in the spectral domain the quantity S(f)=Sσ/c. This gives(5)V=∫-∞∞Scos⁡2πfcx+sin⁡2πfcxdfand rescales the inverse operation like(6)S=1c∫-∞∞Vxcos⁡2πfcx+sin⁡2πfcxdxto match the units [S]=[V]/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 SE=S(f)+S(-f)/2 of S, the sine transform is linked only to the odd part SO=S(f)-S(-f)/2. Thus, an alternative formulation reads(7)V=2∫0∞SEcos⁡2πfcx+SOsin⁡2πfcxdf,and it becomes clear that the even and the odd parts of V and S are connected.Another representation uses the amplitude|S| and the phase α by setting(8)S=2SE2+SO2,tan⁡α=-SOSEwhich gives(9)V=∫0∞Scos⁡2πfcx+αdf.Because the model representations (5) and (7) are linear in S, SE, and SO, 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 thatSE and SO are finite for the range from fL and fU with the bandwidth ΔF=fU-fL and the centre frequency fC=(fL+fU)/2, (7) becomes(10)Vx=2∫fLfUSEcos⁡2πfcx+SOsin⁡2πfcxdf.Here, 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 fC and ΔF the relation ΔF≤2fC 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 orderk (∈N). These coefficients are defined here by the integrals(11)AEO,0AEO,kBEO,k=1ΔF∫fLfUSEO12cos⁡2πkf-fCΔF2sin⁡2πkf-fCΔFdfwhich carry the unit as SE and SO. The coefficients AE(O),0 label the mean values of SE and SO in the spectral domain covered. Then, one can replace the even and odd functions with(12)SEO=AEO,0+∑k=1∞AEO,kcos⁡2πkf-fCΔF+∑k=1∞BEO,ksin⁡2πkf-fCΔFwhich allows performing the spectral integration in (10) analytically. Since the result(13)Vx=2ΔFAE,0cos⁡2πfCcxsincπΔFcx+2ΔFAO,0sin⁡2πfCcxsincπΔFcx+ΔF∑k=1∞AE,kcos⁡2πfCcxsinck-+sinck+-ΔF∑k=1∞BE,ksin⁡2πfCcxsinck--sinck++ΔF∑k=1∞AO,ksin⁡2πfCcxsinck-+sinck++ΔF∑k=1∞BO,kcos⁡2πfCcxsinck--sinck+with(14)sinck±=sincπk±xcΔF.follows, 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 sinck± functions which depend on the order k, the spatial coordinate, and the bandwidth. The latter quantity determines the spatial width of sinc±. Furthermore, the localisation is permitted at xk=±kc/ΔF, where a coefficient for a given k 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 AE(O),0 for which both sinc functions coincide.The second contribution causes a modulation ofsinck± and is given by a sine/cosine with the centre frequency and spatial coordinate in the argument. This dependency makes the basis function for AO,0 vanish at the spatial origin.With respect to the spatial origin, the transformed basis functions of the coefficients forSE and SO are symmetric and antisymmetric, respectively.Some basis functions in the spatial domain are shown in Figures1(a) and 1(b) for fC = 500 GHz and ΔF = 1000 GHz.Figure 1
Basis functions in spatial domain of Fourier coefficients (of orderk = 0, 1, and 12) AE,k, BE,k for even spectral function (a) and AO,k, BO,k for odd spectral function (b). The spectral functions are finite for the interval from 0 GHz to 1000 GHz for which the centre frequency fC = 500 GHz and bandwidth ΔF = 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 x, ΔF, and k. Furthermore, the sinc functions are modulated by sine or cosine with fC and x in the argument. The sinc functions for fixed k become unity at xk=±kc/ΔF, obviously increasing in k.
(a)
(b)


#### 2.3.2. Embedding into Larger Spectral Domain

The functionsSE and SO may be finite in the spectral domain with limits fL and fU or centre frequency fC and bandwidth ΔF. Embedding this domain in a larger one with limits fL′<fL and fU′>fU (ΔF<ΔF′), another set of Fourier coefficients AE(O),k′′ and BE(O),k′′ is obtained with associated basis functions for the spectral domain. Without going into more detail here, these coefficients can be evaluated from AE(O),k and BE(O),k and the scalar products of the basis functions labeled with k and k′. For instance, one finds for the ratio of the means AE(O),0′/AE(O),0=ΔF/ΔF′<1. Basically, AE(O),k and BE(O),k maximise the information per coefficient when fL and fU 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 AE(O),0 is mapped to an infinite number of coefficients AE(O),k′′ and BE(O),k′′ which are mandatory to capture both discontinuities.In the spatial domain, the basis functions forAE(O),k′′ and BE(O),k′′ behave differently than the ones for AE(O),k and BE(O),k. Important to mention is the effect which the larger bandwidth has; sinck′± are spatially narrower than sinck±. In addition, the number of coefficients per spatial domain increases which is expressed by k′=kΔF′/ΔF.ForfC′ = 1873.7 GHz and ΔF′=2fC′ Figure 2 shows some basis functions (k′=0, 1, and 2) in the spatial domain. Indeed, the basis function for AE,0 with fC = 500 GHz and ΔF = 1000 GHz (see Figure 1(a)) is broader.Figure 2
Basis functions in spatial domain of Fourier coefficientsAE,k′′, BE,k′′ with order k′=0, 1, and 2 for even spectral function, being finite in spectral domain with centre frequency fC′ = 1873.7 GHz and bandwidth ΔF′=2fC′. For this domain, the sinc functions, contributing to the basis functions, are narrower than for fC = 500 GHz and ΔF = 1000 GHz (see Figure 1(a)). As an example, the basis function for AE,0 (dashed-black) is given. Since ΔF′≈4ΔF holds, approximately 4 times more Fourier coefficients locate inside a given spatial domain. ΔF′ was chosen to equal the Nyquist frequency fNy=c/(2Δx), having set the spatial sampling increment Δx=40 μ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 lengthlS of the band-limited function S is stated by(15)lS2=c∫fLfUS2df+c∫-fU-fLS2df=2c∫fLfUSE2+SO2df,because only the term SESO remains odd and cancels by the integration. Furthermore, the scaling by the factor c appears which originates in df/dσ=c. Replacing SE(O) with the expression (12) and exploiting that the basis functions are perpendicular in the spectral domain leaves(16)lS2=cΔF∑i=E,O2Ai,02+∑k=1∞Ai,k2+Bi,k2.Thus, the length is given by the sum of the Fourier coefficients squared. According to Parseval’s theorem,(17)lS2=∫-∞∞V2xdx,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 k1≠k2, and the spatial integral yields c/(2ΔF) for k1=k2=0 and c/ΔF for k1=k2>0. Analytically, this is hard to prove; however, this was numerically investigated and is considered to be valid.


#### 2.4.1. Square-Integrable Functions

The functionS is said to be square-integrable, when the condition lS2<∞ holds. Furthermore, if S is square-integrable, then the Fourier series representations in (12) converges towards SE and SO almost everywhere in the spectral domain as the order k grows [10]. Hence, the requirement on S 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 functionVx 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,V(x) is sampled at a finite set of optical path difference locations xi with i∈[1,ND], and ND marks the number of sample points. Usually, the sampling with constant increment Δx=xi+1-xi 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 xi might be unknown. If so, it is mandatory to introduce the zero-path difference x0 which is in the following set that x1=x0 holds.The finite spatial sampling leavesV(x) undetermined between the sampling nodes and outside the limits x1 and xND. Assuming Δx=const. holds, the Nyquist theorem states that the maximum frequency accessible is given by the Nyquist frequency fNy=c/(2Δx). Hence, for the spectral quantities SE and SO to be inferred, a maximum for the upper limit fU,Max=fNy follows from sampling theory. To prevent aliasing, Δx needs to be chosen small enough so that SE and SO vanish below fNy. In case, fU 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 fNy contributes to V(xi). 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 x0 has no profound impact.A diagnostic limitation is that the distancexND-x1 is finite, and, thus, no sampling is achieved below and above these limits. To gain information about SE and SO or the phase α (see (8)), V(x) 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 x1 and, thus, x0 are set to be negative, and the double-sided region is identified for the locations |xi|≤|x0|. This diminishes the maximum optical path difference xND achievable, being positive, and, thus, the length of the single-sided domain is identified by the relation |x0|<xi≤xND. Because xND scales with the order k 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 SE and SO are inferred. Hence, xND 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 asD(xi)=V(xi)+ϵi, and the actual interferometric data is expressed by the vector D→=V→+ϵ→. As spectral quantities are investigated, photons are involved in the measuring principle, and, hence, a part of ϵi 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 ϵi is a sample of a normal distribution with vanishing mean and a certain variance given by the squared noise level σD2. In any case, dedicated diagnostics tests are mandatory to characterise ϵi 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 by(18)Dxi=2∫fLfUSEcos⁡2πfcxi+SOsin⁡2πfcxidf+ϵiwith the additional information: the zero-path difference location x0<0, the spatial limits are restricted by x1=x0, xND>0 and xi<∞, the sampling increment Δx=xi+1-xi is constant and known, SE and SO have the same spectral limits which obey fL≥0 and fU≤fNy=c/(2Δx), and in case of Gaussian noise p(ϵi)=N(0,σD2). To be precise, the unknowns of the model are x0, fL, fU, SE(f) and SO(f), and σD.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 quantitiesSE and SO from an interferometric data set D→, the standard techniques rely on a noiseless model and follow a hierarchical ansatz. After making assumptions on the spectral limits fL and fU, the zero-path difference location x0 is estimated. The next step evaluates a phase α which is a measure of the ratio SO/SE, relying on the data located in the double-sided region. Given the model, the spectral limits, the spatial origin, and the phase, SE and SO 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 D→ 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 D→=V→, 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 interestSE and SO. 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 SE and SO 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 frequencyfNy=c/(2Δx) (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 Δx apart via their basis functions in the spatial domain, and the maximum order probed is artificially blown up to xND/(2Δx) due to the Nyquist assumptions (see Section 2.3). However, if the functions SE and SO are finite in the spectral domain with limits fL>0 and fU<fNy, then the embedding of the smaller domain [fL,fU] into the domain [0,fNy] 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 Δx=40 μm, that is, fL′ = 0 GHz and fU′ = 3747.4 GHz (fC′ = 1873.7 GHz, ΔF′=2fC′) with the wider basis function for the absolute term for fL = 0 GHz and fU = 1000 GHz (fC=500 GHz, ΔF = 1000 GHz).The uncertainty of a Fourier coefficient, relying on the Nyquist assumptions, scales likeσD/fNy (noise level/maximum bandwidth) which follows from the linear uncertainty propagation for (13). But for the band-limited case, the uncertainty would scale like (ΔF/fNy)1/2σD/ΔF, 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 differencex0 is most likely missed by the spatial sampling. One of the standard approaches to estimate x0 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 (AE,0 and AO,0) 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 x0.Withx0 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 SO and a decrease in SE 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 ofSE and SO 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 onSE and SO introduced by nonprobed Fourier coefficients. Hence, the requirement of square-integrability of the spectral functions is not taken into account.


## 2.1. General Definitions

Commonly, the complementary coordinates used for Fourier transformations are the frequencyf and time t, or the wavenumber σ and a spatial coordinate x. For the moment, the latter pair is used to state the basic operations of Fourier transformation. Afterwards, the wavenumber is replaced via σ=f/c for convenience.The real-valued continuous functionsV(x) and Sσ(σ) form a Fourier transformation pair stated by(1)Vx=∫-∞∞Sσcos⁡2πσx+sin⁡2πσxdσ.The inverse operation reads(2)Sσ′=∫-∞∞Vcos⁡2πσ′x+sin⁡2πσ′xdx.Note that so far σ′≠σ holds. To find a relation, one inserts (1) in the above expression. After applying trigonometric identities, the spatial integral becomes(3)∫-∞∞cos⁡2πσ-σ′x+sin⁡2πσ+σ′xdx=δσ-σ′,because the sinusoidal contribution vanishes. Hence,(4)Sσ′=∫-∞∞Sσδσ-σ′dσ,Sσ′=Sσfollow, because the two wavenumber coordinates equal each other. In fact, the delta distribution occurs, because Sσ is a distribution.Replacingσ with f/c gives in the spectral domain the quantity S(f)=Sσ/c. This gives(5)V=∫-∞∞Scos⁡2πfcx+sin⁡2πfcxdfand rescales the inverse operation like(6)S=1c∫-∞∞Vxcos⁡2πfcx+sin⁡2πfcxdxto match the units [S]=[V]/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 SE=S(f)+S(-f)/2 of S, the sine transform is linked only to the odd part SO=S(f)-S(-f)/2. Thus, an alternative formulation reads(7)V=2∫0∞SEcos⁡2πfcx+SOsin⁡2πfcxdf,and it becomes clear that the even and the odd parts of V and S are connected.Another representation uses the amplitude|S| and the phase α by setting(8)S=2SE2+SO2,tan⁡α=-SOSEwhich gives(9)V=∫0∞Scos⁡2πfcx+αdf.Because the model representations (5) and (7) are linear in S, SE, and SO, 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 thatSE and SO are finite for the range from fL and fU with the bandwidth ΔF=fU-fL and the centre frequency fC=(fL+fU)/2, (7) becomes(10)Vx=2∫fLfUSEcos⁡2πfcx+SOsin⁡2πfcxdf.Here, 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 fC and ΔF the relation ΔF≤2fC 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 orderk (∈N). These coefficients are defined here by the integrals(11)AEO,0AEO,kBEO,k=1ΔF∫fLfUSEO12cos⁡2πkf-fCΔF2sin⁡2πkf-fCΔFdfwhich carry the unit as SE and SO. The coefficients AE(O),0 label the mean values of SE and SO in the spectral domain covered. Then, one can replace the even and odd functions with(12)SEO=AEO,0+∑k=1∞AEO,kcos⁡2πkf-fCΔF+∑k=1∞BEO,ksin⁡2πkf-fCΔFwhich allows performing the spectral integration in (10) analytically. Since the result(13)Vx=2ΔFAE,0cos⁡2πfCcxsincπΔFcx+2ΔFAO,0sin⁡2πfCcxsincπΔFcx+ΔF∑k=1∞AE,kcos⁡2πfCcxsinck-+sinck+-ΔF∑k=1∞BE,ksin⁡2πfCcxsinck--sinck++ΔF∑k=1∞AO,ksin⁡2πfCcxsinck-+sinck++ΔF∑k=1∞BO,kcos⁡2πfCcxsinck--sinck+with(14)sinck±=sincπk±xcΔF.follows, 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 sinck± functions which depend on the order k, the spatial coordinate, and the bandwidth. The latter quantity determines the spatial width of sinc±. Furthermore, the localisation is permitted at xk=±kc/ΔF, where a coefficient for a given k 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 AE(O),0 for which both sinc functions coincide.The second contribution causes a modulation ofsinck± and is given by a sine/cosine with the centre frequency and spatial coordinate in the argument. This dependency makes the basis function for AO,0 vanish at the spatial origin.With respect to the spatial origin, the transformed basis functions of the coefficients forSE and SO are symmetric and antisymmetric, respectively.Some basis functions in the spatial domain are shown in Figures1(a) and 1(b) for fC = 500 GHz and ΔF = 1000 GHz.Figure 1
Basis functions in spatial domain of Fourier coefficients (of orderk = 0, 1, and 12) AE,k, BE,k for even spectral function (a) and AO,k, BO,k for odd spectral function (b). The spectral functions are finite for the interval from 0 GHz to 1000 GHz for which the centre frequency fC = 500 GHz and bandwidth ΔF = 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 x, ΔF, and k. Furthermore, the sinc functions are modulated by sine or cosine with fC and x in the argument. The sinc functions for fixed k become unity at xk=±kc/ΔF, obviously increasing in k.
(a)
(b)


### 2.3.2. Embedding into Larger Spectral Domain

The functionsSE and SO may be finite in the spectral domain with limits fL and fU or centre frequency fC and bandwidth ΔF. Embedding this domain in a larger one with limits fL′<fL and fU′>fU (ΔF<ΔF′), another set of Fourier coefficients AE(O),k′′ and BE(O),k′′ is obtained with associated basis functions for the spectral domain. Without going into more detail here, these coefficients can be evaluated from AE(O),k and BE(O),k and the scalar products of the basis functions labeled with k and k′. For instance, one finds for the ratio of the means AE(O),0′/AE(O),0=ΔF/ΔF′<1. Basically, AE(O),k and BE(O),k maximise the information per coefficient when fL and fU 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 AE(O),0 is mapped to an infinite number of coefficients AE(O),k′′ and BE(O),k′′ which are mandatory to capture both discontinuities.In the spatial domain, the basis functions forAE(O),k′′ and BE(O),k′′ behave differently than the ones for AE(O),k and BE(O),k. Important to mention is the effect which the larger bandwidth has; sinck′± are spatially narrower than sinck±. In addition, the number of coefficients per spatial domain increases which is expressed by k′=kΔF′/ΔF.ForfC′ = 1873.7 GHz and ΔF′=2fC′ Figure 2 shows some basis functions (k′=0, 1, and 2) in the spatial domain. Indeed, the basis function for AE,0 with fC = 500 GHz and ΔF = 1000 GHz (see Figure 1(a)) is broader.Figure 2
Basis functions in spatial domain of Fourier coefficientsAE,k′′, BE,k′′ with order k′=0, 1, and 2 for even spectral function, being finite in spectral domain with centre frequency fC′ = 1873.7 GHz and bandwidth ΔF′=2fC′. For this domain, the sinc functions, contributing to the basis functions, are narrower than for fC = 500 GHz and ΔF = 1000 GHz (see Figure 1(a)). As an example, the basis function for AE,0 (dashed-black) is given. Since ΔF′≈4ΔF holds, approximately 4 times more Fourier coefficients locate inside a given spatial domain. ΔF′ was chosen to equal the Nyquist frequency fNy=c/(2Δx), having set the spatial sampling increment Δx=40 μm (black dots).


## 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 orderk (∈N). These coefficients are defined here by the integrals(11)AEO,0AEO,kBEO,k=1ΔF∫fLfUSEO12cos⁡2πkf-fCΔF2sin⁡2πkf-fCΔFdfwhich carry the unit as SE and SO. The coefficients AE(O),0 label the mean values of SE and SO in the spectral domain covered. Then, one can replace the even and odd functions with(12)SEO=AEO,0+∑k=1∞AEO,kcos⁡2πkf-fCΔF+∑k=1∞BEO,ksin⁡2πkf-fCΔFwhich allows performing the spectral integration in (10) analytically. Since the result(13)Vx=2ΔFAE,0cos⁡2πfCcxsincπΔFcx+2ΔFAO,0sin⁡2πfCcxsincπΔFcx+ΔF∑k=1∞AE,kcos⁡2πfCcxsinck-+sinck+-ΔF∑k=1∞BE,ksin⁡2πfCcxsinck--sinck++ΔF∑k=1∞AO,ksin⁡2πfCcxsinck-+sinck++ΔF∑k=1∞BO,kcos⁡2πfCcxsinck--sinck+with(14)sinck±=sincπk±xcΔF.follows, 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 sinck± functions which depend on the order k, the spatial coordinate, and the bandwidth. The latter quantity determines the spatial width of sinc±. Furthermore, the localisation is permitted at xk=±kc/ΔF, where a coefficient for a given k 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 AE(O),0 for which both sinc functions coincide.The second contribution causes a modulation ofsinck± and is given by a sine/cosine with the centre frequency and spatial coordinate in the argument. This dependency makes the basis function for AO,0 vanish at the spatial origin.With respect to the spatial origin, the transformed basis functions of the coefficients forSE and SO are symmetric and antisymmetric, respectively.Some basis functions in the spatial domain are shown in Figures1(a) and 1(b) for fC = 500 GHz and ΔF = 1000 GHz.Figure 1
Basis functions in spatial domain of Fourier coefficients (of orderk = 0, 1, and 12) AE,k, BE,k for even spectral function (a) and AO,k, BO,k for odd spectral function (b). The spectral functions are finite for the interval from 0 GHz to 1000 GHz for which the centre frequency fC = 500 GHz and bandwidth ΔF = 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 x, ΔF, and k. Furthermore, the sinc functions are modulated by sine or cosine with fC and x in the argument. The sinc functions for fixed k become unity at xk=±kc/ΔF, obviously increasing in k.
(a)
(b)


## 2.3.2. Embedding into Larger Spectral Domain

The functionsSE and SO may be finite in the spectral domain with limits fL and fU or centre frequency fC and bandwidth ΔF. Embedding this domain in a larger one with limits fL′<fL and fU′>fU (ΔF<ΔF′), another set of Fourier coefficients AE(O),k′′ and BE(O),k′′ is obtained with associated basis functions for the spectral domain. Without going into more detail here, these coefficients can be evaluated from AE(O),k and BE(O),k and the scalar products of the basis functions labeled with k and k′. For instance, one finds for the ratio of the means AE(O),0′/AE(O),0=ΔF/ΔF′<1. Basically, AE(O),k and BE(O),k maximise the information per coefficient when fL and fU 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 AE(O),0 is mapped to an infinite number of coefficients AE(O),k′′ and BE(O),k′′ which are mandatory to capture both discontinuities.In the spatial domain, the basis functions forAE(O),k′′ and BE(O),k′′ behave differently than the ones for AE(O),k and BE(O),k. Important to mention is the effect which the larger bandwidth has; sinck′± are spatially narrower than sinck±. In addition, the number of coefficients per spatial domain increases which is expressed by k′=kΔF′/ΔF.ForfC′ = 1873.7 GHz and ΔF′=2fC′ Figure 2 shows some basis functions (k′=0, 1, and 2) in the spatial domain. Indeed, the basis function for AE,0 with fC = 500 GHz and ΔF = 1000 GHz (see Figure 1(a)) is broader.Figure 2
Basis functions in spatial domain of Fourier coefficientsAE,k′′, BE,k′′ with order k′=0, 1, and 2 for even spectral function, being finite in spectral domain with centre frequency fC′ = 1873.7 GHz and bandwidth ΔF′=2fC′. For this domain, the sinc functions, contributing to the basis functions, are narrower than for fC = 500 GHz and ΔF = 1000 GHz (see Figure 1(a)). As an example, the basis function for AE,0 (dashed-black) is given. Since ΔF′≈4ΔF holds, approximately 4 times more Fourier coefficients locate inside a given spatial domain. ΔF′ was chosen to equal the Nyquist frequency fNy=c/(2Δx), having set the spatial sampling increment Δx=40 μ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 lengthlS of the band-limited function S is stated by(15)lS2=c∫fLfUS2df+c∫-fU-fLS2df=2c∫fLfUSE2+SO2df,because only the term SESO remains odd and cancels by the integration. Furthermore, the scaling by the factor c appears which originates in df/dσ=c. Replacing SE(O) with the expression (12) and exploiting that the basis functions are perpendicular in the spectral domain leaves(16)lS2=cΔF∑i=E,O2Ai,02+∑k=1∞Ai,k2+Bi,k2.Thus, the length is given by the sum of the Fourier coefficients squared. According to Parseval’s theorem,(17)lS2=∫-∞∞V2xdx,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 k1≠k2, and the spatial integral yields c/(2ΔF) for k1=k2=0 and c/ΔF for k1=k2>0. Analytically, this is hard to prove; however, this was numerically investigated and is considered to be valid.


### 2.4.1. Square-Integrable Functions

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


## 2.4.1. Square-Integrable Functions

The functionS is said to be square-integrable, when the condition lS2<∞ holds. Furthermore, if S is square-integrable, then the Fourier series representations in (12) converges towards SE and SO almost everywhere in the spectral domain as the order k grows [10]. Hence, the requirement on S 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 functionVx 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,V(x) is sampled at a finite set of optical path difference locations xi with i∈[1,ND], and ND marks the number of sample points. Usually, the sampling with constant increment Δx=xi+1-xi 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 xi might be unknown. If so, it is mandatory to introduce the zero-path difference x0 which is in the following set that x1=x0 holds.The finite spatial sampling leavesV(x) undetermined between the sampling nodes and outside the limits x1 and xND. Assuming Δx=const. holds, the Nyquist theorem states that the maximum frequency accessible is given by the Nyquist frequency fNy=c/(2Δx). Hence, for the spectral quantities SE and SO to be inferred, a maximum for the upper limit fU,Max=fNy follows from sampling theory. To prevent aliasing, Δx needs to be chosen small enough so that SE and SO vanish below fNy. In case, fU 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 fNy contributes to V(xi). 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 x0 has no profound impact.A diagnostic limitation is that the distancexND-x1 is finite, and, thus, no sampling is achieved below and above these limits. To gain information about SE and SO or the phase α (see (8)), V(x) 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 x1 and, thus, x0 are set to be negative, and the double-sided region is identified for the locations |xi|≤|x0|. This diminishes the maximum optical path difference xND achievable, being positive, and, thus, the length of the single-sided domain is identified by the relation |x0|<xi≤xND. Because xND scales with the order k 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 SE and SO are inferred. Hence, xND 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 asD(xi)=V(xi)+ϵi, and the actual interferometric data is expressed by the vector D→=V→+ϵ→. As spectral quantities are investigated, photons are involved in the measuring principle, and, hence, a part of ϵi 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 ϵi is a sample of a normal distribution with vanishing mean and a certain variance given by the squared noise level σD2. In any case, dedicated diagnostics tests are mandatory to characterise ϵi 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 by(18)Dxi=2∫fLfUSEcos⁡2πfcxi+SOsin⁡2πfcxidf+ϵiwith the additional information: the zero-path difference location x0<0, the spatial limits are restricted by x1=x0, xND>0 and xi<∞, the sampling increment Δx=xi+1-xi is constant and known, SE and SO have the same spectral limits which obey fL≥0 and fU≤fNy=c/(2Δx), and in case of Gaussian noise p(ϵi)=N(0,σD2). To be precise, the unknowns of the model are x0, fL, fU, SE(f) and SO(f), and σD.The basic model is a starting point and must be amended by diagnostic imperfections and specifics to the interferometer design type.


## 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 functionVx 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,V(x) is sampled at a finite set of optical path difference locations xi with i∈[1,ND], and ND marks the number of sample points. Usually, the sampling with constant increment Δx=xi+1-xi 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 xi might be unknown. If so, it is mandatory to introduce the zero-path difference x0 which is in the following set that x1=x0 holds.The finite spatial sampling leavesV(x) undetermined between the sampling nodes and outside the limits x1 and xND. Assuming Δx=const. holds, the Nyquist theorem states that the maximum frequency accessible is given by the Nyquist frequency fNy=c/(2Δx). Hence, for the spectral quantities SE and SO to be inferred, a maximum for the upper limit fU,Max=fNy follows from sampling theory. To prevent aliasing, Δx needs to be chosen small enough so that SE and SO vanish below fNy. In case, fU 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 fNy contributes to V(xi). 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 x0 has no profound impact.A diagnostic limitation is that the distancexND-x1 is finite, and, thus, no sampling is achieved below and above these limits. To gain information about SE and SO or the phase α (see (8)), V(x) 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 x1 and, thus, x0 are set to be negative, and the double-sided region is identified for the locations |xi|≤|x0|. This diminishes the maximum optical path difference xND achievable, being positive, and, thus, the length of the single-sided domain is identified by the relation |x0|<xi≤xND. Because xND scales with the order k 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 SE and SO are inferred. Hence, xND 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 asD(xi)=V(xi)+ϵi, and the actual interferometric data is expressed by the vector D→=V→+ϵ→. As spectral quantities are investigated, photons are involved in the measuring principle, and, hence, a part of ϵi 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 ϵi is a sample of a normal distribution with vanishing mean and a certain variance given by the squared noise level σD2. In any case, dedicated diagnostics tests are mandatory to characterise ϵi 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 by(18)Dxi=2∫fLfUSEcos⁡2πfcxi+SOsin⁡2πfcxidf+ϵiwith the additional information: the zero-path difference location x0<0, the spatial limits are restricted by x1=x0, xND>0 and xi<∞, the sampling increment Δx=xi+1-xi is constant and known, SE and SO have the same spectral limits which obey fL≥0 and fU≤fNy=c/(2Δx), and in case of Gaussian noise p(ϵi)=N(0,σD2). To be precise, the unknowns of the model are x0, fL, fU, SE(f) and SO(f), and σD.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 quantitiesSE and SO from an interferometric data set D→, the standard techniques rely on a noiseless model and follow a hierarchical ansatz. After making assumptions on the spectral limits fL and fU, the zero-path difference location x0 is estimated. The next step evaluates a phase α which is a measure of the ratio SO/SE, relying on the data located in the double-sided region. Given the model, the spectral limits, the spatial origin, and the phase, SE and SO 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 D→ 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 D→=V→, 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 interestSE and SO. 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 SE and SO 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 frequencyfNy=c/(2Δx) (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 Δx apart via their basis functions in the spatial domain, and the maximum order probed is artificially blown up to xND/(2Δx) due to the Nyquist assumptions (see Section 2.3). However, if the functions SE and SO are finite in the spectral domain with limits fL>0 and fU<fNy, then the embedding of the smaller domain [fL,fU] into the domain [0,fNy] 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 Δx=40 μm, that is, fL′ = 0 GHz and fU′ = 3747.4 GHz (fC′ = 1873.7 GHz, ΔF′=2fC′) with the wider basis function for the absolute term for fL = 0 GHz and fU = 1000 GHz (fC=500 GHz, ΔF = 1000 GHz).The uncertainty of a Fourier coefficient, relying on the Nyquist assumptions, scales likeσD/fNy (noise level/maximum bandwidth) which follows from the linear uncertainty propagation for (13). But for the band-limited case, the uncertainty would scale like (ΔF/fNy)1/2σD/ΔF, 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 differencex0 is most likely missed by the spatial sampling. One of the standard approaches to estimate x0 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 (AE,0 and AO,0) 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 x0.Withx0 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 SO and a decrease in SE 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 ofSE and SO 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 onSE and SO introduced by nonprobed Fourier coefficients. Hence, the requirement of square-integrability of the spectral functions is not taken into account.


## 2.6.1. Noiseless Model

The model used and stated by (18) lacks the noise contribution by definition, and, hence, treats the noisy data D→ 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 D→=V→, 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 interestSE and SO. 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 SE and SO 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 frequencyfNy=c/(2Δx) (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 Δx apart via their basis functions in the spatial domain, and the maximum order probed is artificially blown up to xND/(2Δx) due to the Nyquist assumptions (see Section 2.3). However, if the functions SE and SO are finite in the spectral domain with limits fL>0 and fU<fNy, then the embedding of the smaller domain [fL,fU] into the domain [0,fNy] 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 Δx=40 μm, that is, fL′ = 0 GHz and fU′ = 3747.4 GHz (fC′ = 1873.7 GHz, ΔF′=2fC′) with the wider basis function for the absolute term for fL = 0 GHz and fU = 1000 GHz (fC=500 GHz, ΔF = 1000 GHz).The uncertainty of a Fourier coefficient, relying on the Nyquist assumptions, scales likeσD/fNy (noise level/maximum bandwidth) which follows from the linear uncertainty propagation for (13). But for the band-limited case, the uncertainty would scale like (ΔF/fNy)1/2σD/ΔF, 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 differencex0 is most likely missed by the spatial sampling. One of the standard approaches to estimate x0 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 (AE,0 and AO,0) 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 x0.Withx0 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 SO and a decrease in SE 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 ofSE and SO 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 onSE and SO 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)p(d,v) captures the chance that the outcome d, let us say a data value or set, and the outcome v, a single model parameter or a set, are realised simultaneously.The product rule(19)pd,v=pd∣vpv=pv∣dpdintroduces the conditional probabilities for finding the outcome d, if the outcome v were true and vice versa. By the theorem of Bayes(20)pv∣d=pd∣vpdpv,one conditional probability can be expressed by the other, when the marginal distributions p(d)=∫pv∣dp(d)dv and p(v)=∫pd∣vp(v)dd are known. Hence, Bayes’ theorem captures the information/knowledge gained for v when a certain outcome for d has manifested. For the pdfs occurring in Bayes’ theorem, common names are used, that is, the posterior pv∣d, the likelihood pd∣v, the evidence p(d), and the prior p(v). The link or functional dependence d=f(v)+ϵ enters in the likelihood which takes into account known uncertainties ϵ like, for example, measurement noise. Any knowledge about v before new data is available can be found in the prior pv.Bayes’ rule can be extended to(21)pv,h∣d=pd∣vpdpv∣hph,introducing a set of hyperparameter h which enters per definition solely in the prior pv∣h. The additional pdf p(h) is called hyperprior which allocates trust in h. Apart from having the posterior for the parameters v, the marginalisation with respect to v reveals the posterior(22)ph∣d=∫pv,h∣ddv=pd∣hpdphfor h which measures the plausibility of an outcome of the hyperparameter given the data. Since p(d) does not depend on h, the most likely hyperparameter set is identified by the maximum of pd∣h, assuming p(h) 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 vectory→ (∈RN) be a multivariate normal with mean μ→ (∈RN) and covariance matrix Σ (∈RN×RN); then the pdf becomes(23)py→∣μ→,Σ=Ny→∣μ→,Σ=exp⁡-1/2y→-μ→TΣ-1y→-μ→2πN/2Σ1/2with 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 D→ (∈RND). The parameters of interest l→ (∈RNl) map linearly to the data domain like(24)D→=Ml→+ϵ→,where the ND×Nl dimensional matrixM 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 ϵi=N0,σD2 with vanishing mean and standard deviation σD (noise level), then the Gaussian likelihood(25)pD→∣l→=ND→∣l→,ΣD=exp⁡-1/2D→-Ml→TΣD-1D→-Ml→2πND/2ΣD1/2can be found with the covariance matrix ΣD=σD2δij.(b) Gaussian Prior. The prior information about l→ may be expressed by the multivariate normal(26)pl→=Nl→∣l→Pr,ΣPr=exp⁡-1/2l→-l→PrTΣPr-1l→-l→Pr2πNl/2ΣPr1/2with the prior mean l→Pr and the prior covariance ΣPr.(c) Gaussian Posterior and Evidence. Formally, Bayes’ theorem states the posterior by(27)pl→∣D→=ND→∣l→,ΣDpD→Nl→∣l→Pr,ΣPr.After some algebra, one can show that the posterior(28)pl→∣D→=Nl→∣l→Po,ΣPois a multivariate normal with posterior mean(29)l→Po=l→Pr+ΣPoMTΣD-1D→-Ml→Prand covariance(30)ΣPo=MTΣD-1M+ΣPr-1-1which are both analytically obtained. Furthermore, the evidence reads(31)pD→=exp⁡-1/2D→TΣD-1D→2πND/2ΣD1/2ΣPo1/2ΣPr1/2exp⁡1/2l→PoTΣPo-1l→Poexp⁡1/2l→PrTΣPr-1l→Pr,where 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 M of the parameters l→ to the data domain, as stated by (24), should remain valid. However, the mapping itself may depend on the parameters n→ in a nonlinear way, so that M=M(n→). This leaves the Gaussian likelihood in (25) formally unchanged but is symbolically stated as pD→∣l→,n→.(b) Priors. The Gaussian prior for l→ should be given by pl→∣h→=Nl→Prh→,ΣPrh→, where the prior mean and covariance depend on some of the hyperparameters h→. Similarly, a prior pn→∣h→ follows for the nonlinear parameters. Finally, the hyperparameters have an assigned prior p(h→).(c) Posteriors and Evidence. According to Bayes’ theorem, one can write the joint posterior like(32)pl→∣n→,h→,D→pn→,h→∣D→=pD→∣l→,n→pD→pl→∣h→pn→∣h→ph→,and the conditional amplitude posterior for l→ becomes a multivariate normal(33)pl→∣n→,h→,D→=Nl→∣l→Pon→,h→,ΣPon→,h→.Thus, both, the conditional posterior mean l→Po and covariance ΣPo evaluated by (29) and (30), depend on the nonlinear parameters and hyperparameters. After the trivial marginalisation with respect to l→, the joint posterior for n→ and h→(34)pn→,h→∣D→=exp⁡-1/2D→TΣD-1D→2πND/2ΣD1/21pD→ΣPo1/2ΣPr1/2exp⁡1/2l→PoTΣPo-1l→Poexp⁡1/2l→PrTΣPr-1l→Prpn→∣h→ph→remains. By expressing the posterior, named settings posterior in the following, like(35)pn→,h→∣D→=1K′ΣPo1/2ΣPr1/2exp⁡1/2l→PoTΣPo-1l→Poexp⁡1/2l→PrTΣPr-1l→Prpn→∣h→ph→,the evidence is identified with(36)pD→=K′exp⁡-1/2D→TΣD-1D→2πND/2ΣD1/2.Note, that the dimensionless constant K′, and, thus, the evidence depend on the chosen model, including likelihood and priors. Hence, K′ 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 pn→,h→∣D→, 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 quantity(37)p×K′=pn→,h→∣D→×K′is of interest, because it is numerically accessible. In case, p×K′ has a well distinguishable global maximum, p×K′ can be approximated by a multivariate normal which is estimated by evaluating the Hessian matrix. Thus, one finds the posterior means n→Po and h→Po with the associated posterior covariance. This allows an approximate marginalisation with respect to n→ and/or h→. This can be understood as a propagation of the posterior uncertainties in n→ and h→ to the marginalised posterior pl→∣D→ 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 l→Pr is set to 0, and the priors pn→∣h→ and p(h→) are chosen to be uniform. Then, one can set K′=Kpn→∣h→p(h→) which modifies (37) to(38)p×K=ΣPo1/2ΣPr1/2exp⁡12l→PoTΣPo-1l→Po.


## 3.1. Bayes’ Theorem

The joint probability density function (pdf)p(d,v) captures the chance that the outcome d, let us say a data value or set, and the outcome v, a single model parameter or a set, are realised simultaneously.The product rule(19)pd,v=pd∣vpv=pv∣dpdintroduces the conditional probabilities for finding the outcome d, if the outcome v were true and vice versa. By the theorem of Bayes(20)pv∣d=pd∣vpdpv,one conditional probability can be expressed by the other, when the marginal distributions p(d)=∫pv∣dp(d)dv and p(v)=∫pd∣vp(v)dd are known. Hence, Bayes’ theorem captures the information/knowledge gained for v when a certain outcome for d has manifested. For the pdfs occurring in Bayes’ theorem, common names are used, that is, the posterior pv∣d, the likelihood pd∣v, the evidence p(d), and the prior p(v). The link or functional dependence d=f(v)+ϵ enters in the likelihood which takes into account known uncertainties ϵ like, for example, measurement noise. Any knowledge about v before new data is available can be found in the prior pv.Bayes’ rule can be extended to(21)pv,h∣d=pd∣vpdpv∣hph,introducing a set of hyperparameter h which enters per definition solely in the prior pv∣h. The additional pdf p(h) is called hyperprior which allocates trust in h. Apart from having the posterior for the parameters v, the marginalisation with respect to v reveals the posterior(22)ph∣d=∫pv,h∣ddv=pd∣hpdphfor h which measures the plausibility of an outcome of the hyperparameter given the data. Since p(d) does not depend on h, the most likely hyperparameter set is identified by the maximum of pd∣h, assuming p(h) 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 vectory→ (∈RN) be a multivariate normal with mean μ→ (∈RN) and covariance matrix Σ (∈RN×RN); then the pdf becomes(23)py→∣μ→,Σ=Ny→∣μ→,Σ=exp⁡-1/2y→-μ→TΣ-1y→-μ→2πN/2Σ1/2with 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 D→ (∈RND). The parameters of interest l→ (∈RNl) map linearly to the data domain like(24)D→=Ml→+ϵ→,where the ND×Nl dimensional matrixM 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 ϵi=N0,σD2 with vanishing mean and standard deviation σD (noise level), then the Gaussian likelihood(25)pD→∣l→=ND→∣l→,ΣD=exp⁡-1/2D→-Ml→TΣD-1D→-Ml→2πND/2ΣD1/2can be found with the covariance matrix ΣD=σD2δij.(b) Gaussian Prior. The prior information about l→ may be expressed by the multivariate normal(26)pl→=Nl→∣l→Pr,ΣPr=exp⁡-1/2l→-l→PrTΣPr-1l→-l→Pr2πNl/2ΣPr1/2with the prior mean l→Pr and the prior covariance ΣPr.(c) Gaussian Posterior and Evidence. Formally, Bayes’ theorem states the posterior by(27)pl→∣D→=ND→∣l→,ΣDpD→Nl→∣l→Pr,ΣPr.After some algebra, one can show that the posterior(28)pl→∣D→=Nl→∣l→Po,ΣPois a multivariate normal with posterior mean(29)l→Po=l→Pr+ΣPoMTΣD-1D→-Ml→Prand covariance(30)ΣPo=MTΣD-1M+ΣPr-1-1which are both analytically obtained. Furthermore, the evidence reads(31)pD→=exp⁡-1/2D→TΣD-1D→2πND/2ΣD1/2ΣPo1/2ΣPr1/2exp⁡1/2l→PoTΣPo-1l→Poexp⁡1/2l→PrTΣPr-1l→Pr,where 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 M of the parameters l→ to the data domain, as stated by (24), should remain valid. However, the mapping itself may depend on the parameters n→ in a nonlinear way, so that M=M(n→). This leaves the Gaussian likelihood in (25) formally unchanged but is symbolically stated as pD→∣l→,n→.(b) Priors. The Gaussian prior for l→ should be given by pl→∣h→=Nl→Prh→,ΣPrh→, where the prior mean and covariance depend on some of the hyperparameters h→. Similarly, a prior pn→∣h→ follows for the nonlinear parameters. Finally, the hyperparameters have an assigned prior p(h→).(c) Posteriors and Evidence. According to Bayes’ theorem, one can write the joint posterior like(32)pl→∣n→,h→,D→pn→,h→∣D→=pD→∣l→,n→pD→pl→∣h→pn→∣h→ph→,and the conditional amplitude posterior for l→ becomes a multivariate normal(33)pl→∣n→,h→,D→=Nl→∣l→Pon→,h→,ΣPon→,h→.Thus, both, the conditional posterior mean l→Po and covariance ΣPo evaluated by (29) and (30), depend on the nonlinear parameters and hyperparameters. After the trivial marginalisation with respect to l→, the joint posterior for n→ and h→(34)pn→,h→∣D→=exp⁡-1/2D→TΣD-1D→2πND/2ΣD1/21pD→ΣPo1/2ΣPr1/2exp⁡1/2l→PoTΣPo-1l→Poexp⁡1/2l→PrTΣPr-1l→Prpn→∣h→ph→remains. By expressing the posterior, named settings posterior in the following, like(35)pn→,h→∣D→=1K′ΣPo1/2ΣPr1/2exp⁡1/2l→PoTΣPo-1l→Poexp⁡1/2l→PrTΣPr-1l→Prpn→∣h→ph→,the evidence is identified with(36)pD→=K′exp⁡-1/2D→TΣD-1D→2πND/2ΣD1/2.Note, that the dimensionless constant K′, and, thus, the evidence depend on the chosen model, including likelihood and priors. Hence, K′ 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 pn→,h→∣D→, 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 quantity(37)p×K′=pn→,h→∣D→×K′is of interest, because it is numerically accessible. In case, p×K′ has a well distinguishable global maximum, p×K′ can be approximated by a multivariate normal which is estimated by evaluating the Hessian matrix. Thus, one finds the posterior means n→Po and h→Po with the associated posterior covariance. This allows an approximate marginalisation with respect to n→ and/or h→. This can be understood as a propagation of the posterior uncertainties in n→ and h→ to the marginalised posterior pl→∣D→ 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 l→Pr is set to 0, and the priors pn→∣h→ and p(h→) are chosen to be uniform. Then, one can set K′=Kpn→∣h→p(h→) which modifies (37) to(38)p×K=ΣPo1/2ΣPr1/2exp⁡12l→PoTΣPo-1l→Po.


## 3.2.1. Multivariate Normal

Let the joint pdf for the random vectory→ (∈RN) be a multivariate normal with mean μ→ (∈RN) and covariance matrix Σ (∈RN×RN); then the pdf becomes(23)py→∣μ→,Σ=Ny→∣μ→,Σ=exp⁡-1/2y→-μ→TΣ-1y→-μ→2πN/2Σ1/2with 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 D→ (∈RND). The parameters of interest l→ (∈RNl) map linearly to the data domain like(24)D→=Ml→+ϵ→,where the ND×Nl dimensional matrixM 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 ϵi=N0,σD2 with vanishing mean and standard deviation σD (noise level), then the Gaussian likelihood(25)pD→∣l→=ND→∣l→,ΣD=exp⁡-1/2D→-Ml→TΣD-1D→-Ml→2πND/2ΣD1/2can be found with the covariance matrix ΣD=σD2δij.(b) Gaussian Prior. The prior information about l→ may be expressed by the multivariate normal(26)pl→=Nl→∣l→Pr,ΣPr=exp⁡-1/2l→-l→PrTΣPr-1l→-l→Pr2πNl/2ΣPr1/2with the prior mean l→Pr and the prior covariance ΣPr.(c) Gaussian Posterior and Evidence. Formally, Bayes’ theorem states the posterior by(27)pl→∣D→=ND→∣l→,ΣDpD→Nl→∣l→Pr,ΣPr.After some algebra, one can show that the posterior(28)pl→∣D→=Nl→∣l→Po,ΣPois a multivariate normal with posterior mean(29)l→Po=l→Pr+ΣPoMTΣD-1D→-Ml→Prand covariance(30)ΣPo=MTΣD-1M+ΣPr-1-1which are both analytically obtained. Furthermore, the evidence reads(31)pD→=exp⁡-1/2D→TΣD-1D→2πND/2ΣD1/2ΣPo1/2ΣPr1/2exp⁡1/2l→PoTΣPo-1l→Poexp⁡1/2l→PrTΣPr-1l→Pr,where 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 M of the parameters l→ to the data domain, as stated by (24), should remain valid. However, the mapping itself may depend on the parameters n→ in a nonlinear way, so that M=M(n→). This leaves the Gaussian likelihood in (25) formally unchanged but is symbolically stated as pD→∣l→,n→.(b) Priors. The Gaussian prior for l→ should be given by pl→∣h→=Nl→Prh→,ΣPrh→, where the prior mean and covariance depend on some of the hyperparameters h→. Similarly, a prior pn→∣h→ follows for the nonlinear parameters. Finally, the hyperparameters have an assigned prior p(h→).(c) Posteriors and Evidence. According to Bayes’ theorem, one can write the joint posterior like(32)pl→∣n→,h→,D→pn→,h→∣D→=pD→∣l→,n→pD→pl→∣h→pn→∣h→ph→,and the conditional amplitude posterior for l→ becomes a multivariate normal(33)pl→∣n→,h→,D→=Nl→∣l→Pon→,h→,ΣPon→,h→.Thus, both, the conditional posterior mean l→Po and covariance ΣPo evaluated by (29) and (30), depend on the nonlinear parameters and hyperparameters. After the trivial marginalisation with respect to l→, the joint posterior for n→ and h→(34)pn→,h→∣D→=exp⁡-1/2D→TΣD-1D→2πND/2ΣD1/21pD→ΣPo1/2ΣPr1/2exp⁡1/2l→PoTΣPo-1l→Poexp⁡1/2l→PrTΣPr-1l→Prpn→∣h→ph→remains. By expressing the posterior, named settings posterior in the following, like(35)pn→,h→∣D→=1K′ΣPo1/2ΣPr1/2exp⁡1/2l→PoTΣPo-1l→Poexp⁡1/2l→PrTΣPr-1l→Prpn→∣h→ph→,the evidence is identified with(36)pD→=K′exp⁡-1/2D→TΣD-1D→2πND/2ΣD1/2.Note, that the dimensionless constant K′, and, thus, the evidence depend on the chosen model, including likelihood and priors. Hence, K′ 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 pn→,h→∣D→, 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 quantity(37)p×K′=pn→,h→∣D→×K′is of interest, because it is numerically accessible. In case, p×K′ has a well distinguishable global maximum, p×K′ can be approximated by a multivariate normal which is estimated by evaluating the Hessian matrix. Thus, one finds the posterior means n→Po and h→Po with the associated posterior covariance. This allows an approximate marginalisation with respect to n→ and/or h→. This can be understood as a propagation of the posterior uncertainties in n→ and h→ to the marginalised posterior pl→∣D→ 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 l→Pr is set to 0, and the priors pn→∣h→ and p(h→) are chosen to be uniform. Then, one can set K′=Kpn→∣h→p(h→) which modifies (37) to(38)p×K=ΣPo1/2ΣPr1/2exp⁡12l→PoTΣPo-1l→Po.


## 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 S and, thus, SE and SO 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 SE and SO. 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 SE and SO. 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 toT. 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 by(39)ΣBB′s,t=-stT+s:s≤tt:s>tfor s and t∈[0,T].


### 4.2. Adapted Definition

The even and odd spectral functions, being finite for the interval[fL,fU], 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 fL and the interval length T=fU-fL=ΔF need to be set. Thus, one gets(40)ΣBB′f,f′=-f-fLf′-fLΔF+f-fL:f≤f′f′-fL:f>f′with the unitΣBB′ = Hz. With the normalisation(41)∫fLfU∫fLfUΣBB′f,f′dfdf′=ΔF312,the modified covariance becomes(42)ΣBBf,f′=12ΔF3ΣBB′f,f′,where [ΣBB] = Hz−2 has now the proper unit with respect to the spectral scale.The parametersaE and aO of unit [aE(O)]=V2 are introduced which are defined each as a scaling factor for the associated process. With these scalings the covariances ΣBB,E(O)=aE(O)ΣBB are obtained, so that the units [ΣBB,E(O)]=SEO2 = 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 (43)ΣBB,Ek,k′=aEΔF2∫fLfU∫fLfU12cos⁡2πkf-fCΔF2sin⁡2πkf-fCΔF12cos⁡2πk′f′-fCΔF2sin⁡2πk′f′-fCΔFTΣBBf,f′dfdf′=aEΔF26π2π26-1k′+1k′20→T-1k+1k2δk,k′kk′00→0δk,k′kk′=ΣAE,0,AE,0Σ→TAE,0,AE,k′0→TΣ→AE,k,AE,0ΣAE,k,AE,k′00→0ΣBE,k,BE,k′for all k,k′≥1, and similarly ΣBB,Ok,k′ 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 Σ→(AE(O),k,AE(O),0) and Σ→TAEO,0,AEO,k′, respectively. This is caused by the condition that SE(O) vanishes at the spectral boundaries where the sine vanishes intrinsically for any k 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,ΣBB,EOk,k′ has no off-diagonal elements due to the Kronecker delta δk,k′, meaning that these coefficients are independent on each other. Furthermore, for the infinite-dimensional matrices ΣAEO,k,AEO,k′=ΣBEO,k,BEO,k′=δk,k′/kk′ 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)), lS2 is evaluated by summing the squares of the Fourier coefficients. Because the entries of the main diagonal in the covariances ΣBB,Ek,k′ and ΣBB,Ok,k′ drop with the order squared, the Brownian bridge process ensures square-integrability of SE and SO as long as the scalings aE and aO remain finite.


### 4.5. Signal Envelope

For the even process, the signal level can be estimated by the envelopeσD,E′xk in the data domain. Starting point is the square root of the main diagonal of ΣBB,Ek,k′=k. Since the argument of sinck± (see (14)) localises the even and odd Fourier coefficient at a fixed k=±xkΔF/c in the same data domain, the even and odd contributions of ΣBB,Ek,k′=k must be added for k≥1. 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 ΔF. In total, one finds the envelope as(44)σD,E′xk=±2aE1/22:xk=0121/2πcxkΔF:xk≠0.In the above equation, factor 2 in front of aE was chosen, so that σD,E′(xk) captures most of the signal. An approximation might be convenient, because 121/2/π≈ 1.For the envelopeσD,O′xk 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 finds(45)σD,O′xk=±2aO1/20:xk=0121/2πcxkΔF:xk≠0.Both envelops drop with 1/|xk|, and, thus, most of the signal associated with each process would in the data domain.


## 4.1. Standard Definition

The Brownian bridge is a continuous stochastic process for an interval, say from 0 toT. 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 by(39)ΣBB′s,t=-stT+s:s≤tt:s>tfor s and t∈[0,T].


## 4.2. Adapted Definition

The even and odd spectral functions, being finite for the interval[fL,fU], 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 fL and the interval length T=fU-fL=ΔF need to be set. Thus, one gets(40)ΣBB′f,f′=-f-fLf′-fLΔF+f-fL:f≤f′f′-fL:f>f′with the unitΣBB′ = Hz. With the normalisation(41)∫fLfU∫fLfUΣBB′f,f′dfdf′=ΔF312,the modified covariance becomes(42)ΣBBf,f′=12ΔF3ΣBB′f,f′,where [ΣBB] = Hz−2 has now the proper unit with respect to the spectral scale.The parametersaE and aO of unit [aE(O)]=V2 are introduced which are defined each as a scaling factor for the associated process. With these scalings the covariances ΣBB,E(O)=aE(O)ΣBB are obtained, so that the units [ΣBB,E(O)]=SEO2 = 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 (43)ΣBB,Ek,k′=aEΔF2∫fLfU∫fLfU12cos⁡2πkf-fCΔF2sin⁡2πkf-fCΔF12cos⁡2πk′f′-fCΔF2sin⁡2πk′f′-fCΔFTΣBBf,f′dfdf′=aEΔF26π2π26-1k′+1k′20→T-1k+1k2δk,k′kk′00→0δk,k′kk′=ΣAE,0,AE,0Σ→TAE,0,AE,k′0→TΣ→AE,k,AE,0ΣAE,k,AE,k′00→0ΣBE,k,BE,k′for all k,k′≥1, and similarly ΣBB,Ok,k′ 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 Σ→(AE(O),k,AE(O),0) and Σ→TAEO,0,AEO,k′, respectively. This is caused by the condition that SE(O) vanishes at the spectral boundaries where the sine vanishes intrinsically for any k 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,ΣBB,EOk,k′ has no off-diagonal elements due to the Kronecker delta δk,k′, meaning that these coefficients are independent on each other. Furthermore, for the infinite-dimensional matrices ΣAEO,k,AEO,k′=ΣBEO,k,BEO,k′=δk,k′/kk′ 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)), lS2 is evaluated by summing the squares of the Fourier coefficients. Because the entries of the main diagonal in the covariances ΣBB,Ek,k′ and ΣBB,Ok,k′ drop with the order squared, the Brownian bridge process ensures square-integrability of SE and SO as long as the scalings aE and aO remain finite.


## 4.5. Signal Envelope

For the even process, the signal level can be estimated by the envelopeσD,E′xk in the data domain. Starting point is the square root of the main diagonal of ΣBB,Ek,k′=k. Since the argument of sinck± (see (14)) localises the even and odd Fourier coefficient at a fixed k=±xkΔF/c in the same data domain, the even and odd contributions of ΣBB,Ek,k′=k must be added for k≥1. 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 ΔF. In total, one finds the envelope as(44)σD,E′xk=±2aE1/22:xk=0121/2πcxkΔF:xk≠0.In the above equation, factor 2 in front of aE was chosen, so that σD,E′(xk) captures most of the signal. An approximation might be convenient, because 121/2/π≈ 1.For the envelopeσD,O′xk 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 finds(45)σD,O′xk=±2aO1/20:xk=0121/2πcxkΔF:xk≠0.Both envelops drop with 1/|xk|, 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 D→1 is acquired in terms of Volts dependent on the optical path difference x. 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 D(xi)=D2(xi)-D1(xi) acquired at the spatial grid node xi. Then, the abstract model for the Martin-Puplett interferometer is stated by (46)Dxi=A22∫fLfUSEcos⁡2πfcxi+SOsin⁡2πfcxidf+VOffxi+Ni0,σD2,using the total amplification A of the detection system. Furthermore, the offset VOff marks a diagnostic imperfection which varies with xi. The Gaussian noise contributes to each data sample described by Ni(0,σD2). The unknown quantities in the diagnostic model are the spatial grid xi, the lower and upper spectral boundaries fL and fU of the Fourier transform integral, the even and odd functions SE and SO dependent on frequency f, and the offset.


#### 5.1.2. Interferometric Data

The data setD→, that is, the difference interferogram consists of ND = 788 values (see Figure 3(a)). Merely for graphical presentation a certain x 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 interferogramD→ (black) versus optical path difference x. The data set includes ND = 788 values, and the noise level for each Di reads σD = 132.29 μV. The values of the spatial grid are determined by the standard model xi=i-1Δx+x0 with i∈1,ND, the constant increment Δx = 40 μm, and the origin x0. This origin is estimated as −5.118 mm by a fit of a quadratic polynomial to the maximum max⁡D→ and its neighbouring values to either side. Clearly, D→ deviates from the baseline (red). (b) Difference quotient for different subsets of D→. Since the zig-zag-pattern in the point-to-point variation appears only for the first case (black), the standard model xi+1-xi = const. seems to be inadequate. Instead, a model seems appropriate for which x2i+1-x2i-1 = x2i+2-x2i=const. holds.
(a)
(b)The componentsDi are measured independently on each other, and the noise level for each Di is captured by σD = 132.29 μV, and, thus, the variance of the whole data vector is stated by the matrix ΣD=σD2δij.


#### 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Δx = 40 μm. Hence, the standard model xi=(i-1)Δx+x0 is obvious with i=1,…,ND, and the zero-path difference x0 is a free parameter. However, this model is not accurate. Applying the standard approach which fits a second-order polynomial to the maximum max⁡D→ and its two nearby values, x0=-127.95Δx=-5.118 mm, is inferred for the data set shown in Figure 3(a). Furthermore, the difference quotient evaluated by (Di+1-Di)/Δx is presented versus the optical path difference x=(xi+1+xi)/2 in Figure 3(b). The point-to-point variation of the quotient has a zig-zag pattern which implies that the assumption Δx = const. is incorrect. Indeed, each of two the difference quotients (D2i+1-D2i-1)/(2Δx) and (D2i+2-D2i)/(2Δx) is smoother. Hence, the model(47)xi=i-1Δx+x0+0:ifiis  oddδx:ifiis  evenseems more appropriate, making use of two free parameters: the zero-path difference x0 and a shift δx for every other grid value. The priors p(x0) and p(δx) are set to be uniform.


#### 5.1.4. Offset

The upward trend ofD→ (see Figure 3(a)) is modelled by the offset VOff(xi). Here, a second-order polynomial is chosen to capture the offset by VOff(xi)=b0+b1xi+b2xi2 with the free parameters b0, b1, and b2 being summarised by the vector b→=b0,b1,b2T. The corresponding mapping to the data domain can be expressed as(48)MOffb→=A21x1x121x2x22⋮⋮⋮1xNDxND2b0b1b2with the ND×3 matrix MOff. The joint prior is expressed by the factorisable multivariate normal distribution:(49)pb→=∏i=02N0,σbi,Pr2=N0,Σb→,Pr,where σbi,Pr 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 f  ∈  ]fL,fU[. Thereby, fL and fU are free parameters. The spectral domain is dicretised using the constant increment Δf which is considered as a free parameter as well. Then, the spectral domain is represented by the set fj=fL+(j-0.5)Δf with j=1,…,Nf, covering the band NfΔf=(fU-fL)=ΔF centred at fC=(fU+fL)/2. To be clear, since fL, fU, and Δf are free parameters and will be inferred, the number Nf is variable. This number determines the dimensionality of the vectors S→E and S→O which represent the discretised functions SE(f) and SO(f). The mapping of S→E and S→O to the data domain in (46) is written as(50)A22∫fLfUSEcos⁡2πfcxi+SOsin⁡2πfcxidf≈ME,MOS→ES→Owith the two ND×Nf matrices(51)ME=Acos⁡2πf1cx1⋯cos⁡2πfNfcx1⋮⋱⋮cos⁡2πf1cxND⋯cos⁡2πfNfcxNDΔf,MO=Asin⁡2πf1cx1⋯sin⁡2πfNfcx1⋮⋱⋮sin⁡2πf1cxND⋯sin⁡2πfNfcxNDΔf.The joint prior for the two vectorial quantities S→E and S→O is factorised, and each prior is chosen as a multivariate normal distribution with vanishing mean. Since the Brownian bridge covariance ΣBB (see Section 4) describes functions which vanish at the boundaries fL and fU and are square-integrable, and its signal envelope decays with the optical path difference like the interferometric data at hand, the priors are chosen by(52)pS→EO∣Δf,fL,fU,aEO=N0,aEOΣBBΔf,fL,fUwith the two hyperparameters aE and aO for scaling. For each of these hyperparameters a uniform prior paEO=1/ΔaE(O) is applied.Finally, the joint prior(53)pΔf,fL,fU=1KΔf,fL,fUshould be constant, so that any combination of Δf, fL, and fU has the same probability. Furthermore, the conditions 0≤Δf≤fU-fL, 0≤fL≤fU, and 0≤fU≤ΔfU (global upper limit) must be fulfilled. For example, the upper limit of fU is set here to the Nyquist frequency ΔfU=fNy=c/(2Δx) = 3747.4 GHz.


### 5.2. Bayes’ Theorem

In the following, the linear parameters are summerised by the setl→=S→E,S→O,b→T, the nonlinear parameters by n→=fL,fU,x0,δxT, and the hyperparameters by h→=aE,aO,σb0,Pr,σb1,Pr,σb2,PrT.The matrixM(n→,Δf)=(ME,MO,MOff) maps the parameters l→ to the data domain. Hence, the Gaussian likelihood is written as(54)pD→∣l→,n→,Δf=exp⁡-1/2D→-Ml→TΣD-1D→-Ml→2πND/2ΣD1/2.The prior for the full problem takes the form(55)pl→,Δf,n→,h→=pl→∣Δf,fL,fU,aE,aO,σb0,Pr,σb1,Pr,σb2,PrpΔf,fL,fUpaEpaOpx0pδxpσb0,Prpσb1,Prpσb2,Prwith the multivariate normal prior(56)pl→∣Δf,fL,fU,aE,aO,σb→,Pr=exp⁡-1/2l→TΣPr-1l→2π2Nf+3/2ΣPr1/2for l→, using the (2Nf+3)×(2Nf+3) dimensional covariance matrix(57)ΣPr=aEΣBB000aOΣBB000Σb→,Pr.One gets the joint posterior(58)pl→,n→,h→,Δf∣D→=pl→∣n→,h→,Δf,D→pn→,h→,Δf∣D→with the conditional amplitude posterior given by the multivariate normal(59)pl→∣n→,h→,Δf,D→=exp⁡-1/2l→-l→PoTΣPo-1l→-l→Po2π2Nf+3/2ΣPo1/2,using the posterior covariance matrix ΣPo=MTΣD-1M+ΣPr-1-1 and the mean l→Po=ΣPoMTΣD-1D→. Furthermore, one obtains the settings posterior(60)pn→,h→,Δf∣D→=p×KK=ΣPo1/2ΣPr1/2exp⁡1/2l→PoTΣPo-1l→PoK,where the constant K 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 forx0, δx, fL, fU, Δf, aE and aO, σb0,Pr, σb1,Pr, and σb2,Pr.The values forx0=-127.95Δx=-5.118 mm and δx=0 mm are chosen to form an optical path difference grid xi for which all subsequent nodes are separated by Δx. Then, |xi|≤5.118 mm identifies the double-sided domain (256 data points), and the single-sided domain is marked by 5.118 mm ≤xi≤ 26.362 mm (532 data points).The spectral domain is covered fromfL = 0 GHz to the Nyquist frequency fU=c/(2Δx) = 3747.4 GHz (=ΔF). With the bandwidth and the optical path difference set, the conversion k=|x|ΔF/c (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Δf = 5.68 GHz is used, and the lower and upper spectral limits are set to fL = 0 GHz and fU = 3747.4 GHz (Nyquist assumptions), respectively. The scalings of the processes are chosen by aE=aO=10-8 V2. (a) Prior samples S→E,Pr,Sa and S→O,Pr,Sa 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 D→E(O),Pr,Sa for optical path difference obtained for x0=-5.118 mm and δx = 0. With increasing distance to the origin, the amplitude drops which is well bounded by the signal envelopes σD,E=σD,O for the chosen bandwidth ΔF=fU. Fourier coefficients up to order k = 64 and 330 settle in the double-sided domain x≤x0 and single-sided domain x>x0, respectively.
(a)
(b)To match the maximum order (330) of probed Fourier coefficient, the discretisation is set toΔf=c/(2Δx(256/2+532))≈5.68 GHz which reflects the classical increment. Then, the spectral grid has Nf = 660 elements.The hyperparameters are set to large values likeaE=aO=10-8 V2 to prevent from a prior determined posterior. This unwanted feature would be present, if one would choose aE(O) too small, so that the signal envelops do not include the measured data (see next paragraph).With the chosen spectral priors, sample functions/vectorsS→E,Pr,Sa and S→O,Pr,Sa are drawn and shown in Figure 4(a). The function values are of the order of 10-16 V/Hz decaying towards the spectral boundaries. The mapping of the spectral prior samples to the spatial domain via D→E,Pr,Sa+D→O,Pr,Sa=MES→E,Pr,Sa+MOS→O,Pr,Sa shows that D→E,Pr,Sa+D→O,Pr,Sa 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 D→E,Pr,Sa+D→O,Pr,Sa drops with the distance from the origin. All these data samples are well bounded by the signal envelopes σD,E(O) obtained by σD,E(O)=(A/2)σD,E(O)′ which is the scaled version of the expressions (44) and (45).Forb→ the prior covariance is determined by setting σb0,Pr=1 V, σb1,Pr=1 V/m and σb2,Pr=1 V/m2.(a) Posterior Mean. For the above values, the mean l→Po of the amplitude posterior is given by b→Po=(-61.41 nV, 5.75 μV/m, 2.42 μV/m2)T, S→E,Po and S→O,Po as presented in Figure 5(a). Clearly, S→E,Po peaks below 1000 GHz with the maximum of approximately 6 × 10−18 V/Hz around 400 GHz. Furthermore, S→O,Po is much smaller than S→E,Po below 1000 GHz.Figure 5
Investigation of conditional amplitude posteriorpl→∣Δf,fL,fU,aE,aO,x0,δx,σb0,Pr,σb1,Pr,σb2,Pr,D→ for linear parameters l→=S→E,S→O,b→T. The posterior is obtained for the settings x0=-5.118 mm and δx = 0, fL = 0 GHz, fU = 3747.4 GHz, Δf = 5.68 GHz, aE=aO=10-8 V2, σb0,Pr=1 V, σb1,Pr=1 V/m, and σb2,Pr=1 V/m2. (a) Spectral posterior means S→E(O),Po and one posterior sample. Below 1000 GHz, the means show a finite amplitude. Each sample deviates strongly from its corresponding mean. (b) Posterior mean l→Po and 100 samples l→Po,Sa mapped to give even, odd, and offset contributions in data domain. While in the double-sided region (|x|≤|x0|) most of the data described by S→E,Po, S→E,Po and S→O,Po give same contributions in the single-sided region (x>|x0|). In the single-sided region, the spread of the even (odd) contributions for the posterior samples is large but bounded by the signal envelops σD,E(O). (c) Residuals for l→Po (black) and one sample l→Po,Sa (red) both mapped to data domain. An almost perfect description of the noisy data is achieved by lPo, indicating an overfitting of the data. (d) Measured data D→ and 100 samples of amplitude posterior mapped to data domain. The individual contributions have a much wider spread (see (b)), but the total sum D→Po,Sa agrees well with D→ and its uncertainties. This is caused by the posterior covariance of l→ with suited off-diagonal entries.
(a)
(b)
(c)
(d)The offsetD→Off,Po=MOffb→Po contribution to the data has an upwards trend but is small (see Figure 5(b)). The even and odd contributions D→E(O),Po=ME(O)S→E(O),Po behave differently in the double- and single-sided domains (see Figure 5(b)). While in the double-sided domain D→E,Po is much larger than D→O,Po, D→E,Po and D→O,Po equal each other in the single-sided domain. This indicates an underestimation (overestimation) of the Fourier coefficients above the order 64 for SE,Po (SO,Po). The cause of this finding relies in the choice aE=aO (σD,E=σD,O) which becomes more clear in Section 5.3.3.Since the quantityD→Po=Ml→Po equals almost the data set D→ (see residuals in Figure 5(c)), a nearly perfect match is achieved. Hence, all noise contributions are captured by the mean l→Po, and an overfitting of the data is obvious.(b) Posterior Covariance. The posterior covariance ΣPo is not characterised in detail. From the elements of the main diagonal of ΣPo, one finds the values σb0,Po = 73.30 nV, σb1,Po = 12.58 μV/m, and σb2,Po = 446.67 μV/m2. These posterior uncertainties are quite large when compared to the values of the corresponding mean b→Po, especially for the quadratic coefficient.For the spectral quantities, the square root of the main diagonal elements ofΣPo is of the order of some 10-18 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 l→Po,Sa=S→E,Po,Sa,S→O,Po,Sa,b→Po,SaT are drawn from the amplitude posterior and mapped to the data domain (see Figure 5(b)). Thereby, the contributions D→E,Po,Sa=MES→E,Po,Sa, D→O,Po,Sa=MOS→O,Po,Sa, and D→Off,Po,Sa=MOffb→Po,Sa are split to investigate their interplay. D→Off,Po,Sa has a noticeable width but remains in the vicinity of D→Off,Po for the whole spatial domain.For the double-sided regionD→E,Po,Sa and D→O,Po,Sa form a narrow band around the corresponding posterior mean quantities D→E,Po and D→O,Po. On the contrary, for the single-sided region D→E,Po,Sa and D→O,Po,Sa increase drastically but are restricted by σD,E and σD,O. The deviation is large between the individual and the posterior mean quantities; though, the sum D→Po,Sa=D→E,Po,Sa+D→O,Po,Sa+D→Off,Po,Sa remains well inside the 2-σD band surrounding D→ (see Figures 5(c) and 5(d)). Hence, the posterior covariance ΣPo captures properly the correlations.With the specific values of the settings one gets the numberln⁡p×K≈1.0488×106.


#### 5.3.2. Settings Posterior

As the problem is formulated, the settings posterior is proportional top×K. Its optimisation, that is, finding the global maximum of the settings posterior is a 10-dimensional problem. Since large numbers are involved, one has to investigate ln⁡p×K. In general, if only one parameter is changed, a well distinguishable peak in ln⁡p×K is found. Hence, the optimisation is currently carried out by varying each parameter separately (coordinate descent algorithm). Thereby, ln⁡p×K increases when Δf becomes smaller. This implies an increase of the dimensionality of the involved covariances ΣPr and ΣPo and, thus, a prolonging of the optimisation procedure. To demonstrate this procedure, the parameters fL = 0 GHz, aE=aO=10-8 V2, x0=-5.118 mm, δx=0 mm, σb0,Pr=1 V, σb1,Pr=1 V/m, and σb2,Pr=1 V/m2 are set to the values used for the nonoptimised case (see Section 5.3.1). Scanned roughly in fU, ln⁡p×K shows a peak close to 1000 GHz (see Figure 6(a)). Furthermore, this peak increases for Δf≤5 GHz to ln⁡p×K≈106. The peak is localised at 860 GHz which moves to about 910 GHz for the reduced values aE=aO=10-11 V2 (see Figure 6(b)).Figure 6
Ten-dimensional quantityp×K scanned at different values for fU, Δf, aE, and aO while keeping fL = 0 GHz, x0=-5.118 mm, δx=0 mm, σb0,Pr = 1 V, σb1,Pr = 1 V/m, and σb2,Pr = 1 V/m2. (a) Rough scan in fU and Δf at aE=aO=10-8 V2. A peak for ln⁡p×k with the amplitude of the order of 106 is found close to fU = 1000 GHz. Small Δf are in favour. (b) Scan in fU and Δf close to peaking for two sets of aE and aO. The peak locates in the vicinity of 910 GHz when aE and aO reduce by some orders of magnitude.
(a)
(b)To ease the computational effort but still being able to characteriseln⁡p×K, its maximum is determined dependent on Δf which is scanned in the values 5.68, 4, 3, 2, 1, and 1/2 GHz. Each maximum is captured by the sets fL,Po, fU,Po, aE,Po, aO,Po, x0,Po, δxPo, σb0,Pr,Po, σb1,Pr,Po, and σb2,Pr,Po summarised by Table 1. All maxima locate at very similar sets which seem to converge as Δf becomes smaller. Although, in relative terms, the maximum for a smaller Δf has higher odds (see Figure 7). For example, the odds read 1 : 0.011 when the maxima at Δf = 1/2 GHz and 5.68 GHz are related. Furthermore, taking the values for the maximum at Δf = 1/2 GHz to evaluate ln⁡p×K at Δf = 1/3 GHz and 1/4 GHz gives the odds 1 : 1.013 : 1.018. Thus, the global maximum locates somewhere in the range below Δf≤ 1/4 GHz. But this range cannot be investigated in more detail from numerical point of view. However, the increase in ln⁡p×K when Δf is decreased below 1/2 GHz is interpreted as confirmation that a continuous spectrum is indeed probed.Table 1
Characterisation of maxima ofp×K dependent on Δf. The maxima locate at similar nine-dimensional parameter sets. The odds are evaluated by the ratio p×K(Δf)/p×K(Δf = 1/2 GHz) and relate the maxima globally. Since the odds are rising, the global maximum locates below Δf = 1/2 GHz.
Δ
f (GHz)
≈5.68
4
3
2
1
1/2
f
L
,
P
o (GHz)
33.499
33.243
33.277
32.817
32.445
32.425
f
U
,
P
o (GHz)
913.572
913.243
912.277
912.817
913.445
913.425
a
E
,
P
o (×10-12 V2)
0.9152
1.0127
1.0623
1.1002
1.1224
1.1299
a
O
,
P
o (×10-15 V2)
5.7381
5.7731
5.7776
5.797
5.8126
5.8145
x
0
,
P
o (mm)
−5.115293
−5.115297
−5.115288
−5.115291
−5.115295
−5.115295
δ
x
P
o (μm)
−4.814
−4.814
−4.814
−4.814
−4.814
−4.814
σ
b
0
,
P
r
,
P
o (nV)
62.150
62.148
62.146
62.143
62.142
62.142
σ
b
1
,
P
r
,
P
o (μV/m)
6.031
6.003
6.003
6.000
5.998
5.998
σ
b
2
,
P
r
,
P
o (μV/m2)
11.541
9.884
9.913
9.725
9.575
9.573
ln
⁡
p
×
K
1052035.05
1052037.54
1052038.52
1052039.15
1052039.48
1052039.55
Odds
0.011
0.134
0.357
0.668
0.929
1Figure 7
Odds for maxima ofp×K dependent on spectral increment Δf. The maximum at Δf = 1/2 GHz is used as a reference (unity odd). For Δf = 1/3 GHz and 1/4 GHz, the maximising parameter set (see Table 1) for Δf = 1/2 GHz is used to evaluate the associated odds. The global maximum of p×K is obtained for very small spectral increments.ForΔf = 5.68 GHz, the odds read 1 : 10−1390 when the associated maximum (ln⁡p×K≈1.0520×106) is related to the nonoptimised case (ln⁡p×K≈1.0488×106) investigated in Section 5.3.1. This is caused by having chosen some parameters like aE(O) several orders of magnitude too large with respect to the maximum settings values.ForΔf = 1 GHz, Figures 8(a)–8(f) summarise scans in the parameter pairs (aE,aO), (fL,fU), (x0,δx), (σb0,Pr,σb1,Pr), (σb0,Pr,σb2,Pr) and in σb0,Pr while the remaining parameters are held at the maximum values. In σb0,Pr, σb1,Pr, and σb2,Pr, a skewed distribution is found. The remaining six-dimensional posterior distribution has a high probability in a narrow region for a given Δf. This distribution is well approximated by a multivariate normal. Its mean is given by the maximum values listed in Table 1. The posterior covariance is estimated from the inverse of the Hessian matrix evaluated numerically via the second-order partial derivatives in the vicinity of the maximum. The off-diagonal elements of this covariance are negligible, so that one can factorise the posterior as a product of individual Gaussians. The posterior standard deviations σfL,Po, σfU,Po, σaE,Po, σaO,Po, σx0,Po, and σδx,Po vary little when Δf is changed (see Table 2). The spectral boundaries are well determined within an interval of some GHz. While the uncertainty in the scaling aE is small compared to its posterior mean value, aO is more uncertain. The uncertainty in the zero-path difference is of the order of some μm, and the shift δx is quite certain within some hundreds of nanometers.Table 2
Posterior standard deviations dependent onΔf for spectral boundaries fL and fU, scalings aE and aO of prior covariances, zero-path difference x0, and shift δx. The standard deviations are estimated from the inverse Hessian matrix at the corresponding maximum of p×K.
Δ
f (GHz)
≈5.68
4
3
2
1
1/2
σ
f
L
,
P
o (GHz)
1.742
1.795
1.746
1.762
1.754
1.718
σ
f
U
,
P
o (GHz)
3.101
3.155
2.927
2.869
2.858
2.810
σ
a
E
,
P
o (×10-12 V2)
0.1064
0.1179
0.1242
0.1290
0.1316
0.1327
σ
a
O
,
P
o (×10-15 V2)
1.5912
1.6064
1.6112
1.6164
1.6199
1.6211
σ
x
0
,
P
o (μm)
1.207
1.207
1.207
1.206
1.205
1.205
σ
δ
x
,
P
o (μm)
0.170
0.170
0.170
0.170
0.170
0.170Figure 8
p
×
K scanned in individual nonlinear parameters and hyperparameters in vicinity of maximum of p×K for Δf = 1 GHz (see Table 1). (a) Zero-path difference position x0 and shift δx. (b) Prior standard deviation σb0,Pr for absolute coefficient. (c) Lower and upper limits fL and fU of spectral domain. (d) σb0,Pr and prior standard deviation σb1,Pr for linear coefficient. (e) Scaling aE and scaling aO of Brownian bridge priors for S→E and S→O. (f) σb0,Pr and prior standard deviation σb2,Pr for quadratic coefficient.


#### 5.3.3. Conditional Amplitude Posterior for Maximising Settings

In the following, the amplitude posterior is investigated for the maximising settings givenΔf = 1/2 GHz (see Table 1). For the listed settings, the double-sided domain is identified by |x|≤|x0,Po|≈ 5.115 mm; the single-sided domain is bounded by the lower limit x>|x0,Po| and the upper limit of about 26.36 mm, the centre frequency fC,Po, and the bandwidth ΔF,Po of the spectral domain read 472.925 GHz and 881 GHz, respectively. In addition, the conversion of the spatial coordinate to the order of the Fourier coefficients k=|x|ΔF,Po/c reveals that the double- and single-sided domain contain information up to the 15th and 77th order (see Figures 9(c) and 9(d)), respectively. From this follows that the interferometric part of the 256/788 data values in the double-sided/full domain are modelled best by 62/310 Fourier coefficients and their associated basis functions dependent on fC,Po, ΔF,Po, and x. The signal envelopes σD,E and σD,O are evaluated with aE,Po and aO,Po.Figure 9
(a) Posterior meansS→E,Po and S→O,Po with associated 2-σ uncertainty bands evaluated at maximum of settings posterior for Δf = 1/2 GHz (see Table 1). Most likely, the spectral quantities S→E,Po and S→O,Po are finite only for the domain between fL,Po = 32.425 GHz and fU,Po = 913.425 GHz. This upper limit is approximately 4 smaller than the Nyquist frequency fNy = 3747.4 GHz. The dips in S→E,Po at 557 and 752 GHz are caused by absorption due to atmospheric water vapour. (b) D→E,Po, D→O,Po, and D→Off,Po obtained by mapping posterior means S→E,Po, S→O,Po, b0,Po=-62.142 nV, b1,Po = 5.997 μV/m, and b2,Po=-8.641 μV/m2 to data domain. (c) Cosine transform D→E,Po of S→E,Po dependent on optical path difference x and order k=|x|ΔF,Po/c (bandwidth ΔF,Po=fU,Po-fL,Po) of Fourier coefficients. The envelope σD,E∝aE,Po1/2/k, incorporating the scaling aE,Po = 1.1299 × 10−12 V2, bounds D→E,Po best. (d) Zoom of (b) dependent on x and k. Globally, the sine transform D→O,Po of S→O,Po is best summarised by the envelope σD,O∝aO,Po1/2/k with the scaling aO,Po=5.8145×10-15 V2. Since σD,E≫σD,O holds, the single-sided domain is determined mainly by S→E,Po. (e) Histogram of residuals, that is, difference between measured data D→ and posterior mean prediction D→Po=D→E,Po+D→O,Po+D→Off,Po scaled by noise level σD. The data is well described, because the mean -1.7×10-4 and the variance 0.942 of the residuals are almost vanishing and close to unity, respectively. (f) Residuals versus optical path difference.
(a)
(b)
(c)
(d)
(e)
(f)(a) Posterior Mean. Given the boundaries and increment, the spectral dimension reads Nf = 1762, and, thus, the amplitude posterior mean l→Po=S→E,Po,S→O,Po,b→PoT is 3527-dimensional. The posterior mean values for the three coefficients b0,Po=-62.142 nV, b1,Po = 5.997 μV/m, and b2,Po=-8.641 μV/m2 have magnitudes which are close to the corresponding prior standard deviations at the maximum of the settings posterior (see Table 1). Furthermore, the absolute and the linear mean values are similar to the ones obtained for the nonoptimised case (see Section 5.3.1). The offset D→Off,Po which follows from mapping b→Po is shown in Figures 9(b) and 9(d).The posterior meansS→E,Po and S→O,Po presented in Figure 9(a) approach zero towards the spectral boundaries as expected. Furthermore, S→E,Po is much larger in amplitude than S→O,Po anticipated by the quite different posterior mean scalings aE,Po and aO,Po. Hence, the even process describes most of the interferometric data in the double-sided domain. This is seen by mapping both spectral means to the data domain which gives the even and odd quantities D→E,Po and D→O,Po (see Figures 9(b)–9(d)). In addition, the single-sided domain is described mainly by S→E,Po, because in this domain only the envelope σD,E for the even process is of the order of the interferometric data (see Figure 9(d)). This differs from the findings for the nonoptimised case (aE=aO) for which S→E,Po and S→O,Po determine the single-sided region almost equally.The histogram of the residuals (see Figure9(e)) is approximated very well by the normal distribution N(-1.7×10-4,0.942). Since the mean vanishes almost, and the standard deviation is very close to unity, the data set is well described by the model and the posterior means. Most likely, the data point located at x≈ 7 mm (see Figure 9(f)) is an outlier, because its residual is outside the 3.5-σD band.(b) Posterior Covariance. For the linear parameters the posterior covariance matrix is written like(61)ΣPo=ΣS→E,E,PoΣS→E,O,PoΣS→E,b→,PoΣS→E,O,PoTΣS→O,O,PoΣS→O,b→,PoΣS→E,b→,PoTΣS→O,b→,PoTΣb→,b→,Po,and taking the square root of an element of the main diagonal gives the posterior standard deviation σPi for the ith parameter. The correlation coefficient between two parameters P1 and P2 can be evaluated by ρP1,P2=ΣP1,P2,Po/(σP1σP2).The standard deviationsσb0,Po = 0.501 nV, σb1,Po = 0.093 μV/m, and σb2,Po = 4.119 μV/m2 of the three coefficients are reduced by about two orders of magnitude with respect to the uncertainties obtained for the nonoptimised case. Though, the quadratic coefficient is inferred with small confidence. Coefficients of neighbouring degree have the highest correlation in magnitude and are anticorrelated. For example, one finds ρb0,b1=-0.55 and ρb1,b2=-0.93, while ρb0,b2 = 0.29 remains small.The covariancesΣS→E,E,Po(f,f′) and ΣS→O,O,Po(f,f′) and their associated correlations ρE,E(f,f′) and ρO,O(f,f′) resemble sinc functions centred along the line f=f′ plus a tip at this condition (see Figures 10(a)–10(c) for cross section at f′ = 100.175 GHz). Though, the FWHM widths of 1.5 GHz and 10 GHz are very different for the even and odd covariances/correlations. Regarding the covariances, the tip heights account for 1.3×10-38 V2/Hz2 and 4.4×10-40 V2/Hz2, and towards the spectral boundaries the heights are approaching zero which can be seen by the standard deviations σS→E,Po=ΣS→E,E,Po1/2(f,f) and σS→O,Po=ΣS→O,O,Po1/2(f,f) (see Figure 10(d)).Figure 10
Characterisation of posterior covariance for amplitude parameters obtained at maximum of settings posterior forΔf = 1/2 GHz (see Table 1). (a) Cross section of covariances ΣS→E,E,Po and ΣS→O,O,Po for even and odd function. Each covariance is similar to a sinc function centred at f=f′, where a tip is added on top. (b) Cross section of covariances ΣS→O,O,Po and ΣS→E,O,Po. At a fixed frequency, S→E and S→O are independent, because ΣS→E,O,Po = 0. But finite correlation remains with the neighbouring spectral domain. (c) Correlation. The even (odd) parameters are correlated (ρE(O),E(O)>0.5) over the spectral scale of about 1.5 GHz (10 GHz), while the correlation between even and odd parameters is small and asymmetric with respect to f=f′. (d) Standard deviations σS→E,Po and σS→O,Po. At the spectral limits the uncertainties approach zero. (e) and (f) Correlations of S→E(O) and b→. The correlations are oscillating and decrease generally with the frequency.
(a)
(b)
(c)
(d)
(e)
(f)At the conditionf=f′, the covariance ΣS→E,O,Po and the correlation ρE,O vanish and are asymmetrically close to this line (see Figures 10(b) and 10(c) for cross section at f′ = 100.175 GHz). This means that at a given frequency f the amplitudes S→E and S→O are independent but have a finite and opposite correlation with the neighbouring spectral domains to either side.The correlationsρE,b→ and ρO,b→ are oscillatory and rise towards the spectral origin (see Figures 10(e) and 10(f)).(c) Posterior Samples. From the conditional amplitude posterior samples S→E,Po,Sa, S→O,Po,Sa and b→Po,Sa are drawn. While Figure 11(a) shows one sample for each of the even and odd spectral functions, Figure 11(b) presents 100 samples. The samples form a band around the corresponding posterior mean S→E(O),Po with the width of about twice the standard deviation σS→E(O),Po as shown in Figure 10(d). Hence, the band for the odd function is much smaller.Figure 11
(a) Individual posterior sample functionsS→E,Po,Sa and S→O,Po,Sa. (b) 100 posterior sample functions S→E,Po,Sa and S→O,Po,Sa. The uncertainty at a given frequency is larger for the even function, because for the associated scalings of the two processes aE,Po≫aO,Po holds. (c) 100 spectral samples and 100 offset coefficients mapped to probed data domain. Unlike for the nonoptimised case (see Figure 5(b)), the spread is much smaller, and the quantities connect smoothly at the transition between double- and single-sided regions. (d) 100 samples D→Po,Sa and measured data D→. The samples are in agreement with the uncertainty of the measurements. (e) One sample D→E,Po,Sa and D→O,Po,Sa evaluated inside probed domain (x < 26.36 mm) and outside. With respect to the envelopes σD,E(O), the contributions behave similarly in both domains. (f) 100 samples D→E,Po,Sa and D→O,Po,Sa evaluated inside the probed domain and outside. Outside the probed domain, the area between the envelops σD,E(O) is filled by D→E,Po,Sa and D→O,Po,Sa. Because the spread in the nonprobed domain is wider than the one in the probed domain, the main uncertainty in the even and odd spectra originates in nonprobed Fourier coefficients.
(a)
(b)
(c)
(d)
(e)
(f)The samples mapped to the data domain giveD→E,Po,Sa, D→O,Po,Sa, and D→Off,Po,Sa (see Figure 11(c)) which form much more narrow bands when compared to the nonoptimised case (see Figure 5(b)). In particular, the transition from the double- to the single-sided domain is smooth. This is a consequence of having used the most likely scalings aE,Po and aO,Po and the associated but quite different signal envelops inferred from the data. Furthermore, this explains why the band for D→O,Po,Sa is much smaller than the one for D→E,Po,Sa.For each sample, the sumD→Po,Sa=D→E,Po,Sa+D→O,Po,Sa+D→Off,Po,Sa complies with the measured data D→ and the noise level (see Figure 11(d)).With the samplesS→E,Po,Sa and S→O,Po,Sa at hand, D→E,Po,Sa and D→O,Po,Sa are predicted outside the probed spatial domain. Up to the optical path difference of 70 mm, Figures 11(e) and 11(f) show D→E,Po,Sa and D→O,Po,Sa for one sample and 100 samples, respectively. No difference is obvious in D→E,Po,Sa and D→O,Po,Sa for an individual prediction in the probed and nonprobed domain. However, many samples take values between σD,E(O) outside the probed region. Since the spread in the probed domain is smaller than in the nonprobed domain, the main uncertainty in S→E,Po and S→O,Po originates in Fourier coefficients not probed rather than the noise on measured data.


#### 5.3.4. Marginal Amplitude Posterior

Since only the marginal posteriorpS→E,S→O∣Δf,D→ of the spectral functions is of interest for this specific Fourier spectroscopy problem, formally the marginalisation(62)pS→E,S→O∣Δf,D→=∫pS→E,S→O,b→∣n→,h→,Δf,D→pn→,h→∣Δf,D→db→dn→dh→needs to be carried out. In general, the pdf pS→E,S→O∣Δf,D→ is not Gaussian but can be approximated by a multivariate normal. Performing the marginalisation, that is, the uncertainty propagation, is analytically not possible but achievable numerically described below for Δf = 1/2 GHz.To reduce the numerical efforts, only the most important parametersfL, fU, x0, δx, aE, and aO are taken into consideration for the marginalisation, and the remaining three parameters σb0,Pr, σb1,Pr, and σb2,Pr are held at the maximum values listed in Table 1. The remaining pdf, labeled by p′n→,h→∣Δf,D→ in the following, is well approximated by a factorisable multivariate normal, because the corresponding covariances are negligible (see Section 5.3.2). The posterior mean values and standard deviations are listed in Tables 1 and 2, respectively.From the six-dimensional Gaussian posteriorp′n→,h→∣Δf,D→, 1000 samples are drawn. For each sample, mean and covariance for the conditional amplitude posterior are evaluated. Thereby, a difficulty arises due to the lower and upper spectral limits fL and fU. For given increment Δf, fL and fU determine the spectral grid. To keep this grid fixed, the spectral grid obtained for fL,Po = 32.425 GHz and fU,Po = 913.425 GHz is extended by some tens of GHz on either side. Furthermore, the posterior samples for fL and fU are compared with the extended but fixed spectral grid points and changed to the closest grid point. With the modified boundaries and the remaining samples unchanged, a conditional amplitude posterior follows. For example, Figure 12(a) shows the mean functions S→E,Po and S→O,Po for the 1000 posterior samples of the settings. While S→E,Po changes little but most at the spectral boundaries, S→O,Po varies considerably, mostly in the centre of the spectral domain covered. This is caused mainly by the posterior uncertainty of the zero-path difference x0. However, the spread is symmetric around the most likely mean.Figure 12
(a) MeansS→E,Po (blue) and S→O,Po (red) obtained for 1000 samples of settings posterior at Δf = 1/2 GHz. Only at the spectral limits the means deviate from the most likely mean S→E,Po (cyan) evaluated at the maximum of the settings posterior. The uncertainty in the zero-path difference x0 causes a considerable deviation from the most likely mean S→O,Po (black). (b) Means S→E,Po,Ma (cyan) and S→O,Po,Ma (black) after performing numerical marginalisation of settings parameters at Δf = 1/2 GHz. The marginalised means are similar to the means S→E,Po and S→O,Po (see (a)) at the maximum of the settings posterior except for minor differences towards the spectral boundaries. Having the marginalised covariance at hand, 100 samples S→E,Po,Ma,Sa (blue) and S→O,Po,Ma,Sa (red) are drawn from the marginalised amplitude posterior. Only in the centre of the spectral domain the spread of S→O,Po,Ma,Sa increases to capture the uncertainty in x0 as well (see (a)).
(a)
(b)For each settings sample, 100 samples are drawn frompS→E,S→O,b→∣n→,h→,Δf,D→. Overall, 100000 samples are available for b→, S→E, and S→O at each spectral grid point. From this large sample set, the marginalised means S→E,Po,Ma, S→O,Po,Ma and covariance ΣS→,Po,Ma are evaluated which capture the first and second moments of the pdf pS→E,S→O∣Δf,D→. In doing so, the marginalisation with respect to b→ is achieved implicitly.(a) Posterior Mean. The means S→E,Po,Ma and S→O,Po,Ma (see Figure 12(b)) look similar to S→E,Po, S→O,Po at the maximum of the settings posterior. However, S→E,Po,Ma and S→O,Po,Ma approach zero some GHz below and above fL,Po and fU,Po, respectively.(b) Posterior Covariance. The posterior covariance after marginalisation with respect to the setting parameters is decomposed like(63)ΣS→,Po,Ma=ΣS→E,E,Po,MaΣS→E,O,Po,MaΣS→E,O,Po,MaTΣS→O,O,Po,Ma.ΣS→E,E,Po,Ma (see Figure 13(a) for cross section at f′ = 100.175 GHz) and is very similar to ΣS→E,E,Po evaluated at the maximum of the settings posterior (see Section 5.3.3). The same holds for ΣS→E,O,Po,Ma and ΣS→E,O,Po. However, a significant change is found for ΣS→O,O,Po,Ma (see Figures 13(a) and 13(b) for some cross sections) when compared to ΣS→O,O,Po. Due to the posterior uncertainty in x0, a broadband additive contribution increases ΣS→O,O,Po,Ma up to about one order of magnitude in the centre of the spectral domain. Comparing the standard deviations σS→E,Po,Ma with σS→E,Po, and σS→O,Po,Ma with σS→O,Po (see Figure 13(c)) shows that towards the spectral limits the propagation of the settings posterior uncertainties has little effect. Furthermore, the correlation ρO,O,Ma is high for a broad spectral range in the central region (see Figure 13(d)).Figure 13
Characterisation of posterior covarianceΣS→,Po,Ma numerically marginalised with respect to settings parameters. ΣS→,Po,Ma is decomposed, so that ΣS→E(O),E(O),Po,Ma captures the covariance of the even (odd) spectral quantities. (a) Cross section at frequency f′ = 100.175 GHz. Apart from the extension in the spectral domain, ΣS→E,E,Po,Ma is similar to the covariance ΣS→E,E,Po evaluated at the maximum of the settings posterior (see Figure 10(a)). ΣS→O,O,Po,Ma shows a significant deviation from ΣS→O,O,Po (see also (b)). (b) Cross sections of ΣS→O,O,Po,Ma at frequencies f′ = 100.175 GHz, 400.175 GHz, and 850.175 GHz. Towards the spectral limits ΣS→O,O,Po,Ma behaves like ΣS→O,O,Po. But in the central region a broadband additive contribution is evident originating in the posterior uncertainty of the zero-path difference x0. (c) Standard deviations evaluated from main diagonal of marginalised covariance ΣS→,Po,Ma (red) and covariance ΣS→,Po (black, see Figure 10(d)) evaluated at maximum of settings posterior. The odd quantity increases by a factor of three in the central region after the marginalisation. (d) Correlation ρO,O,Ma at three frequencies. Compared to ρO,O (see Figure 10(c)), a broadband (several 100s GHz) elevation of ρO,O,Ma up to 0.8 follows from the uncertainty in x0 for the central domain.
(a)
(b)
(c)
(d)(c) Posterior Samples. From the marginalised Gaussian posterior pS→E,S→O∣Δf,D→ 100 samples (S→E,Po,Ma,Sa and S→O,Po,Ma,Sa) are drawn and shown in Figure 12(b). As expected, the spread of S→E,Po,Ma,Sa around the mean S→E,Po,Ma is almost equal to the one for S→E,Po,Sa around the mean S→E,Po (see Figure 11(b)). The samples S→O,Po,Ma,Sa have a wider spread in the centre of the spectral domain than S→O,Po,Sa mainly caused by the uncertainty in the zero-path difference.


### 5.4. Figure of Merit for Real-World Interferometer

For the ideal Martin-Puplett interferometer, the odd spectral function, and, hence, the odd process must vanish from theoretical point of view. However, imperfections of a real-world interferometer leave the odd contribution finite in general. This imperfection is captured by the scalingaO and in terms of signal by the envelope σD,O∝aO1/2. By relating the square root of the scalings like(64)FS=aE1/2aE1/2+aO1/2,one can define the figure of merit FS which expresses by a number the signal deviation of a real-world instrument from the ideal case. For an ideal interferometer, FS = 1 holds (aO=0). For the interferometer investigated here, the settings posterior (see Section 5.3.2) carries the information to state the mean as FS = 0.933 with the uncertainty of about 0.01. Hence, about 7% of the signal is converted from the ideal even process to the odd process by the real-world diagnostic.


## 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 D→1 is acquired in terms of Volts dependent on the optical path difference x. 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 D(xi)=D2(xi)-D1(xi) acquired at the spatial grid node xi. Then, the abstract model for the Martin-Puplett interferometer is stated by (46)Dxi=A22∫fLfUSEcos⁡2πfcxi+SOsin⁡2πfcxidf+VOffxi+Ni0,σD2,using the total amplification A of the detection system. Furthermore, the offset VOff marks a diagnostic imperfection which varies with xi. The Gaussian noise contributes to each data sample described by Ni(0,σD2). The unknown quantities in the diagnostic model are the spatial grid xi, the lower and upper spectral boundaries fL and fU of the Fourier transform integral, the even and odd functions SE and SO dependent on frequency f, and the offset.


### 5.1.2. Interferometric Data

The data setD→, that is, the difference interferogram consists of ND = 788 values (see Figure 3(a)). Merely for graphical presentation a certain x 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 interferogramD→ (black) versus optical path difference x. The data set includes ND = 788 values, and the noise level for each Di reads σD = 132.29 μV. The values of the spatial grid are determined by the standard model xi=i-1Δx+x0 with i∈1,ND, the constant increment Δx = 40 μm, and the origin x0. This origin is estimated as −5.118 mm by a fit of a quadratic polynomial to the maximum max⁡D→ and its neighbouring values to either side. Clearly, D→ deviates from the baseline (red). (b) Difference quotient for different subsets of D→. Since the zig-zag-pattern in the point-to-point variation appears only for the first case (black), the standard model xi+1-xi = const. seems to be inadequate. Instead, a model seems appropriate for which x2i+1-x2i-1 = x2i+2-x2i=const. holds.
(a)
(b)The componentsDi are measured independently on each other, and the noise level for each Di is captured by σD = 132.29 μV, and, thus, the variance of the whole data vector is stated by the matrix ΣD=σD2δij.


### 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Δx = 40 μm. Hence, the standard model xi=(i-1)Δx+x0 is obvious with i=1,…,ND, and the zero-path difference x0 is a free parameter. However, this model is not accurate. Applying the standard approach which fits a second-order polynomial to the maximum max⁡D→ and its two nearby values, x0=-127.95Δx=-5.118 mm, is inferred for the data set shown in Figure 3(a). Furthermore, the difference quotient evaluated by (Di+1-Di)/Δx is presented versus the optical path difference x=(xi+1+xi)/2 in Figure 3(b). The point-to-point variation of the quotient has a zig-zag pattern which implies that the assumption Δx = const. is incorrect. Indeed, each of two the difference quotients (D2i+1-D2i-1)/(2Δx) and (D2i+2-D2i)/(2Δx) is smoother. Hence, the model(47)xi=i-1Δx+x0+0:ifiis  oddδx:ifiis  evenseems more appropriate, making use of two free parameters: the zero-path difference x0 and a shift δx for every other grid value. The priors p(x0) and p(δx) are set to be uniform.


### 5.1.4. Offset

The upward trend ofD→ (see Figure 3(a)) is modelled by the offset VOff(xi). Here, a second-order polynomial is chosen to capture the offset by VOff(xi)=b0+b1xi+b2xi2 with the free parameters b0, b1, and b2 being summarised by the vector b→=b0,b1,b2T. The corresponding mapping to the data domain can be expressed as(48)MOffb→=A21x1x121x2x22⋮⋮⋮1xNDxND2b0b1b2with the ND×3 matrix MOff. The joint prior is expressed by the factorisable multivariate normal distribution:(49)pb→=∏i=02N0,σbi,Pr2=N0,Σb→,Pr,where σbi,Pr 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 f  ∈  ]fL,fU[. Thereby, fL and fU are free parameters. The spectral domain is dicretised using the constant increment Δf which is considered as a free parameter as well. Then, the spectral domain is represented by the set fj=fL+(j-0.5)Δf with j=1,…,Nf, covering the band NfΔf=(fU-fL)=ΔF centred at fC=(fU+fL)/2. To be clear, since fL, fU, and Δf are free parameters and will be inferred, the number Nf is variable. This number determines the dimensionality of the vectors S→E and S→O which represent the discretised functions SE(f) and SO(f). The mapping of S→E and S→O to the data domain in (46) is written as(50)A22∫fLfUSEcos⁡2πfcxi+SOsin⁡2πfcxidf≈ME,MOS→ES→Owith the two ND×Nf matrices(51)ME=Acos⁡2πf1cx1⋯cos⁡2πfNfcx1⋮⋱⋮cos⁡2πf1cxND⋯cos⁡2πfNfcxNDΔf,MO=Asin⁡2πf1cx1⋯sin⁡2πfNfcx1⋮⋱⋮sin⁡2πf1cxND⋯sin⁡2πfNfcxNDΔf.The joint prior for the two vectorial quantities S→E and S→O is factorised, and each prior is chosen as a multivariate normal distribution with vanishing mean. Since the Brownian bridge covariance ΣBB (see Section 4) describes functions which vanish at the boundaries fL and fU and are square-integrable, and its signal envelope decays with the optical path difference like the interferometric data at hand, the priors are chosen by(52)pS→EO∣Δf,fL,fU,aEO=N0,aEOΣBBΔf,fL,fUwith the two hyperparameters aE and aO for scaling. For each of these hyperparameters a uniform prior paEO=1/ΔaE(O) is applied.Finally, the joint prior(53)pΔf,fL,fU=1KΔf,fL,fUshould be constant, so that any combination of Δf, fL, and fU has the same probability. Furthermore, the conditions 0≤Δf≤fU-fL, 0≤fL≤fU, and 0≤fU≤ΔfU (global upper limit) must be fulfilled. For example, the upper limit of fU is set here to the Nyquist frequency ΔfU=fNy=c/(2Δx) = 3747.4 GHz.


## 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 D→1 is acquired in terms of Volts dependent on the optical path difference x. 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 D(xi)=D2(xi)-D1(xi) acquired at the spatial grid node xi. Then, the abstract model for the Martin-Puplett interferometer is stated by (46)Dxi=A22∫fLfUSEcos⁡2πfcxi+SOsin⁡2πfcxidf+VOffxi+Ni0,σD2,using the total amplification A of the detection system. Furthermore, the offset VOff marks a diagnostic imperfection which varies with xi. The Gaussian noise contributes to each data sample described by Ni(0,σD2). The unknown quantities in the diagnostic model are the spatial grid xi, the lower and upper spectral boundaries fL and fU of the Fourier transform integral, the even and odd functions SE and SO dependent on frequency f, and the offset.


## 5.1.2. Interferometric Data

The data setD→, that is, the difference interferogram consists of ND = 788 values (see Figure 3(a)). Merely for graphical presentation a certain x 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 interferogramD→ (black) versus optical path difference x. The data set includes ND = 788 values, and the noise level for each Di reads σD = 132.29 μV. The values of the spatial grid are determined by the standard model xi=i-1Δx+x0 with i∈1,ND, the constant increment Δx = 40 μm, and the origin x0. This origin is estimated as −5.118 mm by a fit of a quadratic polynomial to the maximum max⁡D→ and its neighbouring values to either side. Clearly, D→ deviates from the baseline (red). (b) Difference quotient for different subsets of D→. Since the zig-zag-pattern in the point-to-point variation appears only for the first case (black), the standard model xi+1-xi = const. seems to be inadequate. Instead, a model seems appropriate for which x2i+1-x2i-1 = x2i+2-x2i=const. holds.
(a)
(b)The componentsDi are measured independently on each other, and the noise level for each Di is captured by σD = 132.29 μV, and, thus, the variance of the whole data vector is stated by the matrix ΣD=σD2δij.


## 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Δx = 40 μm. Hence, the standard model xi=(i-1)Δx+x0 is obvious with i=1,…,ND, and the zero-path difference x0 is a free parameter. However, this model is not accurate. Applying the standard approach which fits a second-order polynomial to the maximum max⁡D→ and its two nearby values, x0=-127.95Δx=-5.118 mm, is inferred for the data set shown in Figure 3(a). Furthermore, the difference quotient evaluated by (Di+1-Di)/Δx is presented versus the optical path difference x=(xi+1+xi)/2 in Figure 3(b). The point-to-point variation of the quotient has a zig-zag pattern which implies that the assumption Δx = const. is incorrect. Indeed, each of two the difference quotients (D2i+1-D2i-1)/(2Δx) and (D2i+2-D2i)/(2Δx) is smoother. Hence, the model(47)xi=i-1Δx+x0+0:ifiis  oddδx:ifiis  evenseems more appropriate, making use of two free parameters: the zero-path difference x0 and a shift δx for every other grid value. The priors p(x0) and p(δx) are set to be uniform.


## 5.1.4. Offset

The upward trend ofD→ (see Figure 3(a)) is modelled by the offset VOff(xi). Here, a second-order polynomial is chosen to capture the offset by VOff(xi)=b0+b1xi+b2xi2 with the free parameters b0, b1, and b2 being summarised by the vector b→=b0,b1,b2T. The corresponding mapping to the data domain can be expressed as(48)MOffb→=A21x1x121x2x22⋮⋮⋮1xNDxND2b0b1b2with the ND×3 matrix MOff. The joint prior is expressed by the factorisable multivariate normal distribution:(49)pb→=∏i=02N0,σbi,Pr2=N0,Σb→,Pr,where σbi,Pr 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 f  ∈  ]fL,fU[. Thereby, fL and fU are free parameters. The spectral domain is dicretised using the constant increment Δf which is considered as a free parameter as well. Then, the spectral domain is represented by the set fj=fL+(j-0.5)Δf with j=1,…,Nf, covering the band NfΔf=(fU-fL)=ΔF centred at fC=(fU+fL)/2. To be clear, since fL, fU, and Δf are free parameters and will be inferred, the number Nf is variable. This number determines the dimensionality of the vectors S→E and S→O which represent the discretised functions SE(f) and SO(f). The mapping of S→E and S→O to the data domain in (46) is written as(50)A22∫fLfUSEcos⁡2πfcxi+SOsin⁡2πfcxidf≈ME,MOS→ES→Owith the two ND×Nf matrices(51)ME=Acos⁡2πf1cx1⋯cos⁡2πfNfcx1⋮⋱⋮cos⁡2πf1cxND⋯cos⁡2πfNfcxNDΔf,MO=Asin⁡2πf1cx1⋯sin⁡2πfNfcx1⋮⋱⋮sin⁡2πf1cxND⋯sin⁡2πfNfcxNDΔf.The joint prior for the two vectorial quantities S→E and S→O is factorised, and each prior is chosen as a multivariate normal distribution with vanishing mean. Since the Brownian bridge covariance ΣBB (see Section 4) describes functions which vanish at the boundaries fL and fU and are square-integrable, and its signal envelope decays with the optical path difference like the interferometric data at hand, the priors are chosen by(52)pS→EO∣Δf,fL,fU,aEO=N0,aEOΣBBΔf,fL,fUwith the two hyperparameters aE and aO for scaling. For each of these hyperparameters a uniform prior paEO=1/ΔaE(O) is applied.Finally, the joint prior(53)pΔf,fL,fU=1KΔf,fL,fUshould be constant, so that any combination of Δf, fL, and fU has the same probability. Furthermore, the conditions 0≤Δf≤fU-fL, 0≤fL≤fU, and 0≤fU≤ΔfU (global upper limit) must be fulfilled. For example, the upper limit of fU is set here to the Nyquist frequency ΔfU=fNy=c/(2Δx) = 3747.4 GHz.


## 5.2. Bayes’ Theorem

In the following, the linear parameters are summerised by the setl→=S→E,S→O,b→T, the nonlinear parameters by n→=fL,fU,x0,δxT, and the hyperparameters by h→=aE,aO,σb0,Pr,σb1,Pr,σb2,PrT.The matrixM(n→,Δf)=(ME,MO,MOff) maps the parameters l→ to the data domain. Hence, the Gaussian likelihood is written as(54)pD→∣l→,n→,Δf=exp⁡-1/2D→-Ml→TΣD-1D→-Ml→2πND/2ΣD1/2.The prior for the full problem takes the form(55)pl→,Δf,n→,h→=pl→∣Δf,fL,fU,aE,aO,σb0,Pr,σb1,Pr,σb2,PrpΔf,fL,fUpaEpaOpx0pδxpσb0,Prpσb1,Prpσb2,Prwith the multivariate normal prior(56)pl→∣Δf,fL,fU,aE,aO,σb→,Pr=exp⁡-1/2l→TΣPr-1l→2π2Nf+3/2ΣPr1/2for l→, using the (2Nf+3)×(2Nf+3) dimensional covariance matrix(57)ΣPr=aEΣBB000aOΣBB000Σb→,Pr.One gets the joint posterior(58)pl→,n→,h→,Δf∣D→=pl→∣n→,h→,Δf,D→pn→,h→,Δf∣D→with the conditional amplitude posterior given by the multivariate normal(59)pl→∣n→,h→,Δf,D→=exp⁡-1/2l→-l→PoTΣPo-1l→-l→Po2π2Nf+3/2ΣPo1/2,using the posterior covariance matrix ΣPo=MTΣD-1M+ΣPr-1-1 and the mean l→Po=ΣPoMTΣD-1D→. Furthermore, one obtains the settings posterior(60)pn→,h→,Δf∣D→=p×KK=ΣPo1/2ΣPr1/2exp⁡1/2l→PoTΣPo-1l→PoK,where the constant K 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 forx0, δx, fL, fU, Δf, aE and aO, σb0,Pr, σb1,Pr, and σb2,Pr.The values forx0=-127.95Δx=-5.118 mm and δx=0 mm are chosen to form an optical path difference grid xi for which all subsequent nodes are separated by Δx. Then, |xi|≤5.118 mm identifies the double-sided domain (256 data points), and the single-sided domain is marked by 5.118 mm ≤xi≤ 26.362 mm (532 data points).The spectral domain is covered fromfL = 0 GHz to the Nyquist frequency fU=c/(2Δx) = 3747.4 GHz (=ΔF). With the bandwidth and the optical path difference set, the conversion k=|x|ΔF/c (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Δf = 5.68 GHz is used, and the lower and upper spectral limits are set to fL = 0 GHz and fU = 3747.4 GHz (Nyquist assumptions), respectively. The scalings of the processes are chosen by aE=aO=10-8 V2. (a) Prior samples S→E,Pr,Sa and S→O,Pr,Sa 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 D→E(O),Pr,Sa for optical path difference obtained for x0=-5.118 mm and δx = 0. With increasing distance to the origin, the amplitude drops which is well bounded by the signal envelopes σD,E=σD,O for the chosen bandwidth ΔF=fU. Fourier coefficients up to order k = 64 and 330 settle in the double-sided domain x≤x0 and single-sided domain x>x0, respectively.
(a)
(b)To match the maximum order (330) of probed Fourier coefficient, the discretisation is set toΔf=c/(2Δx(256/2+532))≈5.68 GHz which reflects the classical increment. Then, the spectral grid has Nf = 660 elements.The hyperparameters are set to large values likeaE=aO=10-8 V2 to prevent from a prior determined posterior. This unwanted feature would be present, if one would choose aE(O) too small, so that the signal envelops do not include the measured data (see next paragraph).With the chosen spectral priors, sample functions/vectorsS→E,Pr,Sa and S→O,Pr,Sa are drawn and shown in Figure 4(a). The function values are of the order of 10-16 V/Hz decaying towards the spectral boundaries. The mapping of the spectral prior samples to the spatial domain via D→E,Pr,Sa+D→O,Pr,Sa=MES→E,Pr,Sa+MOS→O,Pr,Sa shows that D→E,Pr,Sa+D→O,Pr,Sa 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 D→E,Pr,Sa+D→O,Pr,Sa drops with the distance from the origin. All these data samples are well bounded by the signal envelopes σD,E(O) obtained by σD,E(O)=(A/2)σD,E(O)′ which is the scaled version of the expressions (44) and (45).Forb→ the prior covariance is determined by setting σb0,Pr=1 V, σb1,Pr=1 V/m and σb2,Pr=1 V/m2.(a) Posterior Mean. For the above values, the mean l→Po of the amplitude posterior is given by b→Po=(-61.41 nV, 5.75 μV/m, 2.42 μV/m2)T, S→E,Po and S→O,Po as presented in Figure 5(a). Clearly, S→E,Po peaks below 1000 GHz with the maximum of approximately 6 × 10−18 V/Hz around 400 GHz. Furthermore, S→O,Po is much smaller than S→E,Po below 1000 GHz.Figure 5
Investigation of conditional amplitude posteriorpl→∣Δf,fL,fU,aE,aO,x0,δx,σb0,Pr,σb1,Pr,σb2,Pr,D→ for linear parameters l→=S→E,S→O,b→T. The posterior is obtained for the settings x0=-5.118 mm and δx = 0, fL = 0 GHz, fU = 3747.4 GHz, Δf = 5.68 GHz, aE=aO=10-8 V2, σb0,Pr=1 V, σb1,Pr=1 V/m, and σb2,Pr=1 V/m2. (a) Spectral posterior means S→E(O),Po and one posterior sample. Below 1000 GHz, the means show a finite amplitude. Each sample deviates strongly from its corresponding mean. (b) Posterior mean l→Po and 100 samples l→Po,Sa mapped to give even, odd, and offset contributions in data domain. While in the double-sided region (|x|≤|x0|) most of the data described by S→E,Po, S→E,Po and S→O,Po give same contributions in the single-sided region (x>|x0|). In the single-sided region, the spread of the even (odd) contributions for the posterior samples is large but bounded by the signal envelops σD,E(O). (c) Residuals for l→Po (black) and one sample l→Po,Sa (red) both mapped to data domain. An almost perfect description of the noisy data is achieved by lPo, indicating an overfitting of the data. (d) Measured data D→ and 100 samples of amplitude posterior mapped to data domain. The individual contributions have a much wider spread (see (b)), but the total sum D→Po,Sa agrees well with D→ and its uncertainties. This is caused by the posterior covariance of l→ with suited off-diagonal entries.
(a)
(b)
(c)
(d)The offsetD→Off,Po=MOffb→Po contribution to the data has an upwards trend but is small (see Figure 5(b)). The even and odd contributions D→E(O),Po=ME(O)S→E(O),Po behave differently in the double- and single-sided domains (see Figure 5(b)). While in the double-sided domain D→E,Po is much larger than D→O,Po, D→E,Po and D→O,Po equal each other in the single-sided domain. This indicates an underestimation (overestimation) of the Fourier coefficients above the order 64 for SE,Po (SO,Po). The cause of this finding relies in the choice aE=aO (σD,E=σD,O) which becomes more clear in Section 5.3.3.Since the quantityD→Po=Ml→Po equals almost the data set D→ (see residuals in Figure 5(c)), a nearly perfect match is achieved. Hence, all noise contributions are captured by the mean l→Po, and an overfitting of the data is obvious.(b) Posterior Covariance. The posterior covariance ΣPo is not characterised in detail. From the elements of the main diagonal of ΣPo, one finds the values σb0,Po = 73.30 nV, σb1,Po = 12.58 μV/m, and σb2,Po = 446.67 μV/m2. These posterior uncertainties are quite large when compared to the values of the corresponding mean b→Po, especially for the quadratic coefficient.For the spectral quantities, the square root of the main diagonal elements ofΣPo is of the order of some 10-18 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 l→Po,Sa=S→E,Po,Sa,S→O,Po,Sa,b→Po,SaT are drawn from the amplitude posterior and mapped to the data domain (see Figure 5(b)). Thereby, the contributions D→E,Po,Sa=MES→E,Po,Sa, D→O,Po,Sa=MOS→O,Po,Sa, and D→Off,Po,Sa=MOffb→Po,Sa are split to investigate their interplay. D→Off,Po,Sa has a noticeable width but remains in the vicinity of D→Off,Po for the whole spatial domain.For the double-sided regionD→E,Po,Sa and D→O,Po,Sa form a narrow band around the corresponding posterior mean quantities D→E,Po and D→O,Po. On the contrary, for the single-sided region D→E,Po,Sa and D→O,Po,Sa increase drastically but are restricted by σD,E and σD,O. The deviation is large between the individual and the posterior mean quantities; though, the sum D→Po,Sa=D→E,Po,Sa+D→O,Po,Sa+D→Off,Po,Sa remains well inside the 2-σD band surrounding D→ (see Figures 5(c) and 5(d)). Hence, the posterior covariance ΣPo captures properly the correlations.With the specific values of the settings one gets the numberln⁡p×K≈1.0488×106.


### 5.3.2. Settings Posterior

As the problem is formulated, the settings posterior is proportional top×K. Its optimisation, that is, finding the global maximum of the settings posterior is a 10-dimensional problem. Since large numbers are involved, one has to investigate ln⁡p×K. In general, if only one parameter is changed, a well distinguishable peak in ln⁡p×K is found. Hence, the optimisation is currently carried out by varying each parameter separately (coordinate descent algorithm). Thereby, ln⁡p×K increases when Δf becomes smaller. This implies an increase of the dimensionality of the involved covariances ΣPr and ΣPo and, thus, a prolonging of the optimisation procedure. To demonstrate this procedure, the parameters fL = 0 GHz, aE=aO=10-8 V2, x0=-5.118 mm, δx=0 mm, σb0,Pr=1 V, σb1,Pr=1 V/m, and σb2,Pr=1 V/m2 are set to the values used for the nonoptimised case (see Section 5.3.1). Scanned roughly in fU, ln⁡p×K shows a peak close to 1000 GHz (see Figure 6(a)). Furthermore, this peak increases for Δf≤5 GHz to ln⁡p×K≈106. The peak is localised at 860 GHz which moves to about 910 GHz for the reduced values aE=aO=10-11 V2 (see Figure 6(b)).Figure 6
Ten-dimensional quantityp×K scanned at different values for fU, Δf, aE, and aO while keeping fL = 0 GHz, x0=-5.118 mm, δx=0 mm, σb0,Pr = 1 V, σb1,Pr = 1 V/m, and σb2,Pr = 1 V/m2. (a) Rough scan in fU and Δf at aE=aO=10-8 V2. A peak for ln⁡p×k with the amplitude of the order of 106 is found close to fU = 1000 GHz. Small Δf are in favour. (b) Scan in fU and Δf close to peaking for two sets of aE and aO. The peak locates in the vicinity of 910 GHz when aE and aO reduce by some orders of magnitude.
(a)
(b)To ease the computational effort but still being able to characteriseln⁡p×K, its maximum is determined dependent on Δf which is scanned in the values 5.68, 4, 3, 2, 1, and 1/2 GHz. Each maximum is captured by the sets fL,Po, fU,Po, aE,Po, aO,Po, x0,Po, δxPo, σb0,Pr,Po, σb1,Pr,Po, and σb2,Pr,Po summarised by Table 1. All maxima locate at very similar sets which seem to converge as Δf becomes smaller. Although, in relative terms, the maximum for a smaller Δf has higher odds (see Figure 7). For example, the odds read 1 : 0.011 when the maxima at Δf = 1/2 GHz and 5.68 GHz are related. Furthermore, taking the values for the maximum at Δf = 1/2 GHz to evaluate ln⁡p×K at Δf = 1/3 GHz and 1/4 GHz gives the odds 1 : 1.013 : 1.018. Thus, the global maximum locates somewhere in the range below Δf≤ 1/4 GHz. But this range cannot be investigated in more detail from numerical point of view. However, the increase in ln⁡p×K when Δf is decreased below 1/2 GHz is interpreted as confirmation that a continuous spectrum is indeed probed.Table 1
Characterisation of maxima ofp×K dependent on Δf. The maxima locate at similar nine-dimensional parameter sets. The odds are evaluated by the ratio p×K(Δf)/p×K(Δf = 1/2 GHz) and relate the maxima globally. Since the odds are rising, the global maximum locates below Δf = 1/2 GHz.
Δ
f (GHz)
≈5.68
4
3
2
1
1/2
f
L
,
P
o (GHz)
33.499
33.243
33.277
32.817
32.445
32.425
f
U
,
P
o (GHz)
913.572
913.243
912.277
912.817
913.445
913.425
a
E
,
P
o (×10-12 V2)
0.9152
1.0127
1.0623
1.1002
1.1224
1.1299
a
O
,
P
o (×10-15 V2)
5.7381
5.7731
5.7776
5.797
5.8126
5.8145
x
0
,
P
o (mm)
−5.115293
−5.115297
−5.115288
−5.115291
−5.115295
−5.115295
δ
x
P
o (μm)
−4.814
−4.814
−4.814
−4.814
−4.814
−4.814
σ
b
0
,
P
r
,
P
o (nV)
62.150
62.148
62.146
62.143
62.142
62.142
σ
b
1
,
P
r
,
P
o (μV/m)
6.031
6.003
6.003
6.000
5.998
5.998
σ
b
2
,
P
r
,
P
o (μV/m2)
11.541
9.884
9.913
9.725
9.575
9.573
ln
⁡
p
×
K
1052035.05
1052037.54
1052038.52
1052039.15
1052039.48
1052039.55
Odds
0.011
0.134
0.357
0.668
0.929
1Figure 7
Odds for maxima ofp×K dependent on spectral increment Δf. The maximum at Δf = 1/2 GHz is used as a reference (unity odd). For Δf = 1/3 GHz and 1/4 GHz, the maximising parameter set (see Table 1) for Δf = 1/2 GHz is used to evaluate the associated odds. The global maximum of p×K is obtained for very small spectral increments.ForΔf = 5.68 GHz, the odds read 1 : 10−1390 when the associated maximum (ln⁡p×K≈1.0520×106) is related to the nonoptimised case (ln⁡p×K≈1.0488×106) investigated in Section 5.3.1. This is caused by having chosen some parameters like aE(O) several orders of magnitude too large with respect to the maximum settings values.ForΔf = 1 GHz, Figures 8(a)–8(f) summarise scans in the parameter pairs (aE,aO), (fL,fU), (x0,δx), (σb0,Pr,σb1,Pr), (σb0,Pr,σb2,Pr) and in σb0,Pr while the remaining parameters are held at the maximum values. In σb0,Pr, σb1,Pr, and σb2,Pr, a skewed distribution is found. The remaining six-dimensional posterior distribution has a high probability in a narrow region for a given Δf. This distribution is well approximated by a multivariate normal. Its mean is given by the maximum values listed in Table 1. The posterior covariance is estimated from the inverse of the Hessian matrix evaluated numerically via the second-order partial derivatives in the vicinity of the maximum. The off-diagonal elements of this covariance are negligible, so that one can factorise the posterior as a product of individual Gaussians. The posterior standard deviations σfL,Po, σfU,Po, σaE,Po, σaO,Po, σx0,Po, and σδx,Po vary little when Δf is changed (see Table 2). The spectral boundaries are well determined within an interval of some GHz. While the uncertainty in the scaling aE is small compared to its posterior mean value, aO is more uncertain. The uncertainty in the zero-path difference is of the order of some μm, and the shift δx is quite certain within some hundreds of nanometers.Table 2
Posterior standard deviations dependent onΔf for spectral boundaries fL and fU, scalings aE and aO of prior covariances, zero-path difference x0, and shift δx. The standard deviations are estimated from the inverse Hessian matrix at the corresponding maximum of p×K.
Δ
f (GHz)
≈5.68
4
3
2
1
1/2
σ
f
L
,
P
o (GHz)
1.742
1.795
1.746
1.762
1.754
1.718
σ
f
U
,
P
o (GHz)
3.101
3.155
2.927
2.869
2.858
2.810
σ
a
E
,
P
o (×10-12 V2)
0.1064
0.1179
0.1242
0.1290
0.1316
0.1327
σ
a
O
,
P
o (×10-15 V2)
1.5912
1.6064
1.6112
1.6164
1.6199
1.6211
σ
x
0
,
P
o (μm)
1.207
1.207
1.207
1.206
1.205
1.205
σ
δ
x
,
P
o (μm)
0.170
0.170
0.170
0.170
0.170
0.170Figure 8
p
×
K scanned in individual nonlinear parameters and hyperparameters in vicinity of maximum of p×K for Δf = 1 GHz (see Table 1). (a) Zero-path difference position x0 and shift δx. (b) Prior standard deviation σb0,Pr for absolute coefficient. (c) Lower and upper limits fL and fU of spectral domain. (d) σb0,Pr and prior standard deviation σb1,Pr for linear coefficient. (e) Scaling aE and scaling aO of Brownian bridge priors for S→E and S→O. (f) σb0,Pr and prior standard deviation σb2,Pr for quadratic coefficient.


### 5.3.3. Conditional Amplitude Posterior for Maximising Settings

In the following, the amplitude posterior is investigated for the maximising settings givenΔf = 1/2 GHz (see Table 1). For the listed settings, the double-sided domain is identified by |x|≤|x0,Po|≈ 5.115 mm; the single-sided domain is bounded by the lower limit x>|x0,Po| and the upper limit of about 26.36 mm, the centre frequency fC,Po, and the bandwidth ΔF,Po of the spectral domain read 472.925 GHz and 881 GHz, respectively. In addition, the conversion of the spatial coordinate to the order of the Fourier coefficients k=|x|ΔF,Po/c reveals that the double- and single-sided domain contain information up to the 15th and 77th order (see Figures 9(c) and 9(d)), respectively. From this follows that the interferometric part of the 256/788 data values in the double-sided/full domain are modelled best by 62/310 Fourier coefficients and their associated basis functions dependent on fC,Po, ΔF,Po, and x. The signal envelopes σD,E and σD,O are evaluated with aE,Po and aO,Po.Figure 9
(a) Posterior meansS→E,Po and S→O,Po with associated 2-σ uncertainty bands evaluated at maximum of settings posterior for Δf = 1/2 GHz (see Table 1). Most likely, the spectral quantities S→E,Po and S→O,Po are finite only for the domain between fL,Po = 32.425 GHz and fU,Po = 913.425 GHz. This upper limit is approximately 4 smaller than the Nyquist frequency fNy = 3747.4 GHz. The dips in S→E,Po at 557 and 752 GHz are caused by absorption due to atmospheric water vapour. (b) D→E,Po, D→O,Po, and D→Off,Po obtained by mapping posterior means S→E,Po, S→O,Po, b0,Po=-62.142 nV, b1,Po = 5.997 μV/m, and b2,Po=-8.641 μV/m2 to data domain. (c) Cosine transform D→E,Po of S→E,Po dependent on optical path difference x and order k=|x|ΔF,Po/c (bandwidth ΔF,Po=fU,Po-fL,Po) of Fourier coefficients. The envelope σD,E∝aE,Po1/2/k, incorporating the scaling aE,Po = 1.1299 × 10−12 V2, bounds D→E,Po best. (d) Zoom of (b) dependent on x and k. Globally, the sine transform D→O,Po of S→O,Po is best summarised by the envelope σD,O∝aO,Po1/2/k with the scaling aO,Po=5.8145×10-15 V2. Since σD,E≫σD,O holds, the single-sided domain is determined mainly by S→E,Po. (e) Histogram of residuals, that is, difference between measured data D→ and posterior mean prediction D→Po=D→E,Po+D→O,Po+D→Off,Po scaled by noise level σD. The data is well described, because the mean -1.7×10-4 and the variance 0.942 of the residuals are almost vanishing and close to unity, respectively. (f) Residuals versus optical path difference.
(a)
(b)
(c)
(d)
(e)
(f)(a) Posterior Mean. Given the boundaries and increment, the spectral dimension reads Nf = 1762, and, thus, the amplitude posterior mean l→Po=S→E,Po,S→O,Po,b→PoT is 3527-dimensional. The posterior mean values for the three coefficients b0,Po=-62.142 nV, b1,Po = 5.997 μV/m, and b2,Po=-8.641 μV/m2 have magnitudes which are close to the corresponding prior standard deviations at the maximum of the settings posterior (see Table 1). Furthermore, the absolute and the linear mean values are similar to the ones obtained for the nonoptimised case (see Section 5.3.1). The offset D→Off,Po which follows from mapping b→Po is shown in Figures 9(b) and 9(d).The posterior meansS→E,Po and S→O,Po presented in Figure 9(a) approach zero towards the spectral boundaries as expected. Furthermore, S→E,Po is much larger in amplitude than S→O,Po anticipated by the quite different posterior mean scalings aE,Po and aO,Po. Hence, the even process describes most of the interferometric data in the double-sided domain. This is seen by mapping both spectral means to the data domain which gives the even and odd quantities D→E,Po and D→O,Po (see Figures 9(b)–9(d)). In addition, the single-sided domain is described mainly by S→E,Po, because in this domain only the envelope σD,E for the even process is of the order of the interferometric data (see Figure 9(d)). This differs from the findings for the nonoptimised case (aE=aO) for which S→E,Po and S→O,Po determine the single-sided region almost equally.The histogram of the residuals (see Figure9(e)) is approximated very well by the normal distribution N(-1.7×10-4,0.942). Since the mean vanishes almost, and the standard deviation is very close to unity, the data set is well described by the model and the posterior means. Most likely, the data point located at x≈ 7 mm (see Figure 9(f)) is an outlier, because its residual is outside the 3.5-σD band.(b) Posterior Covariance. For the linear parameters the posterior covariance matrix is written like(61)ΣPo=ΣS→E,E,PoΣS→E,O,PoΣS→E,b→,PoΣS→E,O,PoTΣS→O,O,PoΣS→O,b→,PoΣS→E,b→,PoTΣS→O,b→,PoTΣb→,b→,Po,and taking the square root of an element of the main diagonal gives the posterior standard deviation σPi for the ith parameter. The correlation coefficient between two parameters P1 and P2 can be evaluated by ρP1,P2=ΣP1,P2,Po/(σP1σP2).The standard deviationsσb0,Po = 0.501 nV, σb1,Po = 0.093 μV/m, and σb2,Po = 4.119 μV/m2 of the three coefficients are reduced by about two orders of magnitude with respect to the uncertainties obtained for the nonoptimised case. Though, the quadratic coefficient is inferred with small confidence. Coefficients of neighbouring degree have the highest correlation in magnitude and are anticorrelated. For example, one finds ρb0,b1=-0.55 and ρb1,b2=-0.93, while ρb0,b2 = 0.29 remains small.The covariancesΣS→E,E,Po(f,f′) and ΣS→O,O,Po(f,f′) and their associated correlations ρE,E(f,f′) and ρO,O(f,f′) resemble sinc functions centred along the line f=f′ plus a tip at this condition (see Figures 10(a)–10(c) for cross section at f′ = 100.175 GHz). Though, the FWHM widths of 1.5 GHz and 10 GHz are very different for the even and odd covariances/correlations. Regarding the covariances, the tip heights account for 1.3×10-38 V2/Hz2 and 4.4×10-40 V2/Hz2, and towards the spectral boundaries the heights are approaching zero which can be seen by the standard deviations σS→E,Po=ΣS→E,E,Po1/2(f,f) and σS→O,Po=ΣS→O,O,Po1/2(f,f) (see Figure 10(d)).Figure 10
Characterisation of posterior covariance for amplitude parameters obtained at maximum of settings posterior forΔf = 1/2 GHz (see Table 1). (a) Cross section of covariances ΣS→E,E,Po and ΣS→O,O,Po for even and odd function. Each covariance is similar to a sinc function centred at f=f′, where a tip is added on top. (b) Cross section of covariances ΣS→O,O,Po and ΣS→E,O,Po. At a fixed frequency, S→E and S→O are independent, because ΣS→E,O,Po = 0. But finite correlation remains with the neighbouring spectral domain. (c) Correlation. The even (odd) parameters are correlated (ρE(O),E(O)>0.5) over the spectral scale of about 1.5 GHz (10 GHz), while the correlation between even and odd parameters is small and asymmetric with respect to f=f′. (d) Standard deviations σS→E,Po and σS→O,Po. At the spectral limits the uncertainties approach zero. (e) and (f) Correlations of S→E(O) and b→. The correlations are oscillating and decrease generally with the frequency.
(a)
(b)
(c)
(d)
(e)
(f)At the conditionf=f′, the covariance ΣS→E,O,Po and the correlation ρE,O vanish and are asymmetrically close to this line (see Figures 10(b) and 10(c) for cross section at f′ = 100.175 GHz). This means that at a given frequency f the amplitudes S→E and S→O are independent but have a finite and opposite correlation with the neighbouring spectral domains to either side.The correlationsρE,b→ and ρO,b→ are oscillatory and rise towards the spectral origin (see Figures 10(e) and 10(f)).(c) Posterior Samples. From the conditional amplitude posterior samples S→E,Po,Sa, S→O,Po,Sa and b→Po,Sa are drawn. While Figure 11(a) shows one sample for each of the even and odd spectral functions, Figure 11(b) presents 100 samples. The samples form a band around the corresponding posterior mean S→E(O),Po with the width of about twice the standard deviation σS→E(O),Po as shown in Figure 10(d). Hence, the band for the odd function is much smaller.Figure 11
(a) Individual posterior sample functionsS→E,Po,Sa and S→O,Po,Sa. (b) 100 posterior sample functions S→E,Po,Sa and S→O,Po,Sa. The uncertainty at a given frequency is larger for the even function, because for the associated scalings of the two processes aE,Po≫aO,Po holds. (c) 100 spectral samples and 100 offset coefficients mapped to probed data domain. Unlike for the nonoptimised case (see Figure 5(b)), the spread is much smaller, and the quantities connect smoothly at the transition between double- and single-sided regions. (d) 100 samples D→Po,Sa and measured data D→. The samples are in agreement with the uncertainty of the measurements. (e) One sample D→E,Po,Sa and D→O,Po,Sa evaluated inside probed domain (x < 26.36 mm) and outside. With respect to the envelopes σD,E(O), the contributions behave similarly in both domains. (f) 100 samples D→E,Po,Sa and D→O,Po,Sa evaluated inside the probed domain and outside. Outside the probed domain, the area between the envelops σD,E(O) is filled by D→E,Po,Sa and D→O,Po,Sa. Because the spread in the nonprobed domain is wider than the one in the probed domain, the main uncertainty in the even and odd spectra originates in nonprobed Fourier coefficients.
(a)
(b)
(c)
(d)
(e)
(f)The samples mapped to the data domain giveD→E,Po,Sa, D→O,Po,Sa, and D→Off,Po,Sa (see Figure 11(c)) which form much more narrow bands when compared to the nonoptimised case (see Figure 5(b)). In particular, the transition from the double- to the single-sided domain is smooth. This is a consequence of having used the most likely scalings aE,Po and aO,Po and the associated but quite different signal envelops inferred from the data. Furthermore, this explains why the band for D→O,Po,Sa is much smaller than the one for D→E,Po,Sa.For each sample, the sumD→Po,Sa=D→E,Po,Sa+D→O,Po,Sa+D→Off,Po,Sa complies with the measured data D→ and the noise level (see Figure 11(d)).With the samplesS→E,Po,Sa and S→O,Po,Sa at hand, D→E,Po,Sa and D→O,Po,Sa are predicted outside the probed spatial domain. Up to the optical path difference of 70 mm, Figures 11(e) and 11(f) show D→E,Po,Sa and D→O,Po,Sa for one sample and 100 samples, respectively. No difference is obvious in D→E,Po,Sa and D→O,Po,Sa for an individual prediction in the probed and nonprobed domain. However, many samples take values between σD,E(O) outside the probed region. Since the spread in the probed domain is smaller than in the nonprobed domain, the main uncertainty in S→E,Po and S→O,Po originates in Fourier coefficients not probed rather than the noise on measured data.


### 5.3.4. Marginal Amplitude Posterior

Since only the marginal posteriorpS→E,S→O∣Δf,D→ of the spectral functions is of interest for this specific Fourier spectroscopy problem, formally the marginalisation(62)pS→E,S→O∣Δf,D→=∫pS→E,S→O,b→∣n→,h→,Δf,D→pn→,h→∣Δf,D→db→dn→dh→needs to be carried out. In general, the pdf pS→E,S→O∣Δf,D→ is not Gaussian but can be approximated by a multivariate normal. Performing the marginalisation, that is, the uncertainty propagation, is analytically not possible but achievable numerically described below for Δf = 1/2 GHz.To reduce the numerical efforts, only the most important parametersfL, fU, x0, δx, aE, and aO are taken into consideration for the marginalisation, and the remaining three parameters σb0,Pr, σb1,Pr, and σb2,Pr are held at the maximum values listed in Table 1. The remaining pdf, labeled by p′n→,h→∣Δf,D→ in the following, is well approximated by a factorisable multivariate normal, because the corresponding covariances are negligible (see Section 5.3.2). The posterior mean values and standard deviations are listed in Tables 1 and 2, respectively.From the six-dimensional Gaussian posteriorp′n→,h→∣Δf,D→, 1000 samples are drawn. For each sample, mean and covariance for the conditional amplitude posterior are evaluated. Thereby, a difficulty arises due to the lower and upper spectral limits fL and fU. For given increment Δf, fL and fU determine the spectral grid. To keep this grid fixed, the spectral grid obtained for fL,Po = 32.425 GHz and fU,Po = 913.425 GHz is extended by some tens of GHz on either side. Furthermore, the posterior samples for fL and fU are compared with the extended but fixed spectral grid points and changed to the closest grid point. With the modified boundaries and the remaining samples unchanged, a conditional amplitude posterior follows. For example, Figure 12(a) shows the mean functions S→E,Po and S→O,Po for the 1000 posterior samples of the settings. While S→E,Po changes little but most at the spectral boundaries, S→O,Po varies considerably, mostly in the centre of the spectral domain covered. This is caused mainly by the posterior uncertainty of the zero-path difference x0. However, the spread is symmetric around the most likely mean.Figure 12
(a) MeansS→E,Po (blue) and S→O,Po (red) obtained for 1000 samples of settings posterior at Δf = 1/2 GHz. Only at the spectral limits the means deviate from the most likely mean S→E,Po (cyan) evaluated at the maximum of the settings posterior. The uncertainty in the zero-path difference x0 causes a considerable deviation from the most likely mean S→O,Po (black). (b) Means S→E,Po,Ma (cyan) and S→O,Po,Ma (black) after performing numerical marginalisation of settings parameters at Δf = 1/2 GHz. The marginalised means are similar to the means S→E,Po and S→O,Po (see (a)) at the maximum of the settings posterior except for minor differences towards the spectral boundaries. Having the marginalised covariance at hand, 100 samples S→E,Po,Ma,Sa (blue) and S→O,Po,Ma,Sa (red) are drawn from the marginalised amplitude posterior. Only in the centre of the spectral domain the spread of S→O,Po,Ma,Sa increases to capture the uncertainty in x0 as well (see (a)).
(a)
(b)For each settings sample, 100 samples are drawn frompS→E,S→O,b→∣n→,h→,Δf,D→. Overall, 100000 samples are available for b→, S→E, and S→O at each spectral grid point. From this large sample set, the marginalised means S→E,Po,Ma, S→O,Po,Ma and covariance ΣS→,Po,Ma are evaluated which capture the first and second moments of the pdf pS→E,S→O∣Δf,D→. In doing so, the marginalisation with respect to b→ is achieved implicitly.(a) Posterior Mean. The means S→E,Po,Ma and S→O,Po,Ma (see Figure 12(b)) look similar to S→E,Po, S→O,Po at the maximum of the settings posterior. However, S→E,Po,Ma and S→O,Po,Ma approach zero some GHz below and above fL,Po and fU,Po, respectively.(b) Posterior Covariance. The posterior covariance after marginalisation with respect to the setting parameters is decomposed like(63)ΣS→,Po,Ma=ΣS→E,E,Po,MaΣS→E,O,Po,MaΣS→E,O,Po,MaTΣS→O,O,Po,Ma.ΣS→E,E,Po,Ma (see Figure 13(a) for cross section at f′ = 100.175 GHz) and is very similar to ΣS→E,E,Po evaluated at the maximum of the settings posterior (see Section 5.3.3). The same holds for ΣS→E,O,Po,Ma and ΣS→E,O,Po. However, a significant change is found for ΣS→O,O,Po,Ma (see Figures 13(a) and 13(b) for some cross sections) when compared to ΣS→O,O,Po. Due to the posterior uncertainty in x0, a broadband additive contribution increases ΣS→O,O,Po,Ma up to about one order of magnitude in the centre of the spectral domain. Comparing the standard deviations σS→E,Po,Ma with σS→E,Po, and σS→O,Po,Ma with σS→O,Po (see Figure 13(c)) shows that towards the spectral limits the propagation of the settings posterior uncertainties has little effect. Furthermore, the correlation ρO,O,Ma is high for a broad spectral range in the central region (see Figure 13(d)).Figure 13
Characterisation of posterior covarianceΣS→,Po,Ma numerically marginalised with respect to settings parameters. ΣS→,Po,Ma is decomposed, so that ΣS→E(O),E(O),Po,Ma captures the covariance of the even (odd) spectral quantities. (a) Cross section at frequency f′ = 100.175 GHz. Apart from the extension in the spectral domain, ΣS→E,E,Po,Ma is similar to the covariance ΣS→E,E,Po evaluated at the maximum of the settings posterior (see Figure 10(a)). ΣS→O,O,Po,Ma shows a significant deviation from ΣS→O,O,Po (see also (b)). (b) Cross sections of ΣS→O,O,Po,Ma at frequencies f′ = 100.175 GHz, 400.175 GHz, and 850.175 GHz. Towards the spectral limits ΣS→O,O,Po,Ma behaves like ΣS→O,O,Po. But in the central region a broadband additive contribution is evident originating in the posterior uncertainty of the zero-path difference x0. (c) Standard deviations evaluated from main diagonal of marginalised covariance ΣS→,Po,Ma (red) and covariance ΣS→,Po (black, see Figure 10(d)) evaluated at maximum of settings posterior. The odd quantity increases by a factor of three in the central region after the marginalisation. (d) Correlation ρO,O,Ma at three frequencies. Compared to ρO,O (see Figure 10(c)), a broadband (several 100s GHz) elevation of ρO,O,Ma up to 0.8 follows from the uncertainty in x0 for the central domain.
(a)
(b)
(c)
(d)(c) Posterior Samples. From the marginalised Gaussian posterior pS→E,S→O∣Δf,D→ 100 samples (S→E,Po,Ma,Sa and S→O,Po,Ma,Sa) are drawn and shown in Figure 12(b). As expected, the spread of S→E,Po,Ma,Sa around the mean S→E,Po,Ma is almost equal to the one for S→E,Po,Sa around the mean S→E,Po (see Figure 11(b)). The samples S→O,Po,Ma,Sa have a wider spread in the centre of the spectral domain than S→O,Po,Sa mainly caused by the uncertainty in the zero-path difference.


## 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 forx0, δx, fL, fU, Δf, aE and aO, σb0,Pr, σb1,Pr, and σb2,Pr.The values forx0=-127.95Δx=-5.118 mm and δx=0 mm are chosen to form an optical path difference grid xi for which all subsequent nodes are separated by Δx. Then, |xi|≤5.118 mm identifies the double-sided domain (256 data points), and the single-sided domain is marked by 5.118 mm ≤xi≤ 26.362 mm (532 data points).The spectral domain is covered fromfL = 0 GHz to the Nyquist frequency fU=c/(2Δx) = 3747.4 GHz (=ΔF). With the bandwidth and the optical path difference set, the conversion k=|x|ΔF/c (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Δf = 5.68 GHz is used, and the lower and upper spectral limits are set to fL = 0 GHz and fU = 3747.4 GHz (Nyquist assumptions), respectively. The scalings of the processes are chosen by aE=aO=10-8 V2. (a) Prior samples S→E,Pr,Sa and S→O,Pr,Sa 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 D→E(O),Pr,Sa for optical path difference obtained for x0=-5.118 mm and δx = 0. With increasing distance to the origin, the amplitude drops which is well bounded by the signal envelopes σD,E=σD,O for the chosen bandwidth ΔF=fU. Fourier coefficients up to order k = 64 and 330 settle in the double-sided domain x≤x0 and single-sided domain x>x0, respectively.
(a)
(b)To match the maximum order (330) of probed Fourier coefficient, the discretisation is set toΔf=c/(2Δx(256/2+532))≈5.68 GHz which reflects the classical increment. Then, the spectral grid has Nf = 660 elements.The hyperparameters are set to large values likeaE=aO=10-8 V2 to prevent from a prior determined posterior. This unwanted feature would be present, if one would choose aE(O) too small, so that the signal envelops do not include the measured data (see next paragraph).With the chosen spectral priors, sample functions/vectorsS→E,Pr,Sa and S→O,Pr,Sa are drawn and shown in Figure 4(a). The function values are of the order of 10-16 V/Hz decaying towards the spectral boundaries. The mapping of the spectral prior samples to the spatial domain via D→E,Pr,Sa+D→O,Pr,Sa=MES→E,Pr,Sa+MOS→O,Pr,Sa shows that D→E,Pr,Sa+D→O,Pr,Sa 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 D→E,Pr,Sa+D→O,Pr,Sa drops with the distance from the origin. All these data samples are well bounded by the signal envelopes σD,E(O) obtained by σD,E(O)=(A/2)σD,E(O)′ which is the scaled version of the expressions (44) and (45).Forb→ the prior covariance is determined by setting σb0,Pr=1 V, σb1,Pr=1 V/m and σb2,Pr=1 V/m2.(a) Posterior Mean. For the above values, the mean l→Po of the amplitude posterior is given by b→Po=(-61.41 nV, 5.75 μV/m, 2.42 μV/m2)T, S→E,Po and S→O,Po as presented in Figure 5(a). Clearly, S→E,Po peaks below 1000 GHz with the maximum of approximately 6 × 10−18 V/Hz around 400 GHz. Furthermore, S→O,Po is much smaller than S→E,Po below 1000 GHz.Figure 5
Investigation of conditional amplitude posteriorpl→∣Δf,fL,fU,aE,aO,x0,δx,σb0,Pr,σb1,Pr,σb2,Pr,D→ for linear parameters l→=S→E,S→O,b→T. The posterior is obtained for the settings x0=-5.118 mm and δx = 0, fL = 0 GHz, fU = 3747.4 GHz, Δf = 5.68 GHz, aE=aO=10-8 V2, σb0,Pr=1 V, σb1,Pr=1 V/m, and σb2,Pr=1 V/m2. (a) Spectral posterior means S→E(O),Po and one posterior sample. Below 1000 GHz, the means show a finite amplitude. Each sample deviates strongly from its corresponding mean. (b) Posterior mean l→Po and 100 samples l→Po,Sa mapped to give even, odd, and offset contributions in data domain. While in the double-sided region (|x|≤|x0|) most of the data described by S→E,Po, S→E,Po and S→O,Po give same contributions in the single-sided region (x>|x0|). In the single-sided region, the spread of the even (odd) contributions for the posterior samples is large but bounded by the signal envelops σD,E(O). (c) Residuals for l→Po (black) and one sample l→Po,Sa (red) both mapped to data domain. An almost perfect description of the noisy data is achieved by lPo, indicating an overfitting of the data. (d) Measured data D→ and 100 samples of amplitude posterior mapped to data domain. The individual contributions have a much wider spread (see (b)), but the total sum D→Po,Sa agrees well with D→ and its uncertainties. This is caused by the posterior covariance of l→ with suited off-diagonal entries.
(a)
(b)
(c)
(d)The offsetD→Off,Po=MOffb→Po contribution to the data has an upwards trend but is small (see Figure 5(b)). The even and odd contributions D→E(O),Po=ME(O)S→E(O),Po behave differently in the double- and single-sided domains (see Figure 5(b)). While in the double-sided domain D→E,Po is much larger than D→O,Po, D→E,Po and D→O,Po equal each other in the single-sided domain. This indicates an underestimation (overestimation) of the Fourier coefficients above the order 64 for SE,Po (SO,Po). The cause of this finding relies in the choice aE=aO (σD,E=σD,O) which becomes more clear in Section 5.3.3.Since the quantityD→Po=Ml→Po equals almost the data set D→ (see residuals in Figure 5(c)), a nearly perfect match is achieved. Hence, all noise contributions are captured by the mean l→Po, and an overfitting of the data is obvious.(b) Posterior Covariance. The posterior covariance ΣPo is not characterised in detail. From the elements of the main diagonal of ΣPo, one finds the values σb0,Po = 73.30 nV, σb1,Po = 12.58 μV/m, and σb2,Po = 446.67 μV/m2. These posterior uncertainties are quite large when compared to the values of the corresponding mean b→Po, especially for the quadratic coefficient.For the spectral quantities, the square root of the main diagonal elements ofΣPo is of the order of some 10-18 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 l→Po,Sa=S→E,Po,Sa,S→O,Po,Sa,b→Po,SaT are drawn from the amplitude posterior and mapped to the data domain (see Figure 5(b)). Thereby, the contributions D→E,Po,Sa=MES→E,Po,Sa, D→O,Po,Sa=MOS→O,Po,Sa, and D→Off,Po,Sa=MOffb→Po,Sa are split to investigate their interplay. D→Off,Po,Sa has a noticeable width but remains in the vicinity of D→Off,Po for the whole spatial domain.For the double-sided regionD→E,Po,Sa and D→O,Po,Sa form a narrow band around the corresponding posterior mean quantities D→E,Po and D→O,Po. On the contrary, for the single-sided region D→E,Po,Sa and D→O,Po,Sa increase drastically but are restricted by σD,E and σD,O. The deviation is large between the individual and the posterior mean quantities; though, the sum D→Po,Sa=D→E,Po,Sa+D→O,Po,Sa+D→Off,Po,Sa remains well inside the 2-σD band surrounding D→ (see Figures 5(c) and 5(d)). Hence, the posterior covariance ΣPo captures properly the correlations.With the specific values of the settings one gets the numberln⁡p×K≈1.0488×106.


## 5.3.2. Settings Posterior

As the problem is formulated, the settings posterior is proportional top×K. Its optimisation, that is, finding the global maximum of the settings posterior is a 10-dimensional problem. Since large numbers are involved, one has to investigate ln⁡p×K. In general, if only one parameter is changed, a well distinguishable peak in ln⁡p×K is found. Hence, the optimisation is currently carried out by varying each parameter separately (coordinate descent algorithm). Thereby, ln⁡p×K increases when Δf becomes smaller. This implies an increase of the dimensionality of the involved covariances ΣPr and ΣPo and, thus, a prolonging of the optimisation procedure. To demonstrate this procedure, the parameters fL = 0 GHz, aE=aO=10-8 V2, x0=-5.118 mm, δx=0 mm, σb0,Pr=1 V, σb1,Pr=1 V/m, and σb2,Pr=1 V/m2 are set to the values used for the nonoptimised case (see Section 5.3.1). Scanned roughly in fU, ln⁡p×K shows a peak close to 1000 GHz (see Figure 6(a)). Furthermore, this peak increases for Δf≤5 GHz to ln⁡p×K≈106. The peak is localised at 860 GHz which moves to about 910 GHz for the reduced values aE=aO=10-11 V2 (see Figure 6(b)).Figure 6
Ten-dimensional quantityp×K scanned at different values for fU, Δf, aE, and aO while keeping fL = 0 GHz, x0=-5.118 mm, δx=0 mm, σb0,Pr = 1 V, σb1,Pr = 1 V/m, and σb2,Pr = 1 V/m2. (a) Rough scan in fU and Δf at aE=aO=10-8 V2. A peak for ln⁡p×k with the amplitude of the order of 106 is found close to fU = 1000 GHz. Small Δf are in favour. (b) Scan in fU and Δf close to peaking for two sets of aE and aO. The peak locates in the vicinity of 910 GHz when aE and aO reduce by some orders of magnitude.
(a)
(b)To ease the computational effort but still being able to characteriseln⁡p×K, its maximum is determined dependent on Δf which is scanned in the values 5.68, 4, 3, 2, 1, and 1/2 GHz. Each maximum is captured by the sets fL,Po, fU,Po, aE,Po, aO,Po, x0,Po, δxPo, σb0,Pr,Po, σb1,Pr,Po, and σb2,Pr,Po summarised by Table 1. All maxima locate at very similar sets which seem to converge as Δf becomes smaller. Although, in relative terms, the maximum for a smaller Δf has higher odds (see Figure 7). For example, the odds read 1 : 0.011 when the maxima at Δf = 1/2 GHz and 5.68 GHz are related. Furthermore, taking the values for the maximum at Δf = 1/2 GHz to evaluate ln⁡p×K at Δf = 1/3 GHz and 1/4 GHz gives the odds 1 : 1.013 : 1.018. Thus, the global maximum locates somewhere in the range below Δf≤ 1/4 GHz. But this range cannot be investigated in more detail from numerical point of view. However, the increase in ln⁡p×K when Δf is decreased below 1/2 GHz is interpreted as confirmation that a continuous spectrum is indeed probed.Table 1
Characterisation of maxima ofp×K dependent on Δf. The maxima locate at similar nine-dimensional parameter sets. The odds are evaluated by the ratio p×K(Δf)/p×K(Δf = 1/2 GHz) and relate the maxima globally. Since the odds are rising, the global maximum locates below Δf = 1/2 GHz.
Δ
f (GHz)
≈5.68
4
3
2
1
1/2
f
L
,
P
o (GHz)
33.499
33.243
33.277
32.817
32.445
32.425
f
U
,
P
o (GHz)
913.572
913.243
912.277
912.817
913.445
913.425
a
E
,
P
o (×10-12 V2)
0.9152
1.0127
1.0623
1.1002
1.1224
1.1299
a
O
,
P
o (×10-15 V2)
5.7381
5.7731
5.7776
5.797
5.8126
5.8145
x
0
,
P
o (mm)
−5.115293
−5.115297
−5.115288
−5.115291
−5.115295
−5.115295
δ
x
P
o (μm)
−4.814
−4.814
−4.814
−4.814
−4.814
−4.814
σ
b
0
,
P
r
,
P
o (nV)
62.150
62.148
62.146
62.143
62.142
62.142
σ
b
1
,
P
r
,
P
o (μV/m)
6.031
6.003
6.003
6.000
5.998
5.998
σ
b
2
,
P
r
,
P
o (μV/m2)
11.541
9.884
9.913
9.725
9.575
9.573
ln
⁡
p
×
K
1052035.05
1052037.54
1052038.52
1052039.15
1052039.48
1052039.55
Odds
0.011
0.134
0.357
0.668
0.929
1Figure 7
Odds for maxima ofp×K dependent on spectral increment Δf. The maximum at Δf = 1/2 GHz is used as a reference (unity odd). For Δf = 1/3 GHz and 1/4 GHz, the maximising parameter set (see Table 1) for Δf = 1/2 GHz is used to evaluate the associated odds. The global maximum of p×K is obtained for very small spectral increments.ForΔf = 5.68 GHz, the odds read 1 : 10−1390 when the associated maximum (ln⁡p×K≈1.0520×106) is related to the nonoptimised case (ln⁡p×K≈1.0488×106) investigated in Section 5.3.1. This is caused by having chosen some parameters like aE(O) several orders of magnitude too large with respect to the maximum settings values.ForΔf = 1 GHz, Figures 8(a)–8(f) summarise scans in the parameter pairs (aE,aO), (fL,fU), (x0,δx), (σb0,Pr,σb1,Pr), (σb0,Pr,σb2,Pr) and in σb0,Pr while the remaining parameters are held at the maximum values. In σb0,Pr, σb1,Pr, and σb2,Pr, a skewed distribution is found. The remaining six-dimensional posterior distribution has a high probability in a narrow region for a given Δf. This distribution is well approximated by a multivariate normal. Its mean is given by the maximum values listed in Table 1. The posterior covariance is estimated from the inverse of the Hessian matrix evaluated numerically via the second-order partial derivatives in the vicinity of the maximum. The off-diagonal elements of this covariance are negligible, so that one can factorise the posterior as a product of individual Gaussians. The posterior standard deviations σfL,Po, σfU,Po, σaE,Po, σaO,Po, σx0,Po, and σδx,Po vary little when Δf is changed (see Table 2). The spectral boundaries are well determined within an interval of some GHz. While the uncertainty in the scaling aE is small compared to its posterior mean value, aO is more uncertain. The uncertainty in the zero-path difference is of the order of some μm, and the shift δx is quite certain within some hundreds of nanometers.Table 2
Posterior standard deviations dependent onΔf for spectral boundaries fL and fU, scalings aE and aO of prior covariances, zero-path difference x0, and shift δx. The standard deviations are estimated from the inverse Hessian matrix at the corresponding maximum of p×K.
Δ
f (GHz)
≈5.68
4
3
2
1
1/2
σ
f
L
,
P
o (GHz)
1.742
1.795
1.746
1.762
1.754
1.718
σ
f
U
,
P
o (GHz)
3.101
3.155
2.927
2.869
2.858
2.810
σ
a
E
,
P
o (×10-12 V2)
0.1064
0.1179
0.1242
0.1290
0.1316
0.1327
σ
a
O
,
P
o (×10-15 V2)
1.5912
1.6064
1.6112
1.6164
1.6199
1.6211
σ
x
0
,
P
o (μm)
1.207
1.207
1.207
1.206
1.205
1.205
σ
δ
x
,
P
o (μm)
0.170
0.170
0.170
0.170
0.170
0.170Figure 8
p
×
K scanned in individual nonlinear parameters and hyperparameters in vicinity of maximum of p×K for Δf = 1 GHz (see Table 1). (a) Zero-path difference position x0 and shift δx. (b) Prior standard deviation σb0,Pr for absolute coefficient. (c) Lower and upper limits fL and fU of spectral domain. (d) σb0,Pr and prior standard deviation σb1,Pr for linear coefficient. (e) Scaling aE and scaling aO of Brownian bridge priors for S→E and S→O. (f) σb0,Pr and prior standard deviation σb2,Pr for quadratic coefficient.


## 5.3.3. Conditional Amplitude Posterior for Maximising Settings

In the following, the amplitude posterior is investigated for the maximising settings givenΔf = 1/2 GHz (see Table 1). For the listed settings, the double-sided domain is identified by |x|≤|x0,Po|≈ 5.115 mm; the single-sided domain is bounded by the lower limit x>|x0,Po| and the upper limit of about 26.36 mm, the centre frequency fC,Po, and the bandwidth ΔF,Po of the spectral domain read 472.925 GHz and 881 GHz, respectively. In addition, the conversion of the spatial coordinate to the order of the Fourier coefficients k=|x|ΔF,Po/c reveals that the double- and single-sided domain contain information up to the 15th and 77th order (see Figures 9(c) and 9(d)), respectively. From this follows that the interferometric part of the 256/788 data values in the double-sided/full domain are modelled best by 62/310 Fourier coefficients and their associated basis functions dependent on fC,Po, ΔF,Po, and x. The signal envelopes σD,E and σD,O are evaluated with aE,Po and aO,Po.Figure 9
(a) Posterior meansS→E,Po and S→O,Po with associated 2-σ uncertainty bands evaluated at maximum of settings posterior for Δf = 1/2 GHz (see Table 1). Most likely, the spectral quantities S→E,Po and S→O,Po are finite only for the domain between fL,Po = 32.425 GHz and fU,Po = 913.425 GHz. This upper limit is approximately 4 smaller than the Nyquist frequency fNy = 3747.4 GHz. The dips in S→E,Po at 557 and 752 GHz are caused by absorption due to atmospheric water vapour. (b) D→E,Po, D→O,Po, and D→Off,Po obtained by mapping posterior means S→E,Po, S→O,Po, b0,Po=-62.142 nV, b1,Po = 5.997 μV/m, and b2,Po=-8.641 μV/m2 to data domain. (c) Cosine transform D→E,Po of S→E,Po dependent on optical path difference x and order k=|x|ΔF,Po/c (bandwidth ΔF,Po=fU,Po-fL,Po) of Fourier coefficients. The envelope σD,E∝aE,Po1/2/k, incorporating the scaling aE,Po = 1.1299 × 10−12 V2, bounds D→E,Po best. (d) Zoom of (b) dependent on x and k. Globally, the sine transform D→O,Po of S→O,Po is best summarised by the envelope σD,O∝aO,Po1/2/k with the scaling aO,Po=5.8145×10-15 V2. Since σD,E≫σD,O holds, the single-sided domain is determined mainly by S→E,Po. (e) Histogram of residuals, that is, difference between measured data D→ and posterior mean prediction D→Po=D→E,Po+D→O,Po+D→Off,Po scaled by noise level σD. The data is well described, because the mean -1.7×10-4 and the variance 0.942 of the residuals are almost vanishing and close to unity, respectively. (f) Residuals versus optical path difference.
(a)
(b)
(c)
(d)
(e)
(f)(a) Posterior Mean. Given the boundaries and increment, the spectral dimension reads Nf = 1762, and, thus, the amplitude posterior mean l→Po=S→E,Po,S→O,Po,b→PoT is 3527-dimensional. The posterior mean values for the three coefficients b0,Po=-62.142 nV, b1,Po = 5.997 μV/m, and b2,Po=-8.641 μV/m2 have magnitudes which are close to the corresponding prior standard deviations at the maximum of the settings posterior (see Table 1). Furthermore, the absolute and the linear mean values are similar to the ones obtained for the nonoptimised case (see Section 5.3.1). The offset D→Off,Po which follows from mapping b→Po is shown in Figures 9(b) and 9(d).The posterior meansS→E,Po and S→O,Po presented in Figure 9(a) approach zero towards the spectral boundaries as expected. Furthermore, S→E,Po is much larger in amplitude than S→O,Po anticipated by the quite different posterior mean scalings aE,Po and aO,Po. Hence, the even process describes most of the interferometric data in the double-sided domain. This is seen by mapping both spectral means to the data domain which gives the even and odd quantities D→E,Po and D→O,Po (see Figures 9(b)–9(d)). In addition, the single-sided domain is described mainly by S→E,Po, because in this domain only the envelope σD,E for the even process is of the order of the interferometric data (see Figure 9(d)). This differs from the findings for the nonoptimised case (aE=aO) for which S→E,Po and S→O,Po determine the single-sided region almost equally.The histogram of the residuals (see Figure9(e)) is approximated very well by the normal distribution N(-1.7×10-4,0.942). Since the mean vanishes almost, and the standard deviation is very close to unity, the data set is well described by the model and the posterior means. Most likely, the data point located at x≈ 7 mm (see Figure 9(f)) is an outlier, because its residual is outside the 3.5-σD band.(b) Posterior Covariance. For the linear parameters the posterior covariance matrix is written like(61)ΣPo=ΣS→E,E,PoΣS→E,O,PoΣS→E,b→,PoΣS→E,O,PoTΣS→O,O,PoΣS→O,b→,PoΣS→E,b→,PoTΣS→O,b→,PoTΣb→,b→,Po,and taking the square root of an element of the main diagonal gives the posterior standard deviation σPi for the ith parameter. The correlation coefficient between two parameters P1 and P2 can be evaluated by ρP1,P2=ΣP1,P2,Po/(σP1σP2).The standard deviationsσb0,Po = 0.501 nV, σb1,Po = 0.093 μV/m, and σb2,Po = 4.119 μV/m2 of the three coefficients are reduced by about two orders of magnitude with respect to the uncertainties obtained for the nonoptimised case. Though, the quadratic coefficient is inferred with small confidence. Coefficients of neighbouring degree have the highest correlation in magnitude and are anticorrelated. For example, one finds ρb0,b1=-0.55 and ρb1,b2=-0.93, while ρb0,b2 = 0.29 remains small.The covariancesΣS→E,E,Po(f,f′) and ΣS→O,O,Po(f,f′) and their associated correlations ρE,E(f,f′) and ρO,O(f,f′) resemble sinc functions centred along the line f=f′ plus a tip at this condition (see Figures 10(a)–10(c) for cross section at f′ = 100.175 GHz). Though, the FWHM widths of 1.5 GHz and 10 GHz are very different for the even and odd covariances/correlations. Regarding the covariances, the tip heights account for 1.3×10-38 V2/Hz2 and 4.4×10-40 V2/Hz2, and towards the spectral boundaries the heights are approaching zero which can be seen by the standard deviations σS→E,Po=ΣS→E,E,Po1/2(f,f) and σS→O,Po=ΣS→O,O,Po1/2(f,f) (see Figure 10(d)).Figure 10
Characterisation of posterior covariance for amplitude parameters obtained at maximum of settings posterior forΔf = 1/2 GHz (see Table 1). (a) Cross section of covariances ΣS→E,E,Po and ΣS→O,O,Po for even and odd function. Each covariance is similar to a sinc function centred at f=f′, where a tip is added on top. (b) Cross section of covariances ΣS→O,O,Po and ΣS→E,O,Po. At a fixed frequency, S→E and S→O are independent, because ΣS→E,O,Po = 0. But finite correlation remains with the neighbouring spectral domain. (c) Correlation. The even (odd) parameters are correlated (ρE(O),E(O)>0.5) over the spectral scale of about 1.5 GHz (10 GHz), while the correlation between even and odd parameters is small and asymmetric with respect to f=f′. (d) Standard deviations σS→E,Po and σS→O,Po. At the spectral limits the uncertainties approach zero. (e) and (f) Correlations of S→E(O) and b→. The correlations are oscillating and decrease generally with the frequency.
(a)
(b)
(c)
(d)
(e)
(f)At the conditionf=f′, the covariance ΣS→E,O,Po and the correlation ρE,O vanish and are asymmetrically close to this line (see Figures 10(b) and 10(c) for cross section at f′ = 100.175 GHz). This means that at a given frequency f the amplitudes S→E and S→O are independent but have a finite and opposite correlation with the neighbouring spectral domains to either side.The correlationsρE,b→ and ρO,b→ are oscillatory and rise towards the spectral origin (see Figures 10(e) and 10(f)).(c) Posterior Samples. From the conditional amplitude posterior samples S→E,Po,Sa, S→O,Po,Sa and b→Po,Sa are drawn. While Figure 11(a) shows one sample for each of the even and odd spectral functions, Figure 11(b) presents 100 samples. The samples form a band around the corresponding posterior mean S→E(O),Po with the width of about twice the standard deviation σS→E(O),Po as shown in Figure 10(d). Hence, the band for the odd function is much smaller.Figure 11
(a) Individual posterior sample functionsS→E,Po,Sa and S→O,Po,Sa. (b) 100 posterior sample functions S→E,Po,Sa and S→O,Po,Sa. The uncertainty at a given frequency is larger for the even function, because for the associated scalings of the two processes aE,Po≫aO,Po holds. (c) 100 spectral samples and 100 offset coefficients mapped to probed data domain. Unlike for the nonoptimised case (see Figure 5(b)), the spread is much smaller, and the quantities connect smoothly at the transition between double- and single-sided regions. (d) 100 samples D→Po,Sa and measured data D→. The samples are in agreement with the uncertainty of the measurements. (e) One sample D→E,Po,Sa and D→O,Po,Sa evaluated inside probed domain (x < 26.36 mm) and outside. With respect to the envelopes σD,E(O), the contributions behave similarly in both domains. (f) 100 samples D→E,Po,Sa and D→O,Po,Sa evaluated inside the probed domain and outside. Outside the probed domain, the area between the envelops σD,E(O) is filled by D→E,Po,Sa and D→O,Po,Sa. Because the spread in the nonprobed domain is wider than the one in the probed domain, the main uncertainty in the even and odd spectra originates in nonprobed Fourier coefficients.
(a)
(b)
(c)
(d)
(e)
(f)The samples mapped to the data domain giveD→E,Po,Sa, D→O,Po,Sa, and D→Off,Po,Sa (see Figure 11(c)) which form much more narrow bands when compared to the nonoptimised case (see Figure 5(b)). In particular, the transition from the double- to the single-sided domain is smooth. This is a consequence of having used the most likely scalings aE,Po and aO,Po and the associated but quite different signal envelops inferred from the data. Furthermore, this explains why the band for D→O,Po,Sa is much smaller than the one for D→E,Po,Sa.For each sample, the sumD→Po,Sa=D→E,Po,Sa+D→O,Po,Sa+D→Off,Po,Sa complies with the measured data D→ and the noise level (see Figure 11(d)).With the samplesS→E,Po,Sa and S→O,Po,Sa at hand, D→E,Po,Sa and D→O,Po,Sa are predicted outside the probed spatial domain. Up to the optical path difference of 70 mm, Figures 11(e) and 11(f) show D→E,Po,Sa and D→O,Po,Sa for one sample and 100 samples, respectively. No difference is obvious in D→E,Po,Sa and D→O,Po,Sa for an individual prediction in the probed and nonprobed domain. However, many samples take values between σD,E(O) outside the probed region. Since the spread in the probed domain is smaller than in the nonprobed domain, the main uncertainty in S→E,Po and S→O,Po originates in Fourier coefficients not probed rather than the noise on measured data.


## 5.3.4. Marginal Amplitude Posterior

Since only the marginal posteriorpS→E,S→O∣Δf,D→ of the spectral functions is of interest for this specific Fourier spectroscopy problem, formally the marginalisation(62)pS→E,S→O∣Δf,D→=∫pS→E,S→O,b→∣n→,h→,Δf,D→pn→,h→∣Δf,D→db→dn→dh→needs to be carried out. In general, the pdf pS→E,S→O∣Δf,D→ is not Gaussian but can be approximated by a multivariate normal. Performing the marginalisation, that is, the uncertainty propagation, is analytically not possible but achievable numerically described below for Δf = 1/2 GHz.To reduce the numerical efforts, only the most important parametersfL, fU, x0, δx, aE, and aO are taken into consideration for the marginalisation, and the remaining three parameters σb0,Pr, σb1,Pr, and σb2,Pr are held at the maximum values listed in Table 1. The remaining pdf, labeled by p′n→,h→∣Δf,D→ in the following, is well approximated by a factorisable multivariate normal, because the corresponding covariances are negligible (see Section 5.3.2). The posterior mean values and standard deviations are listed in Tables 1 and 2, respectively.From the six-dimensional Gaussian posteriorp′n→,h→∣Δf,D→, 1000 samples are drawn. For each sample, mean and covariance for the conditional amplitude posterior are evaluated. Thereby, a difficulty arises due to the lower and upper spectral limits fL and fU. For given increment Δf, fL and fU determine the spectral grid. To keep this grid fixed, the spectral grid obtained for fL,Po = 32.425 GHz and fU,Po = 913.425 GHz is extended by some tens of GHz on either side. Furthermore, the posterior samples for fL and fU are compared with the extended but fixed spectral grid points and changed to the closest grid point. With the modified boundaries and the remaining samples unchanged, a conditional amplitude posterior follows. For example, Figure 12(a) shows the mean functions S→E,Po and S→O,Po for the 1000 posterior samples of the settings. While S→E,Po changes little but most at the spectral boundaries, S→O,Po varies considerably, mostly in the centre of the spectral domain covered. This is caused mainly by the posterior uncertainty of the zero-path difference x0. However, the spread is symmetric around the most likely mean.Figure 12
(a) MeansS→E,Po (blue) and S→O,Po (red) obtained for 1000 samples of settings posterior at Δf = 1/2 GHz. Only at the spectral limits the means deviate from the most likely mean S→E,Po (cyan) evaluated at the maximum of the settings posterior. The uncertainty in the zero-path difference x0 causes a considerable deviation from the most likely mean S→O,Po (black). (b) Means S→E,Po,Ma (cyan) and S→O,Po,Ma (black) after performing numerical marginalisation of settings parameters at Δf = 1/2 GHz. The marginalised means are similar to the means S→E,Po and S→O,Po (see (a)) at the maximum of the settings posterior except for minor differences towards the spectral boundaries. Having the marginalised covariance at hand, 100 samples S→E,Po,Ma,Sa (blue) and S→O,Po,Ma,Sa (red) are drawn from the marginalised amplitude posterior. Only in the centre of the spectral domain the spread of S→O,Po,Ma,Sa increases to capture the uncertainty in x0 as well (see (a)).
(a)
(b)For each settings sample, 100 samples are drawn frompS→E,S→O,b→∣n→,h→,Δf,D→. Overall, 100000 samples are available for b→, S→E, and S→O at each spectral grid point. From this large sample set, the marginalised means S→E,Po,Ma, S→O,Po,Ma and covariance ΣS→,Po,Ma are evaluated which capture the first and second moments of the pdf pS→E,S→O∣Δf,D→. In doing so, the marginalisation with respect to b→ is achieved implicitly.(a) Posterior Mean. The means S→E,Po,Ma and S→O,Po,Ma (see Figure 12(b)) look similar to S→E,Po, S→O,Po at the maximum of the settings posterior. However, S→E,Po,Ma and S→O,Po,Ma approach zero some GHz below and above fL,Po and fU,Po, respectively.(b) Posterior Covariance. The posterior covariance after marginalisation with respect to the setting parameters is decomposed like(63)ΣS→,Po,Ma=ΣS→E,E,Po,MaΣS→E,O,Po,MaΣS→E,O,Po,MaTΣS→O,O,Po,Ma.ΣS→E,E,Po,Ma (see Figure 13(a) for cross section at f′ = 100.175 GHz) and is very similar to ΣS→E,E,Po evaluated at the maximum of the settings posterior (see Section 5.3.3). The same holds for ΣS→E,O,Po,Ma and ΣS→E,O,Po. However, a significant change is found for ΣS→O,O,Po,Ma (see Figures 13(a) and 13(b) for some cross sections) when compared to ΣS→O,O,Po. Due to the posterior uncertainty in x0, a broadband additive contribution increases ΣS→O,O,Po,Ma up to about one order of magnitude in the centre of the spectral domain. Comparing the standard deviations σS→E,Po,Ma with σS→E,Po, and σS→O,Po,Ma with σS→O,Po (see Figure 13(c)) shows that towards the spectral limits the propagation of the settings posterior uncertainties has little effect. Furthermore, the correlation ρO,O,Ma is high for a broad spectral range in the central region (see Figure 13(d)).Figure 13
Characterisation of posterior covarianceΣS→,Po,Ma numerically marginalised with respect to settings parameters. ΣS→,Po,Ma is decomposed, so that ΣS→E(O),E(O),Po,Ma captures the covariance of the even (odd) spectral quantities. (a) Cross section at frequency f′ = 100.175 GHz. Apart from the extension in the spectral domain, ΣS→E,E,Po,Ma is similar to the covariance ΣS→E,E,Po evaluated at the maximum of the settings posterior (see Figure 10(a)). ΣS→O,O,Po,Ma shows a significant deviation from ΣS→O,O,Po (see also (b)). (b) Cross sections of ΣS→O,O,Po,Ma at frequencies f′ = 100.175 GHz, 400.175 GHz, and 850.175 GHz. Towards the spectral limits ΣS→O,O,Po,Ma behaves like ΣS→O,O,Po. But in the central region a broadband additive contribution is evident originating in the posterior uncertainty of the zero-path difference x0. (c) Standard deviations evaluated from main diagonal of marginalised covariance ΣS→,Po,Ma (red) and covariance ΣS→,Po (black, see Figure 10(d)) evaluated at maximum of settings posterior. The odd quantity increases by a factor of three in the central region after the marginalisation. (d) Correlation ρO,O,Ma at three frequencies. Compared to ρO,O (see Figure 10(c)), a broadband (several 100s GHz) elevation of ρO,O,Ma up to 0.8 follows from the uncertainty in x0 for the central domain.
(a)
(b)
(c)
(d)(c) Posterior Samples. From the marginalised Gaussian posterior pS→E,S→O∣Δf,D→ 100 samples (S→E,Po,Ma,Sa and S→O,Po,Ma,Sa) are drawn and shown in Figure 12(b). As expected, the spread of S→E,Po,Ma,Sa around the mean S→E,Po,Ma is almost equal to the one for S→E,Po,Sa around the mean S→E,Po (see Figure 11(b)). The samples S→O,Po,Ma,Sa have a wider spread in the centre of the spectral domain than S→O,Po,Sa mainly caused by the uncertainty in the zero-path difference.


## 5.4. Figure of Merit for Real-World Interferometer

For the ideal Martin-Puplett interferometer, the odd spectral function, and, hence, the odd process must vanish from theoretical point of view. However, imperfections of a real-world interferometer leave the odd contribution finite in general. This imperfection is captured by the scalingaO and in terms of signal by the envelope σD,O∝aO1/2. By relating the square root of the scalings like(64)FS=aE1/2aE1/2+aO1/2,one can define the figure of merit FS which expresses by a number the signal deviation of a real-world instrument from the ideal case. For an ideal interferometer, FS = 1 holds (aO=0). For the interferometer investigated here, the settings posterior (see Section 5.3.2) carries the information to state the mean as FS = 0.933 with the uncertainty of about 0.01. Hence, about 7% of the signal is converted from the ideal even process to the odd process by the real-world diagnostic.


## 6. Discussion

### 6.1. Choice of Spectral Priors or Model Plausibility

The results obtained in the previous section rely on the modelMBB presented here with the assumption that the even and odd spectral functions can be described each by a Brownian bridge process and its associated prior covariance. An alternative model, let us say MA, with certain assumptions on the spectral functions, leads to different prior covariances and to a different posterior. In addition, the model MA might be more or less plausible when compared to the model MBB.In principle, the plausibility of a model relative to an alternative can be investigated within the Bayesian framework by rising the abstraction level. Starting from (36), a further factorisation needs to be carried out with respect to the used models MBB or MA. Basically, one can assign the model posterior by pMBB∣D→∝KBB′p(MBB) and pMA∣D→∝KA′p(MA) with the model priors p(MA) and p(MBB). The dimensionless constants KBB′ and KA′ may be obtained by marginalising over all linear and nonlinear parameters and hyperparameters. Then, the model plausibility is captured by the ratio pMBB∣D→/pMA∣D→. This ratio becomes KBB′/KA′, if no model is preferred a priori for which one sets p(MA)=p(MBB)=1/2. Such a model plausibility study was not carried out, because it demands an investigation on its own alongside with a costly numerical treatment. However, some aspects of a plausibility study will be discussed below.The signal envelope, corresponding to the used prior covariance for the even and odd spectral functions, is expected to be an indicator for a competitive model. This envelope should be able to resemble the global trend of the interferometric data (see Figure3(a)). The investigated model MBB seems to have these desired characteristics (see Figure 9).For the even and odd spectral functions, an alternative prior choice could beΣA(f,f′)=g(f)δ(f-f′) which assigns no correlation. For example, if the function g(f) is chosen to remain constant or as a triangular function, centred at fC and approaching 0 towards the spectral limits, then the corresponding covariance ΣA(k,k′) for the Fourier coefficients has a constant value along the main diagonal for all orders. This can be shown analytically by performing the operation stated in (43) on ΣA(f,f′). As a consequence, no decay is imposed on the Fourier coefficients with rising order which is incompatible with square-integrability. Furthermore, the alternative signal envelop is constant in the optical path difference domain, opposing the fall-off in the interferometric data (see Figure 3(a)). Thus, the data should not be described better by either of the two suggested alternative models.A competitive model could use a prior covarianceΣC(k,k′) for the Fourier coefficients which has a dropping amplitude when the order of the coefficients rises. For instance, the main diagonal of ΣC(k,k′) could be chosen like 1/(kk′)e, and a transformation back to the spectral domain would allow further investigations of the properties of ΣC(f,f′). For different but positive exponents e, the model plausibility could be examined.


### 6.2. Comparison with Standard Model

The standard analysis approach for the interferometer investigated here relies on a different modelMS which is set up hierarchically [12], and no covariances of the parameters involved are available. Basically, the quantity of interest, that is, the spectrum, is determined conditionally on the inferred voltage offset VOff, the zero-path difference x0, while setting the shift δx=0 and the phase. Furthermore, spatial window functions are multiplied to the data to retrieve the phase and spectrum, using discrete fast Fourier transformation (DFFT) routines, for the Nyquist assumptions. A consideration of nonprobed Fourier coefficients, locating outside the experimentally accessible spatial domain, is missing completely as well as the influence of the measurement noise. Hence, a comparison of the standard model MS with the Bayesian model MBB is not straightforward. In order to perform a reasonable comparison, the spectral grid with the increment Δf = 3.66 GHz, which follows from the use of DFFT by MS, and the limits fL = 0 GHz, fU = 3747.4 GHz are the same for both models. In addition, settings are determined for MBB which come close to the settings for MS (see second column of Table 3) to resemble the standard model. Given these settings, the even and odd spectra are compared. Then, the plausibility of these settings can be obtained with respect to the maximum of the settings posterior for MBB by evaluating ln⁡p×K and the corresponding odds.Table 3
Comparison of standard modelMS and Bayesian model MBB for spectral increment Δf = 3.66 GHz. MS (second column) is well approximated by MBB when the settings of the third column are chosen. However, an aliasing feature is present in the spectral domain 3000–3500 GHz (see Figure 14(a)). This feature disappears when the spatial settings are changed (fourth column) to the values at the maximum of the settings posterior of MBB (fifth column). The maximum of the settings posterior reveals that the Nyquist assumptions (fL = 0 GHz and fU = 3747.4 GHz) are not plausible, because an overfitting of the data is made.
Model
M
S
M
B
B
M
B
B
M
B
B
f
L (GHz)
0
0
0
33.614
f
U (GHz)
3747.4
3747.4
3747.4
911.912
a
E (×10-12 V2)
N/A
19.8182
19.8182
1.0304
a
O (×10-15 V2)
N/A
129.8670
129.8670
5.7641
x
0 (mm)
-
5.118
-
5.118
-
5.115285
-
5.115285
δ
x (μm)
0
0
-
4.814
-
4.814
σ
b
0
,
P
r (nV)
N/A
61.735
61.735
62.146
σ
b
1
,
P
r (μV/m)
N/A
5.822
5.822
6.007
σ
b
2
,
P
r (μV/m2)
N/A
0.000087
0.000087
10.102
ln
⁡
p
×
K
N/A
1051396.55
1051400.41
1052037.88
Odds
N/A
1
0
-
278.53
1
0
-
276.84
1
Aliasing feature
Yes
Yes
No
No
Overfitting data
Yes
Yes
Yes
NoThe only window function applied here in the standard analysis weighs the single-sided data domain twice as large as the double-sided domain. Furthermore, the settings of the standard approachx0=-5.118 mm and vanishing δx are transferred to the model MBB, and the remaining settings of MBB are optimised (see third column of Table 3). The even and odd quantities S→E,S and S→O,S are inferred with MS, and the means S→E,Po and S→O,Po of the conditional amplitude posterior of the model MBB are evaluated. Figure 14(a) shows the results which are similar in amplitude. For both models, an aliasing feature can be found for the odd spectral quantity in the spectral domain from 3000 GHz to 3500 GHz. This feature originates in the assumed uniform optical path difference grid (δx=0). The differences S→E,S-S→E,Po and S→O,S-S→O,Po settle mostly in the range ±0.1×10-18 V/Hz (see Figure 14(b)).Figure 14
Comparison of even and odd spectral quantities inferred from standard modelMS and Bayesian model MBB. MS exploits discrete fast Fourier transformation routines to evaluate the even and odd spectra S→E,S and S→O,S conditionally on the spectral limits fL = 0 GHz and fU = 3747.4 GHz (Nyquist assumptions), the spectral increment Δf = 3.66 GHz, vanishing shift δx, and the zero-path difference x0=-5.118 mm. Conditional on these settings, the model MBB obtains the means S→E,Po and S→O,Po of the conditional amplitude posterior when the remaining settings are taken at their maximum values (see third column of Table 3). (a) Even and odd spectra. Similar results follow from both models. The elevated amplitude in the odd spectra between 3000 GHz and 3500 GHz is caused by aliasing due to the assumption that the optical path difference grid is uniform. (b) Absolute differences of spectral quantities inferred by both models. Except in the vicinity of the water vapour absorption line at 557 GHz, the difference is small in absolute terms.
(a)
(b)Keeping fixed all settings inMBB but choosing x0 and δx which are present at the maximum of the settings posterior for Δf = 3.66 GHz (see forth column of Table 3), the aliasing feature vanishes completely, because the nonuniformity in the optical path difference grid is taken into account properly. This small change in the spatial settings would make the conditional amplitude posterior about 47 times more likely.The odds 1 :10-278.53 for the values at the maximum of the settings posterior (see fifth column of Table 3) with respect to the settings for the model MBB which mimics MS mark the Nyquist assumptions made by the standard approach very unlikely. This is explained when the corresponding residuals are investigated which show an overfitting of the data if the Nyquist assumptions are used (similar to Figure 5(c)).


### 6.3. Computational Time

The algorithms for the modelMBB and the standard model MS are implemented in Scilab [13]. To obtain the even and odd spectra with MS via DFFT routines, a computational time of about 4 ms is measured. This fast analysis time is exceeded by several orders of magnitude, when the problem is investigated with the Bayesian model.For the modelMBB, the number Nl=2(fU-fL)/Δf+3 of linear parameters, dependent on the spectral domain and increment, gives the dimension of the prior and posterior covariance matrices. Hence, Nl determines the computational times tl, to evaluate the mean and covariance of the conditional amplitude posterior, and tln⁡p×K, to investigate the settings posterior at one point in the parameter space. For the implemented algorithm and the maximising settings listed in Table 1 with decreasing Δf (Nl increases by about one order of magnitude), tl and tln⁡p×K measured increase at least quadratically with Nl (see Table 4). This is caused mainly by the need to invert numerically prior and posterior covariance matrices.Table 4
Measured computational timestl and tln⁡p×K of used algorithm for Bayesian model MBB dependent on number Nl of linear parameters and discretisation increment Δf for maximising settings listed in Table 1. tl and tln⁡p×K rise approximately with Nl2. The investigation of the settings posterior demands a time of tln⁡p×K multiplied with a factor of at least 103 which becomes hours for Δf = 1/2 GHz.
Δ
f (GHz)
≈5.68
4
3
2
1
1/2
N
l
313
443
589
883
1765
3527
t
l (s)
0.072
0.097
0.180
0.421
2.084
11.231
t
ln
⁡
p
×
K (s)
0.079
0.108
0.200
0.471
2.406
13.648The characterisation of the settings posterior (see Section5.3.2) requires a duration of 103tln⁡p×K at least which adds up to hours for small Δf (large Nl) and to much less than one hour for Δf = 4 GHz.The numerical marginalisation described in Section5.3.4 takes about half a day for Δf = 1/2 GHz.


## 6.1. Choice of Spectral Priors or Model Plausibility

The results obtained in the previous section rely on the modelMBB presented here with the assumption that the even and odd spectral functions can be described each by a Brownian bridge process and its associated prior covariance. An alternative model, let us say MA, with certain assumptions on the spectral functions, leads to different prior covariances and to a different posterior. In addition, the model MA might be more or less plausible when compared to the model MBB.In principle, the plausibility of a model relative to an alternative can be investigated within the Bayesian framework by rising the abstraction level. Starting from (36), a further factorisation needs to be carried out with respect to the used models MBB or MA. Basically, one can assign the model posterior by pMBB∣D→∝KBB′p(MBB) and pMA∣D→∝KA′p(MA) with the model priors p(MA) and p(MBB). The dimensionless constants KBB′ and KA′ may be obtained by marginalising over all linear and nonlinear parameters and hyperparameters. Then, the model plausibility is captured by the ratio pMBB∣D→/pMA∣D→. This ratio becomes KBB′/KA′, if no model is preferred a priori for which one sets p(MA)=p(MBB)=1/2. Such a model plausibility study was not carried out, because it demands an investigation on its own alongside with a costly numerical treatment. However, some aspects of a plausibility study will be discussed below.The signal envelope, corresponding to the used prior covariance for the even and odd spectral functions, is expected to be an indicator for a competitive model. This envelope should be able to resemble the global trend of the interferometric data (see Figure3(a)). The investigated model MBB seems to have these desired characteristics (see Figure 9).For the even and odd spectral functions, an alternative prior choice could beΣA(f,f′)=g(f)δ(f-f′) which assigns no correlation. For example, if the function g(f) is chosen to remain constant or as a triangular function, centred at fC and approaching 0 towards the spectral limits, then the corresponding covariance ΣA(k,k′) for the Fourier coefficients has a constant value along the main diagonal for all orders. This can be shown analytically by performing the operation stated in (43) on ΣA(f,f′). As a consequence, no decay is imposed on the Fourier coefficients with rising order which is incompatible with square-integrability. Furthermore, the alternative signal envelop is constant in the optical path difference domain, opposing the fall-off in the interferometric data (see Figure 3(a)). Thus, the data should not be described better by either of the two suggested alternative models.A competitive model could use a prior covarianceΣC(k,k′) for the Fourier coefficients which has a dropping amplitude when the order of the coefficients rises. For instance, the main diagonal of ΣC(k,k′) could be chosen like 1/(kk′)e, and a transformation back to the spectral domain would allow further investigations of the properties of ΣC(f,f′). For different but positive exponents e, the model plausibility could be examined.


## 6.2. Comparison with Standard Model

The standard analysis approach for the interferometer investigated here relies on a different modelMS which is set up hierarchically [12], and no covariances of the parameters involved are available. Basically, the quantity of interest, that is, the spectrum, is determined conditionally on the inferred voltage offset VOff, the zero-path difference x0, while setting the shift δx=0 and the phase. Furthermore, spatial window functions are multiplied to the data to retrieve the phase and spectrum, using discrete fast Fourier transformation (DFFT) routines, for the Nyquist assumptions. A consideration of nonprobed Fourier coefficients, locating outside the experimentally accessible spatial domain, is missing completely as well as the influence of the measurement noise. Hence, a comparison of the standard model MS with the Bayesian model MBB is not straightforward. In order to perform a reasonable comparison, the spectral grid with the increment Δf = 3.66 GHz, which follows from the use of DFFT by MS, and the limits fL = 0 GHz, fU = 3747.4 GHz are the same for both models. In addition, settings are determined for MBB which come close to the settings for MS (see second column of Table 3) to resemble the standard model. Given these settings, the even and odd spectra are compared. Then, the plausibility of these settings can be obtained with respect to the maximum of the settings posterior for MBB by evaluating ln⁡p×K and the corresponding odds.Table 3
Comparison of standard modelMS and Bayesian model MBB for spectral increment Δf = 3.66 GHz. MS (second column) is well approximated by MBB when the settings of the third column are chosen. However, an aliasing feature is present in the spectral domain 3000–3500 GHz (see Figure 14(a)). This feature disappears when the spatial settings are changed (fourth column) to the values at the maximum of the settings posterior of MBB (fifth column). The maximum of the settings posterior reveals that the Nyquist assumptions (fL = 0 GHz and fU = 3747.4 GHz) are not plausible, because an overfitting of the data is made.
Model
M
S
M
B
B
M
B
B
M
B
B
f
L (GHz)
0
0
0
33.614
f
U (GHz)
3747.4
3747.4
3747.4
911.912
a
E (×10-12 V2)
N/A
19.8182
19.8182
1.0304
a
O (×10-15 V2)
N/A
129.8670
129.8670
5.7641
x
0 (mm)
-
5.118
-
5.118
-
5.115285
-
5.115285
δ
x (μm)
0
0
-
4.814
-
4.814
σ
b
0
,
P
r (nV)
N/A
61.735
61.735
62.146
σ
b
1
,
P
r (μV/m)
N/A
5.822
5.822
6.007
σ
b
2
,
P
r (μV/m2)
N/A
0.000087
0.000087
10.102
ln
⁡
p
×
K
N/A
1051396.55
1051400.41
1052037.88
Odds
N/A
1
0
-
278.53
1
0
-
276.84
1
Aliasing feature
Yes
Yes
No
No
Overfitting data
Yes
Yes
Yes
NoThe only window function applied here in the standard analysis weighs the single-sided data domain twice as large as the double-sided domain. Furthermore, the settings of the standard approachx0=-5.118 mm and vanishing δx are transferred to the model MBB, and the remaining settings of MBB are optimised (see third column of Table 3). The even and odd quantities S→E,S and S→O,S are inferred with MS, and the means S→E,Po and S→O,Po of the conditional amplitude posterior of the model MBB are evaluated. Figure 14(a) shows the results which are similar in amplitude. For both models, an aliasing feature can be found for the odd spectral quantity in the spectral domain from 3000 GHz to 3500 GHz. This feature originates in the assumed uniform optical path difference grid (δx=0). The differences S→E,S-S→E,Po and S→O,S-S→O,Po settle mostly in the range ±0.1×10-18 V/Hz (see Figure 14(b)).Figure 14
Comparison of even and odd spectral quantities inferred from standard modelMS and Bayesian model MBB. MS exploits discrete fast Fourier transformation routines to evaluate the even and odd spectra S→E,S and S→O,S conditionally on the spectral limits fL = 0 GHz and fU = 3747.4 GHz (Nyquist assumptions), the spectral increment Δf = 3.66 GHz, vanishing shift δx, and the zero-path difference x0=-5.118 mm. Conditional on these settings, the model MBB obtains the means S→E,Po and S→O,Po of the conditional amplitude posterior when the remaining settings are taken at their maximum values (see third column of Table 3). (a) Even and odd spectra. Similar results follow from both models. The elevated amplitude in the odd spectra between 3000 GHz and 3500 GHz is caused by aliasing due to the assumption that the optical path difference grid is uniform. (b) Absolute differences of spectral quantities inferred by both models. Except in the vicinity of the water vapour absorption line at 557 GHz, the difference is small in absolute terms.
(a)
(b)Keeping fixed all settings inMBB but choosing x0 and δx which are present at the maximum of the settings posterior for Δf = 3.66 GHz (see forth column of Table 3), the aliasing feature vanishes completely, because the nonuniformity in the optical path difference grid is taken into account properly. This small change in the spatial settings would make the conditional amplitude posterior about 47 times more likely.The odds 1 :10-278.53 for the values at the maximum of the settings posterior (see fifth column of Table 3) with respect to the settings for the model MBB which mimics MS mark the Nyquist assumptions made by the standard approach very unlikely. This is explained when the corresponding residuals are investigated which show an overfitting of the data if the Nyquist assumptions are used (similar to Figure 5(c)).


## 6.3. Computational Time

The algorithms for the modelMBB and the standard model MS are implemented in Scilab [13]. To obtain the even and odd spectra with MS via DFFT routines, a computational time of about 4 ms is measured. This fast analysis time is exceeded by several orders of magnitude, when the problem is investigated with the Bayesian model.For the modelMBB, the number Nl=2(fU-fL)/Δf+3 of linear parameters, dependent on the spectral domain and increment, gives the dimension of the prior and posterior covariance matrices. Hence, Nl determines the computational times tl, to evaluate the mean and covariance of the conditional amplitude posterior, and tln⁡p×K, to investigate the settings posterior at one point in the parameter space. For the implemented algorithm and the maximising settings listed in Table 1 with decreasing Δf (Nl increases by about one order of magnitude), tl and tln⁡p×K measured increase at least quadratically with Nl (see Table 4). This is caused mainly by the need to invert numerically prior and posterior covariance matrices.Table 4
Measured computational timestl and tln⁡p×K of used algorithm for Bayesian model MBB dependent on number Nl of linear parameters and discretisation increment Δf for maximising settings listed in Table 1. tl and tln⁡p×K rise approximately with Nl2. The investigation of the settings posterior demands a time of tln⁡p×K multiplied with a factor of at least 103 which becomes hours for Δf = 1/2 GHz.
Δ
f (GHz)
≈5.68
4
3
2
1
1/2
N
l
313
443
589
883
1765
3527
t
l (s)
0.072
0.097
0.180
0.421
2.084
11.231
t
ln
⁡
p
×
K (s)
0.079
0.108
0.200
0.471
2.406
13.648The characterisation of the settings posterior (see Section5.3.2) requires a duration of 103tln⁡p×K at least which adds up to hours for small Δf (large Nl) and to much less than one hour for Δf = 4 GHz.The numerical marginalisation described in Section5.3.4 takes about half a day for Δf = 1/2 GHz.


## 7. Conclusions

The Fourier transform is the heart of Fourier spectroscopy applications. Thereby, the interferometric data has a linear dependence on the even and odd continuous spectra to be inferred. Standard analysis techniques lack appropriate handling of fundamental aspects like noisy measurements, the influence of nonprobed spatial domains linked to Fourier coefficients above a certain order, the estimation of spectral limits, and the propagation of uncertainties of additional parameters like the zero-path difference onto the inferred spectra. For instance, the Nyquist assumption implies the fundamental misconception that the upper spectral limit of spectra to be inferred would depend on the spatial sampling. In addition, a broad spectral bandwidth would follow which increases artificially the number of Fourier coefficients necessary to describe the data. On the contrary, it can be shown analytically that a band-limitation causes spatially extended basis functions (modulated sinc functions) assigned to the Fourier coefficients in the data domain. Thus, several nearby data points are captured sufficiently by less coefficients. This example demonstrates that interferometric data contains more information than usually extracted.As an alternative to the standard analysis techniques, a probabilistic ansatz, relying on Bayes’ theorem, was proposed which is able to capture the fundamental aspects listed above. In general, Bayes’ theorem relates the posterior probability density function of model parameters to the product of the likelihood and the prior probability density function for these parameters. The ansatz presented here uses multivariate normal distributions for the likelihood and the prior for parameters which map linearly to the data domain. This gives straightforwardly an analytical solution for the posterior of these linear parameters in form of a multivariate normal. Though, this amplitude posterior is conditional on the settings parameters, summarising all nonlinear model parameters and hyperparameters. After the trivial marginalisation over the linear parameters, the remaining quantity can be scanned in the settings parameters to investigate their joint posterior. This can be understood as a means of applying Ockham’s razor for the linear problem. With the settings posterior at hand, the marginalisation projects the uncertainties in the settings onto the linear parameters.The example application for the Bayesian approach infers even and odd spectra, which qualify as linear parameters, in the microwave and far-infrared spectral domain and several settings parameters, like the spectral discretisation increment, the spectral limits, the scalings of the even and odd processes, the zero-path difference, and a shift correction to the spatial sampling, given a measured interferometric data set. Each spectrum is modelled by a scaled Brownian bridge process which is able to capture a band-limitation, and the associated covariance is used in the Gaussian prior. This covariance assigns a broadband correlation, but its transform to the domain of Fourier coefficients reveals no correlation (vanishing off-diagonal elements except in connection with the zeroth-order term) between the coefficients. Furthermore, the diagonal elements drop with the square of the order of the coefficients. Hence, the prior information stated by the Brownian bridge covariance considers functions which are square-integrable and, thus, converge globally in the limit when the discretisation increment approaches zero and the order of the Fourier coefficients tends to infinity. In addition, these functions vanish smoothly at the lower and upper spectral boundaries. In the data domain, a signal envelope follows from the Brownian bridge process. This envelope decays with the optical path difference and the spectral bandwidth.For the linear parameters like the even and odd spectra, a conditional amplitude posterior was briefly examined, relying on the Nyquist assumptions. Due to the large upper spectral limit, all noise contributions to the interferometric data are captured by the posterior mean of the linear parameters. This implies an overfitting. Because large and equal values are taken for the two Brownian bridge scalings (large signal envelops), the mapped posterior means of the spectra describe the even and odd parts of the interferometric data to equal parts in the single-sided domain, while the even part dominates in the double-sided domain. This is an indicator that the Fourier coefficients located in the single-sided domain are underestimated and overestimated for the even and odd spectra, respectively. The posterior samples for both spectra show large deviations from the means, and the even and odd contributions, obtained by mapping the samples, form much wider bands than the measurement uncertainty, especially in the single-sided domain. This indicates an unnecessary expanded solution space for the problem. Only by the posterior covariance of the linear parameters, the sum of the mapped samples complies with the data and its uncertainty band. The listed features mark a very unlikely conditional amplitude posterior which is revealed by the settings posterior.The settings posterior for the most important settings is well approximated by the product of individual normal distributions, because no significant correlations could be found. The corresponding posterior means and standard deviations take reasonable values. These values are affected little by the discretisation increment which tends to be small, confirming the proposition that continuous spectral quantities are probed. The upper spectral limit is about a factor four smaller than the Nyquist frequency, and the lower limit is well separated from zero. This reduction of the bandwidth implies that the interferometric data can be described by a number of Fourier coefficients with associated spatial basis functions which is about one-quarter of the amount of data points. The scaling of the even process exceeds the one for the odd process by about two orders of magnitude.For values corresponding to the maximum of the settings posterior at a small discretisation increment, the conditional amplitude posterior was investigated. By the discretisation, the number of the linear parameters exceeds the one of the Fourier coefficients, which is mandatory to describe the data points within the measurement uncertainty, by one order of magnitude. However, the Ockham’s razor principle implemented by the settings posterior limits the solution space, so that the posterior means and samples for the linear spectral parameters, mapped to the data domain, have a smooth transition between the double- and single-sided regions. Due to the much larger scaling with respect to the odd process, the even process and, thus, the even spectrum describe most of the single-sided and double-sided interferogram region. While the probed interferometric data is well described within the uncertainty by the means and samples, the nonprobed data domain, corresponding to Fourier coefficients above a certain order, is filled broadly by these mapped samples. This filling decays with increasing optical path difference and is limited by the signal envelopes which follow from the estimated scalings of the Brownian bridge processes and the bandwidth. Because the spread of the mapped samples in the nonprobed domain exceeds the one in the probed region, the main uncertainties for the even and odd spectra originates in nonprobed Fourier coefficients.The numerically costly marginalisation over some of the settings shows that the zero-path difference changes the covariance for the odd spectral parameters significantly. Basically, a broad increase of the posterior uncertainties and correlation was found.A figure of merit was introduced which states the deviation of a real-world interferometer from an ideal diagnostic. By relating the scalings of the even and odd processes, the used interferometer is characterised as being close to the ideal case for which the odd process must vanish.



---

*Source: 9265084-2017-10-31.xml*