Data acquisition system for optical frequency comb spectroscopy

Data acquisition system for optical frequency comb spectroscopy

Author

Ragnar Seton

Supervisor

Amir Khodabakhsh

Examiners

Aleksandra Foltynowicz-Matyba Ove Axner

2017-09-05

-Richard Moore

ASCII American Standard Code for Information Interchange (text format)

BS Beam Splitter CoW Copy on Write

CUDA Compute Unied Device Architecture (NVIDIAs parallel programming framework)

DAC Data Acquisition Card DAQS Data Acquisition System FFT Fast Fourier Transform FIFO First In, First Out

FLOPS Floating-Point Operations Per Second FTS Fourier Transform Spectrometer FWHM Full Width Half Maximum GUI Graphical User Interface

GUIDE Graphical User Interface Development Environment (part of MATLAB)

HAWKS HITRAN Atmospheric Workstation HITRAN High Resolution Transmission database HLA High Level Assembly

HTML Hyper Text Markup Language HWHM Half-Width Half-Maximum

I/O Input/Output

IPC Inter Process Communication IVI Interchangeable Virtual Instrument JIT Just-In-Time (compilation)

JVM Java Virtual Machine JDBC Java Database Connectivity MI Michelson Interferometer

NI National Instruments corporation NICE Noise Immune Cavity Enhanced

NIR Near Infrared (commonly used for wavelength spectra) ODBC Open Database Connectivity

OFCS Optical Frequency Comb Spectroscopy OOP Object Oriented Programming

OPD Optical Path Dierence

PCI Peripheral Component Interconnect

PHP PHP: Hypertext Preprocessor (recursive backronym) RF Radio Frequency

SLOC Source Lines Of Code UX User Experience

The Optical Frequency Comb Spectroscopy (OFCS) Group at the Department of Physics at Umeå University develops new techniques for extremely high sensitivity trace gas detection, non invasive temperature measurements, and other applications of OFCS.

Their setup used primarily for trace gas detection contains several components that have been developed in-house, including a Fourier Transform Spectrometer (FTS) and an auto-balancing detector. This is the one used in this thesis work and it includes a high frequency data acquisition card (DAC) recording interferograms in excess of 10⁷ double- precision oating point samples per sweep of the FTS's retarder. For acquisition and analysis to be possible in both directions of the retarder the interferograms needs to be analysed in a sub-second timeframe, something not possible with the present software.

The aim of this thesis work has thus been to develop a system with optimized analysis implementations in MATLAB. The latter was a prerequisite from the group to ensure maintainability, as all members are well acquainted with it.

Fullling its primary purpose MATLAB performs vector and matrix computations quite eciently, has mostly fully mutable datatypes, and with recent just-in-time (JIT) compilation optimizations vector resizing performance has improved to what in many instances is perceived as equivalent to preallocated variables. This memory management abstraction, however, also means that explicit control of when arguments are passed by value or by reference to a function is not ocially supported. The following performance ramications naturally increase with the size of the data sets (N) passed as arguments and become quite noticeable even at moderate values of N when dealing with data visualization, a key function in system. To circumvent these problems explicit data references were implemented using some of the undocumented functions of MATLAB's libmxlibrary together with a custom data visualization function.

The main parts of the near real time interferogram analysis are resampling and a Fourier transformation, both of which had functionally complete but not optimized im- plementations. The minimal requirement for the reimplementation of these were simply to improve eciency while maintaining output precision.

On experimentally obtained data the new system's (DAQS) resampling implementation increased sample throughput by a factor of 19 which in the setup used corresponds to 10⁸ samples per second. Memory usage was decreased by 72% or in terms of the theo- retical minimum from a factor 7.1 to 2.0. Due to structural changes in the sequence of execution DAQS has no corresponding implementation of the reference FFT function as the computations performed in it have been parallelized and/or are only executed on de- mand, their combined CPU-time can however in a worst-case scenario reach 75% of that of the reference. The data visualization performance increase (compared to MATLAB's own, as the old system used LabVIEW) depends on the size in pixels of the surface it is visualized on and N, decreasing with the former and increasing with the latter. In the baseline case of a default surface size of 434x342 pixels and N corresponding to one full sweep of the FTS's retarder DAQS oers a 100x speed-up to the Windows 7 version of MATLAB R2014b's^plot.

Secondary to the above was the implementation of acquisition and analysis for both directions of the retarder, a HITRAN reference spectra generator, and functionality to improve the user experience (UX). The rst, though computation time allows for it, has not been implemented due to a delay in the DAC-driver. To provide a generic implementation of the second, the HITRAN database was converted from the text-based format it is distributed in to a MySQL database, a wrapper class providing frequency- span selection and the absorption spectra generation was developed together with a graphical front-end. Finally the improved UX functionality mainly focused on providing easy-access documentation of the properties of the DAC.

In summation, though the primary objectives of optimizing the data analysis func- tions were reached, the end product still requires a new driver for the DAC to provide the full functionality of the reference implementation as the existing one is simply too slow.

Many of DAQS' components can however be used as stand-alone classes and functions until a new driver is available. It is also worth mentioning that National Instruments (NI), the DAC vendor, has according to their technical support no plans to develop native MATLAB drivers as MathWorks will not sell them licenses.

Optical Frequency Comb Spectroscopy (OFCS) Gruppen vid Umeå Universitet utvecklar OFC-baserade tekniker för extremkänslig gasdetektering, beröringsfri temperaturmät- ning, och andra problemspecika applikationer. Flera av de centrala komponenterna i gruppens uppställning för gasdetektion är designade och byggda på plats, bland dessa ingår Fourierspektrometern (FTS) och en auto-balanserande detektor. Det är den, som med sitt högfrekvent datainsamlingskort (DAC) från National Instruments (NI) produc- erar interferogram i storleksordningen 10⁷ yttal med dubbel precision per svep med FTSns retroreektor, har använts i examensarbetet. Interferogrammen behöver anal- yseras i semirealtid (<1 sekund) för att göra kontinuerlig insamling och veriering av signalen möjlig. Målet med examensarbetet har således varit att i MATLAB utveckla en optimerad applikation för just datainsamling och analys (DAQS). Valet av program- språk var givet av gruppen för att säkerställa långsiktigt underhåll av mjukvaran, då alla gruppens medlemmar är väl bekanta med det. Beräkningsmässigt är MATLAB indus- triledande med mycket eektiva vektor- och matrisoperationer, muterbara datatyper, och i senare versioner har just-in-time (JIT) kompileringen eektiviserat iterationsbaserad storleksanpassning av både matriser och vektorer till, för vad en användare ofta upplevs som likvärdigt förallokerade variablers prestanda. Det senare har emellertid medfört en abstraherad minneshanteringsmodell utan ociellt stöd för explicit kontroll av när funk- tionsargument skickas som referenser eller direkta värden. Det här gör att kods prestanda i termer av exekveringstid i värsta fall degraderas linjärt med storleken på dataseten (N) som behandlas, något som vid datavisualisering blir noterbart redan vid moderata vär- den på N. För att kringgå denna problematik utvecklades explicita datareferenser (att jämföra med pekare i C) med hjälp av odokumenterade funktioner i MATLABs ^libmx- bibliotek och en skräddarsydd visualiseringsklass.

För interferogramanalysen tillhandahölls referensimplementationer av omsamplingsru- tinen och den nita, diskreta Fouriertranformationen (FFT), vilka gav en övre gräns för exekveringstid för de nya implementationerna.

På experimentellt erhållen data ökade DAQSs omsamplingsimplementation genom-

ödet med en faktor 19, vilket i den aktuella uppställningen motsvarar omsamplings- frekvens på 10⁸ yttal per sekund. Minnesanvändningen minskade med 72% eller från 7.1 till 2.0 i termer av det teoretiska minimumet. Till följd av en omstrukturerad exekver- ingssekvens och ökad funktionalitet i den nya datavisualiseringsklassen saknar DAQS en Fouriertransformfunktion motsvarande den i referensimplementationen, FFTn sker istäl- let bara vid behov och era delar har parallelliserats. Vid manuellt sekvensierad exekver- ing uppgick dock CPU-tiden för motsvarande funktionalitet i värstafallsscenariot till 75%

av referensen. Som referens till visualiseringsklassen användes MATLAB R2014bs^plot- funktion i Windows 7, där de jämförda tiderna var de från funktionsanrop till färdigritad

gur. Till skillnad från^plotvarierar klassens tider inte bara med N utan bara med stor- leken i pixlar på visualiseringsytan. Vid, av MATLAB given, standardstorlek på ytan och N = 18.5 · 10⁶ (motsvarande ett svep av retroreektorn) visualiserar klassen datan 100x snabbare än ^plot.

I tillägg de primära formulerades även tre sekundära mål; kontinuerlig inhämtning av data för retroreektorns båda färdriktningar, en HITRAN-baserad referensabsorptionssp- ektrumgenerator, och funktionalitet för ökad användarvänlighet. Det förstnämnda visade sig, trots tillräckligt reducerad analystid, inte möjligt att implementera med den bentliga drivrutinen för DACet. För att maximera användbarheten av generatorn konverter- ades HITRAN-databasen till MySQL-format, varefter ett abstraheringslager i form av en MATLAB-klass med metoder för val av absorptionslinjer och spektrumgenerering med ett tillhörande graskt gränssnitt (GUI) utvecklades. Slutligen bestod förbättringarna av användarupplevelsen till stor del av omstrukturering och gruppering av GUI-element samt integration av dokumentation (hjälp-knappar) i GUIt.

Sammanfattningsvis kräver den, för projektet slutgiltiga produkten fortfarande en ny drivrutin för DACet för att uppnå likvärdig funktionalitet med referensimplementatio- nen. Många av DAQSs komponenter kan emellertid användas som fristående klasser och funktioner fram till dess. Nämnvärt är också att NI enligt egen utsago inte har några planer på att utveckla en MATLAB-specik drivrutin då MathWorks inte säljer licenser till dem.

Foltynowicz-Matyba and Ove Axner for their seemingly endless patience, help during the work itself and for letting me stray and explore many of the interesting aspects of this thesis work.

Further I would like to thank Alexandra Johansson for her invaluable input, Olof Rydberg and Edvin Norén for great company during the work, and Elisabeth Wetzer for all the help, encouragement, and proofreading.

Finally this work would never even have begun without the inspiration and help always given by my brother Arvid Åhrberg.

1.1 Background . . . 1

1.2 Scope . . . 1

1.3 Objectives . . . 2

1.4 Thesis outline . . . 2

2 Theory 4 2.1 Fourier Transform Spectrometer . . . 5

2.1.1 Optical Frequency Comb FTS . . . 6

2.2 Resampling . . . 7

2.3 Data reduction/interpolation for visualization . . . 8

2.4 Molecular absorption spectroscopy and HITRAN . . . 10

2.4.1 Lorentzian lineshape . . . 10

2.4.2 Gaussian lineshape . . . 11

2.4.3 Voigt prole . . . 11

2.4.4 van Vleck-Weisskopf lineshape . . . 12

2.4.5 Numerical evaluation . . . 12

3 Method 14 3.1 MATLAB internals . . . 14

3.1.1 Java back end . . . 14

3.1.2 C representation . . . 15

3.2 Design patterns . . . 16

3.2.1 CoW and function arguments . . . 17

3.2.2 Preallocation . . . 17

3.2.3 Boolean object evaluation . . . 19

3.3 Parallelization . . . 20

3.3.1 MATLAB local pool . . . 20

3.3.2 GPU . . . 21

3.4 Notes on optimization . . . 21

3.5 MATLAB GUI programming . . . 22

3.5.1 GUI layout management . . . 22

3.5.2 GUI back end versions . . . 23

3.6 Data visualization . . . 24

4 Implementation 25 4.1 Previous system . . . 25

4.2 New system - DAQS . . . 26

4.2.1 GUI . . . 28

4.2.2 Usage and structure . . . 28

4.2.3 HITRAN . . . 29

4.2.4 MATLAB functionality extensions . . . 29

5.1 Primary objectives . . . 32

5.1.1 Resampling . . . 33

5.1.2 FFT . . . 34

5.2 Secondary objectives . . . 34

5.2.1 HITRAN interface . . . 35

5.3 End product . . . 35

6 Discussion 37 6.1 HITRAN . . . 37

6.2 Code optimization . . . 37

6.3 MATLAB . . . 39

References 40 Appendices 42 A Analytical interferogram derivations 42 A.1 Single frequency interferogram . . . 42

A.2 Multiple frequencies . . . 43

B HITRAN data 45

C mxArray and libmx 48

D Files and descriptions 49

1 Introduction

1.1 Background

Real time data acquisition and analysis have today become a commodity, implemented in everything from augmented reality functions in smart phones to oscilloscopes used by experimental physicists. The tools to implement such functionality rarely oer easy customization though, in consumer markets the vast majority of applications share a common need which is just variations of pattern recognition within image data. In scientic elds applications are as diverse as the term indicates, even when narrowing the scope to spectroscopic ones.

The Erbium-ber OFC setup of the OFCS Group at the Department of Physics at Umeå University includes several unique, in-house developed components that acquire absorption spectra in a fraction of the times required by commercially available systems.

The digitizer used, a National Instruments PCI-5922, is capable of acquiring data at a rate of 15 MS/s per channel with a sample rate dependent resolution ranging from 16 to 24 bit. This requires a real time acquisition and analysis software to be able to perform at least a magnitude more single-precision oating point operations per second, while also visualizing the data continuously.

The initial implementation of this software was based on LabVIEW, using its MATLAB- interconnectivity to implement the computationally intense parts of the analysis. While the interoperability of the two are highly useful it imposes some restrictions on execution time and exchanged data set size (N) as arguments from one to the other are always passed by value. The combination of this, data set sizes in the 10⁷-range, the fact that MATLAB by default uses double precision number representation, and that explicit in- vocations of garbage collection are strongly discouraged in both languages resulted in the need for a new, more optimized implementation.

1.2 Scope

Though the end product is specically designed to t the needs of one spectroscopic setup, the Noise Immune Cavity Enhanced OFCS (NICE-OFCS), it is deliberately developed in a generic fashion to be reusable in similar settings. Hence, only the components of the current setup that are needed to produce the right input for the data acquisition and analysis system (DAQS), and their governing theory, are covered in this report. For a comprehensive description of the setup see [1].

Basic programming skills in MATLAB, a general understanding of object oriented concepts and laser-based spectroscopy as well as familiarity with the more commonly used components therein is expected by the reader. It should be noted that, though this is covered by many text books the nomenclature in this thesis closely resembles that in [2].

As the main body of work consisted of programming in MATLAB, emphasis has been put on the language's strengths and weaknesses in these types of applications and the design choices made to facilitate high cohesion, reusability, and maintainability in

performance critical settings.

1.3 Objectives

The primary objectives for the thesis work have been to implement the following in MATLAB

• A graphical conguration and data acquisition utility for the DAC

• Stand-alone data analysis functions which can be invoked during acquisition

• To-disk saving of raw acquisition data and analysis results

• FTS retarder motor control

The minimum requirement for the above was to perform their individual tasks in less time but with identical numerical results (where applicable) as the reference implementation, for which the end product should work as drop in replacement. This includes not only equivalent functionality but also a user experience (UX) that ts into the workow of the spectroscopy group.

With the primary fullled and given time, improvements in functionality was set as the secondary objectives. Tasks given were

Conguration proles

Saving and loading sets of DAC-conguration parameters to and from named pro-

les.

Two way acquisition

Given sucient data analysis time reduction the acquisition and analysis should run continuously, i.e. in both directions of movement of the retarder.

HITRAN reference spectra

From the parameters of an acquisition a reference absorption-spectrum should be generated from HITRAN-data.

In addition to these, the general goal given for the project was maintainability. In combination with the need for high performance this boiled down to a substantial amount of documentation.

1.4 Thesis outline Theory

The section covers the algorithms and supercially the theory governing the data analysis, data-reduction for visualization, and HITRAN-based numerical evaluation of the basic theory of absorption spectroscopy. The main focus of the last part is to describe it from an implementational perspective, hence constants are generally separated in equations where they in other literature would not be.

Method

Section 3 describes the development environments and tools in terms of their strengths and weaknesses. It also includes subsections covering some of the special considerations which need to be taken when working with MATLAB in performance critical situations.

Implementation

The section includes a brief description of the previous system and the logical struc- ture of the acquisition and analysis before covering the implementation of DAQS.

This is followed by a comparison between the two. A more in-depth description of the intricacies DAQS is available in [3], which should be seen as a complement to the thesis in general and this section in particular.

Results

The most signicant results of performance improvements are presented in the

rst part of this section in terms of performance relative to the reference imple- mentation, with only short examples of the absolute values as these are entirely hardware-dependent. The secondary objectives, end product, and a description of the precision in the calculations performed in DAQS then follows.

Discussion

Thoughts and discussions covering primarily the project itself, its key components and concepts. It includes a short description of the spin-o project for HITRAN- interaction, the eects of MATLAB's high-level programming concepts, and a more general discussion of code optimization in scientic projects.

2 Theory

An OFC can in many absorption based spectroscopic techniques be viewed as the equiva- lent of using thousands of synchronized laser lines in a measurement. The output, a train of extremely short pulses, is produced by locking together the phases of a large number of longitudinal modes of the laser. This technique is ttingly called mode-locking and the details of it can be found in [4]. In gure 2.1, a small subset of pulse train is depicted together with the frequency domain representation of it. The key properties dened in the gure are the oset frequency,

ω₀= −∆φ

tr (1)

where ∆φ is the pulse-to-pulse carrier-envelope phase shift and tr the pulse separation in time, and the angular-frequency spacing, or angular repetition frequency

ωr= 2π tr

. (2)

t

E(t)

t_r

→ ← ∆φ

→ ← τ

/ / / /

ω

I(ω)

→ ← ω0

→ ← ω^r

→ ← δω

/ /

Figure 2.1: Mode-locked pulse train of an OFC in time and frequency domain. The oset frequency, ω0, corresponds to the carrier envelope phase shift, ∆φ, in the time domain and can be seen in the OFCs spectral envelope as well.

Given these two, a frequency

ωn= ω0+ nωr (3)

where n is any positive integer such that ωn falls within the spectral envelope, shown in the frequency domain in gure 2.1, can be obtained. The FWHM of this is dened by the temporal duration of each pulse τ as

δω = 2π/τ (4)

and thus the spectral output of the OFC. Knowing and being able to control these param- eters constitutes the lower limit of what is needed to employ an OFC-based setup together with DAQS, accomplished in the experimental setup with a MenloSystems FC1500-250- WG with a repetition rate of 250 MHz and output in the near infrared (NIR) wavelength range (1.5-1.6 µm).

2.1 Fourier Transform Spectrometer

A FTS is based on the Michelson interferometer (MI), but instead of having one xed and one moving mirror both the reected and transmitted part of the beam are routed to a double-sided retroreector mounted on a rail as shown in gure 2.2a. The displacement of the retarder,

∆L = 4δ_{OP D} (5)

thereby corresponds to four times the change in optical path dierence (OPD), δOP D.

OFC HeNe

T R + RT T T + RR

I (δOP D) ^BS ∆L

(a) Fourier transform spectrometer.

t#λ 4v

$ I(δOPD)[I0]

0 1

(b) Theoretical HeNe inteferogram.

Figure 2.2: Figure (a) shows the layout of an FTS with top output consisting of through the beam splitter (BS) twice transmitted (T T in the gure) and twice reected (RR) OFC beams, middle output of rst transmitted then reected (T R) and vice versa (RT ). (b) shows a HeNe interferogram, lower output, given constant retarder velocity and δOP D⁽⁰⁾ = 0.

As seen in 2.2a the two output beams, not counting the HeNe reference, will have a relative phase-shift of π. The upper (T T + RR in gure 2.2a) with an odd number of reections i.e. ∆φup= (2n + 1)π = π and the lower (middle, T R + RT ) with ∆φdown= 2nπ = 0. When one of these are subtracted from the other, as is done in the current setup by an auto-balancing detector, common-mode noise is considerably reduced while the resulting maximum amplitude corresponds to I0 − (−I₀) = 2I₀, where I0 is the intensity incident on the detector.

The interferogram shown in gure 2.2b has twice the frequency of one from a corre- sponding MI. This can be deduced by equation 5, the constructive interference condition,

δOP D = nλ, (6)

and the destructive one, δ_{OP D} =

 n + 1

2



λ (7)

where n is any positive integer and λ the wavelength. Given that ∆L = v∆t, where v is the retarder velocity, the frequency of the interferogram, fI, will be

fI≡ f_I^{F T S} = 2f_I^{M I} = 4v

λ. (8)

Before moving on to the OFC interferogram it is worth considering the practical impli- cations of fI. As mentioned earlier the DAC operates at discrete sampling-frequencies.

To ensure valid measurement results the Nyquist-Shannon sampling theorem, given in [5], gives an upper limit for the retarders velocity as

v_max= λf_S

8 (9)

where fS is the sample rate.

With equation 8 the interferogram for a single wavelength beam such as the one from a HeNe laser is simply a cosine-wave with a linearly v-dependent frequency and a phase shift, ∆φ, determined by the initial OPD. The eect of a constant versus non-constant retarder velocity can be seen by comparing gure 2.2b and 2.4a, and b (dotted line only).

2.1.1 Optical Frequency Comb FTS

Starting with the expression for a single frequency MI interferogram with a perfect beam- splitter (BS), as described in [2], and an arbitrary initial OPD the expression for the transmitted intensity can be expressed in terms described by equations 6 and 8 as fol- lows,

I (δOP D) = 2RT I0[1 + cos (2πfIt + ∆φ)] (10a)

= 2RT I₀

 1 + cos

 2π4v

λt + ∆φ



(10b)

= 2RT I₀[1 + cos (2π4δ_{OP D}(t) + ∆φ)] (10c)

where R and T are the reection and transmission coecients of the BS. Above, with a perfect one, both are 0.5. Treating an OFC simply as a sum of multiple single frequency lasers results in its FTS interferogram becoming the sum of equation 10, as seen in gure 2.3,

I (δ_{OP D}) = 2RT I₀

N

X

n

[1 + cos (2πf_Int + ∆φ_n)] (11a)

= 2RT I0 N

X

n

 1 + cos

 2π4v

λ_nt + ∆φn



(11b)

= 2RT I0 N

X

n

[1 + cos (2π4δOP Dn(t) + ∆φn)] (11c) where N is the number of frequencies, and ∆φn the pulse-to-pulse carrier envelope phase shift. For a more complete, analytical derivation of equations 10 and 11 refer to appendix A.

No.Frequencies,N

4⁰ 4¹ 4² 4³

Multiples of (dt · frep)⁻¹

0 1 2 3 4 5 6 7 8 9 10

(a) Summation of multiple single frequency signals with normalized amplitudes.

Samples _×10⁶

0 2 4 6 8 10 12 14 16 18

I[V]

-5 -4 -3 -2 -1 0 1 2 3 4 5

(b) Experimental OFC interferogram

Figure 2.3: Comparison of summations of single frequency signals, (a), and experimental OFC interfer- ogram, (a).

2.2 Resampling

Due to mechanical limitations the retarder cannot run at a constant velocity throughout a sweep, which is why a reference beam is needed. With a well dened reference wavelength, λ = 632.9914 nm, for HeNe and a constant sample rate, fS, the velocity of the retarder, v, may be unknown as long as it is less than vmax, as given by equation 9. With samples of a sweep separated by a constant ∆t = 1/fS, and a zero-mean adjusted HeNe interferogram, two adjacent (interpolated) zero-crossings will be separated by λ/2 in the spatial domain. Using these together with the extreme points described by equation

6 and 7 the OFC interferogram can be resampled to a constant sample separation of λ/4with algorithm 2.1¹, thus fullling the equispaced samples requirement of the nite discrete Fourier transform.

Algorithm 2.1: Resampling of OFC interferogram¹ 1. in: H(t), Zero mean HeNe interferogram

in: D(t), OFC interferogram

out: R(λ), Resampled OFC interferogram

2. Fetch indices of zero-crossing pairs (green and black dots in g. 2.4b) i : (sgn (Hi) +sgn (Hi+1) = 0)

i⁰ := i + 1

3. Calculate fractional indices of zero-crossings (red x's) d^z_i := i +_H^Hi

i−H_i0

4. Interpolate OFC interferogram values at zero-crossings R_i := D_i− |d^z_i − i| (D_i− D_i⁰)

5. Calculate fractional indices of local extreme values (magenta x's) d^x := ^d

z i+d^z_i−1

2

6. Determine indices of local extreme pairs, assuming adjacency to d^x k :|H_k| = max 

Hbd^xc

 ,

Hdd^xe



k⁰:= k − 1 + 2 (dd^xe − k)

7. Interpolate values at local extremes R_k:= D_k− |d^x− k| (D_k− D_k⁰)

t v[v0]

(a) Retarder velocity prole.

t I[I0]

(b) HeNe inteferogram.

δ_{OP D}[λ]

I[I0]

(c) Resample points.

Figure 2.4: Example stages from FTS reection to resample points. Notice the HeNe interferograms dependence on the retarder-velocity and the resample points' lack thereof.

2.3 Data reduction/interpolation for visualization

Rendering large datasets using the built in functions in MATLAB, even when hardware- accelerated graphics are supported, is a time consuming task independent of the size

1Here ":" should be read as such that and ":=" as is equal to.

of the surface to visualize it on. In cases where the number of samples in a data set, N, is larger than what can be visually represented with a line plot manually reducing the set is - in the experimental setup - far more ecient than letting MATLAB render it in its entirety.² Given a data visualization implementation expecting surface size- dependent data sizes as input the reverse situation (N smaller than available surface allows) will in certain instances introduce an overhead caused by value interpolation. In the experimental setup however, this overhead is negligible. With the following data set reshaping algorithm the time required from function call to actual rendered data is thus signicantly reduced in most use cases.

Algorithm 2.2: Pre-line plot visualization data set reshaping 1. in: I, data set to be visualized

in: N, length of data set

in: M, desired length of output data out: O, reduced or interpolated input data 2. If M = N:

a. Set output to input O := I

3. Else if M > N:

a. Evenly distribute I in O n := {1, 2, 3, . . . , N } m := n +

j

(n − 1)^{M −N}_{N −1} k O_m := I_n

b. Linearly interpolate values for empty, forward adjacent elements in O (nal step for cases where M ≤ 2N − 1)

i : (mi+1− m_i> 1)

O_mi+1 := O_mi +^O_m^mi+1−Omi

i+1−m_i−1

c. Interpolate remaining values (if any) m⁰_j :=j : j ≤ M, j ∈ N⁺\{m}

O_m⁰

j+1 := O_m⁰

j+ O_m⁰

j− O_m⁰

j−1

 4. Otherwise:

a. Set O to the max and min values of all contiguous subsets of I of size N/M n_ij :=i : i < j_M^N, i ∈ N⁺,j : j_M^N ≤ N, j ∈ N⁺

m :=0, 1, . . . ,^M₂ − 1 O_2m:= max

I_nij O_2m+1:= min

I_nij

2While this is also true for scatter plots only line plots were implemented in DAQS.

Note that the implementation of the algorithm in most programming languages requires checking and calculating values for "residual" indices in both step 3 and 4.

2.4 Molecular absorption spectroscopy and HITRAN The Lambert-Beer law,

I_A= I₀e^−α(ν) (12)

states that the intensity of light transmitted through an absorbing sample, IAdecreases relative to the incident intensity I0 by the exponential factor

α (ν) = ηSpcrelLχ (ν) (13)

η = 10⁻⁶n₀T₀/T

where η is, given a constant temperature, T , a scaling constant [molecules / (cm³·atm)]

containing the Loschmidt constant, n0 [molecules / m³], and T0 = 273.15 K. S, the integrated molecular linestrength [cm⁻¹ / (molecule·cm⁻²)] is given in HITRAN, p is the pressure [atm], crel the relative concentration of the species, L the length of the sample [cm], and χ (ν) is the area-normalized absorption lineshape function [1/cm⁻¹].

The last term describes the broadening of the observed intensity dip, mainly caused by the movement and internal interactions of the molecules.

For gaseous molecules and incident light in the NIR spectrum the absorption causing the amplitude of the dip in transmitted energy is primarily caused by transitions in the rotational and vibrational (ro-vibrational) states of the molecules. The number of such transitions depends on the number of atoms making up the molecule, giving a nonlinear molecule a maximum of 3n degrees of freedom to move, n being the number of atoms.

Three of these represent translational motion in mutually perpendicular directions, an- other three the rotational motion about the x, y, and z axis of the molecule, and the remaining 3n − 6 the number of ways atoms in a nonlinear molecule can vibrate. For detailed information on these vibrational modes and Fourier transform spectrometry in the infrared spectra see [2].

The molecular properties governing ro-vibrational absorption are also documented in HITRAN (a list of all the properties in it is available in table B.1). It currently contains in excess of seven million records, spanning 47 dierent species of gaseous molecules and their most signicant isotopologoues. With these and the exponent from equation 12 given by equation 13 theoretical absorption spectra from any of the 47 species can be generated.

2.4.1 Lorentzian lineshape

When pressure broadening dominates, i.e. when the lifetime of molecular states are shortened mainly by collisions, as is the case at room temperature and pressures above 0.1 [atm] (altitudes below 16 km), the lineshape is described by a Cauchy distribution,

in spectroscopic applications called the Lorentzian lineshape, χ_L(ν) = 1

π

δν_L

(νd)²+ (δνL)² (14)

ν_d= ν − (ν0+ δ_airp)

where νd is the detuning frequency [cm⁻¹], the transition line frequency, ν0 [cm⁻¹], and the pressure-induced line shift, δair are both given in HITRAN. The Lorentzian prole half-width-half-maximum (HWHM) is calculated using equation A12 from [6] as follows

δν_L= T_ref T

n_air

p (γ_air(1 − c_rel) + γ_selfc_rel) (15) where Tref = 296 K, and both the air broadening and self broadening coecients, γair

and γ_selfrespectively, are found in HITRAN. When lacking other data the classical value, as suggested by [7], of the temperature dependence exponent n_air= 0.5can be used.

2.4.2 Gaussian lineshape

In lower pressure situations where the main broadening mechanism is the Doppler eect experienced by the molecules; having a velocity component vk, parallel to the direction of propagation of the light will cause them to interact with a frequency ν0⁰ = ν₀ 1 + v_k/c
Doppler broadening causes a Gaussian lineshape, .

χ_G(ν) = k^χG δν_D exp

"

−

 ν_d δν_D

2#

(16) k^χG = 2p

ln (2) /π

where δνD, the HWHM (in this case) also known as the half Doppler width, is given by δνD = k^δνDν0

pT /M (17)

k^δνD = p2 ln (2) 10³R c

where R is the molar gas constant. The factor 10³comes from R's unit being J / (mol·K)

= m²·kg/(mol·s²·K) and the molar mass, M, being given in units of g/mol.

2.4.3 Voigt prole

In settings where both Doppler and pressure broadening are inuential a convolution of the two named after German physicist Woldemar Voigt is used. Due to the computational cost of convolutions Whitings approximation is used in DAQS, and with the assumption α (ν)  1equation 13 is reduced to

α (ν) = I0− I_T I0

= ∆I I0

= σVL (18)

where σV is given by

σV (ν) = σV (ν0) (1 − x) exp −k^σ₁^Vy²
+ x 1 + y² + k₂^σVx(1 − x)



exp −k^σ₃^Vy^2.25
− 10 10 + k₄^σVy^2.25



(19) x = δν_L

δν_V y = |ν − ν₀|

δν_V = |ν_d| δν_V

σ_V (ν₀) = S

2δν_V (k^σ₅^V + k₆^σVx + k^σ₇^Vx²)

where the constants k_i^σV, i = 1, ..., 7, given in table 2.1, have been adjusted from the original approximation in [8] to account for δνV as HWHM which is given by [9] as such

δν_V = k^δν₁ ^Vδν_L+ q

k₂^δνV (δν_L)²+ (δν_D)² (20)

where, again, k^δν₁ ^V = 0.5346and k₂^δνV = 0.2166have been adjusted. The approximation reduces to the Gaussian prole as δνL→ 0 and is always accurate to within 0.01%.

Table 2.1: k^σ_i^V values

i 1 2 3 4 5 6 7

Value 0.6930 0.0160 0.0841 0.2100 1.065 0.447 0.058

2.4.4 van Vleck-Weisskopf lineshape

The three lineshape functions listed above are by far the most commonly used ones. In the special case of Lorentzian lineshape and extremely long wavelengths, ν0 < 200cm⁻¹, the van Vleck-Weisskopf lineshape oers a slight improvement, as described in [10], and is given by

χ_{V W}(ν) = δν_L π

 ν ν0

2

1

(ν_d)²+ (δνL)² + 1

(ν + ν0)²+ (δνL)²



(21)

2.4.5 Numerical evaluation

To implement the above for a uniformly distributed set of wavenumbers W , of size n and corresponding HITRAN records H, where W0 ≤ H₀^ν0 and Wn≥ H_m^ν0 where m is the total number of records in H, the following algorithm has been used³. Note that for an accurate approximation of the absorbance spectrum A, n ≥ 100m and n  H_m^ν0 − H₀^ν0 needs to be true.

3For brevity only references to equations are included.

Algorithm 2.3: HITRAN based absorption spectra generator³ 1. in: W , uniformly distributed set of wavenumbers [cm⁻¹]

in: H, set of HITRAN records out: A, absorbance spectrum

2. Calculate δνL and δνD as given in equation 15 and 17 from H 3. Calculate δνV as given in equation 20 from table B.2

4. Calculate x and σV (ν₀) as given in equation 19 from H 5. For all records i in H:

a. Set νdi = W − H_i^ν0

b. Calculate yi as given in equation 19 from Hi and δνVi

c. Evaluate and add the results of equation 18 and 19 to A

3 Method

As mentioned in the introduction the DAQS is implemented in the numeric comput- ing environment MATLAB which employs a fourth-generation programming language.

The generational division of programming languages is not necessarily chronological, instead it is based on abstraction. Each generation adds a layer of abstraction onto the hardware, i.e. rst-generation languages are pure machine code, historically pro- grammed with punch cards but rarely used today. Assembly languages are categorized as second-generation, more easily read and written by humans but still architecture- specic. Most general purpose languages are third-generation meaning they are not machine dependent and the host on which it runs or is being compiled generates the machine instructions. Fourth-generation languages do not only abstract hardware but also many concepts such as data representation, storage, mathematical optimization, and are often domain-specic. The layers of abstraction introduced should ideally produce no or very little overhead in terms of execution speed and memory consumption. In MATLAB however this is not true for every aspect, hence the following section describes in general terms the methods used to minimize the overhead introduced in various parts of the implementation of DAQS as well as other measures to optimize MATLAB code.

3.1 MATLAB internals

With execution time being integral to DAQS simply optimizing function calls and code structure proved to be insucient. To increase speed even further the internal structure of MATLAB had to be examined, a not entirely trivial task given the age and maturity of it. The internals have been in development for over 30 years and have during this time provided a relatively consistent interface to many of the best performing and most versatile math-libraries, including MathWorks' own. This requires abstraction not only towards the end user but internally as well. When the dierent back ends of MATLAB are not taken into account the overhead created by the internal abstraction layers can, as noted during the development of DAQS, become quite substantial. One of the most obvious cases of this is when large amounts of data has to be copied due to graphical user interface (GUI) related actions. To avoid this, as the ocial documentation on the subject is rather sparse, the structure of MATLAB's Java back end, which is responsible for GUI-interaction, had to be examined and contrasted with the non-GUI related ones.

3.1.1 Java back end

Since version 6.0, released in November of 2000, MATLAB comes bundled with a pre- congured Java Virtual Machine (JVM). It uses this on hosts where newer versions are available as well, an understandable design choice even with the obvious drawback of missing new features.⁴ With it cross-platform GUI-programming capabilities based on Java's AWT and Swing toolkits are supplied.

4This has rarely been a problem with the notable exception of rendering on new, high resolution/high pixel density monitors.

MATLAB's GUI-components are primarily wrappers for their Java toolkit counter- parts. These, though not ocially available to the programmer are quite easily accessible using built-in functions only. Modifying GUI components by accessing their underlying Java-controls can be useful when only minor changes in the components functionality re- quires major additions to its wrapper's callback. However, since internal structures and toolkit-versions change with MATLAB-releases this is not a practice that should be incor- porated in any distributed project. Instead the still undocumented, but since the initial release supported, Java component constructor family of functions should be used. Given that many toolboxes, including the ones shipped with the base installation of MATLAB rely heavily on them and that they have stayed consistent over the recent major revision update of the graphical back end, the risk of changes in their syntax or functionality can be considered quite small in the foreseeable future. The main advantage of the con- structors are that they, unless explicitly stated otherwise, create Java-components and a generic MATLAB parent component using the toolkits distributed with MATLAB in the structure used by the running instance.

Creating and interacting with Java-components as opposed to their wrappers is in prac- tice usually only necessary when major changes in functionality of an existing or a com- pletely new type of component is needed. An illustrative example of the former is the uiedit-class and its autocomplete-functionality. The lack of control of caret position and text-selection in a uicontrol(’Style’, ’edit’)-control makes inline autocomple- tion quite tedious to implement, with the underlying Swing ^JTextField though it is a fairly straight forward process. With the extra control that using pure Java controls grants any operation that may be native to the MATLAB control can be easily cus- tomized.

3.1.2 C representation

While the technique described above is useful to extend GUI-functionality its use in performance-tuning, primarily by minimizing data copy operations, requires at least a ba- sic understanding of MATLAB's internal data representation. In versions before R2013b when developing mex-les, i.e. compiling C/C++ or Fortran (not covered here) code into functions callable from MATLAB, the C structure-representation, called ^mxArray, of arguments and outputs were, though opaque, accessible to the programmer. The intention of having an opaque datatype represent variables and provide an application programming interface (API) mimicking the MATLAB language is, as in the case of the language itself, to abstract the internal workings and thus have the ability to change it without breaking legacy code support. It does however also impose many of the re- straints present in such a high-level language on an earlier generation language. For this reason, investigative purposes, and pure academic curiosity Peter Boettcher and Gareth Jones documented [11] and supply a header le giving direct access to the elds of the mxArray, a subset of which are listed in table C.1. From their work details on MAT- LAB's copy-on-write (CoW) mechanism can be discerned. The direct eects of it can be examined using the undocumented ^-memoryoption of ^profile as illustrated by listing 3.1 and table 3.1.

Listing 3.1: MATLAB CoW

% cow_profiler.m

function cow_profiler(r, c)

x = rand(r, c); % initial allocation

y = x; % crosslink set, no allocation x(1, 1) = 0; % data copied to y and then modified end

% profiling call

>> profile(’-memory’, ’on’);

>> cow_profiler(1e4, 1e4);

>> profile(’viewer’);

Table 3.1: MATLAB CoW-related proling results

Line Code Total Time Allocated Memory

2 x = rand(r, c); 0.802 s 782780 Kb

3 ^{y = x;} 0 s 0 Kb

4 x(1, 1) = 0; 0.243 s 782780 Kb

As seen the data in ^x is not actually copied until one of the variables value is altered.

Implementationally this corresponds to updating the^crosslinkand^pdataelds of the mxArray; at line 2 the ^mxArray of ^x is created, at line 3 the same happens for ^y and the^crosslink of ^x is updated. Finally a new data segment is allocated, the original data in ^x is copied to it, ^ys ^pdatais updated and both their^crosslinks are reset.

Boettcher and Jones' work covers the MATLAB releases of 2012. In later versions the undocumented functions of the^libmx-library highly suggest that MathWorks has moved to a dynamic ^mxArray structure, a notion supported by tests carried out during this thesis work. Using the undocumented functions of^libmxis comparatively safer than di- rectly accessing the oblique structure and has made custom functionality less hazardous to implement.

Though experimental results strongly indicate the size of an ^mxArray in MATLAB R2014b running on the 64 bit version of Windows 7 to be 112 bytes and can thus be assumed to contain more elds than the 104 byte big one described by [11] the subset listed in table C.1 is sucient to extend mex-functionality to amongst other things imi- tate C-style pointers.

For consistency most of the code in DAQS is designed to avoid unnecessary memory al- locations and while this may not be apparent without knowledge of the CoW mechanism the data-acquisition loops rely heavily on it. The^plot-wrapper class ^screen2also uti- lize one-time allocations in combination with C-style pointers to provide single-allocation data visualization.

3.2 Design patterns

The introduction of just-in-time (JIT) compilation in R2011a has reduced the inherent structural performance problems of MATLAB's object-oriented programming language

but also added a few quirks which can seem incoherent, i.e. seemingly small and some- times insignicant changes in code having a signicant performance impact. A rather unfortunate eect of the object oriented structure which is not aected by JIT, is that implicit object evaluation can in the cases of code branching and variable initialization be considered non-trivial and can easily result in bugs that are hard to locate. Common practice in many similar languages is to always use explicit evaluation and enforce cer- tain code styles that improve readability and are in large considered good practice. Doing the latter often comes with a performance trade-o, this is especially true in the case of MATLAB. Performance critical code may therefore not only deviate from these stan- dards but do in some cases require a structure almost exact opposite of them. Following the guidelines given in [3] the cases in which this style is used have been almost entirely isolated to DAQS-specic les, resulting in what could be described as predictable code structure inconsistency.

3.2.1 CoW and function arguments

As briey described in section 3.1.2 the CoW mechanism is designed to reduce the number of memory allocations. This feature is particularly important to keep in mind when working with large data sets. To avoid unnecessary copies in-place data manipulation should be utilized when possible. This also relates to function arguments though in a slightly more involved way. MATLAB passes arguments to functions by reference when possible, that is when their values are not altered within the function. Hence the general rule would be to avoid assigning new values to input arguments, the most notable exception to this is small functions where in- and output are strongly related. As with anonymous functions it is possible to avoid allocating memory for variables within a function by using the same name for both input and output arguments as illustrated by listing 3.2. The denition of "small functions" in this context is functions where the size of local variables does not depend on N.

Listing 3.2: Minimizing memory usage in functions

function x = myFunction(x, y)

% temporary variable x is not allocated [x, z] = someOperation(x, y); % z is a scalar x = someOtherOperation(x, z);

end

It should be noted that a design optimized for large N often includes global or persistent variables, parent-scoped variables within nested functions and/or heavy use of properties of handle class objects. In small N situations this can result in degraded performance due to the overhead of this structure and it is therefore important to ensure at least an approximate magnitude of N before settling on it.

3.2.2 Preallocation

If there is no way to know the size of a vector at compile time or preallocation for some other reason is not possible the method used to grow said vector aects performance

profoundly, as seen in listing 3.3. In this case JIT has played a major role in performance boosting [12]. Depending on the magnitude of N on-the-y preallocation, i.e. dynami- cally increasing the vector's size by certain a factor when it runs out of space is generally a, performance-wise, safe choice. The other, earlier strongly discouraged ways of either using out of bounds indices or concatenation have in case of indexing, since R2011a increased in speed and are by the 2014-releases on par with most on-the-y methods.

Concatenation, however, is still severely lacking. The small dierence between^end-based and direct indexing is the result of MATLAB having to determine the length of the vector during each iteration in the former and, though the relative dierence in certain situa- tions can be quite big, the actual impact in a loop with any type of computation is in most cases negligible.

Listing 3.3: Dynamically growing vectors

>> data = []; y = 0;

>> tic();

while(y < 1e6)

data = [data, 1];

y = y + rand();

end toc()

Elapsed time is 5687.447569 seconds.

>> length(data) ans =

1999465

>> data = []; y = 0;

>> tic();

while(y < 1e6)

data(end + 1) = 1;

y = y + rand();

end toc()

Elapsed time is 0.438919 seconds.

>> length(data) ans =

1998533

>> data = []; y = 0; i = 1;

>> tic();

while(y < 1e6) data(i) = 1;

y = y + rand();

i = i + 1;

end toc()

Elapsed time is 0.433163 seconds.

>> length(data) ans =

2000989

3.2.3 Boolean object evaluation

As a consequence of being a weakly typed, computation-centric language MATLAB has a special way of evaluating objects as conditions. In many weakly typed languages argument validation are often either forcibly explicit through ^toBoolean-like methods or implicit as built in functions and operations will throw exceptions when supplied arguments of the wrong class. This is not the case for all built-ins in MATLAB, evaluating cells or^structs as booleans will throw exceptions but doing so with vectors or matrices will not. Unfortunately the evaluation of matrices diers depending on if it is empty, as shown in listing 3.4, or contains elements in which case ^if(A), where ^A is the matrix, is equivalent to if(all(A) && ~isempty(A)).

Listing 3.4: Boolean evaluation of empty vector

>> A = [];

>> if(A); disp(’A == true’); end

>> if(~A); disp(’A == false’); end

>> if(~A == A); disp(’A == undefined’); end

>> if(any(A)); disp(’At least one in A == true’); end

>> if(all(A)); disp(’All in A == true’); end All in A == true

Though the behavior of ^all is as documented in [13], the documentation is a bit con- tradictory, rst describing ^all as a natural extension of the logical AND operator but later stating that it will return logical true given an empty vector, a case for which the result of AND is strictly undened.

To circumvent some of the problems of direct object evaluation ^isempty is intended to supply a safe way to test if a variable is empty regardless of class. However, the denition of empty diers between classes as shown in table 3.2. When taking these aspects of MATLAB into account argument validation becomes a more strenuous task than in many other languages. One remedying quality with regards to this is that the language uses short-circuit condition evaluation, i.e. as soon as a condition is logically evaluated the evaluation stops. This way a qualifying condition can be placed in front of a class-specic one (this is the default style of user input based condition evaluation in DAQS).

Table 3.2: ^isempty return values for various arguments content argument

[] {} struct()^∗

Nothing true true false

[] true false false

{} false false true

struct() false false false

∗) Content given as value to key ’field’, except for Nothing.

3.3 Parallelization

Vector and matrix computations are in later version of MATLAB internally parallelized to the extent supported by the CPU and OS. The parallel toolbox expands this with means to perform explicit parallelization. In general this can be implemented using one of two methods, distributed (local or over a network) worker pools or by using a local graphics processing unit (GPU) with multiple cores. Both introduce an overhead in their respec- tive setup and teardown and are therefore mainly appropriate in computation-intense situations, i.e. not only (or necessarily) very large N but when element transformations are complex. The main reasons for these pre- and post-calculation overheads are dis- tribution and gathering of data, respectively. From an end user perspective these are method-agnostic actions, results computed in parallel, called composite variables, need to be merged into one variable per result. Implementationally, this process diers funda- mentally between the two, something that needs to be taken into account when selecting a method of parallelization.

3.3.1 MATLAB local pool

Distributed computations are based on pools of workers, locally Mathworks recommends one worker per logical CPU core, though similar performance can be obtained using nw = ⁿ_2m^cpu workers, where nw is the number of workers, ncpu number of logical cores and m ∈ N⁺. The main benet of running a local pool as opposed to letting MATLAB distribute the computations over the cores is that they will be utilized to a higher degree.

This does not automatically mean shorter execution time, on the contrary, since locally distributed workers in practice can be compared with several user interface-less MATLAB instances running on ⁿ_n^cpu_w cores each, unless the computations in the loop body are at least partially execution order independent and time-consuming the result is often the opposite as exemplied by listing 3.5.

Listing 3.5: ^parfor comparison with execution order dependent computation loop body

>> a = rand(1, 1e6);

>> b = rand(1, 1e6);

>> c = zeros(1, 1e6);

>> tic();

parfor i = 1 : length(a);

c = (a(i) + sum(b)) / (length(b) + 1);

end

d = gather(c);

toc()

% Elapsed time is 370.933364 seconds.

% Average CPU-load core 0: 100%, 1: 100%, 2: 100%, 3: 100%

>> tic();

for i = 1 : length(a);

d = (a(i) + sum(b)) / (length(b) + 1);

end toc()

% Elapsed time is 247.411576 seconds.

% Average CPU-load core 0: 97%, 1: 99%, 2: 97%, 3: 96%

The startup and shutdown of a local pool also introduce additional overheads, the for- mer normally a few seconds though these are one-shot costs. These factors make local pools suitable only in very specic situations and performance gains should, as with all optimization, be conrmed rather than assumed especially in distributed projects where the underlying hardware cannot be controlled.

3.3.2 GPU

Though technically no set up is required for GPU parallelization, checking for the pres- ence of a supported GPU is in any distributed project necessary and adds an overhead time-wise comparable to setting up a local pool of four workers. Implementationally, the parallel toolbox relies on CUDA, a parallel computing platform distributed by video card manufacturer NVIDIA, and is thus limited to their cards. Since even consumer grade video cards contain several hundred streaming cores, compared to eight or twelve logical ones in high-end CPUs, and CUDA is inherently designed for computation, performance gains are obtained in a wider range of situations than with local pools even though GPU cores in most cases operate at lower frequencies and with dierent instruction sets than their CPU counterparts.

Data distribution and gathering are, as it has to be copied to and from the actual video card, quite a bit slower than in-RAM copies. This also imposes the amount of RAM on the video card as the new restriction on address space, i.e. the amount of data that can be used for computations. The restriction diers slightly from the one when working directly in MATLAB as the OS employs several swapping and virtualization mechanisms to accommodate new allocations, making the theoretical maximum memory available larger than that of the physical RAM. For GPU-calculations the amount of RAM available on the video card is (at the moment) the maximum number of bytes available for calculations. However, with proper memory management and code that is designed with parallel execution in mind the GPU parallelization can be distributed with condence.

3.4 Notes on optimization

The general rule of "make it run, make it right, make it fast", while sound princi- pals in wide purpose application development, does not always apply in the same sense when working in domain specic environments and languages. When working with data- analysis in MATLAB a developer will, given a good understanding of the computa- tions performed, have a very good idea of where most of the processing time will be spent. Therefore, to learn and utilize the quirks of the language, and combine them with comments on the related code sections will produce code that is, within reasonable constraints, as maintainable as its sometimes substantially slower but more easily read counterparts.

In special cases with large "black boxes", wide purpose frameworks or toolboxes, the principle should be employed to its fullest extent. An example of this is rendering graph-